DSP

C# FFT的原理、解析及代码实现

2019-07-13 20:11发布

 

一、离散傅里叶变换回顾与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. (一)用迭代法实现
迭代法利用迭代函数运算实现,每次迭代进行一次抽取,直至抽取出的子序列只剩下一点,函数进行回代,并开始进行蝶形运算。先把我写的源代码贴上供大家参考:  
  1. //迭代法求FFT
  2. private void FFT(ref double[] DataIn_r, ref double[] DataIn_i, ref double[] DataOut_r, ref double[] DataOut_i)
  3. {
  4. int len = DataIn_r.Length;
  5. if (len == 1)//抽取到只有一个点时,将数据返回
  6. {
  7. DataOut_r = DataIn_r;
  8. DataOut_i = DataIn_i;
  9. return;
  10. }
  11. len = len / 2;
  12. //申请空间,以存储新抽取的子列
  13. double[] evenData_r = new double[len];
  14. double[] evenData_i = new double[len];
  15. double[] oddData_r = new double[len];
  16. double[] oddData_i = new double[len];
  17. //用以存储下一级函数返回的计算结果
  18. double[] evenResult_r = new double[len];
  19. double[] evenResult_i = new double[len];
  20. double[] oddResult_r = new double[len];
  21. double[] oddResult_i = new double[len];
  22. //按奇偶迭代
  23. for (int i = 0; i < len ; i++)
  24. {
  25. evenData_r[i] = DataIn_r[i * 2];
  26. evenData_i[i] = DataIn_i[i * 2];
  27. oddData_r[i] = DataIn_r[i * 2 + 1];
  28. oddData_i[i] = DataIn_i[i * 2 + 1];
  29. }//迭代
  30. FFT(ref evenData_r, ref evenData_i, ref evenResult_r, ref evenResult_i);
  31. FFT(ref oddData_r, ref oddData_i, ref oddResult_r, ref oddResult_i);
  32. //最下一级函数返回后,进行本级的蝶形计算
  33. double WN_r,WN_i;
  34. for (int i = 0; i < len; i++)
  35. {
  36. WN_r = Math.Cos(2 * Math.PI / (2 * len) * i);
  37. WN_i = -Math.Sin(2 * Math.PI / (2 * len) * i);
  38. DataOut_r[i] = evenResult_r[i] + oddResult_r[i] * WN_r - oddResult_i[i] * WN_i;
  39. DataOut_i[i] = evenResult_i[i] + oddResult_i[i] * WN_r + oddResult_r[i] * WN_i;
  40. DataOut_r[i + len] = evenResult_r[i] - oddResult_r[i] * WN_r + oddResult_i[i] * WN_i;
  41. DataOut_i[i + len] = evenResult_i[i] - oddResult_i[i] * WN_r - oddResult_r[i] * WN_i;
  42. }
  43. evenData_r = evenData_i = evenResult_r = evenResult_i = null;
  44. oddData_i = oddData_r = oddResult_i = oddResult_r = null;
  45. GC.Collect();//释放本级函数占用的资源
  46. }
显而易见的,迭代法对内存的消耗是很大的。下面来介绍下利用旋转因子进行的FFT算法。

(二)用旋转因子法实现

这个方法早就出现了,名字是我自己给起的,这种算法运用一个三重循环,围绕所谓的旋转因子WN来计算,也是比较流行的主流算法。 众所周知,在进行FFT运算之前,必须对时域序列重排。以往的重排算法是基于对下标二进制码的比特反转来实现的,但是这种方法用高级语言实现起来很麻烦。我写了一个简单易懂的算法。先来看下下标二进制码的规律: 我们来看表格第三列,即重排后的下标二进制码。我们竖着看,最高位一列的0、1的反转周期为2,次高位一列的反转周期为4,以此类推。因此我们可以写一个M次的循环,每一个循环对一列比特位进行操作,具体来说就是每隔半周期加上add=Math.Pow(2,M-L)。下面是源码:
  1. //对原数据组进行重排
  2. private void DataSort(ref double[] data_r, ref double[] data_i)
  3. {
  4. if (data_r.Length == 0 || data_i.Length == 0 || data_r.Length != data_i.Length)
  5. return;
  6. int len = data_r.Length;
  7. int[] count = new int[len];
  8. int M = (int)(Math.Log(len) / Math.Log(2));
  9. double[] temp_r = new double[len];
  10. double[] temp_i = new double[len];
  11. for (int i = 0; i < len; i++)
  12. {
  13. temp_r[i] = data_r[i];
  14. temp_i[i] = data_i[i];
  15. }
  16. for (int l = 0; l < M; l++)
  17. {
  18. int space = (int)Math.Pow(2, l);
  19. int add = (int)Math.Pow(2, M - l - 1);
  20. for (int i = 0; i < len; i++)
  21. {
  22. if ((i / space) % 2 != 0)
  23. count[i] += add;
  24. }
  25. }
  26. for (int i = 0; i < len; i++)
  27. {
  28. data_r[i] = temp_r[count[i]];
  29. data_i[i] = temp_i[count[i]];
  30. }
  31. }

再来研究蝶形运算的规律:
  • 我们发现每个蝶形的两个输入点只对本蝶有用,对其他蝶形的计算不产生影响;
  • 同一级中所有蝶形的输入点在同一竖直线上,意味着我们可以按级来运算,对于M级的蝶形,编个M次循环就好了;
  • 所有数据点在运算后不会窜位,即计算后可以将结果存入原来的地址空间。这就免去了像上面的迭代法那样浪费资源的现象;
  • 每级N/2个蝶都需要用到系数WN,这里称它为旋转因子。我们来观察旋转因子WN的规律。以8点的蝶形图为例:
 可见,第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#源码:  
  1. //FFT算法
  2. private void Dit2_FFT(ref double[] data_r, ref double[] data_i,ref double[] result_r,ref double[] result_i)
  3. {
  4. if (data_r.Length == 0 || data_i.Length == 0 || data_r.Length != data_i.Length)
  5. return;
  6. int len = data_r.Length;
  7. double[] X_r = new double[len];
  8. double[] X_i = new double[len];
  9. for (int i = 0; i < len; i++)//将源数据复制副本,避免影响源数据的安全性
  10. {
  11. X_r[i] = data_r[i];
  12. X_i[i] = data_i[i];
  13. }
  14. DataSort(ref X_r, ref X_i);//位置重排
  15. double WN_r,WN_i;//旋转因子
  16. int M = (int)(Math.Log(len) / Math.Log(2));//蝶形图级数
  17. for (int l = 0; l < M; l++)
  18. {
  19. int space = (int)Math.Pow(2, l);
  20. int num = space;//旋转因子个数
  21. double temp1_r, temp1_i, temp2_r, temp2_i;
  22. for (int i = 0; i < num; i++)
  23. {
  24. int p = (int)Math.Pow(2, M - 1 - l);//同一旋转因子有p个蝶
  25. WN_r = Math.Cos(2 * Math.PI / len * p * i);
  26. WN_i = -Math.Sin(2 * Math.PI / len * p * i);
  27. for (int j = 0, n = i; j < p; j++, n += (int)Math.Pow(2, l + 1))
  28. {
  29. temp1_r = X_r[n];
  30. temp1_i = X_i[n];
  31. temp2_r = X_r[n + space];
  32. temp2_i = X_i[n + space];//为蝶形的两个输入数据作副本,对副本进行计算,避免数据被修改后参加下一次计算
  33. X_r[n] = temp1_r + temp2_r * WN_r - temp2_i * WN_i;
  34. X_i[n] = temp1_i + temp2_i * WN_r + temp2_r * WN_i;
  35. X_r[n + space] = temp1_r - temp2_r * WN_r + temp2_i * WN_i;
  36. X_i[n + space] = temp1_i - temp2_i * WN_r - temp2_r * WN_i;
  37. }
  38. }
  39. }
  40. result_r = X_r;
  41. result_i = X_i;
  42. }
我写了一个简单的测试界面,包含了上述的所有算法,效果如下:
我已将源码上传资源,方便大家学习交流:点击打开链接