Trasformata di Fourier discreta che fornisce risultati errati

Question

Sto provando a fare trasformazioni discrete di Fourier in C.

Inizialmente solo il metodo forza bruta. Per prima cosa ho avuto il programma aprire un file di dati (di ampiezze) e mettere i dati in un array (solo uno, dato che mi sto limitando a input di valore reale).

Ma la trasformazione sembrava errata, così invece ho cercato di generare una semplice funzione d'onda e controllare se si trasforma correttamente.

Ecco il mio codice, spogliato delle campane e fischietti:

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

Il programma viene compilato, completa, ma invece di diversi picchi taglienti su un terreno di alimentazione (frequenza), ottengo questo: .

Una singola frequenza o frequenze diverse danno la stessa identica curva della vasca da bagno invertita.

ho controllato diverse fonti circa la DFT e io ancora non so cosa non va, non sembrano essere eventuali errori evidenti con la funzione:

dftreal 

stessa. Mi piacerebbe chiedere aiuto su ciò che potrebbe causare il problema. Sto usando il compilatore MinGW su Windows 7. Grazie!

Answer 1

La buona notizia è che non c'è niente di sbagliato con il tuo dftreal implementazione.

Il problema, se così si può chiamare, è che la forma d'onda di prova che si sta utilizzando include componenti di frequenza che sono di frequenza molto bassa rispetto alla frequenza di campionamento linecount. Corrispondentemente, la DFT mostra che l'energia è concentrata nei primi bin, con qualche dispersione spettrale in bidoni più alti poiché le componenti della frequenza della forma d'onda non sono multipli esatti della frequenza di campionamento.

Se si aumenta le frequenze di prova della forma d'onda, rendendo la frequenza relativa alla frequenza di campionamento con ad esempio:

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

si dovrebbe ottenere una trama come ad esempio:

Answer 2

Caveat: Sono un po 'arrugginito su questo

Se ben ricordo, la cos/parte reale si ottiene lo spettro di frequenza e la parte sin/imag rese lo spettro di potenza. Quindi, quello che penso che tu voglia stampare è lo spettro di frequenza (che è solo outreal[i]) piuttosto che quello che hai fatto (cioè aggiungere outreal e outimag è sbagliato). Potresti tracciare entrambi in colori diversi, ma inizierei semplicemente.

Inoltre, inizierei con dati più semplici.

Ho modificato il theoretical per eseguire un'unica onda coseno [oltre che l'originale]. Questo è prevedibile e dovresti ottenere un singolo picco, cosa che fai, quindi penso che dftreal sia corretto.

Ho anche aggiunto alcune opzioni della riga di comando in modo da poter sperimentare.

Ho anche scoperto che la divisione i/linecount è stata troncata a volte, quindi ho creato una macro FRAG per illustrare/correggere questo. Senza questo, il picco per un ingresso di cos(2 * pi * t) finito in outreal[0] (parte DC?) Invece di outreal[1]

Inoltre, a scopo di test, è possibile ottenere da con molti meno campioni (ad esempio 1000). Questo può anche aiutare a distribuire la trama orizzontalmente per essere più in grado di vedere le cose.

Ad ogni modo, ecco il codice [si prega di perdonare la pulizia stile gratuità]:

#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

si menziona che". .. la trasformazione sembrava sbagliata, ... "

Hai confrontato i risultati con un altro programma o pacchetto software per confermare che i risultati sono, in effetti, sbagliati? Recentemente ho scritto un DFT in C++ e JavaScript e ho confrontato gli output con i risultati di MATLAB, Maple e MathCAD. A volte la differenza è solo una questione di ridimensionamento o formattazione.

Quanto era grande la matrice originale di ampiezze che si desidera immettere? Se pubblichi qui dati di esempio, sono disposto a eseguirlo attraverso il mio programma e vedere se il mio programma produce gli stessi risultati dei tuoi.