快速傅里叶变换（FFT）的原理、实现及代码解析（附C#源码）
来源：互联网 发布：淘宝摄影相机推荐 编辑：程序博客网 时间：2024/05/16 08:38
一、离散傅里叶变换回顾与FFT的引出
对于长度为N点的数字信号序列 ，定义其离散傅里叶变换为：
我们知道，利用系数 的性质可以大大减少DFT的计算量，这种算法就是快速离散傅里叶变换FFT。需要说明的是，FFT不是一种新的变换，而是一种求DFT的快速计算机算法。
对序列 按奇偶分成两列，重写DFT表达式：
他们分别是偶相列和奇数项列的DFT：
。
那么，对于一个 的序列进行不断分解，就可以得出如下所谓的蝶形图：
二、FFT运算的规律及实现
- 观察上面的蝶形图，我们发现，一个N点序列（注意N必须是2的幂次，如果原序列不是，可以填上0。如：1000点的x(n)可以在后面天24个0，使之变为1024点），可以进行M=logN次的抽取，对应蝶形图就有M级。如上图中8点序列的蝶形图就有3级；每一级都是N/2个蝶形运算；
每一个蝶形的计算如下图：
X(k)=X1+X2*WN=X1_r+j(X1_i)+(X2_r+jX2_i)*(WN_r+jWN_i)
展开得到X(k)的实部和虚步分别为：X_r = X1_r + X2_r * WN_r - X2_i * WN_i,
X_i = X1_i + X2_i * WN_r + X2_r * WN_i.
同理得到X2(K)的实部和虚步分别为：X2_r = X1_r - X2_r * WN_r + X2_i * WN_i,
X2_i = X1_i - X2_i * WN_r - X2_r * WN_i.
（一）用迭代法实现
迭代法利用迭代函数运算实现，每次迭代进行一次抽取，直至抽取出的子序列只剩下一点，函数进行回代，并开始进行蝶形运算。先把我写的源代码贴上供大家参考：
//迭代法求FFT private void FFT(ref double[] DataIn_r, ref double[] DataIn_i, ref double[] DataOut_r, ref double[] DataOut_i) { int len = DataIn_r.Length; if (len == 1)//抽取到只有一个点时，将数据返回 { DataOut_r = DataIn_r; DataOut_i = DataIn_i; return; } len = len / 2; //申请空间，以存储新抽取的子列 double[] evenData_r = new double[len]; double[] evenData_i = new double[len]; double[] oddData_r = new double[len]; double[] oddData_i = new double[len]; //用以存储下一级函数返回的计算结果 double[] evenResult_r = new double[len]; double[] evenResult_i = new double[len]; double[] oddResult_r = new double[len]; double[] oddResult_i = new double[len]; //按奇偶迭代 for (int i = 0; i < len ; i++) { evenData_r[i] = DataIn_r[i * 2]; evenData_i[i] = DataIn_i[i * 2]; oddData_r[i] = DataIn_r[i * 2 + 1]; oddData_i[i] = DataIn_i[i * 2 + 1]; }//迭代 FFT(ref evenData_r, ref evenData_i, ref evenResult_r, ref evenResult_i); FFT(ref oddData_r, ref oddData_i, ref oddResult_r, ref oddResult_i); //最下一级函数返回后，进行本级的蝶形计算 double WN_r,WN_i; for (int i = 0; i < len; i++) { WN_r = Math.Cos(2 * Math.PI / (2 * len) * i); WN_i = -Math.Sin(2 * Math.PI / (2 * len) * i); DataOut_r[i] = evenResult_r[i] + oddResult_r[i] * WN_r - oddResult_i[i] * WN_i; DataOut_i[i] = evenResult_i[i] + oddResult_i[i] * WN_r + oddResult_r[i] * WN_i; DataOut_r[i + len] = evenResult_r[i] - oddResult_r[i] * WN_r + oddResult_i[i] * WN_i; DataOut_i[i + len] = evenResult_i[i] - oddResult_i[i] * WN_r - oddResult_r[i] * WN_i; } evenData_r = evenData_i = evenResult_r = evenResult_i = null; oddData_i = oddData_r = oddResult_i = oddResult_r = null; GC.Collect();//释放本级函数占用的资源 }显而易见的，迭代法对内存的消耗是很大的。下面来介绍下利用旋转因子进行的FFT算法。
（二）用旋转因子法实现
这个方法早就出现了，名字是我自己给起的，这种算法运用一个三重循环，围绕所谓的旋转因子WN来计算，也是比较流行的主流算法。
众所周知，在进行FFT运算之前，必须对时域序列重排。以往的重排算法是基于对下标二进制码的比特反转来实现的，但是这种方法用高级语言实现起来很麻烦。我写了一个简单易懂的算法。先来看下下标二进制码的规律：
我们来看表格第三列，即重排后的下标二进制码。我们竖着看，最高位一列的0、1的反转周期为2，次高位一列的反转周期为4，以此类推。因此我们可以写一个M次的循环，每一个循环对一列比特位进行操作，具体来说就是每隔半周期加上add=Math.Pow(2,M-L)。下面是源码：
//对原数据组进行重排 private void DataSort(ref double[] data_r, ref double[] data_i) { if (data_r.Length == 0 || data_i.Length == 0 || data_r.Length != data_i.Length) return; int len = data_r.Length; int[] count = new int[len]; int M = (int)(Math.Log(len) / Math.Log(2)); double[] temp_r = new double[len]; double[] temp_i = new double[len]; for (int i = 0; i < len; i++) { temp_r[i] = data_r[i]; temp_i[i] = data_i[i]; } for (int l = 0; l < M; l++) { int space = (int)Math.Pow(2, l); int add = (int)Math.Pow(2, M - l - 1); for (int i = 0; i < len; i++) { if ((i / space) % 2 != 0) count[i] += add; } } for (int i = 0; i < len; i++) { data_r[i] = temp_r[count[i]]; data_i[i] = temp_i[count[i]]; } }
再来研究蝶形运算的规律：
可见，第L级的旋转因子为：
j代表旋转因子的个数，第L级的旋转因数个数为num=Math.Pow(2, L).(2的幂不会打
)同时，还可以看到，每个蝶的两个输入点下标跨度是不一样的。比如第一级中是相邻两个数据作蝶运算，第二级中是两个输入点隔开一个点，而第三级隔开三个点。不难找到规律：第L级中，第二个输入点的坐标是第一个点的坐标+space，space=Math.Pow(2, L)=num。
利用以上性质，FFT的算法是写一个三重循环，
第一重循环对每一级运算（每级含num=Math.Pow(2, L)个蝶形）；
第二重对每一个旋转因子对应的蝶运算， 那么有几个蝶呢？很简单，每级都应该有N/2个蝶，而每个因子对应N/2 / num个蝶；
第三重循环对每个蝶进行计算，需要注意的一是循环下标起始点的位置，二是每次计算需要申明临时变量来保存输入数据点。下面奉上C#源码：
//FFT算法 private void Dit2_FFT(ref double[] data_r, ref double[] data_i,ref double[] result_r,ref double[] result_i) { if (data_r.Length == 0 || data_i.Length == 0 || data_r.Length != data_i.Length) return; int len = data_r.Length; double[] X_r = new double[len]; double[] X_i = new double[len]; for (int i = 0; i < len; i++)//将源数据复制副本，避免影响源数据的安全性 { X_r[i] = data_r[i]; X_i[i] = data_i[i]; } DataSort(ref X_r, ref X_i);//位置重排 double WN_r,WN_i;//旋转因子 int M = (int)(Math.Log(len) / Math.Log(2));//蝶形图级数 for (int l = 0; l < M; l++) { int space = (int)Math.Pow(2, l); int num = space;//旋转因子个数 double temp1_r, temp1_i, temp2_r, temp2_i; for (int i = 0; i < num; i++) { int p = (int)Math.Pow(2, M - 1 - l);//同一旋转因子有p个蝶 WN_r = Math.Cos(2 * Math.PI / len * p * i); WN_i = -Math.Sin(2 * Math.PI / len * p * i); for (int j = 0, n = i; j < p; j++, n += (int)Math.Pow(2, l + 1)) { temp1_r = X_r[n]; temp1_i = X_i[n]; temp2_r = X_r[n + space]; temp2_i = X_i[n + space];//为蝶形的两个输入数据作副本，对副本进行计算，避免数据被修改后参加下一次计算 X_r[n] = temp1_r + temp2_r * WN_r - temp2_i * WN_i; X_i[n] = temp1_i + temp2_i * WN_r + temp2_r * WN_i; X_r[n + space] = temp1_r - temp2_r * WN_r + temp2_i * WN_i; X_i[n + space] = temp1_i - temp2_i * WN_r - temp2_r * WN_i; } } } result_r = X_r; result_i = X_i; }我写了一个简单的测试界面，包含了上述的所有算法，效果如下：
我已将源码上传资源，方便大家学习交流：点击打开链接
http://lib.csdn.net/article/computernetworks/14327
- 快速傅里叶变换（FFT）的原理、实现及代码解析（附C#源码）
- 快速傅里叶变换（FFT）的原理、实现及代码解析（附C#源码）
- 快速傅里叶变换（FFT）的C#实现及详细注释
- 快速傅里叶变换(FFT)的原理及公式
- 快速傅里叶变换(FFT)的原理及公式
- 【C】快速傅里叶变换（FFT）讲解及实现
- 快速傅里叶变换（FFT）
- 快速傅里叶变换（FFT）
- 快速傅里叶变换（FFT）
- 快速傅里叶变换（FFT）
- 快速傅里叶变换（FFT）
- 快速傅里叶变换（FFT）
- 快速傅里叶变换（FFT）
- 快速离散傅里叶变换（FFT）C++实现
- Python实现快速傅里叶变换（FFT）
- 快速傅里叶变换算法原理简述 基于递归的fft实现
- 快速傅里叶变换(FFT)结果的物理意义是什么？（附Matlab程序）
- 快速傅里叶变换(FFT)算法C++实现代码
- SpringCloud（二）：服务发现与服务注册
- HBase Coprocessor编程案例
- Servlet处理Form表单
- js转换时间戳与日期字符串
- Springboot+React+antd开发问题记（一）
- 快速傅里叶变换（FFT）的原理、实现及代码解析（附C#源码）
- mysql 存储过程 -- 游标的使用
- Vue.js框架练习之购物车
- 【基础中心】笔记二--Java8新特性之Lambda表达式
- UrhoSharp iOS和tvOS 开发入门
- 使用 RMI + ZooKeeper 实现远程调用框架
- vsftpd.conf.5配置文件手册（中文版）
- HierarchicalBeanFactory
- double运算精度问题
