DSP 第4章_0905

1、第4章 快速傅里叶变换(FFT) 第第4章章 快速傅里叶变换快速傅里叶变换(FFT) 4.1 引言引言 4.2 基基2FFT算法算法 4.3 进一步减少运算量的措施进一步减少运算量的措施 4.4 其他快速算法简介其他快速算法简介习题与上机题习题与上机题第4章 快速傅里叶变换(FFT) 4.1 引引 言言 在快速傅里叶变换FFT(Fast Fourier Transform)出现以前，直接用DFT算法进行谱分析和信号的实时处理是不切实际的。算法改进的途径：算法改进的途径：1、把长序列分解成短序列，然后组合出结果。、把长序列分解成短序列，然后组合出结果。2、利用计算中的对称性，减少计算量。、利用计

2、算中的对称性，减少计算量。旋转因子nkNW第4章 快速傅里叶变换(FFT) 根据第3章所讲内容，可知有以下几个特点:(1) 共轭对称性： (2) 周期性：。 (3) 可约性： 。 nkNWnkNWnkNWnkNNnkNnkNnkNWWWW2)(nNkNNnkNnkNWWW)()(nkNnkmNmWWnkNW以N=4为例，利用周期性，利用对称性，54144404)()(,WWWWWWWnNkNNnkNnkN022,)(NNnkNNnkNnkNnkNWWWWWW第4章 快速傅里叶变换(FFT) WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWNNNN

3、NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN101000001010000012302020321000009630642032100000第4章 快速傅里叶变换(FFT) 求出四点DFT实际上只需要一次复数乘法。111111)3() 1 ()2()0() 3()3() 1 ()2()0()2()3() 1 ()2()0() 1 ()3() 1 ()2()0()0() 3()2() 1 ()0(111111111111) 3()2() 1 ()0(WxxxxXxxxxXWxxxxXxxxxXxxxxWWWWXXXX第4章 快速傅里叶变换(FFT) 1

4、、将时域序列按奇偶划分，称之为时间抽取(Decimation in Time，DIT)算法或基2时分算法;2、将频域DFT按奇偶划分，称之为频率抽取(Decimation in Frequency，DIF)算法或基2频分算法。第4章 快速傅里叶变换(FFT) 4.2.2 时域抽取法时域抽取法(基基2FFT)基本原理基本原理设序列x(n)的长度为N，且满足N=2M，M为自然数。按n的奇偶把x(n)分解为两个N/2点的子序列12( )(2 )0112( )(21)0112Nx rxrrNx rxrr， ， ，第4章 快速傅里叶变换(FFT) 则x(n)的DFT为/2 1/2 12(21)00/2

5、1/2 1221200( )( )( )(2 )(21)( )( )knknNNnnNNkrkrNNrrNNkrkkrNNNrrX kx n Wx n Wxr WxrWx r WWx r W偶数奇数第4章 快速傅里叶变换(FFT) (4.2.7) (4.2.8) 运算可用图4.2.1所示的流图符号表示，称为蝶形运算符号。 1210)()()(21NkkXWkXkXkN，1210)()()2(21NkkXWkXNkXkN，第4章 快速傅里叶变换(FFT) 图4.2.1 DIT - FFT蝶式运算信号流图(a) 蝶形流图; (b)蝶形流图简化形式第4章 快速傅里叶变换(FFT) 经过一次奇偶抽取分

6、解后，N点DFT运算图可以用图4.2.2表示。 经过一次分解，可使运算量减少一半。 可以对N/2点DFT再作进一步分解，直到分解成N个1点DFT。共有M级蝶形运算1点DFT就是时域序列本身。第4章 快速傅里叶变换(FFT) 【例1】 将8点序列的DFT用基2时分FFT算法进行分解。解解 (1) 第一次分解。将8点序列x(n)分解为两个4点序列x1(n)和x2(n)，x1(n)为偶数序列，x2(n)为奇数序列，即x1(n)=x(0)，x(2)，x(4)，x(6)x2(n)=x(1)，x(3)，x(5)，x(7)第4章 快速傅里叶变换(FFT) 4点序列x1(n)和x2(n)分别做4点DFT得到X

7、1(k)和X2(k)，由X1(k)和X2(k)通过蝶形构造获得8点DFTX(k)。第一次分解的旋转因子为 (k=0，1，2，3)，运算流图如图4.3.2所示。kW8第4章 快速傅里叶变换(FFT) 图4.2.2 N点DFT基2时分第一次分解运算流图(N=8)第4章 快速傅里叶变换(FFT) (2) 第二次分解。将两个4点序列x1(n)和x2(n)分解为四个2点序列x3(n)、x4(n)、 和。)(3nx)(4nx)4(),0()2(),0()(113xxxxnx)6(),2()3(),1 ()(114xxxxnx)5(),1 ()2(),0()(223xxxxnx)7(),3()3(),1 (

8、)(224xxxxnx第4章 快速傅里叶变换(FFT) 图4.2.3 N点DFT基2时分第二次分解运算流图(N=8)第4章 快速傅里叶变换(FFT) (3) 第三次分解。将四个2点序列分解为8个1点序列，其中)4()1 ()(),0()0()(3635xxnxxxnx)6()1 ()(),2()0()(4847xxnxxxnx)5()1 ()(),1 ()0()(3635xxnxxxnx)7()1 ()(),3()0()(4847xxnxxxnx第4章 快速傅里叶变换(FFT) 图4.2.4 N点DFT基2时分FFT运算流图(N=8)?)4()0(04/xWxN) 7 () 3 (04/xWx

9、aN)5()1 (04/xWxaNaWabN02/第4章 快速傅里叶变换(FFT) 4.2.3 DIT-FFT算法与直接计算算法与直接计算DFT运算量的比较运算量的比较当N=2M 时，其运算流图应有M级蝶形，M级运算总共需要的复数乘次数为lb22MNNCMNlbACN MNN复数加次数为第4章 快速傅里叶变换(FFT) 当N1时，N2(N/2) lbN，所以，DIT-FFT算法比直接计算DFT的运算次数大大减少。例如，N=210=1024时，8 .2045120576 048 1lb22NNN第4章 快速傅里叶变换(FFT) 图4.2.5 DIT-FFT算法与直接计算DFT所需复数乘法次数的比

10、较曲线第4章 快速傅里叶变换(FFT) 4.2.4 DIT-FFT的运算规律及编程思想的运算规律及编程思想1 原位计算原位计算DIT-FFT的运算过程的规律。1、N=2M点的FFT共进行M级运算;2、同一级中，每个蝶形的两个输入数据只对计算本蝶形有用；3、在计算完一个蝶形后，所得输出数据可立即存入原输入数据所占用的存储单元(数组元素)。 利用同一存储单元存储蝶形计算输入、输出数据的方法称为原位(址)计算。第4章 快速傅里叶变换(FFT) 2 旋转因子的变化规律旋转因子的变化规律在N点DIT-FFT运算流图中，每个蝶形都要乘以因子，称其为旋转因子，p为旋转因子的指数。找出旋转因子与运算级数的关系

11、。pNW个不同的旋转因子。级共有数。第表示从左到右的运算级用1 -L2LL第4章 快速傅里叶变换(FFT) 对N=2M的一般情况，第L级的旋转因子为3, 2, 1, 0 31, 0 20 1222/24/JWWWLJWWWLJWWWLJJNpNJJNpNJJNpNLLL时时时第4章 快速傅里叶变换(FFT) 因为所以(4.2.12) (4.2.13)用(4.2.12)和(4.2.13)式确定第L级运算的旋转因子(实际编程序时，L为最外层循环变量)。L-12,0,1,2,21LpJNWWJMLMLMLN222212, 2, 1, 0122LJNJNpNJWWWLMMLLMJp2第4章 快速傅里叶

12、变换(FFT) 3 蝶形运算规律蝶形运算规律设序列x(n)经时域抽选(倒序)后，按图4.2.4所示的次序(倒序)存入数组A中。如果蝶形运算的两个输入数据相距B个点，应用原位计算，则蝶形运算可表示成如下形式：式中下标L表示第L级运算，AL(J)则表示第L级运算后的数组元素A(J)的值（即第L级蝶形的输出数据）。而AL1(J)表示第L级运算前A(J)的值（即第L级蝶形的输入数据）。11( )( )()pLLLNA JAJAJB W11()( )()pLLLNA JBAJAJB W12 0,1,21; 1,2,MLLpJJLM第4章 快速傅里叶变换(FFT) 4. 编程思想及程序框图编程思想及程序框

13、图 观察图4.2.4，归纳出对编程有用的运算规律：第L级中，每个蝶形的两个输入数据相距B=2L1个点；每级有B个不同的旋转因子；同一旋转因子对应着间隔为2L点的2ML个蝶形。三重循环程序实现DIT-FFT运算，程序框图如图4.2.6所示。 第4章 快速傅里叶变换(FFT) 图4.2.6 DIT-FFT运算和程序框图第4章 快速傅里叶变换(FFT) DIT-FFT算法运算流图的输出X(k)为自然顺序。输入序列不是按x(n)的自然顺序排列，排序称为序列x(n)的倒序(倒位)。在运算M级蝶形之前应先对序列x(n)进行倒序。5 序列的倒序序列的倒序规律如图4.2.7所示。第4章 快速傅里叶变换(FFT

14、) 图4.2.7 形成例序的树状图（N=23）第4章 快速傅里叶变换(FFT) 表4.2.1 顺序和倒序二进制数对照表第4章 快速傅里叶变换(FFT) 图4.2.9所示的倒序的程序框图中的虚线框内就是完成计算倒序值的运算流程图。 第4章 快速傅里叶变换(FFT) 图4.2.9 倒序程序框图第4章 快速傅里叶变换(FFT) 第3章介绍的MATLAB函数fft是一个计算DFT的智能程序，如果计算点数N=2M，则自动按DIT-FFT快速算法计算。第4章 快速傅里叶变换(FFT) 4.2.5 频域抽取法频域抽取法FFT(DIF-FFT)在基2FFT算法中，频域抽取法FFT也是一种常用的快速算法，简称D

15、TF- FFT。设序列x(n)长度为N=2M，首先将x(n)前后对半分开，得到两个子序列，其DFT可表示为如下形式：10/2 110/2/2 1/2 1(/2)00/2 1/20( )DFT ( )( )( )( )( )2( )2NknNnNNknknNNnn NNNknk n NNNnnNkNknNNnX kx nx n Wx n Wx n WNx n Wx nWNx nWx nW第4章 快速傅里叶变换(FFT) 式中 将X(k)分解成偶数组与奇数组，当k取偶数(k=2m, m=0, 1, , N/21)时/21( 1)1kNkNkWk 偶数奇数，/220/2 1/20(2 )( )2(

16、)(4.2.14)2NmnNnNmnNnNXmx nx nWNx nx nW第4章 快速傅里叶变换(FFT) 当k取奇数(k=2m+1, m=0, 1, , N/21)时，/2 1(21)0/2 1/20(21)( )2( )2 (4.2.15)NnmNnNnnmNNnNXmx nx nWNx nx nWW令122102)()(2)()(21NnWNnxnxnxNnxnxnxnN，第4章 快速傅里叶变换(FFT) 将x1(n)和x2(n)分别代入(4.2.14)和(4.2.15)式，可得 (4.2.16)(4.2.16)式表明，X(k)按奇偶k值分为两组，其偶数组是x1(n)的N/2点DFT，

17、奇数组则是x2(n)的N/2点DFT。x1(n)、x2(n)和x(n)之间的关系也可用图4.2.10所示的蝶形运算流图符号表示。图4.2.11表示N=8时第一次分解的运算流图。/2 11/20/2 12/20(2 )( )(21)( )NmnNnNmnNnXmx n WXmx n W第4章 快速傅里叶变换(FFT) 图4.2.10 DTFFFT蝶形运算流图符号第4章 快速傅里叶变换(FFT) 图4.2.11 DIF-FFT第一次分解运算流图（N=8）第4章 快速傅里叶变换(FFT) 由于N=2M，N/2仍然是偶数，继续将N/2点DFT分成偶数组合奇数组。图4.2.12表示N=8时第二次分解运算

18、流图。 经过M1次分解，最后分解为2M1个两点DFT，两点DFT就是一个基本蝶形运算流图。 当N=8，经两次分解，便分解为四个两点DFT。 N = 8的完整DIF-FFT运算流图如图4.2.13所示。第4章 快速傅里叶变换(FFT) 图4.2.12 DIF-FFT第二次分解运算流图（N = 8）第4章 快速傅里叶变换(FFT) 图4.2.13 DIF-FFT运算流图（N =8）第4章 快速傅里叶变换(FFT) 这种算法是对X(k)进行奇偶抽取分解的结果，所以称之为频域抽取法FFT。观察图4.2.13可知，DIF-FFT算法与DIT-TTF算法类似。不同的是DIF-FFT算法输入序列为自然顺序，

19、而输出为倒序排列。 要对输出数据进行倒序才能得到自然顺序的X(k)。 另外，蝶形运算略有不同，DIT-FFT蝶形先乘后加(减)，而DIF-FFT蝶形先加(减)后相乘。第4章 快速傅里叶变换(FFT) 4.2.6 IDFT的高效算法的高效算法上述FFT算法流图也可以用于计算IDFT。比较DFT和IDFT的运算公式：1010( )DFT ( )( )1( )IDFT ( )( )NknNnNknNkX kx nx n Wx nx nX k WN第4章 快速傅里叶变换(FFT) 只要将DFT运算式中的系数改变为，最后乘以1/N，就是IDFT运算公式。 在DIT-FFT与DIF-FFT算法中的旋转因子

20、改为，最后的输出再乘以1/N就可以用来计算IDFT。流图的输入是X(k)，输出就是x(n)。 原来的DIT-FFT改为IFFT后，称为DIF-IFFT更合适；DIF-FFT改为IFFT后, 应称为DIT-IFFT。knNWpNWpNWknNW第4章 快速傅里叶变换(FFT) 教材第教材第4章习题与上机题解答章习题与上机题解答快速傅里叶变换（FFT）是DFT的快速算法， 没有新的物理概念。 1 如果某通用单片计算机的速度为平均每次复数乘需要4 s， 每次复数加需要1 s， 用来计算N=1024点DFT， 问直接计算需要多少时间。 用FFT计算呢？照这样计算， 用FFT进行快速卷积对信号进行处理时

21、， 估计可实现实时处理的信号最高频率。 第4章 快速傅里叶变换(FFT) 解解: 当N=1024=210时， 直接计算DFT的复数乘法运算次数为N2=10241024=1 048 576次复数加法运算次数为N（N1）=10241023=1 047 552次直接计算所用计算时间TD为TD=410610242+1 047 552106=5.241 856 s用FFT计算1024点DFT所需计算时间TF为第4章 快速傅里叶变换(FFT) 66F665 10lblb10210245 1010 1024 10 10230.72 msNTNN N 快速卷积时， 需要计算一次N点FFT（考虑到H(k)=DF

22、Th(n)已计算好存入内存）、 N次频域复数乘法和一次N点IFFT。 所以， 计算1024点快速卷积的计算时间Tc约为第4章 快速傅里叶变换(FFT) cF2102471680 s4 1024 s65536 sTT 次复数乘计算时间所以， 每秒钟处理的采样点数（即采样速率）s6102415 625 /65536 10F次 秒由采样定理知， 可实时处理的信号最高频率为smax156257.8125 kHz22Ff第4章 快速傅里叶变换(FFT) 应当说明， 实际实现时， fmax还要小一些。 这是由于实际中要求采样频率高于奈奎斯特速率， 而且在采用重叠相加法时， 重叠部分要计算两次。 重叠部分长

23、度与h(n)长度有关， 而且还有存取数据和指令周期等消耗的时间。 2 如果将通用单片机换成数字信号处理专用单片机TMS320系列， 计算复数乘和复数加各需要10 ns。 请重复做上题。 解解: 与第1题同理。 直接计算1024点DFT所需计算时间TD为TD=1010910242+101091 047 552=20.961 28 ms第4章 快速傅里叶变换(FFT) 用FFT计算1024点DFT所需计算时间TF为99F8810 10lb10 10lb210241010 101024 1020.1536 msNTNNN快速卷积计算时间Tc约为cF3921024 2 0.1536 1010 1010

24、24 0.317 44 msTT次复数乘计算时间第4章 快速傅里叶变换(FFT) 可实时处理的信号最高频率fmax为maxsc1110241 = 3.1158 MHz=1.6129 MHz222fFT由此可见， 用DSP专用单片机可大大提高信号处理速度。 所以， DSP在数字信号处理领域得到广泛应用。 机器周期小于1 ns的DSP产品已上市， 其处理速度更高。 第4章 快速傅里叶变换(FFT) 3 已知X(k)和Y(k)是两个N点实序列x(n)和y(n)的DFT， 希望从X(k)和Y(k)求x(n)和y(n)， 为提高运算效率， 试设计用一次N点IFFT来完成的算法。 解解: 因为x(n)和y

25、(n)均为实序列， 所以， X(k)和Y(k)为共轭对称序列， jY(k)为共轭反对称序列。 可令X(k)和jY(k)分别作为复序列F(k)的共轭对称分量和共轭反对称分量， 即F(k)=X(k)+jY(k)=Fep(k)+Fop(k)计算一次N点IFFT得到f(n)=IFFTF(k)=Ref(n)+j Imf(n)第4章 快速傅里叶变换(FFT) 由DFT的共轭对称性可知Ref(n)=IDFTFep(k)=IDFTX(k)=x(n)j Imf(n)=IDFTFop(k)=IDFTjY(k)=jy(n)故1( ) ( )( )2x nf nfn1( ) ( )( )2jy nf nfn4 设x(

26、n)是长度为2N的有限长实序列， X(k)为x(n)的2N点DFT。 (1) 试设计用一次N点FFT完成计算X(k)的高效算法。 (2) 若已知X(k) ，试设计用一次N点IFFT实现求X(k)的2N点IDFT运算。第4章 快速傅里叶变换(FFT) 解解: 本题的解题思路就是DIT-FFT思想。（1） 在时域分别抽取偶数和奇数点x(n)， 得到两个N点实序列x1(n)和x2(n)： x1(n)=x(2n) n=0, 1, , N1x2(n)=x(2n+1) n=0, 1, , N1根据DIT-FFT的思想， 只要求得x1(n)和x2(n)的N点DFT， 再经过简单的一级蝶形运算就可得到x(n)

27、的2N点DFT。 因为x1(n)和x2(n)均为实序列， 所以根据DFT的共轭对称性， 可用一次N点FFT求得X1(k)和X2(k)。 具体方法如下：第4章 快速傅里叶变换(FFT) 令 y(n)=x1(n)+jx2(n)Y(k)=DFTy(n) k=0, 1, , N1则*11ep*22ep1( )DFT ( )( ) ( )()21j( )DFTj( )( ) ( )()2X kx nYkY kYNkXkx nYkY kYNk2N点DFTx（n）=X（k）可由X1(k)和X2(k)得到122122( )( )( ) 0,1,1()( )( )kNkNX kX kWXkkNX kNX kWX

28、k第4章 快速傅里叶变换(FFT) 这样， 通过一次N点IFFT计算就完成了计算2N点DFT。 当然还要进行由Y(k)求X1(k)、 X2(k)和X(k)的运算(运算量相对很少)。 （2） 与（1）相同， 设x1(n)=x(2n) n=0, 1, , N1x2(n)=x(2n+1) n=0, 1, , N1X1(k)=DFTx1(n) k=0, 1, , N1X2(k)=DFTx2(n) k=0, 1, , N1则应满足关系式122122( )( )( ) 0,1,1()( )( )kNkNX kX kWXkkNX kNX kWXk第4章 快速傅里叶变换(FFT) 由上式可解出1221( )(

29、 )()2 0,1,2,11( )( )()2kNX kX kX kNkNXkX kX kN W由以上分析可得出运算过程如下： 由X(k)计算出X1(k)和X2(k): 1221( )( )()21( )( )()2kNX kX kX kNXkX kX kN W第4章 快速傅里叶变换(FFT) 由X1(k)和X2(k)构成N点频域序列Y(k): Y(k)=X1(k)+jX2(k)=Yep(k)+Yop(k)其中， Yep(k)=X1(k), Yop(k)=jX2(k)， 进行N点IFFT， 得到y(n)=IFFTY(k)=Rey(n)+j Imy(n) n=0, 1, , N1由DFT的共轭对

30、称性知*ep1*op21Re ( ) ( )( )DFT( )( )21jIm ( ) ( )( )DFT( )j( )2y ny ny nYkx ny ny ny nYkx n第4章 快速傅里叶变换(FFT) 由x1(n)和x2(n)合成x(n):12 2( )1 2nxnx nnxn偶数奇数，0n2N1在编程序实现时， 只要将存放x1(n)和x2(n)的两个数组的元素分别依次放入存放x(n)的数组的偶数和奇数数组元素中即可。第4章 快速傅里叶变换(FFT) 5 分别画出16点基2DIT-FFT和DIF-FFT运算流图， 并计算其复数乘次数， 如果考虑三类碟形的乘法计算， 试计算复乘次数。

31、解解： 本题比较简单， 仿照教材中的8点基2DIT-FFT和DIF-FFT运算流图很容易画出16点基2DIT-FFT和DIF-FFT运算流图。 但画图占篇幅较大， 这里省略本题解答， 请读者自己完成。第4章 快速傅里叶变换(FFT) 6* 按照下面的IDFT算法编写MATLAB语言 IFFT程序， 其中的FFT部分不用写出清单， 可调用fft函数。 并分别对单位脉冲序列、 矩形序列、 三角序列和正弦序列进行FFT和IFFT变换， 验证所编程序。 *1( )IDFT( )DFT( )x nX kXkN第4章 快速傅里叶变换(FFT) 解解： 为了使用灵活方便， 将本题所给算法公式作为函数编写ifft46.m如下： %函数ifft46.m%按照所给算法公式计算IFETfunction xn=ifft46(Xk，