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);
}

Threecornered

Don't see an obvious error, but I'd suggest generating test cases and checking that the mathematical properties of the coefficients are sound - E.g. Real value inputs mean symmetric coefficients. – Kling

@Kling i'm sorry, i'm not sure exactly what that means. what are the coefficients in this case? would it be the variable w? And I tried to run this on an A4 at 440 hz, and the DFT returned pretty much 0.00000 for every single frequency for the entire duration – Threecornered

If everything is zero, you have a bigger problem. See answer. – Kling

If the goal is doing the transforms rather than learning to write one, consider using OpenCV. Yes it's a battleship to swap a fly but would do what you need. – Pileum

There is also FFTW library for doing Fourier Transform. – Phylys

@keith I did not go through the OP code, but if it is a "classic" discrete transform of discrete function, you will not calculate all coefficients (because you know them to be simmetrical) nor get simmetry, rather circular simmetry, which can be checked, of course, but it's not as trivial as a plain simmetry check. – Phylys

@Threecornered why avoid complex numbers? When you understand how the complex exponent works, coding is quite simple. Honestly, not.trivial, but still straightforward. – Phylys

Some suggestions not related to your bug. On a modern CPU, use double instead of float if memory is not an issue, since double is much much faster than float, and don't use pow(a,2) but instead use (a*a). pow uses calls to the very expensive exponential function, which is unnecessary if you are just doing squaring. – Expedition

Can you post an image of your input data and output spectrum? It is possible that you do have low frequency components that you don't know about. Even an FFT of what you think is just a pure sine wave will have lots of peaks. The spectrum is not as clean as you imagine. – Expedition

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.

Thetic

Answer 2

With code below - only slightly reorganised to compile and create a fake sample, I do not get all zeroes. I have changed the output call to at the end from:

fprintf(f, "%d\t%f\n", probe[j], mag[j] );

to

if (mag[j] > 1e-7)
    fprintf(f, "%d\t%f\n", probe[j], mag[j] * 10000);

This just makes it easier to see the non-zero data. Maybe the only issue is understanding the scale factor? Note how I faked input to generate a pure tone as a test case.

#include <math.h>

#include <stdio.h>

#define M_PI 3.1415926535

#define SAMPLE_RATE 44100.0f

#define len(array) (sizeof array/sizeof *array)


unsigned psf_sndReadFloatFrames(FILE* inFile,float* frame,int framesToRead)
{
    static float counter = 0;   
    float frequency = 1000;
    float time = counter++;
    float phase = time/SAMPLE_RATE*frequency;
    *frame = (float)sin(phase);
    return counter < SAMPLE_RATE;
}

void discreteFourier(FILE* f)                    
{
    FILE* ifd = 0;
    float frame[1];
    int windowSize = 2205;
    int probe[500];
    float hann[2205];


    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

    int j, n;

    unsigned framesread = 0;
    unsigned totalread = 0;

    for (j=0; j<len(probe);j++) {
        realSum[j] = 0.0;
        imagSum[j] = 0.0;
        mag[j] = 0.0;
    }

    // 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;
    }
    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/SAMPLE_RATE);
            printf("%f breakpoint written\n", totalread/SAMPLE_RATE);
            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;
                if (mag[j] > 1e-7)
                    fprintf(f, "%d\t%f\n", probe[j], mag[j] * 10000);
                realSum[j] = 0.0;
                imagSum[j] = 0.0;
            }
            n=0;
        }
        framesread = psf_sndReadFloatFrames(ifd,frame,1);
    }
}

Kling