星巴克市场营销观念:快速傅立叶变换(FFT)

来源:百度文库 编辑:中财网 时间:2024/03/29 02:36:39

 快速傅立叶变换(FFT)

 

41引言

快速傅立叶变换(FFT)并不是一种新的变换,而是离散傅立叶变换(DFT)的一种快速算法。

DFT的计算在数字信号处理中非常有用。例如在FIR滤波器设计中会遇到从h(n)求H(k)或由H(k)计算h(n),这就要计算DFT;信号的谱分析对通信、图像传输、雷达等都是很重要的,也要计算DFT。因直接计算DFT的计算量与变换区间长度N的平方成正比,当N较大时,计算量太大。

自从1965年图基(J. W. Tukey)和库利(T. W. Coody)在《计算数学》(Math. Computation , Vol. 19, 1965)杂志上发表了著名的《机器计算傅立叶级数的一种算法》论文后,桑德(G. Sand)-图基等快速算法相继出现,又经人们进行改进,很快形成一套高效运算方法,这就是快速傅立叶变换简称FFT(Fast Fourier Transform)。这种算法使DFT的运算效率提高1~2个数量级。

 

42 2 FFT算法

一、直接计算DFT的问题及改进的途径

设x(n)为N点有限长序列,其DFT正变换为

= , k=0,1,…,N-1

其反变换(IDFT)

x(n)= ,n=0,1,…,N-1

二者的差别只在于 的指数符号不同,以及差一个常数乘因子1/N,因而下面我们只讨论DFT正变换的运算量,反变换的运算量是完全相同的。

考虑x(n)为复数序列的一般情况,每计算一个X(k),需要N次复数乘法以及(N-1)次复数加法。因此,对所有N个k值,共需N2次复数乘法及
N(N-1)次复数加法运算。所以直接计算DFT,乘法次数和加法次数都是和N2成正比的,当N很大时,运算量是很可观的,因而需要改进对DFT的计算方法,以减少运算次数。

 

 

       下面讨论减少运算工作量的途径。仔细观察DFT的运算就可看出,利用系数 以下固有特性,就可减小DFT的运算量:

(1) 的对称性 ( *=

(2) 的周期性 = =

(3) 的可约性 = =

由此可得: = = =-1, =- 。这样利用这些特性,可以将长序列的DFT分解为短序列的DFT。快速傅立叶变换算法正是基于这样的思路而发展起来的。它的算法基本上可以分成两大类:时域抽取法FFT(Decimation-In-Time FFT,简称DIT-FFT)和频域抽取法FFT(Decimation-In-Time FFT,简称DIF-FFT)。

 

 

二、时域抽取法基-2 FFT原理

先设序列点数为N=2M,M为整数。如果不满足这个条件,可以人为地加上若干零值点,使之达到这一要求。这种N为2的整数幂的FFT称基-2 FFT。

 

设输入序列长度为N=2M (M为正整数) ,将该序列按时间顺序的奇偶分解为越来越短的子序列,称为按时间抽取(DIT )的FFT算法。也称Cooley - Tukey算法。

将N=2M的序列x(n)先按n的奇偶分成以下两组:

x(2r)=x1(r)      r=0,1,…,N/2-1

x(2r+1)=x2(r)

则x(n)的DFT为:

=DFT[x(n)]= ,k=0,1,…,N-1

 = +

 = +

 = +

= +

 = + ,k=0,1,…,N-1 (4.2.4)

 

式中 分别是x1(r)及x2(r)的N/2点DFT:

      = = ,k=0,1,…,N/2-1

= = ,k=0,1,…,N/2-1

由(4.2.4)式可看出,一个N点DFT已分解成两个N/2点的DFT,它们按(4.2.4)式又组合成一个N点DFT。

 

现讨论 的周期性:

= = =

式中用了 =

 

同理可得: =

再考虑 的以下性质:

       = =-

所以前半部分:

   = + ,k=0,1,…,N/2-1

后半部分:

      = +

= ,(k=0,1,…,N/2-1)

这样,只要求出0到(N/2-1)区间的所有 值,即可求出0到(N-1)区间内的所有 值,这就大大节省了运算。

 

= + ,k=0,1,…,     (4.2.7)

= ,k=0,1,…, (4.2.8)

                     

采用蝶形运算符号表示的图示法,可将上面讨论的分解过程表示于下图中。此图表示N=23=8的情况,其中输出值X(0)到X(3)是由(4.2.7)式给出的,而输出值X(4)到X(7)是由(4.2.8)给出。

 

 

 

既然如此,由于N=2M,因而N/2仍是偶数,可以进一步把每个N/2点子序列再按其奇偶部分分解为两个N/4点的子序列。

     x1(2r) = x3(r)      r=0,1,…,N/4-1

x1(2r+1)=x4(r)

 

 

= +

    = +

    = + ,k=0,1,…,N/4-1

= - ,k=0,1,…,N/4-1

同理, 也可进行同样的分解:

      = + ,k=0,1,…,N/4-1

= - ,k=0,1,…,N/4-1

其中, = =

           = =

 

 

 

进一步具体化:N=8=23

     =

         =x(0)+x(4) ,k=0,1

    = x(0)+ x(4)

         = x(0)+ x(4) = x(0)- x(4)

注意上式中 =

同理, = x(1)+ x(5)

           = x(1)- x(5)

 

 

 

三、DIT-FFT算法与直接计算DFT运算量的比较

由上图得,当N=2M时,其运算流图有M级蝶形,每一级都有N/2个蝶形运算构成。每一个蝶形运算需要1次复数乘和2次复数加。所以每一级运算都需要N/2次复数乘和N次复数加。

M级运算总共需要的复数乘法次数为:

     

M级运算总共需要的复数加法次数为:

     

如直接计算DFT,复数乘法为 次,复数加法 次。

当N=1024时,复数乘之比: ,这样使运算效率提高200多倍。下图为FFT算法和直接DFT算法所需运算量与计算点数N的关系曲线:

显然,N越大,优越性就越明显。

 

 

四、按时间抽选的FFT算法的特点:

1、        原位运算

由图4.2.4可以看出,DIT-FFT的运算过程很有规律。N=2M点的FFT共进行M级运算,每级由N/2个蝶形运算组成。同一级中,每个蝶形的两个输入数据只对计算本蝶形有用,而且每个蝶形的输入、输出数据结点又同在一条水平线上,这就意味着计算完一个蝶形后,所得输出数据可立即存入原输入数据所占用的存储单元。这样,经过M级运算后,原来存放输入序列数据的N个存储单元中便依次存放 的N个值。这种利用同一存储单元存储蝶形计算输入、输出数据的方法称为原位计算。原位计算可节省大量内存,从而使设备成本降低。

 

 

2、        倒序规律

由图4.2.4看出,按原位计算时,FFT的输出 是按正常顺序排列在存储单元中,即按 ,…, 的顺序排列,但是这时输入x(n)却不是按自然顺序存储的,而是按x(0),x(4),…,x(7)的顺序存入存储单元,看起来好象是“混乱无序”的,实际上是有规律的,我们称之为倒序。

        造成倒序的原因是输入x(n)按标号n的偶奇的不断分组而造成。由于N=2M,所以倒序数可用M位二进制数(nM-1nM-2…n02(当N=8=23时,二进制为三位)表示。第一次分组,标为n0。 n为偶数在上半部分,用n0=0表示,n为奇数在下半部分,用n0=1表示;第二次分组,标为n1。偶数部分再分为偶(0)奇(1),奇数部分再分为偶(0)奇(1)…。依次类推,直到M次分组,最后所得二进制倒序数如图示。

 

 

 

下表列出了N=8时以二进制数表示的顺序数和倒序数,由表显而易见,只要将顺序数(n2n1n0)的二进制位倒置,则得对应的二进制倒序值(n0n1n2)。

 

 

 

 

 

 

3、        倒序的实现

设原输入序列x(n)先按自然顺序存入数组A中。例如N=8,A(0),A(1),A(2),A(3),A(4),A(5),A(6),A(7)中依次存放着x(0),x(1),x(2),x(3),x(4),x(5),x(6),x(7)。

 

 

 

 

 

顺序数用I表示,I=1~N-2。倒序数用J表示,与I对应分别为4,2,6,1,5,3。当I=J时不需要交换,I

 

    I=1时,对应的倒序数是4;I=2时,对应的倒序数是2…。倒序数从4到2到…的关系可从表4.2.1得到:每次最高位加1。(注意J用十进制数表示)。如果最高位为0,J直接加N/2,如果最高位为1,则要将最高位归0,次高位加1。但次高位加1时也要判断是否为1或0。程序框图如下图虚线框里所示:

 

 

 

 

 

 

 

 

 

 

 

4、        蝶形运算两个输入数据的“距离”

以图4.2.4的8点FFT为例,其输入是倒位序的,输出是自然顺序的。N=2M,共有第一级蝶形运算,第二级蝶形运算,…,第M级蝶形运算。用L表示第某级运算,每个蝶形的两个输入数据的“距离”为B=2L-1

 

5、        旋转因子的变化规律

仍观察图4.2.4,每级都有N/2个蝶形,每个蝶形都要乘以因子 ,称其为旋转因子,p称为旋转因子的指数。第L级蝶形运算(从左向右数),共有2L-1

L=1,第一级: = ,J=0

L=2,第二级: = ,J=0,1

L=3,第三级: = ,J=0,1,2,3

所以对于N=2M的一般情况,第L级的 为:

,J=0,1,…,2L-1-1

由于2L=2M2L-M=N2L-M ,所以 = = =

所以,第L级的 为: ,J=0,1,…,2L-1-1

 

例如:

L=1,第一级:2L-1=1,J=0,      

L=2,第二级:2L-1=2,J=0,1,   

L=3,第三级:2L-1=4,J=0,1,2,3,

 

6、蝶形运算规律

    设序列x(n)经时域抽选(倒序)后,存入数组X中。如果蝶形运算的两个输入数据相距B个点,应用原位计算,则蝶形运算可表示成如下形式:

   

   

式中 ;J=0,1,…,2L-1-1;

L=1,…,M。

 

DIT-FFT运算和程序框图如下:

         

同一旋转因子对应着间隔为2L点的2M-L个蝶形。

 

五、按频率抽选(DIF)的基2 FFT算法

 

设序列点数为N=2M ,M为整数。

,k=0,1,…,N-1

先把输入按n 的顺序分成前后两半:

,k=0,1,…,N-1

= +

    = +

    = ,k=0,1,…,N-1

= ,k=0,1,…,N-1

=(-1)k=   1,k为偶数

                  -1,k为奇数

分解成偶数组与奇数组,

当k取偶数即k=2r时(r=0,1,…,N/2-1),

= =

 

当k取奇数即k=2r+1时(r=0,1,…,N/2-1),

=

=

        n=0,1,…,N/2-1

 

= ,r=0,1,…,N/2-1

=         (4.2.16)

 

之间可用下图所示的蝶形运算符号表示:

 

(4.2.16)用下图表示,N=8。一次分解流图。

 

 

由于N=2M,N/2仍是偶数,继续将N/2点DFT分成偶数组和奇数组。图4.2.12表示N=8时二次分解运算流图。

  

 

最后完整的分解流图为下图:

  

 

    这种算法是对 进行奇偶抽取分解的结果,所以称之为频域抽取法FFT。

    DIF-FFT算法与DIT-FFT算法类似,不同的是DIF-FFT算法输入序列为自然顺序,而输出为倒序排列。另外,蝶形运算略有不同,DIT-FFT蝶形先乘后加,而DIF-FFT蝶形先加后乘。

    上述两种FFT的算法流图形式不是唯一的。只要保证各节点所连支路及其传输系数不变,改变输入与输出点以及中间结点的排列顺序,就可以得到其他变形的FFT运算流图。

 

六、IDFT的高效算法

   离散傅立叶反变换

x(n)= ,n=0,1,…,N-1

与离散傅立叶正变换( = , k=0,1,…,N-1)相比,只要将DFT中的系数 改变为 ,最后乘以1/N,就是IDFT的运算公式。流图输入为 ,输出为x(n)。因此原来的DIT-FFT改为IFFT后称为DIF-IFFT更合适;原来的DIF-FFT改为IFFT后称为DIT-IFFT更合适。

下图是由DIF-FFT运算流图改成的IFFT运算流图(DIT-IFFT):

    在实际中,有时为了防止运算过程中发生溢出,将1/N分配到每一级蝶形运算中。由于1/N= ,所以每级的每个蝶形输出支路均有一相乘因子1/2。如下图示:

  

如果希望直接调用FFT子程序计算IFFT,则可用下面的方法:

由于    x(n)=

∴      x*(n)=

x(n)=

   = {DFT[ ]}*

这样可以先将 取共轭,然后直接调用FFT子程序进行DFT运算,最后取共轭并乘以1/N得到序列x(n)。

43 进一步减少运算量的措施

研究进一步减少运算量的途径,以程序的复杂度换取计算量的进一步提高

一、   多类蝶形单元运算

 

由图4.2.4已得出结论,N=2M点FFT共需要MN/2次复数乘法。

= ,J=0,1,…,2L-1-1得,当L=1时,只有一种旋转因子 ,所以第一级不需要乘法运算。当L=2时,共有两个旋转因子: =1和 =-j,因此第二级也不需要乘法运算。在DFT中,称其值为±1和±j的旋转因子为无关紧要的旋转因子。

    综上所述,先除去第一、第二两级后,所需复数乘法次数

进一步考虑各级中的无关紧要旋转因子。当L=3时,有两个无关紧要的旋转因子 。因为同一旋转因子对应着2M-L=N/2L个蝶形运算,所以,第三级共有2N/23=N/4个蝶形不需要复数乘法运算。依次类推,当L≥3时,第L级的2个无关紧要的旋转因子减少复数乘法的次数为2N/2L=N/2L-1。这样,L=3至L=M共减少复数乘法次数为

因此

        =

 

在基2 FFT程序中,若包含了所有旋转因子,则称该算法为一类蝶形单元运算;若去掉 =±1的旋转因子,则称之为二类蝶形单元运算;若再去掉 =±j的旋转因子,则称为三类蝶形单元运算;若再处理 = ,则称之为四类蝶形运算。我们将后三种运算称为多类蝶形单元运算。显然蝶形单元越多,编程就越复杂,但当N较大时,乘法运算的减少量是相当可观的。例如,N=4096时,三类蝶形单元运算的乘法次数为一类蝶形单元运算的75%。

 

二、   旋转因子的生成

在FFT运算中,旋转因子

= ,求余弦和正弦函数值的计算量很大,所以编程时,一种方法是在每级运算中直接产生,另一种方法是在FFT程序开始前预先计算好,存放在数组中,作为旋转因子表,在程序执行过程中直接查表得到。

 

三、   实序列的FFT算法

在实际工作中,数据x(n)一般都是实序列。如果直接按FFT运算流图计算,就是把x(n)看成一个虚部为零的复序列进行计算,这就增加了运算时间。处理这个问题有二种方法,一种是早期提出的用一个N点FFT计算N点实序列的FFT。第二种方法是用N/2点FFT计算一个N点实序列的DFT。