DSP

FFT算法的完整DSP实现

2019-07-13 10:00发布

【原文:http://blog.csdn.net/xiahouzuoxin/article/details/9790455】 傅里叶变换或者FFT的理论参考: [1] http://www.dspguide.com/ch12/2.htm       The Scientist and Engineer's Guide to Digital Signal Processing,   By Steven W. Smith, Ph.D. [2] http://blog.csdn.net/v_JULY_v/article/details/6196862,可当作[1]的中文参考 [3] 任意一本数字信号处理教材,上面都有详细的推导DCT求解转换为FFT求解的过程 [4] TI文档:基于TMS320C64x+DSP的FFT实现。 使用baidu/google可以搜索到。 另外,FFT的开源代码可参考: [1] FFTW: http://www.fftw.org/ 最快,最好的开源FFT。 [2] FFTReal: http://ldesoras.free.fr/prod.html#src_fftreal 轻量级FFT算法实现 [3] KISS FFT: http://sourceforge.net/projects/kissfft/ 简单易用的FFT的C语言实现

1. 有关FFT理论的一点小小解释

关于FFT这里只想提到两点: (1)DFT变换对的表达式(必须记住

          —— 称旋转因子
(2)FFT用途——目标只有一个,加速DFT的计算效率。 DFT计算X(k)需要N^2次复数乘法和N(N-1)次复数加法;FFT将N^2的计算量降为
FFT其实是很难的东西,即使常年在这个领域下打拼的科学家也未必能很好的写出FFT的算法。”——摘自参考上面提供的参考文献[1] 因此,我们不必太过纠结于细节,当明白FFT理论后,将已有的算法挪过来用就OK了,不必为闭着教材写不出FFT而郁闷不堪。
FFT的BASIC程序伪代码如下:
IFFT的BASIC程序伪代码如下(IFFT通过调用FFT计算):

FFT算法的流程图如下图,总结为3过程3循环: (1)3过程:单点时域分解(倒位序过程) + 单点时域计算单点频谱 + 频域合成 (2)3循环:外循环——分解次数,中循环——sub-DFT运算,内循环——2点蝶形算法
分解过程或者说倒位序的获得参考下图理解:

2. FFT的DSP实现

下面为本人使用C语言实现的FFT及IFFT算法实例,能计算任意以2为对数底的采样点数的FFT,算法参考上面给的流程图。 [cpp] view plaincopy在CODE上查看代码片派生到我的代码片
  1. /* 
  2.  * zx_fft.h 
  3.  * 
  4.  *  Created on: 2013-8-5 
  5.  *      Author: monkeyzx 
  6.  */  
  7.   
  8. #ifndef ZX_FFT_H_  
  9. #define ZX_FFT_H_  
  10.   
  11. typedef float          FFT_TYPE;  
  12.   
  13. #ifndef PI  
  14. #define PI             (3.14159265f)  
  15. #endif  
  16.   
  17. typedef struct complex_st {  
  18.     FFT_TYPE real;  
  19.     FFT_TYPE img;  
  20. } complex;  
  21.   
  22. int fft(complex *x, int N);  
  23. int ifft(complex *x, int N);  
  24. void zx_fft(void);  
  25.   
  26. #endif /* ZX_FFT_H_ */  

[cpp] view plaincopy在CODE上查看代码片派生到我的代码片
  1. /* 
  2.  * zx_fft.c 
  3.  * 
  4.  * Implementation of Fast Fourier Transform(FFT) 
  5.  * and reversal Fast Fourier Transform(IFFT) 
  6.  * 
  7.  *  Created on: 2013-8-5 
  8.  *      Author: monkeyzx 
  9.  */  
  10.   
  11. #include "zx_fft.h"  
  12. #include   
  13. #include   
  14.   
  15. /* 
  16.  * Bit Reverse 
  17.  * === Input === 
  18.  * x : complex numbers 
  19.  * n : nodes of FFT. @N should be power of 2, that is 2^(*) 
  20.  * l : count by bit of binary format, @l=CEIL{log2(n)} 
  21.  * === Output === 
  22.  * r : results after reversed. 
  23.  * Note: I use a local variable @temp that result @r can be set 
  24.  * to @x and won't overlap. 
  25.  */  
  26. static void BitReverse(complex *x, complex *r, int n, int l)  
  27. {  
  28.     int i = 0;  
  29.     int j = 0;  
  30.     short stk = 0;  
  31.     static complex *temp = 0;  
  32.   
  33.     temp = (complex *)malloc(sizeof(complex) * n);  
  34.     if (!temp) {  
  35.         return;  
  36.     }  
  37.   
  38.     for(i=0; i
  39.         stk = 0;  
  40.         j = 0;  
  41.         do {  
  42.             stk |= (i>>(j++)) & 0x01;  
  43.             if(j
  44.             {  
  45.                 stk <<= 1;  
  46.             }  
  47.         }while(j
  48.   
  49.         if(stk < n) {             /* 满足倒位序输出 */  
  50.             temp[stk] = x[i];  
  51.         }  
  52.     }  
  53.     /* copy @temp to @r */  
  54.     for (i=0; i
  55.         r[i] = temp[i];  
  56.     }  
  57.     free(temp);  
  58. }  
  59.   
  60. /* 
  61.  * FFT Algorithm 
  62.  * === Inputs === 
  63.  * x : complex numbers 
  64.  * N : nodes of FFT. @N should be power of 2, that is 2^(*) 
  65.  * === Output === 
  66.  * the @x contains the result of FFT algorithm, so the original data 
  67.  * in @x is destroyed, please store them before using FFT. 
  68.  */  
  69. int fft(complex *x, int N)  
  70. {  
  71.     int i,j,l,ip;  
  72.     static int M = 0;  
  73.     static int le,le2;