Discrete Fourier Transform yielding incorrect results

Question

I am trying to do discrete Fourier transforms in C.

First, just brute force. First, I had a program that opened a data (amplitudes) file and put the data in an array (only one, since I limit myself to material input).

But the transformation didn’t look like that, so instead I tried to create a simple wave function and check if it is being converted correctly.

Here is my code devoid of 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, terminates, but instead of several sharp peaks in the power (frequency) graph, I get the following: .

The same frequency or different frequencies give exactly the same inverted bath curve.

I checked several sources about DFT and I still don’t know what is going wrong, apparently there are no egregious errors in the function:

 dftreal 

myself. I would like to ask for help on what might cause the problem. I am using the MinGW compiler for Windows 7. Thank you!

Answer 1

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

The problem, so to speak, is that the test waveform you are using includes frequency components that have very low frequencies relative to your linecount sample linecount . Accordingly, the DFT shows that energy is concentrated in the first few bunkers with some spectral leakage to higher bunkers, since the frequency components of the waveform are not exact multiples of the sampling frequency.

If you increase the frequencies of the test signals by setting the frequency relative to the sampling frequency, 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 like:

Answer 2

Caution: I'm a little rusty on this

As I recall, the cos / real part gives the frequency spectrum, and the sin / imag fraction gives the power spectrum. So, I think what you want to print is the frequency spectrum (which is just outreal[i] ), and not what you did (that is, Incorrect addition of outreal and outimag ). You can paint both colors, but I would just start.

In addition, I would start with easier data entry.

I modified theoretical to make one cosine wave [as well as your original]. This is predictable, and you should get one surge that you make, so I think dftreal correct.

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

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

In addition, for testing purposes, you can get much fewer samples (for example, 1000). It can also help spread your story horizontally to see things better.

Anyway, here is the code [please pardon the free style cleaning]:

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

Answer 3

You note that "... the transformation did not look like that ..."

Have you compared the results with another program or software package to confirm that the results are really wrong? I recently wrote DFT in C ++ and JavaScript and compared the results with the results from MATLAB, Maple and MathCAD. Sometimes the only difference is scaling or formatting.

How big was the original amplitude scale you wanted to introduce? If you post a sample of data here, I am ready to run it through my own program and see if my program produces the same results as yours.