Schreiben einer einfachen diskreten Fourier-Transformation für reale Eingaben in C

Question

Also versuche ich, die diskrete Fourier-Transformation in C zu schreiben, um mit echten 32-Bit-Float-WAV-Dateien zu arbeiten. Es liest 2 Frames gleichzeitig ein (eines für jeden Kanal, aber für meine Zwecke nehme ich an, dass sie beide gleich sind und so verwende ich Frame [0]). Dieser Code soll das Amplitudenspektrum für eine Eingabedatei durch Sondieren mit den Frequenzen 20, 40, 60, ..., 10000 ausschreiben. Ich verwende ein Hanning-Fenster für die Eingabefelder. Ich möchte vermeiden, komplexe Zahlen zu verwenden, wenn ich kann. Wenn ich das mache, gibt es mir sehr seltsame Amplituden (von denen die meisten extrem klein sind und nicht mit den richtigen Frequenzen verbunden sind), was mich glauben lässt, dass ich einen fundamentalen Fehler in meiner Berechnung mache. Kann jemand einen Einblick geben, was hier passiert? Hier ist mein 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); 
} 

Answer 1

Mit Code unten - nur leicht reorganisiert zu kompilieren und eine gefälschte Probe zu erstellen, das tue ich nicht Nullen erhalten.

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

zu

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

Dies macht es nur einfacher, das Nicht-Null-Daten zu sehen: Ich habe den Ausgang Aufruf am Ende von geändert. Vielleicht ist das einzige Problem, den Skalenfaktor zu verstehen? Beachten Sie, wie ich Eingaben eingegeben habe, um einen reinen Ton als Testfall zu erzeugen.

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

Answer 2

Ich denke, der Fehler ist in der Berechnung des Winkels. Das Inkrement des Winkels für jede Probe ist abhängig von der Abtastfrequenz. So etwas (Sie scheinen 44100 Hz zu haben):

angle = (2.0 * M_PI * probe[j] * n)/44100; 

Ihr Probenfenster einen vollen Zyklus für die niedrigste sondierten Frequenz 20Hz enthalten wird. Wenn Sie eine Schleife bis zu 2205 durchführen, wäre dieser Winkel 2 * M_PI. Was Sie gesehen haben, war vermutlich Aliasing, weil Ihre Referenz die Frequenz 2205Hz hatte und alle Frequenzen oberhalb von 1102Hz mit niedrigeren Frequenzen frequentiert wurden.