程序博客网 > 社交网络电影观后感

FFT算法的完整DSP实现

来源：互联网 发布：社交网络电影观后感 编辑：程序博客网 时间：2024/04/29 07:10

【原文：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 <math.h>  
  13. #include <stdlib.h>  
  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<n; i++) {  
  39.         stk = 0;  
  40.         j = 0;  
  41.         do {  
  42.             stk |= (i>>(j++)) & 0x01;  
  43.             if(j<l)  
  44.             {  
  45.                 stk <<= 1;  
  46.             }  
  47.         }while(j<l);  
  48.   
  49.         if(stk < n) {             /* 满足倒位序输出 */  
  50.             temp[stk] = x[i];  
  51.         }  
  52.     }  
  53.     /* copy @temp to @r */  
  54.     for (i=0; i<n; 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<SAMPLE_NODES;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<SAMPLE_NODES;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. }  

程序在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 plaincopy在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]);  


（2）使用IFFT验证：输入信号的FFT获得的信号再IFFT，则的到的信号与原信号相同

可能大家发现输入信号上面的最后IFFT的信号似乎不同，这是因为FFT和IFFT存在精度截断误差（也叫数据截断噪声，意思就是说，我们使用的float数据类型数据位数有限，没法完全保留原始信号的信息）。因此，IFFT之后是复数（数据截断噪声引入了虚数域，只不过值很小），所以在绘图时使用了计算幅值的方法，

C代码中：

[cpp] view plaincopy在CODE上查看代码片派生到我的代码片
  1. OUTPUT[i] = sqrt(x[i].real*x[i].real + x[i].img*x[i].img)*1024;  

matlab代码中：

[cpp] view plaincopy在CODE上查看代码片派生到我的代码片
  1. 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 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. #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_ */  

[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.  * TEST OK 2014.01.14 
  11.  * == 2014.01.14 
  12.  *   Replace @BitReverse(x,x,N,M) by refrence to  
  13.  *   <The Scientist and Engineer's Guide to Digital Signal Processing> 
  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. }  

另外，可能还需要您在其它头文件中定义COMPLEX的复数类型

[cpp] view plaincopy在CODE上查看代码片派生到我的代码片
  1. typedef struct {  
  2.     float real;  
  3.     float imag;  
  4. } COMPLEX;  

=====================

增加：FFT频谱结果显示

=====================

[plain] view plaincopy在CODE上查看代码片派生到我的代码片
  1. clc;  
  2. clear;  
  3.   
  4. % Read a real audio signal  
  5. [fname,pname]=uigetfile(...  
  6.     {'*.wav';'*.*'},...  
  7.     'Input wav File');  
  8. [x fs nbits] = wavread(fullfile(pname,fname));  
  9.   
  10. % Window  
  11. % N = 1024;  
  12. N = size(x,1);  
  13. x = x(1:N, 1);  
  14.   
  15. % FFT  
  16. y = fft(x);  
  17. % 频率对称，取N/2  
  18. y = y(1:N/2);  
  19.   
  20. % FFT频率与物理频率转换  
  21. x1 = 1:N/2;  
  22. x2 = (1:N/2)*fs/(N/2);  % /1000表示KHz  
  23. log_x2 = log2(x2);  
  24.   
  25. % plot  
  26. figure,  
  27. subplot(2,2,1);plot(x);  
  28. xlabel('Time/N');ylabel('Mag');grid on  
  29. title('原始信号');  
  30. subplot(2,2,2);plot(x1, abs(y));  
  31. xlabel('Freq/N');ylabel('Mag');grid on  
  32. title('FFT信号/横坐标为点数');  
  33. subplot(2,2,3);plot(x2,abs(y));  
  34. xlabel('Freq/Hz');ylabel('Mag');grid on  
  35. title('FFT信号/横坐标为物理频率');  
  36. subplot(2,2,4);plot(log_x2,abs(y));  
  37. xlabel('Freq/log2(Hz)');ylabel('Mag');grid on  
  38. title('FFT信号/横坐标为物理频率取log');  
  39.   
  40. % 更常见是将幅值取log  
  41. y = log2(abs(y));  
  42. figure,  
  43. plot(x2,y);  
  44. xlabel('Freq/Hz');ylabel('Mag/log2');grid on  
  45. title('幅值取log2');  
0 0
社交网络电影观后感 社交网络电影观后感
原创粉丝点击
热门IT博客
ios的福利软件ios的福利软件
张佳玮 詹姆斯 知乎张佳玮 詹姆斯 知乎
网桥mac地址表
如何选购床 知乎
散兵坑公式源码
淘宝的链接怎么发微信
心酸的故事 知乎
c语言窗口程序
淘宝买烟怎么搜索2017
网络贷款有多少人不还
windows dpi 修改工具
pdg转pdf mac
狗万万博软件
windows资源管理器重启
java 99乘法表
决战武林坐骑晋级数据
慢性前列腺炎 知乎
淘宝双11宣传
淘宝商家电话怎么查
阿里云邮 免费邮箱
热门问题 老师的惩罚 人脸识别 我在镇武司摸鱼那些年 重生之率土为王 我在大康的咸鱼生活 盘龙之生命进化 天生仙种 凡人之先天五行 春回大明朝 姑娘不必设防，我是瞎子 巡航导殚 陆地巡航 自适应巡航害死人 全速自适应巡航 自适应巡航有什么弊端 全速自适应巡航系统 自适应巡航有用吗 趋势巡航的博客 自动巡航怎么用 春风650巡航版 巡航什么意思 巡航是什么意思 巡航控制指示灯亮怎么回事 豪爵250巡航太子2018年产 自适应巡航车型 巡航导殚弹道道弹区别 r&b巡茶 巡视员 巡视通 巡视员级别 一级巡视员 巡视组 巡视 巡视条例 二级巡视员 省委巡视组 省委巡视 副巡视员 巡视组电话 一级巡视员相当于什么级别 巡视工作条例 2018年各省巡视组电话 巡视员什么级别 2019省巡视组时间安排 巡视组主要查什么 中央巡视组举报电话 2019全国巡视组投诉方式 2019年各省巡视组电话 中央巡视组电话 2019全国巡视组投诉电话 康军被巡视组带走

程序博客网，程序员的互联网技术博客家园。csdn论坛精品 msdn技术资料都在这里