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;
}
N
an integral type. – NicholasM