反向 fft 产生接近但错误的输出答案

【问题标题】：inverse fft producing close but wrong output反向 fft 产生接近但错误的输出
【发布时间】：2013-12-11 16:50:54
【问题描述】：

我正在执行 radix-2 dif 逆 fft。我正在使用共轭和缩放的属性来返回结果。我对输入向量进行共轭，执行常规 radix-2 fft（不是 ifft），对结果进行共轭，然后按 1.0/N 缩放。但是，我没有得到正确的结果：

int main(){ 
    const int n = 4;
    complex<double> x[n];

        // Test signal

    x[0] = complex<double>(10,0);
    x[1] = complex<double>(-2,0);
    x[2] = complex<double>(-2,2);
    x[3] = complex<double>(-2,-2);

        print(x,n);

    fft_inverse(x,n); 

    print(x,n); 

}
//dif fft. works
void fft(complex<double> X[], int N){
        if(N == 1){return;} 

    complex<double> *temp = new complex<double>[N]; 
    for(int i=0; i<N; i++){
        temp[i]=X[i];
        }   
        for(int i = 0; i<N/2; i++){
        complex<double> tw(cos(-2*M_PI*i/N),sin(-2*M_PI*i/N)); 
            X[i] = temp[i] + temp[i+N/2];
            X[i+N/2] = temp[i]-temp[i+N/2];
            X[i+N/2] = X[i+N/2]*tw;   
        }

        fft(X,N/2);
        fft(X+N/2,N/2);
}
void fft_inverse(complex<double> X[], int N){
    //conjugate
    for(int i = 0; i<N;i++){
        X[i] = conj(X[i]);
    }
    //perform fft
    fft(X,N);
    //conjugate again
    for(int i = 0; i<N;i++){
        X[i] = conj(X[i]);
    }
    //scale by 1.0/N
    double norm_N = 1.0/N;
    for(int i = 0; i<N;i++){
        X[i] *= norm_N;
    }
}

这是我的结果：输入：

(10,0) (-2,0) (-2,2) (-2,-2)

输出：

(1,-0) (3,1) (2.5,-0.5) (3.5,-0.5)

输出应该是：

(1,0) (2,0) (3,0) (4,0)

发生了什么事？我已经测试了我的fft 的输出应该是什么并收到了正确的结果，所以我不确定问题是什么。

【问题讨论】：

标签： c++ math fft

【解决方案1】：

看起来您的代码给出了正确的输出，但 bin 的顺序错误：

octave> X = [ 10, -2, -2+2i, -2-2i ] 
X =

   10 +  0i   -2 +  0i   -2 +  2i   -2 -  2i

octave> x = ifft(X)
x =

   1.00000 + 0.00000i   2.50000 - 0.50000i   3.00000 + 1.00000i   3.50000 - 0.50000i

fft_inverse() 看起来不错，所以我怀疑fft() 可以使用一些进一步的测试/调试 - 可能您需要为索引做一些位反转。

【讨论】：

我同意fft 可能是罪魁祸首。如果时域信号确实是真实的1, 2, 3, 4，那么它的 FFT 应该是对称的。它不是。
我倾向于同意，但是，我已经在输入 (1,0), (2,0), (3,0), (4,0) 上测试了 fft 并收到了正确的结果..
尝试不同和/或更长的 FFT - 将结果与已知良好的 FFT 进行比较，例如MATLAB/Octave/Wolfram Alpha。
好吧，我们都是对的。实函数的 FFT 是对称的，你的结果也是对称的 :) 只是它被移动了，当你只有四个值时很难看到。
好的。感谢您的帮助

【解决方案2】：

让我们回顾一下基础知识：长度为 2N 的逆 fft 在单位根 q^i 处计算 2N 次多项式 p(z)，其中 q^N = -1。为了快速做到这一点，多项式分为偶数和奇数部分，p(z)=pe(z^2)+z*po(z^2)。那么

p(q^i)    =pe(q^2i)+q^i*po(q^2i)
p(q^(N+i))=pe(q^2i)-q^i*po(q^2i)

在实现中，pe和po的值是通过递归调用长度为N的ifft获得的。

正向 fft 现在是它的逆运算，根据多项式的值计算多项式的系数。为此，通过颠倒上述公式，将p(q^i) 的值转换为pe(q^2i) 和po(q^2i) 的值，

pe(q^2i)=0.5*(p(q^i)+p(q^(N+i)))
po(q^2i)=0.5*(p(q^i)-p(q^(N+i)))*q^(-i)

你应该从你的代码中识别出来。然后 fft 被递归调用以从它们现在计算的值中确定 pe 和 po 的系数。因子 0.5 通常被省略，并在 fft 例程之外作为数组长度除以收集。

您现在出错的地方在于，pe 和 po 的系数在 p 的系数序列中交替出现、交错排列。著名的蝴蝶图案。

void fft(complex<double> X[], int N){
    fft_internal(X,N,1);
}

void fft_internal(complex<double> X[], int N, int step){
    if(N == 1){return;} 

    complex<double> *temp = new complex<double>[N]; 
    for(int i=0; i<N; i++){
        temp[i]=X[step*i];
    }   
    for(int i = 0; i<N/2; i++){
    complex<double> tw(cos(-2*M_PI*i/N),sin(-2*M_PI*i/N)); 
        X[2*step*i     ] =  temp[i] + temp[i+N/2];
        X[2*step*i+step] = (temp[i] - temp[i+N/2])*tw;   
    }

    fft_internal(X     ,N/2, 2*step);
    fft_internal(X+step,N/2, 2*step);
}

使用正确重新排序的输入 ([2][1])，

// Test signal

x[0] = complex<double>(10,0);
x[1] = complex<double>(-2,2);
x[2] = complex<double>(-2,0);
x[3] = complex<double>(-2,-2);

这会返回预期的结果。

【讨论】：