Question

I am wondering if it is possible to store the transposed matrix of a 1D FFT call from FFTW. Consider my matrix nrows_1 x w_size.At the moment it is stored in chunks of size w_size

for (ix = 0 ; ix < nrows_1 ; ix++)
    {
      plan =  fftw_plan_dft_1d(w_size, &source_data[ix*w_size], &transposed_data[ix*w_size],
                   FFTW_FORWARD, FFTW_ESTIMATE);
      fftw_execute(plan);
    }

So I would like to transpose the result matrix using the FFTW call.

Was it helpful?

Solution

You can use the Advanced Complex DFT of FFTW. All the dft may be performed at once thanks to parameter howmany and the output is transposed, if parameters ostride and odist are properly set.

The following code demonstrates how it works. It is compiled by gcc main.c -o main -lfftw3 -lm :

#include <fftw3.h>
#include <stdio.h>

int main(int argc, char* argv[]){

    int n=9;
    int m=4;

    fftw_complex *datasource=fftw_malloc(sizeof(fftw_complex)*m*n);
    fftw_complex *dataoutput=fftw_malloc(sizeof(fftw_complex)*m*n);

    int i,j;
    for(i=0;i<n;i++){
        for(j=0;j<m;j++){
            datasource[i*m+j][0]=(i+j)%4+j*j+i*i*i;
            datasource[i*m+j][1]=0;
        }
    }

    printf("input :\n");
    for(i=0;i<n;i++){
        for(j=0;j<m;j++){
            printf("%g ",datasource[i*m+j][0]);
        }
        printf("\n");
    }



    fftw_plan plan;
    for(i=0;i<n;i++){
        plan =  fftw_plan_dft_1d(m, &datasource[i*m], &dataoutput[i*m],
                FFTW_FORWARD, FFTW_ESTIMATE);
        fftw_execute(plan);
        fftw_destroy_plan(plan);
    }

    printf("expected output, real part, not transposed :\n");
    for(i=0;i<n;i++){
        for(j=0;j<m;j++){
            printf("%g ",dataoutput[i*m+j][0]);
        }
        printf("\n");
    }

    plan=fftw_plan_many_dft(1, &m, n,
            datasource, NULL, 1, m,
            dataoutput, NULL,n, 1,
            FFTW_FORWARD, FFTW_ESTIMATE);

    fftw_execute(plan);
    fftw_destroy_plan(plan);

    printf("output, real part, already transposed :\n");
    for(j=0;j<m;j++){
        for(i=0;i<n;i++){
            printf("%g ",dataoutput[j*n+i][0]);
        }
        printf("\n");
    }

    fftw_free(datasource);
    fftw_free(dataoutput);

    return 0;
}

I guess my answer comes too late to be useful...

Licensed under: CC-BY-SA with attribution
Not affiliated with StackOverflow
scroll top