C++实现FFT算法（迭代版本）

上次我们用递归实现了FFT的算法（可以参考我的另一篇文章快速傅里叶变换学习（超详细，附代码实现）_Patarw_Li的博客-CSDN博客），这次我们实现FFT算法的迭代版本。

下面是8点FFT的蝶形图，如果看懂了这张图，那么写出迭代版本的代码也是轻而易举，推荐看这个老师的视频：快速傅里叶变换FFT_哔哩哔哩_bilibili

下面是DFT的公式，也可以帮我们理解这个图：

下面是旋转因子Wn，我们的代码里面也会涉及到这个：

将离散傅里叶变换公式拆分成奇偶项，则前 N/2 个点可以表示为：

同理，后 N/2 个点可以表示为：

 

由此可知，后 N/2 个点的值完全可以通过计算前 N/2 个点时的中间过程值确定。对 A[k] 与 B[k] 继续进行奇偶分解，直至变成 2 点的 DFT，这样就可以避免很多的重复计算，实现了快速离散傅里叶变换（FFT）的过程。

下面介绍一下蝶形图（如果已经理解了蝶形图则可以直接跳过）。

二点FFT的蝶形图

中间的为蝶形运算符，往上加，往下减（x(1)应该是x(1)*Wn，要乘旋转因子）。

四点FFT的蝶形图

可以看到四点FFT的蝶形图就是在两点的基础上又加了一层： 

 八点FFT的蝶形图

同样的，八点的FFT蝶形图也是在四点的基础上加了一层：

所以很容易直到我们要迭代的次数time=，如果是八点的话就要迭代三次。

还有一个就是上面蝶形图左边的 x(i) 的排序不是按照升序或降序拍的，而是将 0−7 的二进制写出来后，将二进制的高位、低位互换后得到的

C++代码：

#include <iostream>
#include <vector>
#include <complex>
#include <cmath>
using namespace std;

//定义圆周率派
#define PI acos(-1)

//FFT，输入P，返回y
void FFT(vector<complex<double>> P, vector<complex<double>> &y, int n){
    int time = log(n) / log(2);//迭代次数
    complex<double> Wn(cos(-2*PI/n), sin(-2*PI/n));//旋转因子
    //交换输入位置
    for(int i = 0;i < n/2;i++){
        if(i%2 == 1){
            complex<double> temp = P[i];
            P[i] = P[i+n/2-1];
            P[i+n/2-1] = temp;
        }
    }
    
    //要进行time次循环
    for(int i = 0;i < time;i++){
        int m = pow(2,i);
        //求第i次迭代的输出y
        for(int k = 0;k < n/(m*2);k++){
            for(int j = 0;j < m;j++){
                y[j+k*2] = P[j+k*2] + pow(Wn, j)*P[j+k*2+m];
                y[j+k*2+m] = P[j+k*2] - pow(Wn, j)*P[j+k*2+m];
            }
        }
        //把y的值赋为P，作为下一次迭代的输入
        for(int j = 0;j < n;j++)
            P[j] = y[j];
    }
    
}


int main() {
    //输入
    vector<complex<double>> P = {1, 1,1,1};
    //vector<complex<double>> P = {0, 1,0,3};

    //n的长度必须为2的幂次
    int n = P.size();
    //输出
    vector<complex<double>> y(n);
    
    FFT(P, y, n);
    
    //输出j结果
    for(int i = 0;i < n;i++){
        cout<<y[i]<<endl;
    }   
}

测试：

（1）输入：1,1,1,1；输出：4,0,0,0

 （2）输入：0,1,0,3；输出：4,2j,-4,-2j