第四章 快速傅里叶变换n

1、2n直接计算直接计算DFT的问题及改进的途径的问题及改进的途径 3 4 设复序列设复序列x(n) 长度为长度为N点，其点，其DFT为为10( )( )NnkNnX kx n Wk=0，N-1 （1）计算一个）计算一个X(k) 值的运算量值的运算量复数乘法次数：复数乘法次数： N复数加法次数：复数加法次数： N15（2）计算全部）计算全部N个个X(k) 值的运算量值的运算量复数乘法次数：复数乘法次数： N2复数加法次数：复数加法次数： N(N1)（3）对应的实数运算量）对应的实数运算量1100( )( )Re ( )Im ( )ReImNNnknknkNNNnnX kx n Wx njx nWj

2、W10Re ( ) ReIm ( ) ImNnknkNNnx nWx nWRe ( ) ImIm ( ) RenknkNNjx nWx nW6一次复数乘法一次复数乘法： 4次实数乘法次实数乘法 2次实数加法次实数加法 一个一个X(k)数值数值 ： 4N次实数乘法次实数乘法2N+2(N-1)= 2(2N-1)次实数加法次实数加法 所以所以 整个整个N点点DFT运算共需要：运算共需要：N2(2N-1)= 2N(2N-1)实数乘法次数：实数乘法次数：4 N2实数加法次数：实数加法次数：7N点点DFT的复数乘法次数举例的复数乘法次数举例NN2NN22464404941612816384864256 6

3、5 536 16256512 262 144 3210281024 1 048 576 结论结论：当：当N很大时，其运算量很大，对实时性很强的信号处理很大时，其运算量很大，对实时性很强的信号处理来说，要求计算速度快，因此需要改进来说，要求计算速度快，因此需要改进DFT的计算方法，以的计算方法，以大大减少运算次数。大大减少运算次数。 8 nkNW主要原理是利用系数主要原理是利用系数 的以下特性对的以下特性对DFT进行分解：进行分解： nkNW（1）对称性）对称性 ()nkNW()k N nNW（2）周期性）周期性 ()()nN kn kNnkNNNWWW（3）可约性）可约性 mnknkmNNWW

4、/nknk mNNmWW另外，另外，12/NNWkNNkNWW)2/(9101(2 )( )xrx r设设N2L，将，将x(n)按按 n 的奇偶分为两组：的奇偶分为两组： 2(21)( )xrx rr =0，1， 12N10( ) ( )( )NnkNnX kDFT x nx n W则则1010)()(NnnnkNNnnnkNWnxWnx为奇数为偶数11120)12(1202) 12()2(NrkrNNrrkNWrxWrx1202212021)()(NrrkNkNNrrkNWrxWWrx)()(21kXWkXkN1010)()(NnnnkNNnnnkNWnxWnx为奇数为偶数式中，式中，X1(

5、k)和和X2(k)分别是分别是x1(n)和和x2(n)的的N/2的的DFT。另外，式中另外，式中k的取值范围是：的取值范围是：0，1， ，N/21 。12因此，因此， 只能计算出只能计算出X(k)的前一半值。的前一半值。12( )( )( )kNX kXkW Xk后一半后一半X(k) 值，值， N/2 ， N/2 1， ，N ？rkNW2(2)2r NkNW利用利用可得到可得到 1()2NXk2 1(2)120( )Nr NkNrx r W2 1120( )NrkNrx r W)(1kX同理可得同理可得22()( )2NXkXk13考虑到考虑到 kNkNNNkNNWWWW2)2(因此可得后半部

6、分因此可得后半部分X(k) )2()2()2(221NkXWNkXNkXNkN12( )( )( )kNX kX kW Xk由前半部分由前半部分X(k) )()(21kXWkXkNk=0，1， ，N/21k=0，1， ，N/2114222 N2N蝶形运算式蝶形运算式蝶形运算信号蝶形运算信号流图符号流图符号 因此，只要求出因此，只要求出2个个N/2点的点的DFT，即，即X1(k)和和X2(k)，再经过，再经过蝶形运算就可求出全部蝶形运算就可求出全部X(k)的值。运算量由次复数乘减为的值。运算量由次复数乘减为运算量减少一半。运算量减少一半。)()()2()()()(2121kXWkXNkXkXWk

7、XkXkNkN150NW1NW2NW3NW以以N=8为例，为例，分解为分解为2个个4点的点的DFT，然后做然后做8/2=4次蝶形运次蝶形运算即可求出所有算即可求出所有8点点X(k)的值。的值。16复数乘法次数： N2复数加法次数： N(N1)复数乘法次数： 2*(N/2)2+N/2=N2/2+N/2复数加法次数： 2*(N/2)(N/21)+2*N/2=N2/2nN点 17 由于N2L，因而N/2仍是偶数 ，可以进一步把每个N/2点子序列再按其奇偶部分分解为两个N/4点的子序列。 以N/2点序列x1(r)为例 1314(2 )( )0,1,1(21)( )4xlx lNlxlxl则有 rkNN

8、rWrxkX212011)()(klNNllkNNlWlxWlx)12(21401221401) 12()2(lkNNlkNlkNNlWlxWWlx41404241403)()()()(42/3kXWkXkNk=0,1, 14N18且且13/24( )( )4kNNXkXkWXkk=0,1, 14N由此可见，一个由此可见，一个N/2点点DFT可分解成两个可分解成两个N/4点点DFT。 同理，也可对同理，也可对x2(n)进行同样的分解，求出进行同样的分解，求出X2(k)。 192013/40( )lkNlx l W02(0)(4)xW x0(0)(4)NxW x 对此例对此例N=8，最后剩下的是

9、，最后剩下的是4个个N/4= 2点的点的DFT，2点点DFT也可以由蝶形运算来完成。以也可以由蝶形运算来完成。以X3(k)为例。为例。 /4 133/40( )( )NlkNlXkx l Wk=0, 1即即03323(0)(0)(1)XxW x13323(1)(0)(1)XxW x12(0)(4)xW x0(0)(4)NxW x这说明，这说明，N=2M的的DFT可全部由蝶形运算来完成。可全部由蝶形运算来完成。21N=8按时间抽取法按时间抽取法FFT信号流图信号流图 22由按时间抽取法FFT的信号流图可知，当N=2L时，共有 级蝶形运算；每级都由 个蝶形运算组成，而每个蝶形有 次复乘、 次复加，

10、因此每级运算都需 次复乘和 次复加。 LN/2 N/2 12N23这样这样 级运算总共需要：级运算总共需要： L复数乘法： NNLN2log22复数加法： NNLN2log直接直接DFT算法运算量算法运算量 复数乘法： 复数加法： N2N(N1)直接计算直接计算DFT与与FFT算法的计算量之比为算法的计算量之比为MNNNNNM222log2log224NN2计算量之比M NN2计算量之比M 2414.012816 38444836.641644.025665 5361 02464.0864125.4512262 1442 304113.816256328.010241 048 5765 120

11、204.83210288012.820484 194 30411 264372.464404919221.4NN2log2NN2log225n序列的逆序排列n同址运算（原位运算）n蝶形运算两节点间的距离n 的确定rNW26返回返回27)(01221)()(BINMMDECnnnnnn 由于由于 x(n) 被反复地按奇、偶分组，所以流图输被反复地按奇、偶分组，所以流图输入端的入端的排列不再是顺序的，但仍有规律可循：排列不再是顺序的，但仍有规律可循： 因为因为 N=2M ， 对于任意对于任意 n（0n N-1)，可以用可以用M个二进制个二进制码表示为：码表示为：10,01221nnnnnMM n

12、反复按奇、偶分解时，即按二进制码的反复按奇、偶分解时，即按二进制码的“0” “1” 分解。分解。n序列的逆序排列2829自然顺序自然顺序 n二进制数二进制数倒位序二进制数倒位序二进制数倒位序顺序数倒位序顺序数0000000010011004201001023011110641000011510110156110011371111117 n30为目前计算机已具有的功能为目前计算机已具有的功能31返回返回32 某一列任何两个节点某一列任何两个节点k 和和j 的节点变量进行蝶形运算后，的节点变量进行蝶形运算后，得到结果为下一列得到结果为下一列k、j两节点的节点变量，而和其他节点变量两节点的节点变量，

13、而和其他节点变量无关。这种原位运算结构可以节省存储单元，降低设备成本。无关。这种原位运算结构可以节省存储单元，降低设备成本。)(kX)2(NkX)(1kX)(2kXkNW运算前运算前运算后运算后)2(NkA)(kA)2(NkA)(kA例例n同址运算（原位运算）33返回返回34以以N=8为例：为例：第一级蝶形，距离为：第一级蝶形，距离为：第二级蝶形，距离为：第二级蝶形，距离为：第三级蝶形，距离为：第三级蝶形，距离为：规律规律：对于共：对于共L级的蝶形而言，其级的蝶形而言，其m级蝶形运算的节级蝶形运算的节 点间的距离为点间的距离为12412mn蝶形运算两节点间的距离 35rNWrNWrNW返回返回

14、rNWrNW36rNW以N=8为例：0,10224/lWWWWmllNrNm时，1 , 0,2422/lWWWWmlllNrNm时，3 , 2 , 1 , 0,382lWWWWmlllNrNm时，级：第LNM,212 , 2 , 1 , 0,12LlrNlWWLMLMLMLN2222LMLMMLMLlNlNjlNjlNrNWeeWW222222rNWn 的确定 37n算法原理算法原理 再把输出再把输出X(k)按按k的奇偶分组的奇偶分组先把输入按先把输入按n的顺序分成前后两半的顺序分成前后两半设序列长度为设序列长度为N=2L，L为整数为整数 前半子序列前半子序列x(n) 后半子序列后半子序列)2

15、(Nnx 0n12N0n12N3810( )( )NnkNnX kx n W由由DFT定义得定义得/2 110/2( )( )NNnknkNNnn Nx n Wx n W12/0)2(12/0)2()(NnkNnNNnnkNWNnxWnx12/02)2()(NnnkNkNNWWNnxnxk=0，1， ，N39由于由于 1222jNNjNNeeW/2 120( )( )()2NNknkNNnNX kx nx nWW所以所以kkNNW) 1(2则则 12/0)2() 1()()(NnnkNkWNnxnxkXk=0，1， ，N40然后按然后按k的奇偶可将的奇偶可将X(k)分为两部分分为两部分 221

16、krkrr=0，1， ，12N则式则式 12/0)2() 1()()(NnnkNkWNnxnxkX可转化为可转化为 nrNNnWNnxnxrX212/0)2()()2(12/02/)2()(NnnrNWNnxnx)12(12/0)2()() 12(rnNNnWNnxnxrXnrNnNNnWWNnxnx212/0)2()(41/2 1/20(2 )( )()2NnrNnNXrx nx nW令令 nNWNnxnxnxNnxnxnx)2()()()2()()(21n=0，1， ，12N代入代入 /2 120(21) ( )()2NnnrNNnNXrx nx nWWnrNNnnrNNnWnxrWnxr

17、2120221201)() 12()()2(r=0，1， ，12N可得可得为为2个个N/2点的点的DFT，合起来正好是，合起来正好是N点点X(k)的值。的值。42nNWNnxnxnxNnxnxnx)2()()()2()()(21将将称为蝶形运算称为蝶形运算与时间抽选基与时间抽选基2FFT算法中的蝶形运算符号略有不同。算法中的蝶形运算符号略有不同。43例例 按频率抽取，将按频率抽取，将N点点DFT分解为两个分解为两个N/2点点DFT的组合的组合(N=8) 44 与时间抽取法的推导过程一样，由于与时间抽取法的推导过程一样，由于 N=2L，N/2仍然是仍然是一个偶数，因而可以将每个一个偶数，因而可以

18、将每个N/2点点DFT的输出再分解为偶数组的输出再分解为偶数组与奇数组，这就将与奇数组，这就将N/2点点DFT进一步分解为两个进一步分解为两个N/4点点DFT。 N=845按频率抽取的按频率抽取的FFT（N=8）信号流图）信号流图 返回返回 1 10NW2NWx(0)x(1)x(2)x(3) 1 1x(4)x(5)x(6)x(7)0NW1NW2NW3NWX(0)X(4)X(2)X(6)X(1)X(5)X(3)X(7)0NW2NW 1 1 1 1 1 1 1 10NW0NW0NW0NW46 由由P31图图 与与P42图图 相比较，初看起来，相比较，初看起来，DIF法与法与DIT法的区别是：法的区

19、别是： DIF输入是自然顺序，输入是自然顺序， 输出是倒位序的，输出是倒位序的， 这与这与DIT法正好相反。法正好相反。 但但这不是实质性的区别，因为这不是实质性的区别，因为DIF法与法与DIT法一样，都可将输入或输出法一样，都可将输入或输出进行重排，使二者的输入或输出顺序变成自然顺序或倒位序顺序。进行重排，使二者的输入或输出顺序变成自然顺序或倒位序顺序。 DIF的基本蝶形（的基本蝶形（ P42图图 ）与）与DIT的基本蝶形（的基本蝶形（ P31 ）则有所不同，）则有所不同，这才是实质的不同，这才是实质的不同，DIF的复数乘法只出现在减法之后，的复数乘法只出现在减法之后，DIT则是先则是先作复

20、乘后再作加减法。作复乘后再作加减法。 47 但是，但是，DIF与与DIT就运算量来说则是相同的，即都有就运算量来说则是相同的，即都有M级级（列）运算，（列）运算， 每级运算需每级运算需N/2个蝶形运算来完成，总共需要个蝶形运算来完成，总共需要 次次复乘与复乘与 次复加次复加 ，DIF法与法与DIT法都可进行原位运法都可进行原位运算。算。 频率抽取频率抽取FFT算法的输入是自然顺序，输出是倒位序的。算法的输入是自然顺序，输出是倒位序的。 因此因此运算完毕后，运算完毕后， 要通过变址计算将倒位序转换成自然位序，然后再输要通过变址计算将倒位序转换成自然位序，然后再输出。出。 转换方法与时间抽取法相同

21、。转换方法与时间抽取法相同。 NN2log2NNLN2log48 如果将如果将DIT的基本蝶形加以转置，就得到的基本蝶形加以转置，就得到DIF的基本蝶形的基本蝶形; 反过来反过来, 将将DIF的基本蝶形加以转置，的基本蝶形加以转置， 就得到就得到DIT的基本蝶形，的基本蝶形， 因而因而DIT法与法与DIF法的基本蝶形是互为转置的。法的基本蝶形是互为转置的。 按照转置定按照转置定理，理， 两个流图的输入两个流图的输入-输出特性必然相同。转置就是将流图输出特性必然相同。转置就是将流图的所有支路方向都反向，并且交换输入与输出，但节点变量的所有支路方向都反向，并且交换输入与输出，但节点变量值不交换。值

22、不交换。 因此可以说，有多少种按时间抽取的因此可以说，有多少种按时间抽取的FFT流图流图就存在多少种按频率抽取的就存在多少种按频率抽取的FFT流图。频率抽取法与时间抽流图。频率抽取法与时间抽取法是两种等价的取法是两种等价的FFT运算。运算。 1 10NW2NWx(0)x(1)x(2)x(3) 1 1x(4)x(5)x(6)x(7)0NW1NW2NW3NWX(0)X(4)X(2)X(6)X(1)X(5)X(3)X(7)0NW2NW 1 1 1 1 1 1 1 10NW0NW0NW0NW49n频率抽取法输入是自然顺序，输出是倒位序的；时间抽取法正好相反。n频率抽取法的基本蝶形与时间抽取法的基本蝶形

23、有所不同。n频率抽取法运算量与时间抽取法相同。n频率抽取法与时间抽取法的基本蝶形是互为转置的。 504.5 N为复合数的为复合数的FFT算法算法 上面讨论的是序列的点数上面讨论的是序列的点数N为为2的幂次的幂次(即即N=2M)情况下，按时间抽情况下，按时间抽取和按频率抽取的基取和按频率抽取的基 -2FFT算法的基本原理。这种基算法的基本原理。这种基-2FFT算法在实算法在实际中使用得最多，因为它的程序简单，效率高，使用方便。但实际上际中使用得最多，因为它的程序简单，效率高，使用方便。但实际上无法保证总是处理长度为无法保证总是处理长度为2的整数幂次的序列。若不满足的整数幂次的序列。若不满足N=2

24、M，可将，可将x(n)增补一些零值点，以使增补一些零值点，以使N增长到最邻近的一个增长到最邻近的一个2M数值。数值。 有限长序有限长序列补零之后，并不影响其频谱列补零之后，并不影响其频谱X(ej)，只不过其频谱的抽样点数增加，只不过其频谱的抽样点数增加了，所造成的结果是增加了计算量而已。了，所造成的结果是增加了计算量而已。 但是，有时计算量增加太多，但是，有时计算量增加太多，浪费较大。浪费较大。 51 例如例如, x(n)的点数的点数N=300，则须补到，则须补到N=29=512，要补，要补212个零值点，个零值点，因而人们才研究因而人们才研究N2M时的时的FFT算法。算法。 若若N是一个复合

25、数，即它可以分解成一些因子的乘积，则可以用是一个复合数，即它可以分解成一些因子的乘积，则可以用FFT的一般算法，即混合基的一般算法，即混合基FFT算法，如库利算法，如库利-图基（图基（CooleyTukey）算法，而基）算法，而基 -2 算法只是这种一般算法的特例。算法只是这种一般算法的特例。 这里，这里， 我们我们不作详细介绍，不作详细介绍， 感兴趣的话可参考文献感兴趣的话可参考文献1。 总之，不管采用什么方法，总之，不管采用什么方法，计算计算DFT的高效算法是把计算长度为的高效算法是把计算长度为N的序列的的序列的DFT逐次分解成计算长度较短序列的逐次分解成计算长度较短序列的DFT。这是很多

26、高效。这是很多高效算法的标准方法和基本原理。算法的标准方法和基本原理。 52 例例4-5 已知信号已知信号x(t)=0.15 sin(2f1t)+sin(2f2t)-0.1sin(2f3t)，其中，其中, f1=1Hz， f2=2 Hz，f3=3Hz。 x(t)包含三个正弦波，但从时域波形图包含三个正弦波，但从时域波形图3-26（a）来看，似乎是一个正弦信号，很难看到小信号的存在，）来看，似乎是一个正弦信号，很难看到小信号的存在， 因为因为它被大信号所掩盖。它被大信号所掩盖。 取取fs=32 Hz作频谱分析。作频谱分析。 解解 因因fs=32 Hz，故，故 nnnnTxnx326sin1 .

27、0324sin322sin15. 0)()(53图图 已知信号的时域图和幅频图已知信号的时域图和幅频图 2 10120102030(a)(b)0051015816nkx(n)|X(k)|54 该信号为周期信号，其周期为该信号为周期信号，其周期为N=32。现对。现对x(n)作作32点的离散傅点的离散傅里叶变换里叶变换(DFT)，其幅度特性，其幅度特性|X(k)|如如下页下页图所示。图中仅给出了图所示。图中仅给出了k=0, 1, , 15的结果。的结果。k=16, 17, , 32的结果可由的结果可由|X(N-k)|=|X(k)|得出。因得出。因N=32, 故频率分辨率故频率分辨率F=fs/N=1

28、Hz。 由由下页下页图图 (b) 可知可知, k=1, 2, 3所对应所对应的频谱即为频率的频谱即为频率f1=1 Hz，f2=2 Hz, f3=3Hz的正弦波所对应的频谱。的正弦波所对应的频谱。幅频图如幅频图如下页下页图图(b), 该图中小信号成分清楚显示出来。该图中小信号成分清楚显示出来。 可见小信号可见小信号成分在时域中很难辨识而在频域中容易识别。成分在时域中很难辨识而在频域中容易识别。 幅度：幅度：2Hz，16*(2/fs)=1 1Hz，2.4*(2/fs)=0.15 3Hz，1.6*(2/fs)=0.155 如果我们事先不知道信号的最高频率，可以根据信号的时域波形如果我们事先不知道信号

29、的最高频率，可以根据信号的时域波形图来估计它。例如，图来估计它。例如， 某信号的波形如图某信号的波形如图 3-23 所示。所示。 先找出相邻的波先找出相邻的波峰与波谷之间的距离，如图中峰与波谷之间的距离，如图中t1，t2，t3，t4。 然后，选出其中最小的然后，选出其中最小的一个如一个如t4。这里。这里, t4可能就是由信号的最高频率分量形成的。可能就是由信号的最高频率分量形成的。 峰与谷之峰与谷之间的距离就是周期的一半。间的距离就是周期的一半。 因此，最高频率为因此，最高频率为 )(214Hztfh知道道fh后就能确定采样频率后就能确定采样频率 hsff256估算信号最高频率估算信号最高频率

30、fh tt3t4ot1t2x(t)57MN 2IDFT公式公式 10)(1NknkNWkXNkXIDFTnxDFT公式公式 nkNNnWnxnxDFTkX10)()()(比较可以看出，比较可以看出， nkNWnkNWMN211IDFT多出多出M个个1/2可分解到可分解到M级蝶形运算中。级蝶形运算中。585910)(1)(NknkNWkXNkXIDFT10)(1NknkNWkXN10)(1NknkNWkXN1( )( )IFFT X kFFT XkN( )X k 求共轭( )XkFFT 求( )FFT Xk( )FFT XkN 除以( )x n 求共轭60n用用FFT进行谱分析的进行谱分析的Ma

31、tlab实现实现n用用CZT进行谱分析的进行谱分析的Matlab实现实现n在在Matlab中使用的线性调频中使用的线性调频z变换函数为变换函数为czt，其，其调用格式为调用格式为nX= czt(x, M, W, A)n其中，其中，x是待变换的时域信号是待变换的时域信号x(n)，其长度为，其长度为N，M是变是变换的长度，换的长度，W确定变换的步长，确定变换的步长，A确定变换的起点。若确定变换的起点。若M= N，A= 1，则，则CZT变成变成DFT。61例例4.6 设模拟信号设模拟信号 ，以，以 t= 0.01n (n=0: N-1) 进行取样，试用进行取样，试用fft函数对其做频谱分析。函数对其

32、做频谱分析。N分别为：分别为：(1) N=45；(2) N=50；(3) N=55；(2) N=60。 ( )2sin(4)5cos(8)x ttt程序清单如下程序清单如下 %计算计算N=45的的FFT并绘出其幅频曲线并绘出其幅频曲线N=45;n=0:N-1;t=0.01*n;q=n*2*pi/N;x=2*sin(4*pi*t)+5*cos(8*pi*t);y=fft(x,N);figure(1)subplot(2,2,1)plot(q,abs(y)title(FFT N=45)62%计算计算N=50的的FFT并绘出其幅频曲线并绘出其幅频曲线N=50;n=0:N-1;t=0.01*n;q=n*

33、2*pi/N;x=2*sin(4*pi*t)+5*cos(8*pi*t);y=fft(x,N);figure(1)subplot(2,2,2)plot(q,abs(y)title(FFT N=50)63%计算计算N=55的的FFT并绘出其幅频曲线并绘出其幅频曲线N=55;n=0:N-1;t=0.01*n;q=n*2*pi/N;x=2*sin(4*pi*t)+5*cos(8*pi*t);y=fft(x,N);figure(1)subplot(2,2,3)plot(q,abs(y)title(FFT N=55)64%计算计算N=60的的FFT并绘出其幅频曲线并绘出其幅频曲线N=60;n=0:N-1

34、;t=0.01*n;q=n*2*pi/N;x=2*sin(4*pi*t)+5*cos(8*pi*t);y=fft(x,N);figure(1)subplot(2,2,4)plot(q,abs(y)title(FFT N=60)6502468050100150FFT N=4502468050100150FFT N=5002468050100150FFT N=5502468050100150FFT N=60从图中可以看出，这几种情况下均有较好的精度。从图中可以看出，这几种情况下均有较好的精度。 66分析：由分析：由t=0.01n进行取样可得，采样频率进行取样可得，采样频率fs=100Hz。而连续。

35、而连续信号的最高模拟角频率为信号的最高模拟角频率为8 ，由，由2 f可得，最高频可得，最高频率为率为8 /2 =4Hz。因此，满足采样定理的要求。因此，满足采样定理的要求。 采样序列为采样序列为( )2cos(4)5cos(8)x nTnTn48( )2cos5cos100100 x nnn即即为周期序列，周期为周期序列，周期N=50。将程序中将程序中plot改为改为stem函数，则可以更清楚地看出频谱。函数，则可以更清楚地看出频谱。6702468050100150FFT N=4502468050100150FFT N=5002468050100150FFT N=550246805010015

36、0FFT N=6068 信号消噪信号消噪 假设信号在传输过程中，受到噪声的干扰，则在接收假设信号在传输过程中，受到噪声的干扰，则在接收端得到的信号由于受到噪声的干扰，信号将难以辨识。用端得到的信号由于受到噪声的干扰，信号将难以辨识。用FFT方法消噪就是对含噪信号的频谱进行处理，将噪声所方法消噪就是对含噪信号的频谱进行处理，将噪声所在频段的在频段的X(k)值全部置零后，再对处理后的值全部置零后，再对处理后的X(k)进行离散傅进行离散傅里叶反变换（里叶反变换（IFFT）可得原信号的近似结果。这种方法要）可得原信号的近似结果。这种方法要求噪声与信号的频谱不在同一频段，求噪声与信号的频谱不在同一频段，

37、 否则，否则， 将很难处理。将很难处理。 69 例例 将上述消噪原理用于语音消噪。图将上述消噪原理用于语音消噪。图4-30是一个实际例子，语音是一个实际例子，语音信号受到了强烈的啸叫噪声干扰，无法听清语音意思，信号受到了强烈的啸叫噪声干扰，无法听清语音意思， 如图如图4-30(a)， 信号淹没在噪声中（信噪比只有信号淹没在噪声中（信噪比只有-10dB）。试用）。试用FFT方法消噪。方法消噪。 先作先作FFT分析，得到其功率谱如图分析，得到其功率谱如图4-30(b)，可见在频率，可见在频率2.5 kHz附近有一极强分量。附近有一极强分量。 这就是啸叫噪声干扰。图这就是啸叫噪声干扰。图4-30(b

38、)中频率在中频率在30800Hz范围是语音信号。对频谱进行修正，去除噪声频段，即将大于范围是语音信号。对频谱进行修正，去除噪声频段，即将大于2.5 kHz部分的部分的X(k)值全部置为零，值全部置为零， 图图4-30(c)是去噪后的功率谱。再是去噪后的功率谱。再由反变换（由反变换（IFFT）重构信号得到原语音信号如图）重构信号得到原语音信号如图4-30(d)。这时信噪。这时信噪比为比为14dB，提高了，提高了24dB。 这就是早期的数字式录音音乐中所采用的这就是早期的数字式录音音乐中所采用的消噪方法。消噪方法。 70图 4-30 语音信号消噪过程信号淹没在啸叫噪声中； (b) 信号与噪声的功率

39、谱； (a) (c) 去噪后的功率谱； (d) 重构原语音信号71n离散哈特莱变换定义离散哈特莱变换定义 n 设设x(n),n=0,1,N-1，为一实序列，其，为一实序列，其DHT定义为定义为1, 2 , 1 , 0)2()()()(10NkknNcasnxnxDHTkXNnH式中，cas()=cos+sin。其逆变换(IDHT)为1, 2 , 1 , 0)2()(1)()(10NnknNcaskXNkXIDHTnxNkHH72已知已知x(n)的的DHT，则，则DFT可用下式求得可用下式求得： 11( )( )()( )()22HHHHX kXkXNkj XkXNk DHT与与DFT之间的关系

40、之间的关系 731）实数序列的）实数序列的FFT以上讨论的FFT算法都是复数运算,包括序列x(n)也认为是复数,但大多数场合,信号是实数序列,任何实数都可看成虚部为零的复数,例如,求某实信号x(n)的复谱,可认为是将实信号加上数值为零的虚部变成复信号(x(n)+j0),再用FFT求其离散傅里叶变换。这种作法很不经济,因为把实序列变成复序列,存储器要增加一倍,且计算机运行时,即使虚部为零,也要进行涉及虚部的运算,浪费了运算量。合理的解决方法是利用复数据FFT对实数据进行有效计算,下面介绍两种方法。74 （1）用 一个N点FFT同时计算两个N点实序列的DFT 设x (n)、y (n)是彼此独立的两

41、个N点实序列,且 X (k)=DFTx (n) ， Y (k)=DFTy(n) 则X (k)、Y(k)可通过一次FFT运算同时获得。首先将x (n)、y(n)分别当作一复序列的实部及虚部,令 g(n)=x (n)+jy(n)75 kjGkGkjYkjXkYkXkYkjYkjXkXkjYkXkGirriiririr kNGkGjkNGkGkNGkGkXngiirr212121ReDFT*通过通过g(n)的的FFT运算结果运算结果G(k),由上式可得到由上式可得到X (k) 的值。的值。 通过通过FFT运算可获得运算可获得g(n)的的DFT值值利用离散傅里叶变换的共轭对称性利用离散傅里叶变换的共轭

42、对称性76 kNGkGjkNGkGkNGkGkjYngjiirr212121ImDFT* kNGkGjkNGkGkYrrii2121作一次作一次点复序列的点复序列的FFT，再通过加、减法运算就可以，再通过加、减法运算就可以将将X(k)与与Y(k)分离出来。显然，这将使运算效率提高一倍。分离出来。显然，这将使运算效率提高一倍。77 设设x(n)是是2N点的实序列点的实序列,现人为地将现人为地将x(n)分为偶数组分为偶数组x1(n)和奇数组和奇数组x2(n) x1(n)=x(2n) n=0,1,N-1 x2(n)=x(2n+1) n=0,1,N-1然后将然后将x1(n)及及x2(n)组成一个复序列

43、：组成一个复序列： y(n)=x1(n)+jx2(n)（2） 用一个用一个N点的点的FFT运算获得一个运算获得一个2N点实序列的点实序列的DFT78通过通过N点点FFT运算可得到：运算可得到： Y(k)=X1(k)+jX2(k) ，N点点根据前面的讨论，得到根据前面的讨论，得到)()(2)()()(21)(*2*1kNYkYjkXkNYkYkX79 为求为求 2N 点点 x（n）所对应）所对应 X（k），需求出），需求出 X（k）与）与 X1（k）、）、X2（k）的关系）的关系 ： 10)12(21022120) 12()2()()(NnknNNnnkNNnnkWnxWnxWnxkX10210

44、) 12()2(NnnkNkNNnnkNWnxWWnx 1011022101011)12()()()2()()(NnnkNNnnkNNnNnnkNnkNWnxWnxkXWnxWnxkX801）由由x1(n)及及x2(n)组成复序列，经组成复序列，经FFT运算求得运算求得 Y(k)，2）利用共轭对称性求出）利用共轭对称性求出 X1(k)、X2(k)，3）最后利用上式求出）最后利用上式求出 X（k）, 达到用一个达到用一个N点的点的FFT计算计算一个一个2N点的实序列的点的实序列的DFT的目的。的目的。 X(k)=X1(k)+W2Nk X2(k)所以所以81 2) 线性卷积的线性卷积的FFT算法算

45、法 线性卷积是求离散系统响应的主要方法之一线性卷积是求离散系统响应的主要方法之一,许多重要应许多重要应用都建立在这一理论基础上用都建立在这一理论基础上,如卷积滤波等。如卷积滤波等。 以前曾讨论了用圆周卷积计算线性卷积的方法归纳如下以前曾讨论了用圆周卷积计算线性卷积的方法归纳如下: 将长为将长为N2的序列的序列x(n)延长到延长到L,补补L-N2个零，个零， 将长为将长为N1的序列的序列h(n)延长到延长到L,补补L-N1个零，个零， 如果如果LN1+N2-1,则圆周卷积与线性卷积相等则圆周卷积与线性卷积相等,此时此时,可用可用FFT计算线性卷积，方法如下计算线性卷积，方法如下: 82a.计算计

46、算X(k)=FFTx(n)b. 求求H(k)=FFTh(n)c. 求求Y(k)=H(k)X(k) k=0L-1d. 求求y(n)=IFFTY(k) n=0L-1 可见，只要进行二次可见，只要进行二次FFT,一次一次IFFT就可就可完成线性卷积计算。完成线性卷积计算。 计算表明计算表明,L32时时,上述计算线性卷积的上述计算线性卷积的方法比直接计算线卷积有明显的优越性方法比直接计算线卷积有明显的优越性,因此因此,也称上述循环卷积方法为快速卷积法。也称上述循环卷积方法为快速卷积法。83 上述结论适用于上述结论适用于 x(n)、h(n) 两序列长度比较接近两序列长度比较接近或相等的情况或相等的情况,

47、如果如果x(n)、h(n) 长度相差较多长度相差较多,例如例如, h(n) 为某滤波器的单位脉冲响应为某滤波器的单位脉冲响应,长度有限长度有限,用来处理用来处理一个很长的输入信号一个很长的输入信号 x(n),或者处理一个连续不断的或者处理一个连续不断的信号，按上述方法信号，按上述方法, h(n) 要补许多零再进行计算要补许多零再进行计算,计算计算量有很大的浪费，或者根本不能实现。量有很大的浪费，或者根本不能实现。 为了保持快速卷积法的优越性为了保持快速卷积法的优越性,可将可将 x(n) 分为许分为许多段多段,每段的长度与每段的长度与 h(n) 接近接近 , 处理方法有两种：处理方法有两种：84

48、h(n)x(n)85 假定假定 xi(n) 表示表示 x（n）序列的第）序列的第i段段 ：01) 1()()(22NiniNnxnxiiinxnx)()(iiiinynhnxnhnxny)()(*)()(*)()()(*)()(nhnxnyii其中其中 于是输出可分解为：于是输出可分解为： 则输入序列可表为：则输入序列可表为： 86 1）先对）先对 h（n）及）及 xi（n）补零，补到具有）补零，补到具有N点长度，点长度，N=N1+N2-1。 一般选一般选 N=2M。 2）用基）用基2 FFT计算计算 yi（n）=xi（n）*h（n）。）。 3）重叠部分相加构成最后的输出序列。）重叠部分相加构

49、成最后的输出序列。 iinyny)()(由于由于 yi（n）的长度为）的长度为N，而，而xi（n）的长度为）的长度为N2，因此相邻，因此相邻两段两段 yi（n）序列必然有）序列必然有N-N2=N1-1点发生重叠。点发生重叠。87计算步骤：计算步骤：a. 事先准备好滤波器参数事先准备好滤波器参数 H（k）=DFTh(n)，N点点b. 用用N点点FFT计算计算Xi（k）=DFTxi(n)c. Yi(k)=Xi(k)H(k)d. 用用N点点IFFT求求 yi（n）=IDFTYi(k)e. 将重叠部分相加将重叠部分相加8889 这种方法和第一种方法稍有不同，即将上面分段这种方法和第一种方法稍有不同，即

50、将上面分段序列中补零的部分不是补零，而是保留原来的输入序列中补零的部分不是补零，而是保留原来的输入序列值，这时，如利用序列值，这时，如利用DFT实现实现h(n)和和xi(n)的循环的循环卷积，则每段卷积结果中有卷积，则每段卷积结果中有N1-1个点不等于线性卷个点不等于线性卷积值需舍去。积值需舍去。 重叠保留法与重叠相加法的计算量差不多，但省重叠保留法与重叠相加法的计算量差不多，但省去了重叠相加法最后的相加运算。去了重叠相加法最后的相加运算。 9091 4）用）用FFT计算相关函数计算相关函数 相关的概念很重要，互相关运算广泛应用于信号相关的概念很重要，互相关运算广泛应用于信号分析与统计分析，如通过相关函数峰值的检测测量两个分析与统计分析，如通过相关函数峰值的检测测量两个信号的时延差等。信号的时延差等。 两个长为两个长为N的实离散时间序列的实离散时间序列 x（n）与）与y（n）的）的互相关函数定义为互相关函数定义为 1