Pisanie prostej, dyskretnej transformaty Fouriera dla rzeczywistych danych wejściowych w C

Question

Próbuję więc zapisać dyskretną transformatę Fouriera w C, aby pracować z prawdziwymi 32-bitowymi plikami WAV typu float. Czyta się w 2 klatkach na raz (po jednym dla każdego kanału, ale dla moich celów zakładam, że oba są takie same i dlatego używam ramki [0]). Ten kod ma wypisać widmo amplitudy dla pliku wejściowego poprzez sondowanie go częstotliwościami 20,40,60, ..., 10000. Używam okna Hanning na wejściowych klatkach. Chcę, jeśli mogę, unikać stosowania liczb zespolonych. Kiedy to robię, daje mi to bardzo dziwne amplitudy (z których większość jest bardzo mała i nie są związane z prawidłowymi częstotliwościami), co powoduje, że wierzę, że popełniam podstawowy błąd w moich obliczeniach. Czy ktoś może dać wgląd w to, co się tutaj dzieje? Tu jest mój kodu:

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

Z kodu poniżej - tylko nieznacznie zreorganizowana do kompilacji i tworzyć fałszywą próbkę, robię nie uzyskać wszystkie zera. Zmieniłem wezwanie do wyjścia na końcu od:

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

do

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

To właśnie sprawia, że ​​łatwiej zobaczyć dane niezerową. Może jedynym problemem jest zrozumienie współczynnika skali? Zwróć uwagę, jak sfałszowałem dane wejściowe, aby wygenerować czysty ton jako przypadek testowy.

#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

Myślę, że błąd jest w obliczaniu kąta. Wzrost kąta dla każdej próbki zależy od częstotliwości próbkowania. Coś takiego (to wydaje się mieć 44100Hz):

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

Twoje okno próbka będzie zawierać jeden pełny cykl dla najniższych częstotliwości 20Hz sondowane. Jeśli zapętlisz n do 2205, wówczas ten kąt będzie wynosił 2 * M_PI. To, co widziałaś, było zapewne aliasingowe, ponieważ twoja referencja miała częstotliwość 2205Hz, a wszystkie częstotliwości powyżej 1102Hz były aliasingowane do niższych częstotliwości.