Discrete Fourier transform giving incorrect results

Question

I am trying to do discrete Fourier transforms in C.

Initially just the brute-force method. First I had the program open a data file (of amplitudes) and put the data into an array (just one, since I'm limiting myself to real-valued input).

But the transform looked wrong, so instead I tried to generate a simple wavefunction and check whether it transforms properly.

Here's my code, stripped of the bells and whistles:

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#define M_PI 3.14159265358979323846

//the test wavefunction
double theoretical(double t)
{
    double a = sin(M_PI * t) + 2 * sin(2 * M_PI * t) + 4 * sin(4 * M_PI * t); 
    return a;
}
//-------------------------------------------------------------------------
void dftreal(double inreal[], double outreal[], double outimag[], int linecount) 
{
    int n, k;
    for (k = 0; k < linecount; k++)
    {
        double sumreal = 0;
        double sumimag = 0;

        for (n = 0; n < linecount; n++) 
        {
            double angle = 2 * M_PI * n * ( k / (double) linecount);
            sumreal +=  inreal[n] * cos(angle);
            sumimag +=  inreal[n] * sin(angle);
        }
        outreal[k] = sumreal;
        outimag[k] = sumimag;
    }
}
//=========================================================================
int main(void)
{
     int linecount = 44100;
     //creates all necessary arrays
     double inreal[linecount], outreal[linecount], outimag[linecount], p[linecount]; 

     FILE *fout = fopen("Output.txt", "w");

     for (int i = 0 ; i < linecount ; ++i)
     {
         inreal[i] = theoretical( i / (double) linecount);
     }

     //actually computes the transform
     dftreal(inreal, outreal, outimag, linecount); 

     for (int i = 0 ; i < linecount ; ++i)
     {
          p[i] = 2*(outreal[i] * outreal[i] + outimag[i] * outimag[i]);
          fprintf(fout, "%f %f \n", (i / (double) linecount), p[i]);
     }
     fclose(fout);
     printf("\nEnd of program");
     getchar();
     return 0;
}

The program compiles, completes, but instead of several sharp peaks on a power(frequency) plot, I get this: .

A single frequency or different frequencies give the exact same inverted-bathtub-curve.

I checked several sources about the DFT and I still don't know what is going wrong, there do not seem to be any glaring errors with the function:

dftreal

itself. I'd like to ask for help on what could be causing the issue. I'm using the MinGW compiler on Windows 7. Thanks!

Theodor

The DFT requieres high decimal precision. Have you tried changing your fprintf statements so they can include more decimals? Like %.10f or more... – Mechelle

Answer 1

The good news is that there is nothing wrong with your dftreal implementation.

The problem, if one could call it that, is that the test waveform you are using includes frequency components which are of very low frequencies relative to your sampling rate linecount. Correspondingly, the DFT shows that the energy is concentrated in the first few bins, with some spectral leakage into higher bins since the waveform frequency components aren't exact multiples of the sampling rate.

If you increase the test waveform frequencies by making the frequency relative to the sampling frequency with for example:

//the test wavefunction
double theoretical(double t)
{
  double f = 0.1*44100;
  double a = sin(2 * M_PI * f * t) + 2 * sin(4 * M_PI * f * t) + 4 * sin(8 * M_PI * f * t); 
  return a;
}

You should get a plot such as:

Retral

You are, of course, right. I don't know how I managed to mangle the math so badly, the program works as intended now. Thank you very much! – Theodor

Answer 2

Caveat: I'm a bit rusty on this

As I remember it, the cos/real part yields the frequency spectrum and the sin/imag part yields the power spectrum. So, what I think you want to print is the frequency spectrum (which is just outreal[i]) rather than what you did (i.e. adding outreal and outimag is wrong). You could plot both in different colors but I'd start simply.

Also, I'd start with simpler data input.

I changed the theoretical to do a single cosine wave [as well as your original]. This is predictable and you should get a single spike, which you do, so I think dftreal is correct.

I also added some command line options so you can experiment.

I also found that dividing i / linecount got truncated sometimes, so I created a FRAG macro to illustrate/correct this. Without this, the spike for an input of cos(2 * pi * t) ended up in outreal[0] (DC part?) instead of outreal[1]

Also, for testing purposes, you can get by with far fewer samples (e.g. 1000). This may also help spread out your plot horizontally to be better able to see things.

Anyway, here's the code [please pardon the gratuitous style cleanup]:

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>

//#define M_PI 3.14159265358979323846
#define M_2PI (M_PI * 2.0)

int opt_A;                              // algorithm to use
int linecount;                          // number of samples
int opt_a;                              // use absolute value on output
int opt_s;                              // use squared value on output
int opt_i;                              // use integer output index
int opt_j;                              // output imaginary part
int opt_c;                              // .csv compatibility

time_t todlast;

// the first was your original and would truncate
#if 0
#define FRAG(_i)    ((_i) / linecount)
#else
#define FRAG(_i)    ((double) (_i) / linecount)
#endif

void
pgr(int k,int n,int count)
{
    time_t todnow = time(NULL);

    if ((todnow - todlast) >= 1) {
        todlast = todnow;
        printf("\r%d %d ",count,k);
        fflush(stdout);
    }
}

// the test wavefunction
double
theoretical(double t)
{
    double a;

    switch (opt_A) {
    case 0:
        a = 0.0;
        a += sin(M_PI * t);
        a += 2.0 * sin(2.0 * M_PI * t);
        a += 4.0 * sin(4.0 * M_PI * t);
        break;
    default:
        a = cos(M_2PI * t * opt_A);
        break;
    }

    return a;
}

//-------------------------------------------------------------------------
void
dftreal(double inreal[], double outreal[], double outimag[], int linecount)
{
    int n;
    int k;

    for (k = 0; k < linecount; k++) {
        double sumreal = 0;
        double sumimag = 0;
        double omega_k = (M_2PI * k) / (double) linecount;

        for (n = 0; n < linecount; n++) {
            // omega k:
#if 0
            double angle = M_2PI * n * FRAG(k);
#else
            double angle = omega_k * n;
#endif

            sumreal += inreal[n] * cos(angle);
            sumimag += inreal[n] * sin(angle);

            pgr(k,n,linecount);
        }

        outreal[k] = sumreal;
        outimag[k] = sumimag;
    }
}

//=========================================================================
int
main(int argc,char **argv)
{
    char *cp;

    --argc;
    ++argv;

    for (;  argc > 0;  --argc, ++argv) {
        cp = *argv;
        if (*cp != '-')
            break;

        switch (cp[1]) {
        case 'a':  // absolute value
            opt_a = ! opt_a;
            break;
        case 's':  // square the output value
            opt_s = ! opt_s;
            break;
        case 'c':  // .csv output
            opt_c = ! opt_c;
            break;
        case 'i':  // integer output
            opt_i = ! opt_i;
            break;
        case 'j':  // imaginary output
            opt_j = ! opt_j;
            break;
        case 'A':
            opt_A = atoi(cp + 2);
            break;
        case 'N':
            linecount = atoi(cp + 2);
            break;
        }
    }

    if (linecount <= 0)
        linecount = 44100 / 2;

    // creates all necessary arrays
    double inreal[linecount];
    double outreal[linecount];
    double outimag[linecount];
#if 0
    double p[linecount];
#endif

    FILE *fout = fopen("Output.txt", "w");

    for (int i = 0; i < linecount; ++i)
        inreal[i] = theoretical(FRAG(i));

    todlast = time(NULL);

    // actually computes the transform
    dftreal(inreal, outreal, outimag, linecount);

    for (int i = 0; i < linecount; ++i) {
#if 0
        p[i] = 2 * (outreal[i] * outreal[i] + outimag[i] * outimag[i]);
        fprintf(fout, "%f %f\n", (i / (double) linecount), p[i]);
#endif

#if 0
        fprintf(fout, "%f %f %f\n", (i / (double) linecount),
            outreal[i] * outreal[i],
            outimag[i] * outimag[i]);
#endif

#if 0
        fprintf(fout, "%f %f\n",
            (i / (double) linecount),outreal[i] * outreal[i]);
#endif

        if (opt_a) {
            outreal[i] = fabs(outreal[i]);
            outimag[i] = fabs(outimag[i]);
        }

        if (opt_s) {
            outreal[i] *= outreal[i];
            outimag[i] *= outimag[i];
        }

        if (! opt_c) {
            if (opt_i)
                fprintf(fout, "%d ",i);
            else
                fprintf(fout, "%f ",FRAG(i));
        }

        fprintf(fout, "%f",outreal[i]);
        if (opt_j)
            fprintf(fout, "%f",outimag[i]);
        fprintf(fout, "\n");
    }
    fclose(fout);

    printf("\nEnd of program\n");

    //getchar();
    return 0;
}

Mears

Answer 3

You mention that "... the transform looked wrong, ..."

Did you compare the results with another program or software package to confirm that the results are, indeed, wrong? I recently wrote a DFT in C++ and JavaScript and compared the outputs to results from MATLAB, Maple, and MathCAD. Sometimes the difference is just a matter of scaling or formatting.

How big was the original array of amplitudes you wanted to input? If you post sample data here, I am willing to run it through my own program and see if my program outputs the same results as yours.

Lamentable