Writing a simple Discrete Fourier Transform for real inputs in C

Question

So I'm trying to write the Discrete Fourier Transform in C to work with real 32-bit float wav files. It reads in 2 frames at a time (one for each channel, but for my purposes I'm assuming they are both the same and so I use frame[0]). This code is supposed to write out the amplitude spectrum for an input file by probing it with frequencies 20,40,60,...,10000. I am using a Hanning window on the input frames. I want to avoid using complex numbers if I can. When I run this, it gives me some very strange amplitudes (most of which are extremely small, and are not associated with the correct frequencies), which makes me believe I am making a fundamental mistake in my computation. Can somebody offer some insight into what is happening here? Here is my code:

int windowSize = 2205;
int probe[500];
float hann[2205];
int j, n;
// initialize probes to 20,40,60,...,10000
for (j=0; j< len(probe); j++) {
    probe[j] = j*20 + 20;
    fprintf(f, "%d\n", probe[j]);
}
fprintf(f, "-1\n");
// setup the Hann window
for (n=0; n< len(hann); n++) {
    hann[n] = 0.5*(cos((2*M_PI*n/(float)windowSize) + M_PI))+0.5;
}

float angle = 0.0;
float w = 0.0; // windowed sample
float realSum[len(probe)]; // stores the real part of the probe[j] within a window
float imagSum[len(probe)]; // stores the imaginary part of probe[j] within window
float mag[len(probe)]; // stores the calculated amplitude of probe[j] within a window
for (j=0; j<len(probe);j++) {
    realSum[j] = 0.0;
    imagSum[j] = 0.0;
    mag[j] = 0.0;
}

n=0; //count number of samples within current window
framesread = psf_sndReadFloatFrames(ifd,frame,1);
totalread = 0;
while (framesread == 1){
    totalread++;

    // window the frame with hann value at current sample
    w = frame[0]*hann[n];

    // determine both real and imag product values at sample n for all probe freqs times the windowed signal
    for (j=0; j<len(probe);j++) {
        angle = (2.0 * M_PI * probe[j] * n) / windowSize;
        realSum[j] = realSum[j] + (w * cos(angle));
        imagSum[j] = imagSum[j] + (w * sin(angle));
    }
    n++;
    // checks to see if current window has ended
    if (totalread % windowSize == 0) {
        fprintf(f, "B(%f)\n", totalread/44100.0);
        printf("%f breakpoint written\n", totalread/44100.0);
        for (j=0; j < len(mag); j++) { // print out the amplitudes 
            realSum[j] = realSum[j]/windowSize;
            imagSum[j] = imagSum[j]/windowSize;
            mag[j] = sqrt(pow((double)realSum[j],2)+pow((double)imagSum[j],2))/windowSize;
            fprintf(f, "%d\t%f\n", probe[j], mag[j]);
            realSum[j] = 0.0;
            imagSum[j] = 0.0;
        }
        n=0;
    }
    framesread = psf_sndReadFloatFrames(ifd,frame,1);
}

Answer 1

I think the error is in the calculation of the angle. The increment of the angle for each sample is dependent on the sampling frequency. Something like this (you seem to have 44100Hz):

angle = (2.0 * M_PI * probe[j] * n) / 44100;

Your sample window will contain one full cycle for your lowest probed frequency 20Hz. If you loop n up to 2205 then that angle would be 2*M_PI. What you saw was probably aliasing because your reference had the frequency 2205Hz and all frequencies above 1102Hz was aliased to lower frequencies.