1. 前言
2. 背景知识
这里说明一下公式中各字母的含义:
s(n)代表时域经过采样后(离散化)的信号, n即时域的采样点。N代表时域中总的采样点个数。这里打个比方,假设有一个能量信号,其持续了1秒钟,我们将这1秒钟的能量信号每隔0.1秒采样一次,则总共可以得到11个采样点(0秒的时候也采了一次)那么N=11,n∈(0,10),那么s(5)就可以代表在0.5秒时信号的幅度或者能量,在下面的分析中我们都认为N是一个偶数,并且满足=N,其中k是一个整数(以底为2为条件进行的FFT又叫基2-FFT,当然根据不同需求还有基6-FFT,基8-FFT等等,根据数据量选用相应的FFT会让效率提升许多)。
u可以理解为频域中的采样点,F(u)则为频域采样点所对应的频域信号强度。所以上式可以理解为频域中某点的信号强度是由时域中各点的信号强度乘上一个系数后进行累加而得到。
这里的系数显然指的就是:
这个系数有个别名叫旋转因子,因为其包含了频率谱中的相位信息(后面如果有时间写图像处理内容的话我们会发现一幅图像的轮廓,细节等主要特征信息都包含在相位谱中,幅度谱只携带了各区域的亮度信息)。这个旋转因子一般用W来表示:
3.1 FFT理论
由于旋转因子可以分解成cos与sin的形式,大家知道三角函数是有周期性与对称性的,故经过简单的推导不难得出旋转因子有如下一些特性:
1.周期性:
2.共轭对称性:
3.其他特性:
其中周期性是一个比较重要的特性,从这个特性中我们可以知道如果有办法将总周期长度N缩短为原先的1/2(最后频域采样点数的长度不能变,仍为N)那么后半段的DFT则可以由前半段来表示。
这里我们将DFT序列按照奇偶进行重新排列:
然后根据其他特性中的第三个特性,将奇偶序列的旋转因子上面的因数2放到下面去,并且将奇序列中不含n的旋转因子分离出来:
这样我们只需要计算式(5)中的奇序列和偶序列(2个N/2点的DFT)便可以得到式(6)的结果,计算量基本是减少了一半。这样,如果我们对奇偶序列分别继续做分解,直到只剩N/2个2点DFT,计算量便会大大减少。
3.2 FFT分解与整合
对于即使是刚接触FFT的同学上面的论都不会很难理解,自己推导几遍就能掌握,很多教程往往到这边就会直接给出FFT的蝶形图,很多人刚看到一般会有这样的感觉:蝶形图与上面的公式关联性很强,但实际推导又推不出来,那我们下面就从分解看起,最后再到蝶形图的整合,让大家有一个比较清晰的流程。
3.2.1FFT分解
用例子来说明会便于理解,这边用很有代表性的8点DFT举例,首先要搞清楚一件事,上面的所有理论分析只是FFT的分解,最终我们要反过来将N/2个2点的DFT组合成频域各点的值就是蝶形图了。
所以下面我们先按照上面的理论画一张8点DFT的分解图,这样大家会看的更明白。
Fig.1
上图便是8点DFT的图,由于在时域是8个点所以在频域也是8个点,我们用8个层来代表8个频域点,后面的每一层都与第一层(u=0层)有着相同的结构,然后我们在每一层中将旋转因子与每个时域点相乘也列出来,他们一起求和便是一个频域点的结果(就是每层最上面那个DFT等式,也就是式1), 下面我们根据式(5)将前4层旋转因子与时域点的乘积按照奇偶进行分组 (为什么只分解前4层?因为用前4层就可以计算出后四层的值了,参考式5,式6):
Fig.2
这样就完成了一次分解,每分解一次,层数就会减少一半,注意,在第一次分解中间有一个加减号,这个加减号只能是相同层进行加减,比如分解后第一阶(Stage 1)中上面的u=0层加上下面的u=0层得到的是前一阶(Stage 0) u=0那一层的结果(注意奇数分解得到的结果前面还要乘上一个旋转因子,这里的N等于前面一阶的N,参考式5,式6,这里的u是这一阶所要计算的层的u)。那么Stage 1中上面的u=0层减去下面的u=0层的得到的就是前一阶中u=0+N/2=4那一层的结果,上面的u=1层加上下面的u=1层的得到的就是前一阶中u=1那一层的结果,上面的u=1层减去下面的u=1层的得到的就是前一阶中u=1+N/2=5那一层的结果,依此类推。
那么我们照着步骤继续分解,便会最终得到如下的图:
3.2.1 FFT整合
Fig.6的四幅图合起来就形成如下的形式:
Fig.7
Fig.8
3.3 FFT计算量分析
分析到这里我们不妨看一下FFT与DFT计算量上面的差距,这个网上已经有很多介绍了,我这边也没有什么新见解,只是想整合一下FFT涉及的知识。先看DFT,我们假设有一个N点的时域序列,根据公式(1),这样计算一个频域点就需要N次复数乘法与N-1次复数加法(注意,根据公式2,旋转因子可以分解成一个复数,所以这边计算都是复数。)那么N点的时域序列经过DFT就会变成N点的频域序列,所以计算N个频域点就需要N*N次复数乘法与N*(N-1)次复数加法。
再来看FFT,对于每一个蝴蝶运算基元(见Fig.5)我们有1次复数乘法,2次复数加法,FFT每一级都有N/2个蝴蝶运算基元,共有L=层(N=
)所以FFT共有M次复数乘法与A次复数加法:
所以如果不考虑常数项系数以及常数项,那么FFT与DFT的计算量可以用如下公式表示:
转换成对应于FFT级数,L=,N=
有:
可以看出FFT是根据层数按线性增长,而DFT却是以指数形式的增长,这样数据越多,FFT的优势就越明显。
3.4 总结
到这里我们的FFT算法部分就分析完了,FFT根据旋转因子的一些特性对DFT进行了优化,这样能够让多数据的时频转换能够在计算机上运行起来,之后我准备再写一篇用C/C++来实现FFT的文章,并给出一些建议。
这是本人第一次开博,目前对博客的规划是专注于图像处理的算法以及代码实现,但由于本人的工作性质不固定(什么都搞,大体是嵌入式方向),也有可能中间插一些硬件电路或者Linux内核之类的文章。有人会问,现在有OpenCv等移植性好的图像处理库,也有像Matlab这样强大的模拟软件,为什么还要知道这些算法呢,但我认为要在图像处理方面做的好,光会用那些库或者函数是远远不够的,对于不同的项目,不同的环境,往往不是一种方法或几种方法能解决的,有时候还需要自己的灵感,或者叫经验,这时候你可能需要在原算法上加些东西,这些东西可能是另一个算法,可能是一些参数优化,还有可能是你自己的见解,但不管怎样这都要建立在对原算法了解的基础上,不然你将无从下手。
好了,这篇文章到此结束,欢迎大家讨论。