C++实现FFT算法(迭代版本)

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

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

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

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

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

二点FFT的蝶形图

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

四点FFT的蝶形图

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

 八点FFT的蝶形图

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

所以很容易直到我们要迭代的次数time=log_{2}(n),如果是八点的话就要迭代三次。

还有一个就是上面蝶形图左边的 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