C++ 快速傅里叶变换

快速傅里叶变换

1 快速傅立换变换的简介
- 1.1 傅里叶变换的不足
- 1.2 快速傅里叶变换
2 快速傅里叶变换的实现
- 2.1一些变量的说明：
- 2.2代码实现
- 2.3程序的输出
3 小结

之前我们已经在文章里对傅里叶变换给出了自己角度的阐述，那么接下来我们就进一步讲述快速傅里叶变换，并且给出代码实现。

1 快速傅立换变换的简介

1.1 傅里叶变换的不足

对于一个长度为 M M M 的信号序列来讲，如果我们要进行傅里叶变换，根据公式：
F ( u ) = ∑ x = 0 M − 1 f ( x ) e − j 2 π ( u x / M ) ) (1) F(u)= \sum_{x=0}^{M-1}f(x)e^{-j2\pi(ux/M)}\tag{1}) F(u)=x=0∑M−1f(x)e−j2π(ux/M))(1)
由上面公式可以看出，计算 N 2 N^2 N2 的复数乘法和 N 2 N^2 N2 的复数加法，在信号序列长度较长时，会导致计算量巨大，使得计算速度非常慢。值得庆幸的是，在后人的不断探索中，发现了傅里叶变换中变换因子 W N u x W_{N}^{ux} WNux 中的一些规律，总结如下：

W N = e − j 2 π / N W_N=e^{-j2\pi/N} WN=e−j2π/N

W 0 = 1 W^{0}=1 W0=1, W N / 2 = − 1 W^{N/2}=-1 WN/2=−1

W N N + r = W N r W_{N}^{N+r}=W_{N}^{r} WNN+r=WNr, W N N / 2 + r = − W N r W_{N}^{N/2+r}=-W_{N}^{r} WNN/2+r=−WNr, W 2 N 2 r = W N r W_{2N}^{2r}=W_{N}^{r} W2N2r=WNr

1.2 快速傅里叶变换

矩阵形式下的傅里叶变换：
[ F ( 0 ) F ( 1 ) F ( 2 ) F ( 3 ) ] = [ 1 1 1 1 1 W 1 − 1 − W 1 1 − 1 1 − 1 1 − W 1 − 1 W 1 ] [ f ( 0 ) f ( 1 ) f ( 2 ) f ( 3 ) ] (2) \left[ \begin{matrix} F(0) \\ F(1) \\ F(2) \\ F(3) \end{matrix} \right] = \left[ \begin{matrix} 1 & 1 & 1 & 1\\ 1 & W^1 & -1 & -W^1\\ 1 & -1 & 1 & -1\\ 1 & -W^1 & -1 & W^1 \end{matrix} \right] \left[ \begin{matrix} f(0) \\ f(1)\\ f(2) \\ f(3) \end{matrix} \right] \tag{2} ⎣⎢⎢⎡F(0)F(1)F(2)F(3)⎦⎥⎥⎤=⎣⎢⎢⎡11111W1−1−W11−11−11−W1−1W1⎦⎥⎥⎤⎣⎢⎢⎡f(0)f(1)f(2)f(3)⎦⎥⎥⎤(2)
将该矩阵的第二列和第三列交换，可得：
[ F ( 0 ) F ( 1 ) F ( 2 ) F ( 3 ) ] = [ 1 1 1 1 1 − 1 W 1 − W 1 1 1 − 1 − 1 1 − 1 − W 1 W 1 ] [ f ( 0 ) f ( 2 ) f ( 1 ) f ( 3 ) ] (3) \left[ \begin{matrix} F(0) \\ F(1) \\ F(2) \\ F(3) \end{matrix} \right] = \left[ \begin{matrix} 1 & 1 & 1 & 1\\ 1 & -1 & W^1 & -W^1\\ 1 & 1 & -1 & -1\\ 1 & -1 & -W^1 & W^1 \end{matrix} \right] \left[ \begin{matrix} f(0) \\ f(2)\\ f(1) \\ f(3) \end{matrix} \right] \tag{3} ⎣⎢⎢⎡F(0)F(1)F(2)F(3)⎦⎥⎥⎤=⎣⎢⎢⎡11111−11−11W1−1−W11−W1−1W1⎦⎥⎥⎤⎣⎢⎢⎡f(0)f(2)f(1)f(3)⎦⎥⎥⎤(3)
4点的FFT快速算法信号流图如下所示：

我们可以从信号流图的左侧观察到原序列发生了变换，即变化后的序列索引对应的元素与变化前不一致，要想实现此变换也是比较简单的，只需要将原位置元素的索引的二进制左右调换后重新赋予新索引对应的元素即可，例如：

f ( 0 ) f(0) f(0)排序后为 f ( 0 ) f(0) f(0)， 0 0 0的二进制为 00 00 00，左右调换后为 00 00 00，所以不变。
f ( 1 ) f(1) f(1)排序后为 f ( 2 ) f(2) f(2)， 1 1 1的二进制为 01 01 01，左右调换后为 10 10 10，所以变为2。
f ( 2 ) f(2) f(2)排序后为 f ( 1 ) f(1) f(1)， 2 2 2的二进制为 10 10 10，左右调换后为 01 01 01，所以变为1。
f ( 3 ) f(3) f(3)排序}后为 f ( 3 ) f(3) f(3)， 3 3 3的二进制为 11 11 11，左右调换后为 11 11 11，所以不变。

W r W^{r} Wr因子分布的规律：

m = 0 m=0 m=0级时， W 2 r W_{2}^{r} W2r， r = 0 r=0 r=0
m = 1 m=1 m=1级时， W 4 r W_{4}^{r} W4r， r = 0 , 1 r=0, 1 r=0,1
m = 2 m=2 m=2级时， W 8 r W_{8}^{r} W8r， r = 0 , 1 , 2 , 3 r=0, 1, 2, 3 r=0,1,2,3
……
m m m级时， W 2 m + 1 r W_{2^{m+1}}^{r} W2m+1r， r = 0 , 1 , 2 , … … , 2 m − 1 r=0, 1, 2,……, 2^{m}-1 r=0,1,2,……,2m−1

2 快速傅里叶变换的实现

声明：本代码的主体是从一位博主处copy来的，本想注明出处，但是因未及时收藏导致后续竟找不到此博主的相关信息，对此我深表遗憾。本人在原博主代码的基础上加以修改（增加了反傅里叶变换功能、序列长度不为2的幂次方时对序列进行拓展的功能），并伴以详细的注释，以飨大家。此外，由于本人能力问题，代码还存在诸多不完美之处，例如：1、将序列填充至快速傅里叶变换长度之后，需要手动定义后续复数数组的长度以进行傅里叶变换；2、在实现功能的过程中，函数有些繁杂，且某些函数内部没有进行优化，还望诸位看客多多见谅。

2.1一些变量的说明：

2.2代码实现

#include <list>
#include <cmath>
#include <string>
#include <vector>
#include <iostream>const int PI = 3.1415926;
using namespace std;//定义复数结构
struct Complex
{double imaginary;double real;
};//进行FFT的数据，此处默认长度为2的幂次方
//double Datapre[] = {1, 2, 6, 4, 48, 6, 7, 8, 58, 10, 11, 96, 13, 2, 75, 16};
double Datapre[] = {100, 2, 56, 4, 48, 6, 45, 8, 58, 10};//复数乘法函数
Complex ComplexMulti(Complex One, Complex Two);//将输入的任意数组填充，以满足快速傅里叶变换
void Data_Expand(double *input, vector<double> &output, const int length);//将输入的实数数组转换为复数数组
void Real_Complex(vector<double> &input, Complex *output, const int length);//快速傅里叶变换函数
void FFT(Complex *input, Complex *StoreResult, const int length);//快速傅里叶反变换
void FFT_Inverse(Complex *arrayinput, Complex *arrayoutput, const int length);int main()
{//获取填充后的Data;vector<double> Data;Data_Expand(Datapre, Data, 10);//打印填充后的序列for (auto data : Data){cout << data << endl;}//因为下面的复数结构未采用动态数组结构，所以需要根据实际情况制定数组的大小//将输入数组转化为复数数组Complex array1D[16];//存储傅里叶变换的结果Complex Result[16];//存储反傅里叶变换的结果Complex Result_Inverse[16];Real_Complex(Data, array1D, 16);FFT(array1D, Result, 16);FFT_Inverse(Result, Result_Inverse, 16);//基于范围的循环，利用auto自动判断数组内元素的数据类型for (auto data : Result_Inverse){//输出傅里叶变换后的幅值cout << data.real << "\t" << data.imaginary << endl;}return 0;
}Complex ComplexMulti(Complex One, Complex Two)
{//Temp用来存储复数相乘的结果Complex Temp;Temp.imaginary = One.imaginary * Two.real + One.real * Two.imaginary;Temp.real = One.real * Two.real - One.imaginary * Two.imaginary;return Temp;
}//此处的length为原序列的长度
void Data_Expand(double *input, vector<double> &output, const int length)
{int i = 1, flag = 1;while (i < length){i = i * 2;}//获取符合快速傅里叶变换的长度int length_Best = i;if (length_Best != length){flag = 0;}if (flag){//把原序列填充到Vector中for (int j = 0; j < length; ++j){output.push_back(input[j]);}}else{//把原序列填充到Vector中for (int j = 0; j < length; ++j){output.push_back(input[j]);}//把需要拓展的部分进行填充for (int j = 0; j < length_Best - length; j++){output.push_back(0);}}
}//此处的length为填充后序列的长度
void Real_Complex(vector<double> &input, Complex *output, const int length)
{for (int i = 0; i < length; i++){output[i].real = input[i];output[i].imaginary = 0;}
}//FFT变换函数
//input: input array
//StoreResult: Complex array of output
//length: the length of input array
void FFT(Complex *input, Complex *StoreResult, const int length)
{//获取序列长度在二进制下的位数int BitNum = log2(length);//存放每个索引的二进制数，重复使用，每次用完需清零list<int> Bit;//暂时存放重新排列过后的输入序列vector<double> DataTemp1;vector<double> DataTemp2;//遍历序列的每个索引for (int i = 0; i < length; ++i){//flag用来将索引转化为二进制//index用来存放倒叙索引//flag2用来将二值化的索引倒序int flag = 1, index = 0, flag2 = 0;//遍历某个索引二进制下的每一位for (int j = 0; j < BitNum; ++j){//十进制转化为长度一致的二进制数，&位运算符作用于位，并逐位执行操作//从最低位（右侧）开始比较//example://a = 011, b1 = 001, b2 = 010 ,b3 = 100//a & b1 = 001, a & b2 = 010, a & b3 = 000int x = (i & flag) > 0 ? 1 : 0;//把x从Bit的前端压进Bit.push_front(x);//等价于flag = flag << 1，把flag的值左移一位的值赋给flagflag <<= 1;}//将原数组的索引倒序，通过it访问Bit的每一位for (auto it : Bit){//example:其相当于把二进制数从左到右设为2^0，2^1，2^2//Bit = 011//1: index = 0 + 0 * pow(2, 0) = 0//2: index = 0 + 1 * pow(2, 1) = 2//3: index = 2 + 1 * pow(2, 2) = 6 = 110index += it * pow(2, flag2++);}//把Data[index]从DataTemp的后端压进，其实现了原序列的数据的位置调换//变换前：f(0), f(1), f(2), f(3), f(4), f(5), f(6), f(7)//变换后：f(0), f(4), f(2), f(6), f(1), f(5), f(3), f(7)DataTemp1.push_back(input[index].real);DataTemp2.push_back(input[index].imaginary);//清空BitBit.clear();}for (int i = 0; i < length; i++){//将数据转移到复数结构里，StoreResult[i]的索引与原序列的索引是不一样的，一定要注意StoreResult[i].real = DataTemp1[i];StoreResult[i].imaginary = DataTemp2[i];}//需要运算的级数int Level = log2(length);Complex Temp, up, down;//进行蝶形运算for (int i = 1; i <= Level; i++){//定义旋转因子Complex Factor;//没有交叉的蝶形结的距离(不要去想索引)//其距离表示的是两个彼此不交叉的蝶型结在数组内的距离int BowknotDis = 2 << (i - 1);//同一蝶形计算中两个数字的距离（旋转因子的个数）//此处距离也表示的是两个数据在数组内的距离(不要去想索引)int CalDis = BowknotDis / 2;for (int j = 0; j < CalDis; j++){//每一级蝶形运算中有CalDis个不同旋转因子//计算每一级（i）所需要的旋转因子Factor.real = cos(2 * PI / pow(2, i) * j);Factor.imaginary = -sin(2 * PI / pow(2, i) * j);//遍历每一级的每个结for (int k = j; k < length - 1; k += BowknotDis){//StoreResult[k]表示蝶形运算左上的元素//StoreResult[k + CalDis]表示蝶形运算左下的元素//Temp表示蝶形运算右侧输出结构的后半部分Temp = ComplexMulti(Factor, StoreResult[k + CalDis]);//up表示蝶形运算右上的元素up.imaginary = StoreResult[k].imaginary + Temp.imaginary;up.real = StoreResult[k].real + Temp.real;//down表示蝶形运算右下的元素down.imaginary = StoreResult[k].imaginary - Temp.imaginary;down.real = StoreResult[k].real - Temp.real;//将蝶形运算输出的结果装载入StoreResultStoreResult[k] = up;StoreResult[k + CalDis] = down;}}}
}void FFT_Inverse(Complex *arrayinput, Complex *arrayoutput, const int length)
{//对傅里叶变换后的复数数组求共轭for (int i = 0; i < length; i++){arrayinput[i].imaginary = -arrayinput[i].imaginary;}//一维快速傅里叶变换FFT(arrayinput, arrayoutput, length);//时域信号求共轭，并进行归一化for (int i = 0; i < length; i++){arrayoutput[i].real = arrayoutput[i].real / length;arrayoutput[i].imaginary = -arrayoutput[i].imaginary / length;}
}

2.3程序的输出

3 小结

通过快速傅里叶变换，算法的时间复杂度由 O ( n 2 ) O(n^{2}) O(n2)变化为 O ( n l o g ( n ) ) O(nlog(n)) O(nlog(n))，FFT在信号长度较大时加速效果较为明显。

至此，本文内容已全部结束，有道是”独乐乐不如众乐乐“，欢迎大家在评论区发表自己意见、交流想法，谢谢！！！