CODEWITHSUNDEEP
cfftw
asdfsdfdfd
asdfsdfdfd

Reputation: 23

FFTW guru interface

I don't understand the guru interface of FFTW. Let me explain how I thought it worked based on the manual and this question How to use fftw Guru interface and maybe someone can clear up my misunderstanding.

fftw_plan fftw_plan_guru64_dft(
 int rank, const fftw_iodim64 *dims,
 int howmany_rank, const fftw_iodim64 *howmany_dims,
 fftw_complex *in, fftw_complex *out,
 int sign, unsigned flags);

Suppose we want to calculate the DFT of interleaved multidimensional arrays, such as the six 2x2 arrays (each with a different colour) in this picture.

interleaved dfts

Because the dfts have stride 3 in the vertical direction, and stride 2 in the horizontal direction, I thought we would need rank = 2 and dims = {(2, 3, 3), (2, 2, 2)}. The starting points are a 3 x 2 subarray, so I thought howmany_rank = 2, howmany_dims = {(3, 1, 1), (2, 1, 1)}.

However, this is not actually what FFTW does. I made a smaller example that is easy to calculate by hand, consisting of 4 DFTs of size 2x1 (indicated by colours). Each dft is of the form (+-1, 0) which has as output (+-1, +-1), but that is not what FFTW calculates.

small example

Here is the code I used to calculate the DFT.

#include <stdio.h>
#include <stdlib.h>
#include <complex.h>
#include <math.h>
#include <fftw3.h>

int main()
{
    fftw_complex* X = fftw_malloc(8 * sizeof(fftw_complex));
    fftw_iodim* sizes = malloc(2 * sizeof(fftw_iodim));
    fftw_iodim* startingPoints = malloc(2 * sizeof(fftw_iodim));
    sizes[0].n = 2; sizes[0].is = 2; sizes[0].os = 2;
    sizes[1].n = 1; sizes[1].is = 2; sizes[1].os = 2;
    startingPoints[0].n = 2; startingPoints[0].is = 1; startingPoints[0].os = 1;
    startingPoints[1].n = 2; startingPoints[1].is = 1; startingPoints[1].os = 1;

    fftw_plan plan = fftw_plan_guru_dft(2, sizes, 2, startingPoints, X, X, FFTW_FORWARD, FFTW_ESTIMATE);

    X[0] = 1.0; X[1] = -1.0;
    X[2] = 1.0; X[3] = -1.0;
    X[4] = 0.0; X[5] = 0.0;
    X[6] = 0.0; X[7] = 0.0;

    fftw_execute(plan);

    printf("\nOutput in row-major order:\n");
    for (int i = 0; i < 8; i++) {
        printf("%lf + %lfi, ", creal(X[i]), cimag(X[i]));
    }

    return 0;
}

Upvotes: 2

Views: 219

Answers (1)

smitsyn
smitsyn

Reputation: 722

Strides even for major axes are in "units", i.e. doubles or fftw_complexes, not number of rows: https://www.fftw.org/fftw3_doc/Guru-vector-and-transform-sizes.html#Guru-vector-and-transform-sizes

My guess is that in major axis strides have to be multiplied by the distance between consecutive rows, also in units. So for the arrays their iodims.is and iodims.os strides should be 4*3 == 12.

Upvotes: 1

Related Questions