Manolete
Manolete

Reputation: 3517

Result of FFTW 1D stored in transposed way

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.

Upvotes: 2

Views: 166

Answers (1)

francis
francis

Reputation: 9817

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...

Upvotes: 2

Related Questions