数字信号处理丁玉美PPT课件

1 4.1 引言 DFT 是数字信号处理的基础，是 一种重要变换。 快速算法的引出，使数字信号处 理技术得到广泛应用。 各种快速算法不断被采用第1页/共77页2T)(txa)( txaDA/)(nx)(nxDFT)()(fXkXa数字计算机N足够大.1.1 DFT提供了计算机处理的可能性 模拟信号的频谱分析 4.1.2 DFT的运算量 10)()()(NnnkNWnxnxDFTkXk=0,1,2,N1计算机运算时： 第2页/共77页0k0) 1(0100) 1() 1 ()0()0(NNNNWNxWxWxX1k 0111(1)1(1)(0)(1)(1)NNNNXxWxWx NW 2k 0 212(1) 2(2)(0)(1)(1)NNNNXxWxWx NW 1 Nk0111(1)1(1)(0)(1)(1)NNNNNNNX NxWxWx NW N项 N个 计算一个 N点DFT ，共需次复乘。2N以做一次复乘1s计，若N =4096，所需时间为由于计算量大，且要求相当大的内存，难以实现实时处理，限制了DFT的应用。s16777216)4096(2s17第3页/共77页4 长期以来，人们一直在寻求一种能提高DFT运算速度的方法。FFT便是Cooley&Tukey 在1965 年提出的的快速算法，它可以使运算速度提高几百倍，从而使数字信号处理学科成为一个新兴的应用学科。4.1.3 FFT算法的设计思想 1利用nkNW的特点NjNeW2具有 1）周期性 kNnNnkNWW)( )(NknNW第4页/共77页2）共轭性 knNNnkNWW)()(3）对称性 kNNkNWW)2(4）恒等性122NNNNWW5）可约性nNnNWW22ReImj8N0NW1NW2NW3NW4NW5NW6NW7NW)(kNnNW2把N 点DFT化为几组点数较少的DFT运算 N点DFT运算的复乘次数为2N次，若将N点DFT化为2 组，则复乘次数约为22222NN次。第5页/共77页64.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）第6页/共77页7 4.2.1 直接计算DFT的问题及改进的途径 k=0, 1, , N-110)()(NnknNWnxkX10)(1)(NkknNWkXNnxn=0, 1, , N-1DFT及IDFT的定义第7页/共77页8可见，DFT 与 IDFT 的计算成本基本相同。直接计算N点DFT 时：对应一个k需要N次复数乘和（N-1）次复数加；对所有N个k值，则需要 N复数乘和N (N-1)次复数加 。其中:一次复数乘需要4次实数乘和2次实数加方能完成。当N较大时，运算量很大。第8页/共77页9 例如：当N=8时，DFT需要64次复数乘；当N=1024时，DFT所需复数乘为1048576次，即一百多万次复数乘。为了减少运算次数，改进计算的方法：1）利用旋转因子的对称性、周期性和可约性； 旋转因子（twiddle factor）2）使长序列变短。NjNeW/2第9页/共77页104.2.2 基2时域抽取法原理 （库利-图基算法）设序列x（n）的长度为N，且M为自然数N-point DFT,N is evenMN2)(log2NM 1,.,1 , 0,)()(10NkWnxkXNnknN第10页/共77页11将其一分为二，分成偶数和奇数序列项（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-1第11页/共77页12偶数序列 the range: 02nN-2 (N is even) 0n(N/2)-1奇数序列 the range: 12n+1N-1 (N is even) 02nN-2 0n(N/2)-1第12页/共77页13 奇数偶数nknNnknNWnxWnxkX)()()(120)12(120)2() 12()2(NrrkNNrrkNWrxWrx1202212021)()(NrkrNkNNrkrNWrxWWrx则x（n）的DFT：第13页/共77页14由于所以krNkrNjkrNjkrNWeeW2/2/2222/2 1/2 112/2/20012( )( )( )( )( )NNkrkkrNNNrrkNX kx r WWx r WX kW Xkk=0,1,N-1第14页/共77页15上式说明,按n的奇偶性将 分解为两个N/2长的序列 和 ,则N点DFT可分解为两个N/2点的DFT来计算.)(nx)(1nx)(2nx第15页/共77页16 其中 N/2-point DFT: k=0,1,N/2-1 )()()(112/02/11rxDFTWrxkXNrkrN12/022/22)()()(NrkrNrxDFTWrxkXk=0,1,N/2-1第16页/共77页17因此，可写出两个N/2点的方程:)()()(21kXWkXkXkN)2()2()2(2)2(1NkXWNkXNkXNkN12,.,1 , 0Nk第17页/共77页18 而120)2(2)()2(11NrrNkNWrxNkXkNNkNWW222)()()2(1112021kXWrxNkXNrkrN第18页/共77页19同理而1222jNNjNNeeW)()2(22kXNkX第19页/共77页20所以)()()2()()()(2121kXWkXNkXkXWkXkXkNkN12,.1 , 0Nk第20页/共77页21表示上述算法可用蝶形结（ butterfly）kNWX X1 1( (k k) )X X2 2( (k k) )X(k)X(k+N/2)()(21kXWkXkN)()(21kXWkXkN第21页/共77页22 Example : 如N=4 x(n)=x0, x1, x2, x3 even: x0, x2 even: x0 odd: x2 odd: x1, x3 even: x1 odd: x3第22页/共77页23 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点的序列 第23页/共77页24x(0)x(1)x(2)x(3)X1(0)X1(1)X2(0)X2(1)X(0)X(1)X(2)X(3)4点点DFT2点点DFT the butterfly（蝶形运算）04W14W02W02W1点序列的DFT就是序列本身，即不用计算 -1-1-1-1第24页/共77页25 如N4，则将 x1(r) 再按r的奇偶进一步分解成两个N/4点长的子序列：) 12()()2()(1413lxlxlxlx14,.,1 , 0Nl第25页/共77页2614/0)12(2/114/022/11) 12()2()(NllkNNlklNWlxWlxkX12,.,1 , 0Nk)()(42/3kXWkXkN14/04/42/14/04/3)()(NlklNkNNlklNWlxWWlx第26页/共77页27 其中)()()(314/04/33lxDFTWlxkXNlklN)()()(414/04/44lxDFTWlxkXNlklN14,.,1 , 0Nk第27页/共77页28由X3（k）和X4（k）的周期性和 的对称性，得kNNkNWW2/4/2/)()()4()()()(42/3142/31kXWkXNkXkXWkXkXkNkN14,.,1 , 0Nk第28页/共77页29同理，得)()()4()()()(62/5262/52kXWkXNkXkXWkXkXkNkN14,.,1 , 0Nk第29页/共77页30 其中) 12()()2()(2625lxlxlxlx)()()(514/04/55lxDFTWlXkXNlklN)()()(614/04/66lxDFTWlXkXNlklN第30页/共77页31 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 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)02W02W02W02W04W04W14W14W08W28W18W38W第31页/共77页32 4.2.3 IDFT的高效算法10)(1)()(NkknNWkXNkXIDFTnx这样IFFT可与FFT用同一子程序实现。*10*)(1NkknNWkXN*)(1kXDFTN第32页/共77页33IDFT的运算规律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 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)*02W02W02W02W04W04W14W14W08W28W18W38W1/N第33页/共77页34 求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/N第34页/共77页35 Example:Determine the x(n) by IDFT.X=5, -1, 1, -1 Solution: n=0,1,2,3*304*)(41)(1)(kknkXWXDFTNnx第35页/共77页361)(41)(41)0(3210*3004XXXXXWxkk1)15(41)(41) 1 (*304jjXWxkkk2)(41)2(30*24kkkXWx第36页/共77页37 所以 x(n)=1, 1, 2, 1 1)(41)3(*3034kkkXWx第37页/共77页38%the program in MATLAB: X=input(Please input X=); N=length(X); x=fft(conj(X),N); x=conj(x/N)第38页/共77页39 Example :用FFT算法计算下面信号的 8点DFT； x(n)=4 3 2 0 1 2 3 1然后，再用 IFFT恢复原信号。第39页/共77页40 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+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 2第40页/共77页41 IFFT(X)=1/NFFT(X*)*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 merging4 4-2-j6-2-j61212-2+j6-2+j65-j-j5-j-j5+j-j5+j-j5-j+j5-j+j5+j+j5+j+j4 41212-2-j6-2-j6-2+j6-2+j65-j-j5-j-j5-j+j5-j+j5+j-j5+j-j5+j+j5+j+j1616-8-8-4-4-12j-12j10-j210-j2-j2-j210+j210+j2-j2-j21212-20-2020204 42020-2 (1+j)-2 (1+j)-4j-4j2 (1-j)2 (1-j)3232-24-2416160 0-8-8-16-1624248 8 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 204W04W14W14W08W18W28W38W乘以1/N=1/8得原序列第41页/共77页42 4.2.4 FFT的计算成本最简单的FFT 是Cooley-Tukey法因为NMNM2log2N点DFT有M级蝶形运算，每一级都有N/2个蝶形运算结构；每个蝶形运算结构都有1次复数乘和2次复数加；第42页/共77页 8442222111111118-DFT8-DFT4-DFT4-DFT2-DFT2-DFT1-DFT1-DFTstage 1stage 1stage 2stage 2stage 3stage 3第43页/共77页44所以，M级蝶形运算有复数乘M级蝶形运算有复数加NNMNCM2log22)2(NNMNCA2log)2(第44页/共77页45直接DFT的复数乘和复数加均近似为 。当N1时，所以DIT-FFT大大减少了运算次数。NNN22log)2(2N第45页/共77页46 Example : for N=8 Solution : M= =3 stages（三级） for FFT, the total cost is （FFT的总计算成本是） MN/2=3 8/2=12 complex multiplications （复数乘）8log2第46页/共77页47 for directly DFT, the total cost is（直接计算DFT的成本是） N=8=64 (complex multiplications) The ratio is: （比值是） 实际上，当N=2048时，直接运算与 FFT算法的计算量之比为372.433. 51264第47页/共77页48 4.2.5 DIT-FFT的运算规律及编程思想1、原位运算：利用同一单元存储蝶形计算的输入、输出数据。每个蝶形的输入和输出均为相同位数。原位运算可节省大量内存，因而硬件成本低。2、旋转因子的变化规律： N点DFT，共有M级蝶形运算，每级有N/2个蝶形。第48页/共77页49 称为旋转因子，p称为旋转因子的指数。 为了编写程序，应找出旋转因子 与运算 级数的关系。 设L=1，2，M，表示从左到右的运算 级数，第L级有 个不同的旋转因子，pNW12LpNW第49页/共77页50对于 ，各级的旋转因子表示为L=1时，L=2时，L=3时，823N0,24/JWWWJJNpNL1 , 0,22/JWWWJJNpNL3 , 2 , 1 , 0,2JWWWJJNpNL第50页/共77页51对于 的一般情况，第L级的旋转因子为MN2JpNLWW212,.,2 , 1 , 01LJ，由于MLMLMLN2222第51页/共77页52所以通过上式，可以计算第L级运算的旋转因子。LMMLJNJNpNWWW2212,.,1 , 01LJLMJp2第52页/共77页53 3、蝶形运算规律设序列x（n）经过时域倒序存入数组X如蝶形运算的两个输入数据相距B个点，应用原位运算，可得：pNLLLpNLLLWBJXJXBJXWBJXJXJX)()()()()()(1111MLJJpLLM,.,1; 12,.,1 , 0;21第53页/共77页54 4、编程思想及程序框图其它运算规律：第L级中，每个蝶形的两个输入数据相距 个点；同一旋转因子对应着间隔为 点的 个蝶形（对应第几组蝶形）12LBLM2L2第54页/共77页55 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 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)02W02W02W02W04W04W14W14W08W28W18W38W第55页/共77页56编程的运算方法：从输入端（第1级）开始，逐级进行，共进行M级运算。在进行第L级运算时，依次求出 个不同的旋转因子，然后计算其对应相同的旋转因子的 个蝶形。可用三重循环程序实现DIT-FFT运算。12LLM2第56页/共77页57开始开始输入输入x(n),M倒序倒序L=1.MJ=0,B-1输出输出结束结束MN212LBLMJP2LNJk2 , 1,pNpNWBkXkXBkXWBkXkXkX)()()()()()(（1）（1）（2）（2）(3)第57页/共77页585、序列的倒置（bit reversal）倒序是有规律的。由于 ，所以顺序数可用M位二进制数（ ）表示。MN20121. nnnnMM第58页/共77页59用硬件电路和汇编语言程序产生倒序很容易，用高级语言倒序的规律为：倒序数是在M位二进制数最高位加1，逢2向右进位。第59页/共77页604.2.6 频率抽取法FFT(DIF-FFT)设序列x（n）长度为 ，将其前后对半分开，得：MN212/12/010)()()()()(NNnknNNnknNNnknNWnxWnxWnxnxDFTkX第60页/共77页61式中12/0)2/(12/0)2()(NnNnkNNnknNWNnxWnx奇数，偶数kkWkkNN1, 1) 1(2/knNkNNNnWNnxWnx)2()(2/12/0第61页/共77页62再将X（k）分解成偶数组和奇数组 k为偶数时：/2 120/2 1/20(2 ) ( )()2 ( )()2NrnNnNrnNnNXrx nx nWNx nx nW第62页/共77页63k为奇数时：rnNnNNnnrNNnWWNnxnxWNnxnxrX2/12/0)12(12/0)2()()2()() 12(第63页/共77页64令12,.,1 , 0,)2()()()2()()(21NnWNnxnxnxNnxnxnxnN第64页/共77页65得12/02/212/02/1)() 12()()2(NnrnNNnrnNWnxrXWnxrX12,.,1 , 0Nr第65页/共77页66 DIF-FFT蝶形运算流图x x(n)(n)x x(n+N/2)(n+N/2)x x1 1(n)=x(n)+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-+第66页/共77页67N=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)0NW1NW2NW2NW2NW3NW0NW0NW0NW0NW0NW0NW第67页/共77页68注意：DIT-FFT和DIF-FFT的算法流图不 是唯一的。 其变形运算流图见P108。第68页/共77页69 4.3 进一步减少运算量的措施以程序的复杂度换取计算量的进一步提高。4.3.1 多类蝶形单元运算第一级旋转因子可简化：第二级旋转因子可简化：称为无关紧要的旋转因子。10NWjWWNNN4/01和第69页/共77页70其复数乘法次数可减少为： CM（2）=（M-2）*N/2当L=3时，第三级蝶形有两个无关紧要旋转因子，同一因子对应2(M-L)=N/2L级蝴蝶结,所以第三级共有(书110页)32*24NN第70页/共77页71依次类推，从L=3到L=M共可减少复数乘法次数为：DIT-FFT的复数乘法次数为：22231MLLNN2) 3(2)22()2(2)2(MNNMNCM第71页/共77页72另外，也可用实数乘法减少计算量。包含所有旋转因子称为一类蝶形单元运算；去掉 为二类；去掉 为三类；依次类推，称为多类蝶形单元运算。N=4096时，三类与一类比，仅75%。2/2)1 (8/jWNN1rNWjWrN第72页/共77页734.3.2 旋转因子的生成 直接查表，提高速度，多占内存。4.3.3 实序列的FFT算法 两个实序列，构造序列y（n） 1212( )( )( ) ( )( )( )( )epopy nx njx nDFT x nYkDFT jx nYk第73页/共77页74由于x（n）为实序列，所以X（k）具有共轭对称性：)()()(21kXWkXkXkN12,.,1 , 0),()(*NkkXkNX第74页/共77页75 4.4 分裂基FFT算法任意基：基较大，则程序或硬件将很复杂，基大于8无意义。分裂基：采用基2和基4的混合算法，可有效提高运算速度。第75页/共77页76第四章作业 1.课本127页第1题,第4题 2.请用学号的后8位做N=8点的DIT-FFT,DIF-FFT,要求画出求解的各步骤.第76页/共77页77感谢您的观看。第77页/共77页