博客
关于我
强烈建议你试试无所不能的chatGPT,快点击我
FFT算法的完整DSP实现
阅读量:7121 次
发布时间:2019-06-28

本文共 2904 字,大约阅读时间需要 9 分钟。

傅里叶变换或者FFT的理论参考:

[1] 

      The Scientist and Engineer's Guide to Digital Signal Processing,   By Steven W. Smith, Ph.D.

[2] ,可当作[1]的中文参考

[3] 任意一本数字信号处理教材,上面都有详细的推导DCT求解转换为FFT求解的过程

[4] TI文档:基于TMS320C64x+DSP的FFT实现。 使用baidu/google可以搜索到。

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,算法参考上面给的流程图。

 

/* * 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)#endiftypedef 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_ */

 

/* * 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图形验证代码

 

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代码中:

OUTPUT[i] = sqrt(x[i].real*x[i].real + x[i].img*x[i].img)*1024;

 

matlab代码中:

 

subplot(2,2,3); plot(abs(z));title('IFFT');

 

所以IFFT的结果将sin函数的负y轴数据翻到了正y轴。另外,在CCS v5的图形中我们将显示信号的幅度放大了1024倍便于观察,而matlab中没有放大。

 

你可能感兴趣的文章
小程序navigationStyle: custom 下web-view高度异常填坑
查看>>
性能与可伸缩性
查看>>
这一篇就够了——APP瘦身总结
查看>>
写篇文章不难
查看>>
解决TextView上的文字与其他view的间距问题
查看>>
从机器智能得到的启示:对年轻人第一份工作的建议
查看>>
object c 一个改变图标颜色的方法
查看>>
Pinpoint-java性能分析最佳实践_开源PaaS Rainbond
查看>>
Kafka 消息序列化和反序列化(上)
查看>>
箭头函数
查看>>
npm--依赖管理 2
查看>>
设计模式 开闭原则
查看>>
UI设计小白怎样学才能快速入门?
查看>>
用Redis轻松实现秒杀系统
查看>>
埋在MySQL数据库应用中的17个关键问题!
查看>>
APP瘦身这一篇就够了
查看>>
Spark in action on Kubernetes - 存储篇(一)
查看>>
Java springboot B2B2C o2o多用户商城 springcloud架构:服务消费(Feign)
查看>>
15. C#数据结构与算法 -- 栈
查看>>
2015年8月27日【文件权限管理及grep正则和扩展正则表达式】-JY1506402-19+liuhui880818...
查看>>