数字信号处理(丁玉美版)教案第4章VIP

1、1第四章第四章 快速傅立叶变换快速傅立叶变换（fft） fast fourier transform 数字信号处理2 4.1 引言引言 dft 是数字信号处理的基础，是数字信号处理的基础，是是 一种重要变换。一种重要变换。 快速算法的引出，使数字信号处快速算法的引出，使数字信号处 理技术得到广泛应用。理技术得到广泛应用。 各种快速算法不断被采用各种快速算法不断被采用3t)(txa)( txada/)(nx)(nxdft)()(fxkxa数字计算机n足够大.1.1 dft提供了计算机处理的可能性 模拟信号的频谱分析 4.1.2 dft的运算量 10)()()(nnnknwnxnxdftkxk=0

2、,1,2,n1计算机运算时： 0k0) 1(0100) 1() 1 () 0() 0(nnnnwnxwxwxx1k 0111(1)1(1)(0)(1)(1)nnnnxxwxwx nw 2k 0212(1)2(2)(0)(1)(1)nnnnxxwxwx nw 1nk0111(1)1(1)(0)(1)(1)nnnnnnnx nxwx wx nw n项 n个 计算一个 n点dft ，共需次复乘。2n以做一次复乘1s计，若n =4096，所需时间为由于计算量大，且要求相当大的内存，难以实现实时处理，限制了dft的应用。s16777216)4096(2s175 长期以来，人们一直在寻求一种能提高dft运

3、算速度的方法。fft便是cooley&tukey 在1965 年提出的的快速算法，它可以使运算速度提高几百倍，从而使数字信号处理学科成为一个新兴的应用学科。4.1.3 fft算法的设计思想 1利用nknw的特点njnew2具有 1）周期性 knnnnknww)( )(nknnw2）共轭性 knnnnknww)()(3）对称性 knnknww)2(4）恒等性122nnnnww5）可约性nnnnww22reimj8n0nw1nw2nw3nw4nw5nw6nw7nw)(knnnw2把n 点dft化为几组点数较少的dft运算 n点dft运算的复乘次数为2n次，若将n点dft化为2 组，则复乘次数约为2

4、2222nn次。74.2 基基 2 抽取抽取fft 算法算法（the decimation-in-time radix-2 fft algorithm）fft分为两大类：分为两大类：时域抽取时域抽取fft（decimation-in-time fft，简称，简称dit-fft）频域抽取频域抽取fft（decimation-in-frequency fft， 简称简称dif-fft）8 4.2.1 直接计算直接计算dft的问题及改进的途径的问题及改进的途径 k=0, 1, , n-110)()(nnknnwnxkx10)(1)(nkknnwkxnnxn=0, 1, , n-1dft及及idft的

5、定义的定义9可见，可见，dft 与与 idft 的计算成本基本相同。的计算成本基本相同。直接计算直接计算n点点dft 时：时：对应一个对应一个k需要需要n次复数乘和（次复数乘和（n-1）次）次复数加；复数加；对所有对所有n个个k值，则需要值，则需要 n复数乘和复数乘和n (n-1)次复数加次复数加 。其中其中:一次复数乘需要一次复数乘需要4次实数乘和次实数乘和2次实数加方能次实数加方能完成。当完成。当n较大时，运算量很大。较大时，运算量很大。10 例如：当例如：当n=8时，时，dft需要需要64次复数乘；次复数乘；当当n=1024时，时，dft所需复数乘为所需复数乘为1048576次，即一百多

6、万次复数乘。次，即一百多万次复数乘。为了减少运算次数，为了减少运算次数，改进计算的方法：改进计算的方法：1）利用旋转因子的对称性、周期性和可约性；）利用旋转因子的对称性、周期性和可约性； 旋转因子旋转因子（twiddle factor）2）使长序列变短。）使长序列变短。njnew/2114.2.2 基基2时域抽取法原理时域抽取法原理 （库利（库利-图基算法）图基算法）设序列设序列x（n）的长度为）的长度为n，且，且m为自然数为自然数n-point dft,n is evenmn2)(log2nm 1,.,1 , 0,)()(10nkwnxkxnnknn12将其一分为二，分成偶数和奇数序列项将其

7、一分为二，分成偶数和奇数序列项（the even-indexed and odd-indexed terms）则则n/2点的序列为：点的序列为： even: x1(r)=x(2r) , r=0, 1, , n/2-1 odd: x2(r)=x(2r+1) , r=0, 1, , n/2-113偶数序列偶数序列 the range: 02nn-2 (n is even) 0n(n/2)-1奇数序列奇数序列 the range: 12n+1n-1 (n is even) 02nn-2 0n(n/2)-114 奇数偶数nknnnknnwnxwnxkx)()()(120)12(120)2() 12()

8、2(nrrknnrrknwrxwrx1202212021)()(nrkrnknnrkrnwrxwwrx则则x（n）的）的dft：15由于由于所以所以krnkrnjkrnjkrnweew2/2/2222/2 1/2 112/2/20012( )( )( )( )( )nnkrkkrnnnrrknx kx r wwx r wx kw xkk=0,1,n-116n上式说明,按n的奇偶性将 分解为两个n/2长的序列 和 ,则n点dft可分解为两个n/2点的dft来计算.)(nx)(1nx)(2nx17 其中其中 n/2-point dft: k=0,1,n/2-1 )()()(112/02/11rxd

9、ftwrxkxnrkrn12/022/22)()()(nrkrnrxdftwrxkxk=0,1,n/2-118因此，可写出两个因此，可写出两个n/2点的方程点的方程:)()()(21kxwkxkxkn)2()2()2(2)2(1nkxwnkxnkxnkn12,.,1 ,0nk19 而而120)2(2)()2(11nrrnknwrxnkxknnknww222)()()2(1112021kxwrxnkxnrkrn20同理同理而而1222jnnjnneew)()2(22kxnkx21所以所以)()()2()()()(2121kxwkxnkxkxwkxkxknkn12,.1 , 0nk22表示上述算法

10、可用蝶形结（表示上述算法可用蝶形结（ butterfly）knw)()(21kxwkxkn)()(21kxwkxkn23 example : 如如n=4 x(n)=x0, x1, x2, x3 even: x0, x2 even: x0 odd: x2 odd: x1, x3 even: x1 odd: x324 bit reversal/shuffling process x x0 0 x x2 2x x1 1x x3 3x x0 0 x x2 2x x1 1x x2 2x x0 0 x x1 1x x2 2x x3 3混序或反混序或反序序码位倒置码位倒置分成四个分成四个1点的序列点的序列

11、25 the butterfly（蝶形运算）（蝶形运算）04w14w02w02w1点序列的点序列的dft就是序列本身，即不用计算就是序列本身，即不用计算 -1-1-1-126 如如n4，则，则将将 x1(r) 再按再按r的奇偶进一步分解成两个的奇偶进一步分解成两个n/4点长的子序列：点长的子序列：) 12()()2()(1413lxlxlxlx14,.,1 , 0nl2714/0)12(2/114/022/11) 12()2()(nllknnlklnwlxwlxkx12,.,1 , 0nk)()(42/3kxwkxkn14/04/42/14/04/3)()(nlklnknnlklnwlxwwl

12、x28 其中其中)()()(314/04/33lxdftwlxkxnlkln)()()(414/04/44lxdftwlxkxnlkln14,.,1 , 0nk29由由x3（k）和）和x4（k）的周期性和）的周期性和 的对称性，得的对称性，得knnknww2/4/2/)()()4()()()(42/3142/31kxwkxnkxkxwkxkxknkn14,.,1 , 0nk30同理，得同理，得)()()4()()()(62/5262/52kxwkxnkxkxwkxkxknkn14,.,1 , 0nk31 其中其中) 12()()2()(2625lxlxlxlx)()()(514/04/55lx

13、dftwlxkxnlkln)()()(614/04/66lxdftwlxkxnlkln32 8点点dit-fft运算流图运算流图x(0)x(0)x(4)x(4)x(2)x(2)x(6)x(6)x(1)x(1)x(5)x(5)x(3)x(3)x(7)x(7)x x3 3(0)(0)x x3 3(1)(1)x x4 4(0)(0)x x4 4(1)(1)x x5 5(0)(0)x x5 5(1)(1)x x6 6(0)(0)x x6 6(1)(1)x x1 1(0)(0)x x1 1(1)(1)x x1 1(2)(2)x x1 1(3)(3)x x2 2(0)(0)x x2 2(1)(1)x x2

14、 2(2)(2)x x2 2(3)(3)x(0)x(0)x(1)x(1)x(2)x(2)x(3)x(3)x(4)x(4)x(5)x(5)x(6)x(6)x(7)x(7)02w02w02w02w04w04w14w14w08w28w18w38w33 4.2.3 idft的高效算法的高效算法10)(1)()(nkknnwkxnkxidftnx这样这样ifft可与可与fft用同一子程序实现。用同一子程序实现。*10*)(1nkknnwkxn*)(1kxdftn34idft的运算规律的运算规律x*(0)x*(0)x*(4)x*(4)x*(2)x*(2)x*(6)x*(6)x*(1)x*(1)x*(5)x

15、*(5)x*(3)x*(3)x*(7)x*(7)x x3 3(0)(0)x x3 3(1)(1)x x4 4(0)(0)x x4 4(1)(1)x x5 5(0)(0)x x5 5(1)(1)x x6 6(0)(0)x x6 6(1)(1)x x1 1(0)(0)x x1 1(1)(1)x x1 1(2)(2)x x1 1(3)(3)x x2 2(0)(0)x x2 2(1)(1)x x2 2(2)(2)x x2 2(3)(3)x(0)*x(0)*x(1)*x(1)*x(2)*x(2)*x(3)*x(3)*x(4)*x(4)*x(5)*x(5)*x(6)*x(6)*x(7)*x(7)*02w0

16、2w02w02w04w04w14w14w08w28w18w38w1/n35 求求ifft，也可用，也可用dit-fft的流程来实现。的流程来实现。x(0)x(0)x(4)x(4)x(2)x(2)x(6)x(6)x(1)x(1)x(5)x(5)x(3)x(3)x(7)x(7)x(0)x(0)x(1)x(1)x(2)x(2)x(3)x(3)x(4)x(4)x(5)x(5)x(6)x(6)x(7)x(7)0nw1nw2nw3nw0nw0nw0nw0nw0nw0nw2nw2nw1/n1/n1/n1/n1/n1/n1/n1/n1/n1/n1/n1/n1/n1/n1/n1/n36 example:dete

17、rmine the x(n) by idft.x=5, -1, 1, -1 solution: n=0,1,2,3*304*)(41)(1)(kknkxwxdftnnx371)(41)(41)0(3210*3004xxxxxwxkk1)15(41)(41) 1 (*304jjxwxkkk2)(41)2(30*24kkkxwx38 所以所以 x(n)=1, 1, 2, 1 1)(41)3(*3034kkkxwx39%the program in matlab: x=input(please input x=); n=length(x); x=fft(conj(x),n); x=conj(x/n)

18、40 example :用用fft算法计算下面信号的算法计算下面信号的 8点点dft； x(n)=4 3 2 0 1 2 3 1 然后，再用然后，再用 ifft恢复原信号。恢复原信号。41 solution:4 4-3-32 20 0-1-1-2-23 31 14 42 2-1-13 3-3-30 0-2-21 14 4-1-12 23 3-3-3-2-20 01 13 35 55 5-1-1-5-5-1-11 1-1-11 1-j-j1 1-j-j8 85+j5+j-2-25-j5-j-4-4-1+j-1+j-6-6-1-j-1-j1 1(1-j)/(1-j)/-j-j-(1+j)/-(1+

19、j)/4 45+j+j5+j+j-2+j6-2+j65-j+j5-j+j12125+j-j5+j-j-2-j6-2-j65-j-j5-j-jshufflingshufflingdft mergingdft mergingx x0 0x x1 1x x2 2x x3 3x x4 4x x5 5x x6 6x x7 7 2 2 2 2 2 242 ifft(x)=1/nfft(x*)*4 45 5- -j j- -j j- -2 2- -j j6 65 5+ +j j- -j j1 12 25 5- -j j+ +j j- -2 2+ +j j6 65 5+ +j j+ +j js sh hu u

20、f ff fl li in ng gd df ft t m me er rg gi in ng g4 4- -2 2- -j j6 61 12 2- -2 2+ +j j6 65 5- -j j- -j j5 5+ +j j- -j j5 5- -j j+ +j j5 5+ +j j+ +j j4 41 12 2- -2 2- -j j6 6- -2 2+ +j j6 65 5- -j j- -j j5 5- -j j+ +j j5 5+ +j j- -j j5 5+ +j j+ +j j1 16 6- -8 8- -4 4- -1 12 2j j1 10 0- -j j2 2- -j j2

21、21 10 0+ +j j2 2- -j j2 21 12 2- -2 20 02 20 04 42 20 0- -2 2 ( (1 1+ +j j) )- -4 4j j2 2 ( (1 1- -j j) )3 32 2- -2 24 41 16 60 0- -8 8- -1 16 62 24 48 8 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 204w04w14w14w08w18w28w38w乘以1/n=1/8得原序列43 4.2.4 fft的计算成本的计算成本最简单的最简单的fft 是是cooley-tukey法法因为因为nmnm2log2n点点dft有有m级蝶形运算，

22、每一级都有级蝶形运算，每一级都有n/2个个蝶形运算结构；蝶形运算结构；每个蝶形运算结构都有每个蝶形运算结构都有1次复数乘和次复数乘和2次复数次复数加；加； 8442222111111118-dft8-dft4-dft4-dft2-dft2-dft1-dft1-dftstage 1stage 1stage 2stage 2stage 3stage 345所以，所以，m级蝶形运算有复数乘级蝶形运算有复数乘m级蝶形运算有复数加级蝶形运算有复数加nnmncm2log22)2(nnmnca2log)2(46直接直接dft的复数乘和复数加均近似为的复数乘和复数加均近似为 。当。当n1时，时，所以所以dit

23、-fft大大减少了运算次数。大大减少了运算次数。nnn22log)2(2n47 example : for n=8 solution : m= =3 stages（三级）（三级） for fft, the total cost is （fft的总计算成本是）的总计算成本是） mn/2=3 8/2=12 complex multiplications （复数乘）（复数乘）8log248 for directly dft, the total cost is（直接计算（直接计算dft的成本是）的成本是） n=8=64 (complex multiplications) the ratio is:

24、（比值是）（比值是） 实际上，当实际上，当n=2048时，直接运算与时，直接运算与 fft算法的计算量之比为算法的计算量之比为372.433. 5126449 4.2.5 dit-fft的运算规律及编程思想的运算规律及编程思想1、原位运算：、原位运算：利用同一单元存储蝶形计算的输入、输出利用同一单元存储蝶形计算的输入、输出数据。数据。每个蝶形的输入和输出均为相同位数。每个蝶形的输入和输出均为相同位数。原位运算可节省大量内存，因而硬件成本低。原位运算可节省大量内存，因而硬件成本低。2、旋转因子的变化规律：、旋转因子的变化规律： n点点dft，共有，共有m级蝶形运算，每级有级蝶形运算，每级有n/2

25、个蝶形。个蝶形。50 称为旋转因子，称为旋转因子，p称为旋转因子的指数。称为旋转因子的指数。 为了编写程序，应找出旋转因子为了编写程序，应找出旋转因子 与运算与运算 级数的关系。级数的关系。 设设l=1，2，m，表示从左到右的运算，表示从左到右的运算 级数，第级数，第l级有级有 个不同的旋转因子，个不同的旋转因子，pnw12lpnw51对于对于 ，各级的旋转因子表示为，各级的旋转因子表示为l=1时，时，l=2时，时，l=3时，时，823n0,24/jwwwjjnpnl1 , 0,22/jwwwjjnpnl3 , 2 , 1 , 0,2jwwwjjnpnl52对于对于 的一般情况，第的一般情况，

26、第l级的旋转因子为级的旋转因子为mn2jpnlww212,.,2 , 1 , 01lj，由于由于mlmlmln222253所以所以通过上式，可以计算第通过上式，可以计算第l级运算的旋转因子。级运算的旋转因子。lmmljnjnpnwww2212,.,1 , 01ljlmjp254 3、蝶形运算规律、蝶形运算规律设序列设序列x（n）经过时域倒序存入数组）经过时域倒序存入数组x如蝶形运算的两个输入数据相距如蝶形运算的两个输入数据相距b个点，个点，应用原位运算，可得：应用原位运算，可得：pnlllpnlllwbjxjxbjxwbjxjxjx)()()()()()(1111mljjpllm,.,1; 1

27、2,.,1 , 0;2155 4、编程思想及程序框图、编程思想及程序框图其它运算规律：其它运算规律：第第l级中，每个蝶形的两个输入数据相距级中，每个蝶形的两个输入数据相距 个点；个点；同一旋转因子对应着间隔为同一旋转因子对应着间隔为 点的点的 个蝶形（对应第几组蝶形）个蝶形（对应第几组蝶形）12lblm2l256 8点点dit-fft运算流图运算流图x(0)x(0)x(4)x(4)x(2)x(2)x(6)x(6)x(1)x(1)x(5)x(5)x(3)x(3)x(7)x(7)x x3 3(0)(0)x x3 3(1)(1)x x4 4(0)(0)x x4 4(1)(1)x x5 5(0)(0)

28、x x5 5(1)(1)x x6 6(0)(0)x x6 6(1)(1)x x1 1(0)(0)x x1 1(1)(1)x x1 1(2)(2)x x1 1(3)(3)x x2 2(0)(0)x x2 2(1)(1)x x2 2(2)(2)x x2 2(3)(3)x(0)x(0)x(1)x(1)x(2)x(2)x(3)x(3)x(4)x(4)x(5)x(5)x(6)x(6)x(7)x(7)02w02w02w02w04w04w14w14w08w28w18w38w57编程的运算方法：编程的运算方法：从输入端（第从输入端（第1级）开始，逐级进行，共级）开始，逐级进行，共进行进行m级运算。级运算。在进

29、行第在进行第l级运算时，依次求出级运算时，依次求出 个不同的个不同的旋转因子，然后计算其对应相同的旋转因子的旋转因子，然后计算其对应相同的旋转因子的 个蝶形。个蝶形。可用三重循环程序实现可用三重循环程序实现dit-fft运算。运算。12llm258(3)595、序列的倒置（、序列的倒置（bit reversal）倒序是有规律的。倒序是有规律的。由于由于 ，所以顺序数可用，所以顺序数可用m位二进制位二进制数（数（ ）表示。）表示。mn20121. nnnnmm60用硬件电路和汇编语言程序产生倒序很容易，用硬件电路和汇编语言程序产生倒序很容易，用高级语言倒序的规律为：用高级语言倒序的规律为：倒序数

30、是在倒序数是在m位二进制数最高位加位二进制数最高位加1，逢逢2向右进位。向右进位。614.2.6 频率抽取法频率抽取法fft(dif-fft)设序列设序列x（n）长度为）长度为 ，将其前后，将其前后对半分开，得：对半分开，得：mn212/12/010)()()()()(nnnknnnnknnnnknnwnxwnxwnxnxdftkx62式中12/0)2/(12/0)2()(nnnnknnnknnwnnxwnx奇数，偶数kkwkknn1, 1) 1(2/knnknnnnwnnxwnx)2()(2/12/063再将再将x（k）分解成偶数组和奇数组）分解成偶数组和奇数组 k为偶数时：为偶数时：/2

31、120/2 1/20(2 ) ( )()2 ( )()2nrnnnnrnnnnxrx nx nwnx nx nw64k为奇数时：为奇数时：rnnnnnnnrnnnwwnnxnxwnnxnxrx2/12/0)12(12/0)2()()2()() 12(65令令12,.,1 , 0,)2()()()2()()(21nnwnnxnxnxnnxnxnxnn66得得12/02/212/02/1)() 12()()2(nnrnnnnrnnwnxrxwnxrx12,.,1 , 0nr67 dif-fft蝶形运算流图蝶形运算流图x x(n)(n)x x(n+n/2)(n+n/2)x x1 1(n)=x(n)+

32、x(n/2)(n)=x(n)+x(n/2)x x2 2(n)=x(n)-x(n+n/2)(n)=x(n)-x(n+n/2)nnwnnw-+68n=8时，时，dif-fft蝶形运算流图蝶形运算流图x(0)x(0)x(4)x(4)x(2)x(2)x(6)x(6)x(1)x(1)x(5)x(5)x(3)x(3)x(7)x(7)x(0)x(0)x(1)x(1)x(2)x(2)x(3)x(3)x(4)x(4)x(5)x(5)x(6)x(6)x(7)x(7)0nw1nw2nw2nw2nw3nw0nw0nw0nw0nw0nw0nw69注意：注意：dit-fft和和dif-fft的算法流图不的算法流图不 是唯一的。是