FFTW of double samples saved in file

Question

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;

}

Answer 1

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.)

Answer 2

Take a good long look at what is wrong with the following code:

#include <stdio.h>

void prn (float *a, int n);

int main (void)
{
    float floatarray[] = { 0.0, 5.0, 10.0, 15.0, 20.0, 25.0 };

    prn (floatarray, sizeof floatarray/sizeof *floatarray);

    return 0;
}

void prn (float *a, int n)
{
    int i;
    double *d = (double *)a;

    for (i = 0; i < n; i++)
        printf (" a[%d] : %01.0lf\n", i, d[i]);
}

Run it, what happens? Why?

What is the sizeof a float? What is the sizeof a double? What happens when the cast double *d = (double *)a; takes place? How many bytes are expected in each element of d? How many are there? Your

double* doubleBuffer = castToDouble(buffer, bufferLen);

causes the exact same problem.

Element-by-Element Cast/Assignment Required

Whether you use a new double array allocated with malloc (or calloc), you must do an element-by-element cast/assignment to complete the cast of the array. Here is the same code with a few additions using a variable length array to avoid allocating memory dynamically:

#include <stdio.h>

void castfloatdouble (double *d, float *f, int n);
void prn (float *a, int n);
void prndouble (double *d, int n);

int main (void)
{
    float floatarray[] = { 0.0, 5.0, 10.0, 15.0, 20.0, 25.0 };
    int n = (int)(sizeof floatarray/sizeof *floatarray);
    double doublearray[n];  /* memset is a good idea to zero a VLA */

    prn (floatarray, n);  /* incorrect results due to typesize mismatch */
    putchar ('\n');

    castfloatdouble (doublearray, floatarray, n);

    prndouble (doublearray, n);  /* correct results after deep cast/assign */
    putchar ('\n');

    return 0;
}

void castfloatdouble (double *d, float *f, int n)
{
    int i;

    for (i = 0; i < n; i++)
        d[i] = (double)f[i];
}

void prn (float *a, int n)
{
    int i;
    double *d = (double *)a;

    for (i = 0; i < n; i++)
        printf (" a[%d] : %01.0lf\n", i, d[i]);
}

void prndouble (double *d, int n)
{
    int i;

    for (i = 0; i < n; i++)
        printf (" a[%d] : %01.0lf\n", i, d[i]);
}

Output Showing Incorrect & Correct Results

$ ./bin/castfloatdouble2
 a[0] : 2048
 a[1] : 16777220
 a[2] : 805306499
 a[3] : 0
 a[4] : 0
 a[5] : 0

 a[0] : 0
 a[1] : 5
 a[2] : 10
 a[3] : 15
 a[4] : 20
 a[5] : 25