I am trying to calculate FFT on 53k double samples using FFTW library and on it's basis guess what's the fundamental frequency of the signal. Samples are generated by sndfile library on the basis of wav input file (program loads in wav file, generates samples of double data and saves them to a text file). Each sample's order of magnitude is ranging from -0.0009 to +0.0009. After calculating FFT using function presented below, I receive in[0][0] = -2.142. Either this is an incorrect method of retrieving data or the input file is incorrect. What am I doing wrong? Is it the wrong data passing to FFTW function, is the data stored in file incorrect or the function I have written is incorrect? Buffer is a float array containing samples.
static fftw_complex* calculateFourier(float* buffer, const unsigned int bufferLen)
{
unsigned int i;
//const unsigned short windowSize = 1024;
fftw_plan result;
fftw_complex *out, *in;
in = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * bufferLen);
out = (fftw_complex*) fftw_malloc( sizeof(fftw_complex) * bufferLen);
double* doubleBuffer = castToDouble(buffer, bufferLen);
for(i = 0; i < bufferLen; i++)
{
in[i][0] = doubleBuffer[i];
in[i][1] = 0.000000000000;
}
result = fftw_plan_dft_1d(bufferLen, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
//result = fftw_plan_dft_r2c_1d(fileSize, castToDouble(buffer, bufferLen), out, FFTW_MEASURE);
fftw_execute(result);
fftw_destroy_plan(result);
return out;
}
static double* castToDouble(float* buffer, const unsigned int bufferLen)
{
unsigned int i;
double* doubleBuffer = (double*)malloc( sizeof(double) * bufferLen);
for(i = 0; i < bufferLen; i++)
doubleBuffer[i] = (double)buffer[i];
return doubleBuffer;
}
I think the problem is that you haven't applied a window function to your raw data. As a result, the Fourier transform has a strong frequency component that corresponds to the length of your sample buffer.
Try applying a Hanning window by replacing this line:
in[i][0] = doubleBuffer[i];
with this:
in[i][0] = doubleBuffer[i] * 0.5 * (1.0 - cos(2.0 * M_PI * i / (bufferLen - 1)));
(Obviously this window function can be cached if you need to use it repeatedly.)