相信网上现在有很多关于FFT的教程,我曾经也参阅了很多网上的教程,感觉都不怎么通俗易懂。在基本上的研究FFT,并且通过编程的形式实现之后。我决定写一篇通俗易懂的关于FFT的讲解。因此我在接下来的叙述中尽量非常通俗细致的讲解。 本人最早知道傅里叶变换的时候是沉迷于音乐的频谱跳动无法自拔,当时就很想做一个音乐频谱显示器。搜阅了很多资料之后,才了解到傅里叶变换,和FFT。当然这都是以前的事情了,经过了系统的学习+2个星期的研究,自制了一个FFT的算法,不敢说速度上的优势,但是个人认为是一种通俗易懂的实现方法。经过实际的VC++模拟实验、和STM32跑的也很成功。 首先,要会FFT,就必须要对DFT有所了解,因为两者之间本质上是一样的。在此之前,先列出离散傅里叶变换对(DFT): 其中: 但是FFT之所以称之为快速傅里叶变换,就利用了以下的几个性质(重中之重!) 周期性: 对称性: 可约性: 先把这仨公式放到这,接下来会用到。 根据这几个特性,就可以将一个长的DFT运算分解为若干短序列的DFT运算的组合,从而减少运算量。在这里,为了方便理解,我就用了一个按时间抽取的快速傅里叶变换(DIT-FFT)的方法。 首先,将一个序列x(n)一分为二 对于 将x(n)按照奇偶分组 x(2r)=x1(r) x(2r+1)=x2(r) r=0,1,…..(N/2)-1 那么将DFT也分为两组来预算 这个时候,我们上边提的三个性质中的可约性就起到作用了: 回顾一下: 那么这个式子就可以化为: 则X1(k)和X2(k)都是长度为N/2的序列 x1(k)和x2(k)的N/2点的离散傅里叶变换 其中: K=0,1,2…N/2-1 至此,一个N点的DFT就被分解为2个N/2的DFT。但是X1(k),和X2(k)只有N/2个点,也就是说式子(1-1)只是DFT前半部分。要求DFT的后半部分可以利用其对称性求出后半部分为: (式1-2) 那么式(1-1)和(1-2)就可以专用一个蝶形信号流图符号来表示。如图: 以N=8为例,可以用下图表示: 通过这样的分解,每一个N/2点DFT只需(N^2)/4次复数相乘计算,明显的节省了运算量。 以此类推,继续将已经得出的X1(k)和X2(k)按照奇偶分解,还是和上边一样的方法。那么就得出了百度上都可以找到的一大推的这个图片了。 (笑) 对于这张图片我想强调的一点就是第二阶蝶形运算的时候,WN0之后不应该是WN1吗,为什么是W2N了,这个问题之前困扰了我一段时间,但是不要忘了,前者的 W0N的展开是W0N/2 因为此时N已经按照奇偶分开了,所以是N/2 而W2N/2也就是 W2N是根据可约性得出的,在这里不能混淆.。 对于运算效率就不用多提了 以上就是FFT算法的理论内容了,接下来就是用C语言对这个算法的实现了,对于FFT算法C语言的实现,网上的方法层出不穷,介于本人比较懒(懒得看别人的程序),再加上自给自足丰衣足食的原则,我自己也写了一个个人认为比较通俗易懂的程序,并且为了帮助读者理解,我特意尽量减少了库函数的使用,一些基本的函数都是自己写的(难免有很多BUG),但是作为FFT算法已经够用了,目前这个程序只能处理2^N的数据,理论上来讲如果不够2^N的话可以在后面数列补0这种操作为了简约我也就没有实现,但是主要的功能是具备的,下面是代码:
这个程序中input这个数组是输入信号,在这里只模拟抽样了8次,输出的数据也是input,如果想看其它序列的话,可以改变FFT_LENGTH的值以及 input里的内容,程序输出的是实部和虚部的模,如果单纯想看实部或者虚部的幅度的话,请自行修改程序~ (本文摘自21ic论坛,感谢原作者辛苦分享) |
|
来自: 西北望msm66g9f > 《培训》