在C语言中编写一个简单的离散傅里叶变换，用于处理实数输入。

Question

我正在尝试用C语言编写离散傅里叶变换，以便处理32位浮点wav文件。每次读取2帧数据（每个声道一帧，但为了我的目的，我假设它们都相同，因此我只使用frame [0]）。该代码旨在通过用频率20、40、60、......、10000来探测输入文件并输出幅度谱。我在输入帧上使用Hanning窗口。如果可以，我想避免使用复数。当我运行此代码时，它给我一些非常奇怪的振幅（其中大多数非常小，并且与正确频率不相关），这使我相信我在计算中犯了一个基本错误。有人能否提供关于这里发生了什么的见解？以下是我的代码：

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

- maniciam

---

看不出明显的错误，但建议生成测试用例并检查系数的数学属性是否正确。例如，实数输入意味着对称系数。 - Keith

还有一个用于执行傅里叶变换的 FFTW 库。 - Stefano Sanfilippo

@keith 我没有仔细查看 OP 代码，但如果它是离散函数的“经典”离散变换，你不会计算所有系数（因为你知道它们是对称的），也不会得到对称性，而是循环对称性，这可以被检查，当然，但它并不像简单的对称性检查那么平凡。 - Stefano Sanfilippo

一些与您的错误无关的建议。在现代 CPU 上，如果内存不是问题，请使用 double 而不是 float，因为 double 比 float 快得多得多，而且不要使用 pow(a,2)，而是使用 (a*a)。 pow 使用非常昂贵的指数函数调用，如果您只需要进行平方运算，则这是不必要的。 - Mark Lakata

你能否发布一张输入数据和输出频谱的图片？有可能存在你不知道的低频分量。即使是你认为只有一个纯正正弦波的FFT也会有很多峰值。频谱并不像你想象的那么干净。 - Mark Lakata

Answer 1

我认为错误在于角度的计算。每个样本的角度增量取决于采样频率。 类似这样（你似乎有44100Hz）：

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

你的样本窗口将包含最低探测频率20Hz的一个完整周期。如果你将n循环增加到2205，那么该角度将为2*M_PI。你看到的可能是混叠，因为你的参考频率为2205Hz，所有高于1102Hz的频率都被混叠到了较低的频率。

- rdrmntn

Answer 2

使用以下代码 - 只略微重新组织以编译并创建一个假样本，我没有得到全部零。 我已将输出调用从结尾更改为：

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

- Keith