3
votes

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;

}

2
Your output function is only outputting the first 50 of the 500 results . It seems like you were attempting to output every 5th result or something.M.M
Please make that N an integral type.NicholasM
I was only outputting the first complete cycle of the input. Now the code outputs all the cycles and I created a gist of the output showing the input, output, and corresponding hz value of the cycle.sugarmuff

2 Answers

6
votes

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.

2
votes

Your input are 10 periods of a pure real sine wave. This means:

  • Your output is pure imaginary caused by the odd sine function (sin(x) = -sin(-x)) In the case of a cosine function the output were pure real because cosines is an even function (cos(x) = cos(-x)
  • There are only two frequency lines, one in the positive frequency band and one in the negative band.
  • Because you have 10 periods in your input, your lines should occur at out[10][1] and out[N-10][1].
  • Because you didn't normalize the output and the number of samples is 500, the amplitudes should be out[10][1] = -250 and out[490][1] = +250.
  • All other values should be zero.