上次我们用递归实现了FFT的算法(可以参考我的另一篇文章快速傅里叶变换学习(超详细,附代码实现)_Patarw_Li的博客-CSDN博客),这次我们实现FFT算法的迭代版本。
下面是8点FFT的蝶形图,如果看懂了这张图,那么写出迭代版本的代码也是轻而易举,推荐看这个老师的视频:快速傅里叶变换FFT_哔哩哔哩_bilibili
下面是DFT的公式,也可以帮我们理解这个图:
下面是旋转因子Wn,我们的代码里面也会涉及到这个:
下面介绍一下蝶形图(如果已经理解了蝶形图则可以直接跳过)。
二点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