FFT 实施产生毛刺

Question

我在白噪声的 FFT 图中遇到了一个奇怪的故障：

我用参考程序检查过，虽然噪声文件似乎没问题。是实现中的错误吗？

void four1(float data[], int nn, int isign) {
    int n, mmax, m, j, istep, i;
    float wtemp, wr, wpr, wpi, wi, theta;
    float tempr, tempi;

    n = nn << 1;
    j = 1;
    for (int i = 1; i < n; i += 2) {
        if (j > i) {
            tempr = data[j];
            data[j] = data[i];
            data[i] = tempr;
            tempr = data[j + 1];
            data[j + 1] = data[i + 1];
            data[i + 1] = tempr;
        }
        m = n >> 1;
        while (m >= 2 && j > m) {
            j -= m;
            m >>= 1;
        }
        j += m;
    }
    mmax = 2;
    while (n > mmax) {
        istep = 2 * mmax;
        theta = TWOPI / (isign * mmax);
        wtemp = sin(0.5 * theta);
        wpr = -2.0 * wtemp * wtemp;
        wpi = sin(theta);
        wr = 1.0;
        wi = 0.0;
        for (m = 1; m < mmax; m += 2) {
            for (i = m; i <= n; i += istep) {
                j = i + mmax;
                tempr = wr * data[j] - wi * data[j + 1];
                tempi = wr * data[j + 1] + wi * data[j];
                data[j] = data[i] - tempr;
                data[j + 1] = data[i + 1] - tempi;
                data[i] += tempr;
                data[i + 1] += tempi;
            }
            wr = (wtemp = wr) * wpr - wi * wpi + wr;
            wi = wi * wpr + wtemp * wpi + wi;
        }
        mmax = istep;
    }
}

Answer 1

除了一些小改动外，这段代码似乎是从 C 中的数字食谱第 2 版中取出来的。此函数的文档（摘自书中）指出：

Replaces data[1..2*nn] by its discrete Fourier transform, if isign is input as 1; or replaces data[1..2*nn] by nn times its inverse discrete Fourier transform, if isign is input as −1. data is a complex array of length nn or, equivalently, a real array of length 2*nn. nn MUST be an integer power of 2 (this is not checked for!).

给定一个从 1 开始索引的输入数组，此实现会产生正确的结果。您可以通过分配大小为 2*nn+1 的 C 数组并从索引 1 开始填充数组来选择使用相同的索引约定。或者，您可以传递一个大小为 2*nn 的数组，该数组已从索引开始填充0，但调用 four1(data-1, nn, isign)（注意 data 数组上的 -1 偏移量）。

FFT 实施产生毛刺

FFT implementation produces a glitch

c

c++

android

fft