FFT详解

1. 前言

        由于本人最近在做一些图像处理的项目,在看到频域分析的时候发现FFT是一个难点,但其又是频域分析里面基础中的基础,在网上找了一些资料后发现大多数文章分析的不是很透彻,故准备写一篇文章,一是为了传播一下相关知识,第二也是为了巩固一下自己的基础。当然此文章的描述也可能有错误,欢迎大家指正。

2. 背景知识

        FFT(Fast Fourier Transformation)故名意思,是傅里叶变换的快速计算方法,这里的傅里叶变换特指离散域的傅里叶变换(DFT),这篇文章假设读者已经对DFT有了一定的了解,故不会再叙述DFT的原理,推导以及含义(涉及DFT的相关书籍有很多,DSP,语音信号处理,图像处理都会有介绍,建议不懂的读者去看一下相关书籍),下面只给出DFT计算公式:

      

       这里说明一下公式中各字母的含义:

       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的旋转因子分离出来:

       写成这样的形式,傅里叶变换的周期性便一目了然了,由于旋转因子是以N/2为底,所以可以认为旋转因子的周期是N/2,然而由于频域分量u的范围还是从0~N,所以u的范围从0~N/2-1与u从N/2~N-1对于周期为N/2的旋转因子来说结果是一样的。再说, 因为是以N为底,所以其周期为N,根据其他特性中的第二个特性,我们知道其前半段的值只需乘上一个-1便可以得到后半段的值,即:

        在这样的思想下,我们固然会想到将F(u)一分为二,因为这样只需要计算前半段的旋转因子即可:

         这样我们只需要计算式(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那一层的结果,依此类推。

        那么我们照着步骤继续分解,便会最终得到如下的图:


Fig.3
        这样分解就完成了,从上图可以发现每分解一次层数便少掉一半,后面我们要做的工作就是从上图的最右边开始往左进行整合,这样整合起来的图就是蝶形图了,大家有兴趣可以自己先整合一下,我在下一节中会给出整合的步骤。

3.2.1 FFT整合

        到这一部分网上介绍的就很多了,我们这边按照Fig.3再画一遍,如果读者看了这篇上面的内容相信画出蝶形图是没有问题了。从Fig.3的最右边开始,这里注意一下时域的序列经过分解顺序已经给重排了,重排的顺序是(0-4),(2-6),(1-5),(3-7)。然后由于分解到最后u只等于0了,所以时域前面的旋转因子都等于1了。然后Stage 3每两个方块相加就得到前一层(Stage 2)的u=0那一层,相减就得到Stage 2的u=1那一层,注意奇数的块要乘上一个旋转系数,由于第三层u=0,所以也可以写成,当然很多教程为了让旋转因子的基数N与原始序列长度保持一致,都会将其写成(对于Stage 3就是)形式所以Stage 2 到Stage 3的蝴蝶图就可以写成如下形式:
Fig.4
         在上图中每一对都是一个蝴蝶运算基元,表示成:
Fig.5
         其中上面一行是做加法运算:s1+s2,下面一行做的是减法运算:s1-W*s2。
         接着我们看Stage 2的上面两个方块的u=0层,按照原理,Stage 2的两个u=0层相加(也就是 s'0+s'2+*(s'4+s'6),这里的s'x指的是每个方块中时域序列与前面的旋转系数相乘的结果,对于Stage 2,s'0= ,参考Fig.3 Stage 2上面的方块,s'0+s'4也即是从前一层Stage 3计算得到的结果 )得到的是Stage 1上方方块的u=0层(Fig.6 a)。 Stage 2的两个u=0层相减( s'0+s'2-(s'4+s'6) )得到的是Stage 1上方方块的u=2层 (注意是u=2层,不是u=1层,因为在Stage 1中N=4相减得到后一半的层),见(Fig.6 b),由于u=0所以s'4+s'6前面的旋转系数就可以不要了。Stage 2的两个u=1层相加( s'0+s'2+(s'4+s'6) )得到的是Stage 1 u=1层的结果,见(Fig.6 c)。Stage 2的两个u=1层相减( s'0+s'2-(s'4+s'6) )得到的是Stage 1 u=3层的结果,见(Fig.6 d)。


Fig.6 a,b,c,d


        Fig.6的四幅图合起来就形成如下的形式:

Fig.7

       之后的步骤都是按照Fig.6的模式来的,相信大家看懂Fig.6的图,之后的蝴蝶图应该可以很容易的就画出来,这边不再赘述,直接给出最终的图:

Fig.8

        这样一个8点的FFT蝴蝶图便完成了,相信大家看了这篇文章之后会对FFT有更进一步的了解。


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这样强大的模拟软件,为什么还要知道这些算法呢,但我认为要在图像处理方面做的好,光会用那些库或者函数是远远不够的,对于不同的项目,不同的环境,往往不是一种方法或几种方法能解决的,有时候还需要自己的灵感,或者叫经验,这时候你可能需要在原算法上加些东西,这些东西可能是另一个算法,可能是一些参数优化,还有可能是你自己的见解,但不管怎样这都要建立在对原算法了解的基础上,不然你将无从下手。

       好了,这篇文章到此结束,欢迎大家讨论。






      



评论 4
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值