c - 在 C 中为实数输入编写一个简单的离散傅里叶变换

Question

所以我尝试用 C 语言编写离散傅里叶变换来处理真正的 32 位浮点 wav 文件。它一次读取 2 帧(每个 channel 一个，但出于我的目的，我假设它们是相同的，所以我使用帧 [0])。这段代码应该通过用频率 20,40,60,...,10000 探测输入文件来写出幅度谱。我在输入帧上使用汉宁窗。如果可以的话，我想避免使用复数。当我运行它时，它给了我一些非常奇怪的振幅(其中大部分非常小，并且与正确的频率无关)，这让我相信我在计算中犯了一个根本性的错误。有人可以深入了解这里发生的事情吗？这是我的代码:

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

我认为错误在于角度的计算。每个样本的角度增量取决于采样频率。像这样(你似乎有 44100Hz):

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

您的样本窗口将包含一个完整的周期，用于您的最低探测频率 20Hz。如果将 n 循环到 2205，则该角度将为 2*M_PI。您看到的可能是混叠，因为您的引用频率为 2205Hz，而所有高于 1102Hz 的频率都混叠为较低的频率。