C에서 실제 입력에 대한 간단한 이산 푸리에 변환 작성

Question

그래서 C에서 실수 푸리에 변환을 작성하여 실제 32 비트 부동 WAV 파일로 작업하려고합니다. 한 번에 2 프레임 씩 (각 채널마다 하나씩 읽습니다. 그러나 제 목적으로는 둘 다 같아서 프레임 [0]을 사용한다고 가정합니다). 이 코드는 주파수 20, 40, 60, ..., 10000으로 입력 파일을 조사하여 진폭 스펙트럼을 작성합니다. 입력 프레임에 해닝 창을 사용하고 있습니다. 가능한 경우 복소수 사용을 피하고 싶습니다. 이것을 실행할 때, 그것은 아주 이상한 진폭을 제공합니다. (대부분이 매우 작고 정확한 주파수와 관련이 없습니다.) 이것은 제가 계산에서 근본적인 실수를 저지르고 있다고 믿게합니다. 누군가 여기서 일어나는 일에 대해 통찰력을 줄 수 있습니까? 약간 컴파일하고 가짜 샘플을 만들 개편, 나는 하지이 모두 0을받을 수 있나요 - 아래 코드와

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

가 : 여기 내 코드입니다. 나는에서 마지막에 출력 호출을 변경 :

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

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

에이 그냥 쉽게 비 - 제로 데이터를 볼 수 있습니다. 어쩌면 유일한 문제는 스케일 팩터를 이해하는 것일까 요? 테스트 케이스로 순수 톤을 생성하기 위해 입력을 위장한 방법에 유의하십시오.

#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

각도의 계산에 오차가 있다고 생각합니다. 각 샘플의 각도 증가는 샘플링 빈도에 따라 다릅니다. 이 같은 뭔가 (당신이에서 44100Hz을 갖고있는 것 같다) :

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

샘플 윈도우가 가장 낮은 프로브 된 주파수 20Hz에서에 대해 하나의 전체주기를 포함합니다. n을 2205까지 반복하면 그 각도는 2 * M_PI가됩니다. 레퍼런스에 2205Hz의 주파수가 있고 1102Hz 이상의 모든 주파수에 주파수가 더 낮기 때문에 앨리어싱을 보았습니다.