FFTW guru interface

Viewed 127

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;
}
1 Answers
Related