CODEWITHSUNDEEP
cfftw
TimeFall
TimeFall

Reputation: 97

In-place FFT with FFTW3

I'm trying to understand in-place vs. out-of-place FFTs with FFTW3 in c. For starters, I decided to use a one dimensional transform. I did the out-of-place real to complex transform first and then I used the complex output from that transform to do an in-place complex to real transform, hoping to get the original real input data back. I do not, however, and I was just wondering why? My sample code to do this is below.

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

void one_dimensional_transform(void)
{
    float *in;
    fftwf_complex *out, *inplace;
    int N = 8;
    fftwf_plan plan;
    int i;
    FILE *fd;

    // Do out of place transform. For this, in has N elements and out has N/2 + 1.
    in  = fftwf_alloc_real(N);
    out = fftwf_alloc_complex(N/2 + 1); 
    inplace = fftwf_alloc_complex(2 * ((N/2) + 1));

    plan = fftwf_plan_dft_r2c_1d(N, in, out, FFTW_ESTIMATE);

    fd = fopen("input_data.txt", "w");

    for(i = 0; i < N; i++)
    {   
        in[i] = sin(2.0 * M_PI * (2.0 * M_PI * i / N));
        fprintf(fd, "%f\n", in[i]);
    }   
    fclose(fd);

    fftwf_execute(plan);
    fftwf_destroy_plan(plan);
    fftwf_cleanup();

    // Do in-place transform
    for(i = 0; i < N/2 + 1; i++)
    {   
        inplace[i] = out[i];
    }   

    plan = fftwf_plan_dft_c2r_1d(N, (fftwf_complex *)inplace, (float *)inplace, FFTW_ESTIMATE);

    // I think that I need to include the complex conjugates to the complex version of
    // inplace (as this is what Mesinger does), and that's why we have the size of
    // inplace being 2 * (N/2 + 1). 
    for(i = 1; i < N/2; i++)
    {   
        inplace[N/2 + i] = conjf(inplace[N/2 - i]); 
    }   

    for(i = 0; i < 2 * (N/2 + 1); i++)
    {   
        inplace[i] /= (float)N;
        fftwf_execute(plan);
        fftwf_destroy_plan(plan);
        fftwf_cleanup();

        fd = fopen("inplace_data.txt", "w");

        for(i = 0; i < N; i++)
        {
            fprintf(fd, "%f\n", crealf(inplace[i]));
        }
        fclose(fd);

        fftwf_free(in);
        fftwf_free(out);
        fftwf_free(inplace);
}

From the FFTW3 docs, they say that the real data has to be of size 2(N/2 + 1). Also, since the real to complex transform cuts off the second half of the complex array due to Hermitian symmetry, I decided to explicitly put those elements back in, but I'm not sure if that's needed. I also rescaled the output of the in-place transform by N, because the docs said that the transforms are not normalized.

Upvotes: 2

Views: 1599

Answers (1)

SleuthEye
SleuthEye

Reputation: 14577

The main issue is that the posted code actually runs fftwf_execute (and fftw_destroy!) multiple times due to the loop

for(i = 0; i < 2 * (N/2 + 1); i++)
{   
    inplace[i] /= (float)N;
    fftwf_execute(plan);
    fftwf_destroy_plan(plan);
    fftwf_cleanup();
    ...
}

You probably meant to execute the following:

for(i = 0; i < 2 * (N/2 + 1); i++)
{   
    inplace[i] /= (float)N;
}
fftwf_execute(plan);
fftwf_destroy_plan(plan);
fftwf_cleanup();

Also, the result is stored as it were in a float array, so you should read back and save the data with:

  for(i = 0; i < N; i++)
  {
    fprintf(fd, "%f\n", ((float*) inplace)[i]));
  }

instead of only extracting the real parts with crealf(inplace[i]) (which would take only every second value).

Finally with respect to your comment

Also, since the real to complex transform cuts off the second half of the complex array due to Hermitian symmetry, I decided to explicitly put those elements back in, but I'm not sure if that's needed.

I can assure you that it's not needed. The complex to real transform doesn't need the second half of the spectrum since it is based on the knowledge that it has Hermitian symmetry.

Upvotes: 2

Related Questions