DSP

FFT算法的完整DSP实现

2019-07-13 11:54发布

站内文章 / DSP
11483 0 1188
文章转载自：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 plaincopyprint?在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_ */  
/* * zx_fft.h * * Created on: 2013-8-5 * Author: monkeyzx */ #ifndef ZX_FFT_H_ #define ZX_FFT_H_ typedef float FFT_TYPE; #ifndef PI #define PI (3.14159265f) #endif typedef struct complex_st { FFT_TYPE real; FFT_TYPE img; } complex; int fft(complex *x, int N); int ifft(complex *x, int N); void zx_fft(void); #endif /* ZX_FFT_H_ */
[cpp] view plaincopyprint?在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;  
  74.     static FFT_TYPE sR,sI,tR,tI,uR,uI;  
  75.   
  76.     M = (int)(log(N) / log(2));  
  77.   
  78.     /* 
  79.      * bit reversal sorting 
  80.      */  
  81.     BitReverse(x,x,N,M);  
  82.   
  83.     /* 
  84.      * For Loops 
  85.      */  
  86.     for (l=1; l<=M; l++) {   /* loop for ceil{log2(N)} */  
  87.         le = (int)pow(2,l);  
  88.         le2 = (int)(le / 2);  
  89.         uR = 1;  
  90.         uI = 0;  
  91.         sR = cos(PI / le2);  
  92.         sI = -sin(PI / le2);  
  93.         for (j=1; j<=le2; j++) {   /* loop for each sub DFT */  
  94.             //jm1 = j - 1;   
  95.             for (i=j-1; i<=N-1; i+=le) {  /* loop for each butterfly */  
  96.                 ip = i + le2;  
  97.                 tR = x[ip].real * uR - x[ip].img * uI;  
  98.                 tI = x[ip].real * uI + x[ip].img * uR;  
  99.                 x[ip].real = x[i].real - tR;  
  100.                 x[ip].img = x[i].img - tI;  
  101.                 x[i].real += tR;  
  102.                 x[i].img += tI;  
  103.             }  /* Next i */  
  104.             tR = uR;  
  105.             uR = tR * sR - uI * sI;  
  106.             uI = tR * sI + uI *sR;  
  107.         } /* Next j */  
  108.     } /* Next l */  
  109.   
  110.     return 0;  
  111. }  
  112.   
  113. /* 
  114.  * Inverse FFT Algorithm 
  115.  * === Inputs === 
  116.  * x : complex numbers 
  117.  * N : nodes of FFT. @N should be power of 2, that is 2^(*) 
  118.  * === Output === 
  119.  * the @x contains the result of FFT algorithm, so the original data 
  120.  * in @x is destroyed, please store them before using FFT. 
  121.  */  
  122. int ifft(complex *x, int N)  
  123. {  
  124.     int k = 0;  
  125.   
  126.     for (k=0; k<=N-1; k++) {  
  127.         x[k].img = -x[k].img;  
  128.     }  
  129.   
  130.     fft(x, N);    /* using FFT */  
  131.   
  132.     for (k=0; k<=N-1; k++) {  
  133.         x[k].real = x[k].real / N;  
  134.         x[k].img = -x[k].img / N;  
  135.     }  
  136.   
  137.     return 0;  
  138. }  
  139.   
  140. /* 
  141.  * Code below is an example of using FFT and IFFT. 
  142.  */  
  143. #define  SAMPLE_NODES              (128)   
  144. complex x[SAMPLE_NODES];  
  145. int INPUT[SAMPLE_NODES];  
  146. int OUTPUT[SAMPLE_NODES];  
  147.   
  148. static void MakeInput()  
  149. {  
  150.     int i;  
  151.   
  152.     for ( i=0;i
  153.     {  
  154.         x[i].real = sin(PI*2*i/SAMPLE_NODES);  
  155.         x[i].img = 0.0f;  
  156.         INPUT[i]=sin(PI*2*i/SAMPLE_NODES)*1024;  
  157.     }  
  158. }  
  159.   
  160. static void MakeOutput()  
  161. {  
  162.     int i;  
  163.   
  164.     for ( i=0;i
  165.     {  
  166.         OUTPUT[i] = sqrt(x[i].real*x[i].real + x[i].img*x[i].img)*1024;  
  167.     }  
  168. }  
  169.   
  170. void zx_fft(void)  
  171. {  
  172.     MakeInput();  
  173.   
  174.     fft(x,128);  
  175.     MakeOutput();  
  176.   
  177.     ifft(x,128);  
  178.     MakeOutput();  
  179. }  
/* * zx_fft.c * * Implementation of Fast Fourier Transform(FFT) * and reversal Fast Fourier Transform(IFFT) * * Created on: 2013-8-5 * Author: monkeyzx */ #include "zx_fft.h" #include #include /* * Bit Reverse * === Input === * x : complex numbers * n : nodes of FFT. @N should be power of 2, that is 2^(*) * l : count by bit of binary format, @l=CEIL{log2(n)} * === Output === * r : results after reversed. * Note: I use a local variable @temp that result @r can be set * to @x and won't overlap. */ static void BitReverse(complex *x, complex *r, int n, int l) { int i = 0; int j = 0; short stk = 0; static complex *temp = 0; temp = (complex *)malloc(sizeof(complex) * n); if (!temp) { return; } for(i=0; i>(j++)) & 0x01; if(j
程序在TMS320C6713上实验，主函数中调用zx_fft()函数即可。 FFT的采样点数为128，输入信号的实数域为正弦信号，虚数域为0，数据精度定义FFT_TYPE为float类型，MakeInput和MakeOutput函数分别用于产生输入数据INPUT和输出数据OUTPUT的函数，便于使用CCS 的Graph功能绘制波形图。这里调试时使用CCS v5中的Tools -> Graph功能得到下面的波形图（怎么用自己琢磨，不会的使用CCS 的Help）。 输入波形
输入信号的频域幅值表示
FFT运算结果
对FFT运算结果逆变换(IFFT)
如何检验运算结果是否正确呢？有几种方法： （1）使用matlab验证，下面为相同情况的matlab图形验证代码 [cpp] view plaincopyprint?在CODE上查看代码片派生到我的代码片
  1. SAMPLE_NODES = 128;  
  2. i = 1:SAMPLE_NODES;  
  3. x = sin(pi*2*i / SAMPLE_NODES);  
  4.   
  5. subplot(2,2,1); plot(x);title('Inputs');  
  6. axis([0 128 -1 1]);  
  7.   
  8. y = fft(x, SAMPLE_NODES);  
  9. subplot(2,2,2); plot(abs(y));title('FFT');  
  10. axis([0 128 0 80]);  
  11.   
  12. z = ifft(y, SAMPLE_NODES);  
  13. subplot(2,2,3); plot(abs(z));title('IFFT');  
  14. axis([0 128 0 1]);  
SAMPLE_NODES = 128; i = 1:SAMPLE_NODES; x = sin(pi*2*i / SAMPLE_NODES); subplot(2,2,1); plot(x);title('Inputs'); axis([0 128 -1 1]); y = fft(x, SAMPLE_NODES); subplot(2,2,2); plot(abs(y));title('FFT'); axis([0 128 0 80]); z = ifft(y, SAMPLE_NODES); subplot(2,2,3); plot(abs(z));title('IFFT'); axis([0 128 0 1]);

（2）使用IFFT验证：输入信号的FFT获得的信号再IFFT，则的到的信号与原信号相同 可能大家发现输入信号上面的最后IFFT的信号似乎不同，这是因为FFT和IFFT存在精度截断误差（也叫数据截断噪声，意思就是说，我们使用的float数据类型数据位数有限，没法完全保留原始信号的信息）。因此，IFFT之后是复数（数据截断噪声引入了虚数域，只不过值很小），所以在绘图时使用了计算幅值的方法， C代码中： [cpp] view plaincopyprint?在CODE上查看代码片派生到我的代码片
  1. OUTPUT[i] = sqrt(x[i].real*x[i].real + x[i].img*x[i].img)*1024;  
OUTPUT[i] = sqrt(x[i].real*x[i].real + x[i].img*x[i].img)*1024; matlab代码中： [cpp] view plaincopyprint?在CODE上查看代码片派生到我的代码片
  1. subplot(2,2,3); plot(abs(z));title('IFFT');  
subplot(2,2,3); plot(abs(z));title('IFFT'); 所以IFFT的结果将sin函数的负y轴数据翻到了正y轴。另外，在CCS v5的图形中我们将显示信号的幅度放大了1024倍便于观察，而matlab中没有放大。
=================
 更正 更正 。。。 ================= 上面程序中的BitReverse函数由于使用了malloc函数，当要分配的n比较大时，在DSP上运行会出现一定的问题，因此改用伪代码中提供的倒位序方法更可靠。 修正后的完整FFT代码文件粘贴如下，在实际的DSP项目中测试通过，可直接拷贝复用。 [cpp] view plaincopyprint?在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. #include "../Headers/Global.h"   
  12.   
  13. #define TYPE_FFT_E     float    /* Type is the same with COMPLEX member */        
  14.   
  15. #ifndef PI   
  16. #define PI             (3.14159265f)   
  17. #endif   
  18.   
  19. typedef COMPLEX TYPE_FFT;  /* Define COMPLEX in Config.h */  
  20.   
  21. extern int fft(TYPE_FFT *x, int N);  
  22. extern int ifft(TYPE_FFT *x, int N);  
  23.   
  24. #endif /* ZX_FFT_H_ */  
/* * zx_fft.h * * Created on: 2013-8-5 * Author: monkeyzx */ #ifndef _ZX_FFT_H #define _ZX_FFT_H #include "../Headers/Global.h" #define TYPE_FFT_E float /* Type is the same with COMPLEX member */ #ifndef PI #define PI (3.14159265f) #endif typedef COMPLEX TYPE_FFT; /* Define COMPLEX in Config.h */ extern int fft(TYPE_FFT *x, int N); extern int ifft(TYPE_FFT *x, int N); #endif /* ZX_FFT_H_ */
[cpp] view plaincopyprint?在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.  * TEST OK 2014.01.14 
  11.  * == 2014.01.14 
  12.  *   Replace @BitReverse(x,x,N,M) by refrence to  
  13.  *    
  14.  */  
  15.   
  16. #include "zx_fft.h"   
  17.   
  18. /* 
  19.  * FFT Algorithm 
  20.  * === Inputs === 
  21.  * x : complex numbers 
  22.  * N : nodes of FFT. @N should be power of 2, that is 2^(*) 
  23.  * === Output === 
  24.  * the @x contains the result of FFT algorithm, so the original data 
  25.  * in @x is destroyed, please store them before using FFT. 
  26.  */  
  27. int fft(TYPE_FFT *x, int N)  
  28. {  
  29.     int i,j,l,k,ip;  
  30.     static int M = 0;  
  31.     static int le,le2;  
  32.     static TYPE_FFT_E sR,sI,tR,tI,uR,uI;  
  33.   
  34.     M = (int)(log(N) / log(2));  
  35.   
  36.     /* 
  37.      * bit reversal sorting 
  38.      */  
  39.     l = N / 2;  
  40.     j = l;  
  41.     //BitReverse(x,x,N,M);   
  42.     for (i=1; i<=N-2; i++) {  
  43.         if (i < j) {  
  44.             tR = x[j].real;  
  45.             tI = x[j].imag;  
  46.             x[j].real = x[i].real;  
  47.             x[j].imag = x[i].imag;  
  48.             x[i].real = tR;  
  49.             x[i].imag = tI;  
  50.         }  
  51.         k = l;  
  52.         while (k <= j) {  
  53.             j = j - k;  
  54.             k = k / 2;  
  55.         }  
  56.         j = j + k;  
  57.     }  
  58.   
  59.     /* 
  60.      * For Loops 
  61.      */  
  62.     for (l=1; l<=M; l++) {   /* loop for ceil{log2(N)} */  
  63.         le = (int)pow(2,l);  
  64.         le2 = (int)(le / 2);  
  65.         uR = 1;  
  66.         uI = 0;  
  67.         sR = cos(PI / le2);  
  68.         sI = -sin(PI / le2);  
  69.         for (j=1; j<=le2; j++) {   /* loop for each sub DFT */  
  70.             //jm1 = j - 1;   
  71.             for (i=j-1; i<=N-1; i+=le) {  /* loop for each butterfly */  
  72.                 ip = i + le2;  
  73.                 tR = x[ip].real * uR - x[ip].imag * uI;  
  74.                 tI = x[ip].real * uI + x[ip].imag * uR;  
  75.                 x[ip].real = x[i].real - tR;  
  76.                 x[ip].imag = x[i].imag - tI;  
  77.                 x[i].real += tR;  
  78.                 x[i].imag += tI;  
  79.             }  /* Next i */  
  80.             tR = uR;  
  81.             uR = tR * sR - uI * sI;  
  82.             uI = tR * sI + uI *sR;  
  83.         } /* Next j */  
  84.     } /* Next l */  
  85.   
  86.     return 0;  
  87. }  
  88.   
  89. /* 
  90.  * Inverse FFT Algorithm 
  91.  * === Inputs === 
  92.  * x : complex numbers 
  93.  * N : nodes of FFT. @N should be power of 2, that is 2^(*) 
  94.  * === Output === 
  95.  * the @x contains the result of FFT algorithm, so the original data 
  96.  * in @x is destroyed, please store them before using FFT. 
  97.  */  
  98. int ifft(TYPE_FFT *x, int N)  
  99. {  
  100.     int k = 0;  
  101.   
  102.     for (k=0; k<=N-1; k++) {  
  103.         x[k].imag = -x[k].imag;  
  104.     }  
  105.   
  106.     fft(x, N);    /* using FFT */  
  107.   
  108.     for (k=0; k<=N-1; k++) {  
  109.         x[k].real = x[k].real / N;  
  110.         x[k].imag = -x[k].imag / N;  
  111.     }  
  112.   
  113.     return 0;  
  114. }  
/* * zx_fft.c * * Implementation of Fast Fourier Transform(FFT) * and reversal Fast Fourier Transform(IFFT) * * Created on: 2013-8-5 * Author: monkeyzx * * TEST OK 2014.01.14 * == 2014.01.14 * Replace @BitReverse(x,x,N,M) by refrence to * */ #include "zx_fft.h" /* * FFT Algorithm * === Inputs === * x : complex numbers * N : nodes of FFT. @N should be power of 2, that is 2^(*) * === Output === * the @x contains the result of FFT algorithm, so the original data * in @x is destroyed, please store them before using FFT. */ int fft(TYPE_FFT *x, int N) { int i,j,l,k,ip; static int M = 0; static int le,le2; static TYPE_FFT_E sR,sI,tR,tI,uR,uI; M = (int)(log(N) / log(2)); /* * bit reversal sorting */ l = N / 2; j = l; //BitReverse(x,x,N,M); for (i=1; i<=N-2; i++) { if (i < j) { tR = x[j].real; tI = x[j].imag; x[j].real = x[i].real; x[j].imag = x[i].imag; x[i].real = tR; x[i].imag = tI; } k = l; while (k <= j) { j = j - k; k = k / 2; } j = j + k; } /* * For Loops */ for (l=1; l<=M; l++) { /* loop for ceil{log2(N)} */ le = (int)pow(2,l); le2 = (int)(le / 2); uR = 1; uI = 0; sR = cos(PI / le2); sI = -sin(PI / le2); for (j=1; j<=le2; j++) { /* loop for each sub DFT */ //jm1 = j - 1; for (i=j-1; i<=N-1; i+=le) { /* loop for each butterfly */ ip = i + le2; tR = x[ip].real * uR - x[ip].imag * uI; tI = x[ip].real * uI + x[ip].imag * uR; x[ip].real = x[i].real - tR; x[ip].imag = x[i].imag - tI; x[i].real += tR; x[i].imag += tI; } /* Next i */ tR = uR; uR = tR * sR - uI * sI; uI = tR * sI + uI *sR; } /* Next j */ } /* Next l */ return 0; } /* * Inverse FFT Algorithm * === Inputs === * x : complex numbers * N : nodes of FFT. @N should be power of 2, that is 2^(*) * === Output === * the @x contains the result of FFT algorithm, so the original data * in @x is destroyed, please store them

热门文章