Diskrete Fourier-Transformation mit falschen Ergebnissen

Question

Ich versuche, diskrete Fourier-Transformationen in C durchzuführen.

Zunächst nur die Brute-Force-Methode. Zuerst ließ ich das Programm eine Datendatei (von Amplituden) öffnen und legte die Daten in ein Array (nur eins, da ich mich auf reellwertige Eingaben beschränke).

Aber die Transformation sah falsch aus, also versuchte ich stattdessen, eine einfache Wellenfunktion zu erzeugen und zu überprüfen, ob sie sich richtig transformiert.

Hier ist mein Code, beraubt Drum und Dran:

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

Das Programm erstellt, abgeschlossen ist, sondern von mehreren scharfen Spitzen auf einer Leistung (Frequenz) Plot, erhalte ich diese: .

Eine einzelne Frequenz oder verschiedene Frequenzen ergeben genau die umgekehrte Badewannenkurve.

Ich habe mehrere Quellen über die DFT und ich weiß immer noch nicht, was falsch läuft, es nicht mit der Funktion keine eklatanten Fehler zu sein scheint:

dftreal 

selbst. Ich möchte um Hilfe bitten, was das Problem verursachen könnte. Ich verwende den MinGW-Compiler unter Windows 7. Danke!

Answer 1

Die gute Nachricht ist, dass es nichts falsch mit Ihrer dftreal Implementierung ist.

Das Problem, wenn man es so nennen könnte, ist, dass die von Ihnen verwendete Testwellenform Frequenzkomponenten enthält, die relativ zu Ihrer Abtastrate sehr niedrige Frequenzen haben . Entsprechend zeigt die DFT, dass die Energie in den ersten paar Bins konzentriert ist, mit einem gewissen spektralen Verlust in höhere Bins, da die Wellenformfrequenzkomponenten keine exakten Vielfachen der Abtastrate sind.

Wenn Sie die Prüfwellenform Frequenzen erhöhen, indem die Frequenz in Bezug auf die Abtastfrequenz mit zum Beispiel:

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

Sie ein Grundstück erhalten sollen wie:

Answer 2

Caveat: Ich bin ein bisschen rostig auf dieser

Wie ich mich erinnere, die cos/Realteil liefert das Frequenzspektrum und die sin/imag Teil des Leistungsspektrums ergibt. Also, was ich denke, dass Sie drucken möchten, ist das Frequenzspektrum (das ist nur outreal[i]), anstatt was Sie getan haben (d. H. outreal und outimag ist falsch). Sie könnten beide in verschiedenen Farben zeichnen, aber ich würde einfach anfangen.

Auch würde ich mit einfacheren Dateneingabe beginnen.

Ich habe die theoretical geändert, um eine einzelne Kosinuswelle [sowie Ihr Original] zu machen. Dies ist vorhersehbar und Sie sollten einen einzigen Spike bekommen, was Sie tun, also denke ich, dftreal ist richtig.

Ich habe auch einige Befehlszeilenoptionen hinzugefügt, so dass Sie experimentieren können.

Ich fand auch, dass die Teilung i/linecount wurde manchmal abgeschnitten, so dass ich ein FRAG Makro erstellt, um dies zu veranschaulichen/zu korrigieren. Ohne diesen endete die Spitze für eine Eingabe von cos(2 * pi * t) bis in outreal[0] (DC Teil?) Statt outreal[1]

Auch für Testzwecke können Sie mit deutlich weniger Proben erhalten, indem (beispielsweise 1000). Dies kann auch helfen, Ihre Handlung horizontal zu verteilen, um besser in der Lage zu sein, Dinge zu sehen.

Wie auch immer, hier ist der Code [der unentgeltliche Stil Bereinigungs bitte verzeihen]:

#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

Sie erwähnen, dass“. .. die Transformation sah falsch aus, ... "

Haben Sie die Ergebnisse mit einem anderen Programm oder Softwarepaket verglichen, um zu bestätigen, dass die Ergebnisse tatsächlich falsch sind? Ich schrieb vor kurzem eine DFT in C++ und JavaScript und verglich die Ausgaben mit Ergebnissen von MATLAB, Maple und MathCAD. Manchmal ist der Unterschied nur eine Frage der Skalierung oder Formatierung.

Wie groß war die ursprüngliche Amplitude, die Sie eingeben wollten? Wenn Sie Beispieldaten hier posten, bin ich bereit, es durch mein eigenes Programm zu führen und zu sehen, ob mein Programm die gleichen Ergebnisse wie Ihre ausgibt.