FFTW simply will not return values other than infinity, values that approach zero, or negative infinity

Question

I am trying to calculate the DFT of a sine wave at 20hz.

First I fill a std::vector with 10 cycles of the sine function at 20 hz:

std::vector<double> sinx;
double samplerate = 1000.0;
double frequency = 20.0;
double num_cycles = 10.0;
for (int i=0; i<samplerate/frequency*num_cycles; i++){
    sinx.push_back(sin(2.0*M_PI*((double)i)*frequency/samplerate));
}
double N = sinx.size();

Then I initiate the FFTW input and output arrays as well as the FFTW Plan:

fftw_complex *in, *out;
fftw_plan p;

in = (fftw_complex*) fftw_malloc(sizeof(fftw_complex)* N);
out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);
p = fftw_plan_dft_1d(N, in, out, FFTW_FORWARD, FFTW_ESTIMATE);

Then I fill the input array with the sine wave function values:

for (int i=0; i<N; i++){
    in[i][0] = sinx[i];
}

Then I calcuate the FFT:

fftw_execute(p);

After calculating the FFT, I check the results and clean up.

    double hz_per_index = samplerate/2.0/N*num_cycles;
for (int i=0; i<N/num_cycles; i++){
    std::cout<<"hz: "<<hz_per_index*i<<" out: "<<fabs(out[i][0])<<" "<<out[i][1]<<" in: "<<in[i][0]<<std::endl;
}

I would expect a spike around 20hz, however there is no spike at all. All the values that come back are always either approacting zero, infinity, or negative infinity. My first thought was that the input was overwriting itself, but check the input after calculating and it is correct. I've tried running FFTW in different modes, FFTW_FORWARD, FFTW_BACKWARD, FFTW_ESTIMATE, FFTW_RELATIVE, nothing is helping. I've tried calculating different sine functions at different frequencies, sample rates, number of cycles.

The sad part is that I was able to get this to work once. Then I continued working, saw that it wasn't working, then re-opened the file that worked, and now it doesn't work anymore!

Anyone else come across this?

EDIT: NEW CODE AS PER SUGGESTIONS

It's still not working. I tried initializing the imaginary portion of the complex numbers to 0, and I still had the same problem, so then here's code utilizing a real in/complex out function of the fftw library.

I'm basically getting 0 for everything. Here's the full output of all 5 cycles of the input that shows both the input, output, and hz value:

https://gist.github.com/anonymous/ed419f4cfb43887c810e

double *in;
fftw_complex* out;
fftw_plan p;

std::vector<double> sinx;
double samplerate = 1000.0;
double frequency = 20.0;
double num_cycles = 5.0;
for (int i=0; i<samplerate/frequency*num_cycles; i++){
    sinx.push_back(sin(2.0*M_PI*((double)i)*frequency/samplerate));
}
int N = sinx.size();

in = (double*) fftw_malloc(sizeof(double)* N);
out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);
p = fftw_plan_dft_r2c_1d(N, in, out, FFTW_ESTIMATE);

for (int i=0; i<N; i++){
    in[i] = sinx[i];
}

fftw_execute(p);

double hz_per_index = samplerate/2.0/N*num_cycles;
for (int i=0; i<N; i++){
    std::cout<<"hz: "<<hz_per_index * (i % (int)(N/num_cycles))<<" out: "<<out[i][0]<<" in: "<<in[i]<<std::endl;
}

fftw_destroy_plan(p);
fftw_free(in); fftw_free(out);

return 0;

}

Answer 1

You're allocating an array of fftw_complexes, which consist of two elements — the real and complex components — but you're only initializing the real component of each sample. This is probably leaving the complex components containing random data, causing unexpected crazy results!

If you don't need to deal with complex samples — which is likely — you may want to use one of FFTW's real-data DFT functions, which takes a double array as input. Otherwise, initialize all the complex components (in[i][1]) to zero.

Additionally, you're not looking at the complex component of the output. There may be something significant there you're missing out on.