中北大学数字信号处理原理及应用快速傅里叶变换

1、第第4 4章章 快速傅里叶变换快速傅里叶变换1/86第第4 4章章 快速傅里叶变换快速傅里叶变换第第4 4章章 快速傅里叶变换快速傅里叶变换2/86 DFT是信号分析与处理中的一种重要变是信号分析与处理中的一种重要变换。直接计算换。直接计算DFT的计算量与变换区间长度的计算量与变换区间长度N的平方成正比，计算量太大，所以在快速的平方成正比，计算量太大，所以在快速傅里叶变换傅里叶变换(简称简称FFT)出现以前，直接用出现以前，直接用DFT算法进行谱分析和信号的实时处理是不切实算法进行谱分析和信号的实时处理是不切实际的。直到际的。直到1965年发现了年发现了DFT的一种快速算的一种快速算法以后，情

2、况才发生了根本的变化。法以后，情况才发生了根本的变化。 第第4 4章章 快速傅里叶变换快速傅里叶变换3/864.1 DFT4.1 DFT的运算量分析的运算量分析4.1.1 4.1.1 直接计算直接计算DFTDFT的运算量的运算量时域频域都是离散的有限长序列时域频域都是离散的有限长序列( )x n( )X k1010( )( ) (01)1( )( ) (01)NnkNnNnkNnX kx n WkNx nX k WnNN DFT计算量太大计算量太大第第4 4章章 快速傅里叶变换快速傅里叶变换4/86 复数乘法：复数乘法：2MCNNN 2(1) (1)ACN NNN 直接计算直接计算DFTDFT

3、需要：需要： 复数加法：复数加法：实数乘法：实数乘法：实数加法：实数加法：24MRN 2222(1)2 424ARN NNNNN直接计算直接计算DFTDFT所需计算量与所需计算量与 成正比！成正比！2N第第4 4章章 快速傅里叶变换快速傅里叶变换5/86 从上面的分析看到，在从上面的分析看到，在DFT计算中，不论计算中，不论是乘法和加法，运算量均与是乘法和加法，运算量均与N2成正比。因此，成正比。因此，N较大时，运算量十分可观。例，计算较大时，运算量十分可观。例，计算N=10点点的的DFT，需要，需要100次复数相乘，而次复数相乘，而N=1024点时点时，需要，需要1048576（一百多万）次

4、复数乘法，如（一百多万）次复数乘法，如果要求实时处理，则要求有很高的计算速度才果要求实时处理，则要求有很高的计算速度才能完成上述计算量。能完成上述计算量。 反变换反变换IDFT与与DFT的运算结构相同，只是的运算结构相同，只是多乘一个常数多乘一个常数1/N，所以二者的计算量相同。，所以二者的计算量相同。第第4 4章章 快速傅里叶变换快速傅里叶变换6/864.1.2 4.1.2 改善改善DFTDFT运算效率的基本途径运算效率的基本途径()()nkk N nnkNNNWWW1、 的对称性：的对称性：nkNWnknk lNNNWW 2、 的周期性：的周期性：nkNW（ 为整数为整数）lnkmnkNm

5、NWW 3、 的可约性：的可约性：nkNW()22NNnknknkNNNNWWWW nknk mNN mWW 第第4 4章章 快速傅里叶变换快速傅里叶变换7/86 例如，对例如，对4点点DFT，直接计算需要，直接计算需要42=16次复数乘法。根据上述次复数乘法。根据上述 的性质，可写成的性质，可写成如下的矩阵形式如下的矩阵形式 knNW0 00 10 (1)1 01 11 (1)2 02 12 (1)(1) 0(1) 1(1) (1)(0)(0)(1)(1)(1)(0)(1)(1)(2)(0)(1)(1)(1)(0)(1)(1)NNNNNNNNNNNNNNNNNNNXxWxWx NWXxWxW

6、x NWXxWxWx NWX NxWxWx NW DFT写成如下的矩阵形式写成如下的矩阵形式 第第4 4章章 快速傅里叶变换快速傅里叶变换8/8600004444012344440246444403694444(0)(0)(1)(1)(2)(2)(3)(3)XxWWWWXxWWWWXxWWWWXxWWWW 令令00004444012344440246444403694444WWWWWWWWWWWWWWWWW 00004444012344440202444403214444WWWWWWWWWWWWWWWWnkNW利利用用的的周周期期性性第第4 4章章 快速傅里叶变换快速傅里叶变换9/86利利用用

7、的的共共轭轭对对称称性性nkNW00004444010144440000444401014444WWWWWWWWWWWWWWWW 114411441111(0)(0)11(1)(1)1111(2)(2)11(3)(3)XxWWXxXxWWXx 则，则，4 4点点DFTDFT矩阵公式变为矩阵公式变为00004444012344440202444403214444WWWWWWWWWWWWWWWWW 第第4 4章章 快速傅里叶变换快速傅里叶变换10/86将该矩阵的第二列和第三列交换，得到将该矩阵的第二列和第三列交换，得到114411441111(0)(0)11(1)(2)1111(2)(1)11(3

8、)(3)XxWWXxXxWWXx 由此得到由此得到 1414(0)(0)(2)(1)(3)(1)(0)(2)(1)(3)(2)(0)(2)(1)(3)(3)(0)(2)(1)(3)XxxxxXxxxxWXxxxxXxxxxW 第第4 4章章 快速傅里叶变换快速傅里叶变换11/86 快速傅里叶算法的基本思想快速傅里叶算法的基本思想（1）利用利用 的性质的性质减少计算量。减少计算量。（2）把长序列的）把长序列的DFT分解成短序列的分解成短序列的DFT，也可以有，也可以有效的减少效的减少DFT运算中复数乘法和复数加法的次数。运算中复数乘法和复数加法的次数。 如果信号长度为如果信号长度为 N，它可表示

9、成：，它可表示成： 当当 时，上式可写成时，上式可写成 ，因子，因子 称称为基（为基（radix）。）。 当当FFT算法中进行序列分解时是采用算法中进行序列分解时是采用则称为基则称为基-2（radix-2）FFT。 knNWmrrrN21rrrrm21mrN rmN2第第4 4章章 快速傅里叶变换快速傅里叶变换12/86 N=2M，4.2 4.2 时间抽取的基时间抽取的基-2FFT-2FFT算法算法.2.1.算法的的基本原理算法的的基本原理1212( )( )( ) 012()( )( ) 2kNkNX kX kW XkNkNX kX kW Xk （）则：则： ( )DFT( )

10、 , 0,1,X kx nn kN1( )(2 ),x rxr 2( )(21),01,2Nx rxrr设：设：若若奇数序号奇数序号偶数序号偶数序号 点序列点序列2N 和和 的的 点点DFT1( )x r2( )xr2N1122( )DFT( ),01,( )DFT( )2X kx rNkXkx r 第第4 4章章 快速傅里叶变换快速傅里叶变换13/86 10)()(NnknNWnxkX/2 1/2 12(21)00(2 )(21)NNkrkrNNrrxr WxrW /2 1/2 1221200( )( )NNkrkrkNNNrrx r Wxr WW/2 1/2 11/22/20012( )(

11、 )( )( )0,1,.,/ 21NNkrkkrNNNrrkNx r WWx r WXkW XkkN第第4 4章章 快速傅里叶变换快速傅里叶变换14/86对于对于 后后N/2的的DFT ：由于由于 （周期性）（周期性）因此因此)(kXkNNkNWW 2/)()()2(21kXWkXNkXkN 12,.,1 , 0 Nk第第4 4章章 快速傅里叶变换快速傅里叶变换15/86例：画出例：画出N=8N=8的时间抽取基的时间抽取基2FFT2FFT算法流图。算法流图。( ) (0), (1), (2), (3), (4), (5), (6), (7)x nxxxxxxxx ( )DFT ( )X kx

12、 n 解：解： 将将 按时间抽取基按时间抽取基2FFT2FFT算法进行分解：算法进行分解：( )x n1( ) (0), (2), (4), (6)x nxxxx 2( ) (1), (3), (5), (7)x nxxxx 1111(0),(1),(2),(3)xxxx 2222(0),(1),(2),(3)xxxx 4点点序列序列第第4 4章章 快速傅里叶变换快速傅里叶变换16/86仍然是偶数仍然是偶数第第4 4章章 快速傅里叶变换快速傅里叶变换17/86 由于由于 , ,因而因而N/2N/2仍是偶数仍是偶数, ,可以进一可以进一步把每个步把每个N/2N/2点的序列再按其奇偶部分分解点的序

13、列再按其奇偶部分分解为两个为两个N/4N/4的子序列的子序列从而从而 可表示为可表示为 MN2 3141221( )()( )()x lxlx lxl0,1,(1)4Nl )(1kX4 14 12(21)1121200( )(2 )(21)NNklklNNllX kxl WxlW 324( )( )kNXkWXk0,1,14Nk 第第4 4章章 快速傅里叶变换快速傅里叶变换18/86因而有因而有对对 也可进行同样的分解：也可进行同样的分解： )(2kX25/2625/26( )( )( ) 014()( )( ) 4kNkNXkXkWXkNkNXkXkWXk （）13/2413/24( )(

14、)( ) 014()( )( ) 4kNkNX kXkWXkNkNX kXkWXk （）21342134( )( )( ) 014()( )( ) 4kNkNX kXkWXkNkNX kXkWXk （）22562256( )( )( ) 014()( )( ) 4kNkNXkXkWXkNkNXkXkWXk （）2/2kkNNWW统统一一系系数数：第第4 4章章 快速傅里叶变换快速傅里叶变换19/86第一次分解：第一次分解：第二次分解：第二次分解：1( ) (0), (2), (4), (6)x nxxxx 2( ) (1), (3), (5), (7)x nxxxx 3( ) (0), (4)

15、x nxx 33(0),(1)xx 4( ) (2), (6)x nxx 44(0),(1)xx 5( ) (1), (5)x nxx 55(0),(1)xx 6( ) (3), (7)x nxx 66(0),(1)xx 1111(0),(1),(2),(3)xxxx2222(0),(1),(2),(3)xxxx2点点序列序列第第4 4章章 快速傅里叶变换快速傅里叶变换20/86 如下图所示：如下图所示： 那么依次类推，经过那么依次类推，经过M-1次分解后，将次分解后，将N点点DFT分解分解成成N/2个两点个两点DFT第第4 4章章 快速傅里叶变换快速傅里叶变换21/86第三次分解：第三次分解

16、：3( ) (0), (4)x nxx 33(0),(1)xx 4( ) (2), (6)x nxx 44(0),(1)xx 5( ) (1), (5)x nxx 55(0),(1)xx 6( ) (3), (7)x nxx 66(0),(1)xx 7( ) (0)x nx 8( ) (4)x nx 9( ) (2)x nx 10( ) (6)xnx 11( ) (1)xnx 12( ) (5)xnx 13( ) (3)xnx 14( ) (7)xnx 第第4 4章章 快速傅里叶变换快速傅里叶变换22/86 例如例如N=8时，分解为时，分解为4个两点个两点DFT，其输出分别为，其输出分别为例如

17、：例如：即即上式中用到了上式中用到了3456( ),( ),( ),( )XkXkXkXk0,1k 4 116646400( )( )( ),0,1NklklNNllXkx l Wx l Wk 00066262(0)(0)(1)(3)(7)(3)(7)NXxW xxW xxW x 11066262(1)(0)(1)(3)(7)(3)(7)NXxW xxW xxW x 2j1j0221NWeeW 第第4 4章章 快速傅里叶变换快速傅里叶变换23/86 下图为下图为N=8时的一个完整基时的一个完整基-2DIT-FFT运算流图运算流图第第4 4章章 快速傅里叶变换快速傅里叶变换24/86对对 点序列

18、：点序列： 4.2.2 4.2.2 运算量运算量N N点点N/2N/2点点N/4N/4点点1 1点点2MN 1 1点序列的点序列的DFTDFT等于本身序列值等于本身序列值最后按时间抽取基最后按时间抽取基2FFT2FFT算法的计算法的计算量仅由蝶行运算产生。算量仅由蝶行运算产生。第第4 4章章 快速傅里叶变换快速傅里叶变换25/86一次蝶式运算需：一次蝶式运算需：1 1次复数乘法和次复数乘法和2 2次复数加法次复数加法按时间抽取基按时间抽取基-2FFT-2FFT共有：共有： 个运算蝶个运算蝶2log22NNMN 基基-2FFT-2FFT 复数乘法：复数乘法： 复数加法：复数加法：2log2NN2

19、2log22logNNNN 直接直接DFTDFT 2N2N1212 ( )( )( ) 012()( )( )2kNkNX kX kW XkNkNX kX kW Xk （）第第4 4章章 快速傅里叶变换快速傅里叶变换26/860100200300400500600700800900100001234567891011x 105直接计直接计算算DFTDFTFFTFFT第第4 4章章 快速傅里叶变换快速傅里叶变换27/86DFT和和FFT运算量比较运算量比较function Xk=FourierTran(xn,N)if nargin 2 N = length(xn);end n = 0 : N -

20、 1;k = 0 : N - 1;WN = exp( - j * 2 * pi / N);nk = n * k;WNnk = WN . nk;Xk = xn * WNnk;第第4 4章章 快速傅里叶变换快速傅里叶变换28/86Matlab程序演示程序演示1 1、直接利用定义计算、直接利用定义计算DFTDFTx=ones(1,1024);f=() FourierTran(x,1024);timeit(f)运行结果：运行结果：ans = 0.7239ans = 0.72392 2、利用、利用FFTFFT算法计算算法计算DFTDFTf1=() fft(x,1024);timeit(f1)运行结果：运

21、行结果：ans = 1.5766e-005ans = 1.5766e-0050.7239 / (1.5766e-005) = 4.5915e+004第第4 4章章 快速傅里叶变换快速傅里叶变换29/864.2.3.FFT4.2.3.FFT算法的特点算法的特点1 1）原位计算（同址运算）原位计算（同址运算）第第4 4章章 快速傅里叶变换快速傅里叶变换30/86m表示第表示第m级迭代，级迭代，i，j表示数据所在的行数表示数据所在的行数1111( )( )( )( )( )( )rmmmNrmmmNAiAiAj WAjAiAj W第第4 4章章 快速傅里叶变换快速傅里叶变换31/862 2）输入序列

22、的序号及整序规律）输入序列的序号及整序规律DITFFT算法的输入序列的排序看起来似乎算法的输入序列的排序看起来似乎很乱，但仔细分析就会发现这种排序是很有很乱，但仔细分析就会发现这种排序是很有规律的。由于规律的。由于N=2M，所以顺序数可用，所以顺序数可用M位二位二进制数进制数(nM-1nM-2n1n0)表示。表示。 第第4 4章章 快速傅里叶变换快速傅里叶变换32/86(0) (1) (2) (3) (4) (5) (6) (7)xxxxxxxx(0) (2) (4) (6)xxxx(1) (3) (5) (7)xxxx(0) (4)xx(2) (6)xx(1) (5)xx(3) (7)xx(

23、0) (4) (2) (6) (1) (5) (3) (7)xxxxxxxx000 010 100 110001 011 101 111000 100010 110001 101011 111000100010110001101011111000 001 010 011 100 101 110 111第第4 4章章 快速傅里叶变换快速傅里叶变换33/86输入的混序输入的混序是通过输入正序序列按是通过输入正序序列按码位倒置码位倒置实现的实现的。自然顺序自然顺序( )n二进制数二进制数 倒位序二进制数倒位序二进制数 倒位顺序倒位顺序 ( )n0123456700000101001110010111

24、011100010001011000110101111104261537第第4 4章章 快速傅里叶变换快速傅里叶变换34/86x(0)x(1)x(2)x(3)x(4)x(5)x(6)x(7)A(0)A(1)A(2)A(3)A(4)A(5)A(6)A(7)A(0)A(1)A(2)A(3)A(4)A(5)A(6)A(7)x(0)x(4)x(2)x(6)x(1)x(5)x(3)x(7)输入数据的变址处理第第4 4章章 快速傅里叶变换快速傅里叶变换35/86 3 3）各类蝶形运算两节点的）各类蝶形运算两节点的“距离距离”及及 的变化规律的变化规律rNW对对N = 2M点点FFT，输入倒位序，输出自然序

25、，第，输入倒位序，输出自然序，第m级级运算每个蝶形的两节点距离为运算每个蝶形的两节点距离为 2m1，第，第m级运算：级运算：1111111( )( )(2)(2)( )(2)mrmmmNmmrmmmNAiAiXiWAiAiXiW 1111( )( )( )( )( )( )rmmmNrmmmNAiAiAj WAjAiAj W第第4 4章章 快速傅里叶变换快速傅里叶变换36/86蝶形运算两节点的第一个节点为蝶形运算两节点的第一个节点为i值，表示值，表示成成M位二进制数，左移位二进制数，左移M m位，把右边空位，把右边空出的位置补零，结果为出的位置补零，结果为r的二进制数。的二进制数。rNW 的的

26、确确定定2( )2Mmri （1 1）直接计算法）直接计算法第第4 4章章 快速傅里叶变换快速傅里叶变换37/86第第4 4章章 快速傅里叶变换快速傅里叶变换38/86级数 的取值范围 重复组数 第一级 第二级 第三级 第m级 第M级 1rNW0NW2/N0NW4/NNW4/N0NW0NW0NW8/NNW8/2NNW8/3NNWmNNW2/1NWmNNW2/22NW mmNNW2/)12(112/NNW8/NmN 2/ （2 2）逆推法）逆推法第第4 4章章 快速傅里叶变换快速傅里叶变换39/86 4.2.4 4.2.4、DITDIT算法的其他形式流图算法的其他形式流图n输入自然序输出倒位序输

27、入自然序输出倒位序n输入输出均自然序输入输出均自然序n相同几何形状相同几何形状n输入自然序输出倒位序输入自然序输出倒位序第第4 4章章 快速傅里叶变换快速傅里叶变换40/86 第第4 4章章 快速傅里叶变换快速傅里叶变换41/86第第4 4章章 快速傅里叶变换快速傅里叶变换42/86第第4 4章章 快速傅里叶变换快速傅里叶变换43/864.3.1 算法的基本原理算法的基本原理 设序列设序列x(n)长度为长度为N=2M，首先将，首先将x(n)前后对半分开，前后对半分开，得到两个子序列，其得到两个子序列，其DFT可表示为如下形式：可表示为如下形式： 12/02/12/0)2/(12/012/12/

28、010)2/()()2/()()()()()(DFT)(NnnkNNkNNnkNnNNnnkNNNnnkNNnnkNNnnkNWNnxWnxWNnxWnxWnxWnxWnxnxkX4.3 频域抽取基频域抽取基-2 FFT算法算法第第4 4章章 快速傅里叶变换快速傅里叶变换44/86由于由于N点点DFT按按k的奇偶分组可分为两个的奇偶分组可分为两个N/2的的DFTk取偶数时（取偶数时（k=2r，r=0,1,.,N/2-1） /21( 1)1kNkNkWk 为为偶偶数数为为奇奇数数/2 120/2 1/20(2 )( )2( )2NrnNnNrnNnNXrx nx nWNx nx nW 第第4 4

29、章章 快速傅里叶变换快速傅里叶变换45/86设设2221(21)010(21)( )()2( )()2NNNnrNnnnrNnNXrx nx nWNx nx nWW 12( )( )+ ()20,1,12( )( )()2nNNx nx nx nNnNx nx nx nW 第第4 4章章 快速傅里叶变换快速傅里叶变换46/86 将将x1(n)和和x2(n)分别代入分别代入 和和 式，可得式，可得/2 11/20/2 12/20(2 )( )(21)( )NrnNnNrnNnXrx nWXrx nW (2 )Xr(21)Xr则则X(2r)和和X(2r+1)分别是分别是x1(n)和和x2(n)的的

30、 N / 2点点DFT，记为记为X1(k)和和X2(k)第第4 4章章 快速傅里叶变换快速傅里叶变换47/86DIF-FFT的一次分解运算流图（N=8）：第第4 4章章 快速傅里叶变换快速傅里叶变换48/86 由于由于N/2仍然为仍然为2的整数幂，继续将的整数幂，继续将N/2点点DFT分成偶数组和奇数组，这样每个分成偶数组和奇数组，这样每个N/2点点DFT又可分解成两个又可分解成两个N/4点点DFT，其输，其输入序列分别是按上下对半分图开后通过蝶入序列分别是按上下对半分图开后通过蝶式运算构成的式运算构成的4个子序列个子序列 ,如下图如下图第第4 4章章 快速傅里叶变换快速傅里叶变换49/86第

31、第4 4章章 快速傅里叶变换快速傅里叶变换50/86 按照以上方法继续分解下去，经过按照以上方法继续分解下去，经过M-1次次分解，最后分解为分解，最后分解为N/2个两点个两点DFT，这，这N/2个个2点点DFT的输出就是的输出就是N点点DFT的结果的结果X(k) ，如下图，如下图第第4 4章章 快速傅里叶变换快速傅里叶变换51/86第第4 4章章 快速傅里叶变换快速傅里叶变换52/86第第4 4章章 快速傅里叶变换快速傅里叶变换53/86蝶形运算两节点的第一个节点为值蝶形运算两节点的第一个节点为值k，表，表示成示成M位二进制数，左移位二进制数，左移m-1位，把右边位，把右边空出的位置补零，结果

32、为空出的位置补零，结果为r的二进制数。的二进制数。rNW 的的确确定定12( )2mrk第第4 4章章 快速傅里叶变换快速傅里叶变换54/86 以上所讨论的以上所讨论的FFT的运算方法同样可用的运算方法同样可用于于IDFT的运算，简称为的运算，简称为IFFT。即快速。即快速傅傅里叶里叶反变换。从反变换。从IDFT的定义出发，可以的定义出发，可以导出下列二种利用导出下列二种利用FFT来计算来计算IFFT的方的方法法。4.4 快速傅里叶反变换快速傅里叶反变换第第4 4章章 快速傅里叶变换快速傅里叶变换55/8.1.1.稍微变动稍微变动FFTFFT程序和参数可实现程序和参数可实现IF

33、FTIFFT1010DFT1IDFT( )( )( )( )( )( )NnkNnNnkNkX kx nx n Wx nX kX k WN只要只要DFT的每个系数的每个系数 换成换成 ,最后再乘以最后再乘以常数常数1/N就可以得到就可以得到IDFT的快速算法的快速算法-IFFT。另外另外,可以将常数可以将常数1/N分配到每级运算中分配到每级运算中, ,也就是每级也就是每级蝶形运算均乘以蝶形运算均乘以1/2。 MN)21(1 nkNWnkNW第第4 4章章 快速傅里叶变换快速傅里叶变换56/86第第4 4章章 快速傅里叶变换快速傅里叶变换57/8.2不改变不改变FFTFFT的程

34、序直接实现的程序直接实现IFFTIFFT因为因为所以所以 101,.,0,)(1)(NknkNNnWkXNnx 10)(1)(NknkNWkXNnx两边取共轭有：两边取共轭有：101( )( ),0,.,1NnkNkx nXk WnNN 1DFT( )XkN 第第4 4章章 快速傅里叶变换快速傅里叶变换58/864.6 实序列的实序列的FFT算法算法n根据序列根据序列DFT的共轭对称性知道，任意的共轭对称性知道，任意复数序列的实部的复数序列的实部的DFT对应于其对应于其DFT共共轭对称分量，而其虚部的轭对称分量，而其虚部的DFT对应于其对应于其DFT共轭反对称分量，故可用一次共轭反对称分量，故

35、可用一次N点点DFT计算两个计算两个N点实序列的点实序列的DFT。第第4 4章章 快速傅里叶变换快速傅里叶变换59/86设设 和和 为两个长度为为两个长度为N的实序列，的实序列，按如下方式构造新序列按如下方式构造新序列 ：1( )x n2( )x n( )y n12( )( )( )y nx njx nN点点DFT为为( )DFT ( )( )( )epopY ky nYkYk 4.6.1利用频谱对称性求实信号的利用频谱对称性求实信号的FFT第第4 4章章 快速傅里叶变换快速傅里叶变换60/86根据共轭对称性，则根据共轭对称性，则*1*21( )DFT( ) ( )()( )21( )DFT(

36、 ) ( )()( )2epNNopNNYkx nY kYNkRnYkjx nY kYNkRn 所以所以*11*221( )DFT( )( ) ( )()( )2( )DFT( )( ) ( )()( )2epNNopNNXkx nYkY kYNkRnjXkxnjYkY kYNkRn 第第4 4章章 快速傅里叶变换快速傅里叶变换61/86设一个设一个2N点的序列点的序列 ，现按偶、奇进行，现按偶、奇进行分解得：分解得：( )x n12( )(2 )01( )(21)x nxnnNx nxn 再构造复数序列再构造复数序列y(n)12( )( )( )y nx njx n第第4 4章章 快速傅里叶

37、变换快速傅里叶变换62/86122122( )( )( )01()( )( )kNkNY kXkWXkkNY kNXkWXk 这相当于一个这相当于一个N点点DFT运算加上运算加上DIT-FFT蝶蝶形运算，当形运算，当N较大时，可节省近一半的计算较大时，可节省近一半的计算量。量。然后先求出然后先求出x1(n)和和x2(n)的的N点点DFTX1(k)和和X2(k)，因为，因为x1(n)和和x2(n)分别是原序列分别是原序列x(n)的的偶、奇序列，这与按时间抽取的偶、奇序列，这与按时间抽取的FFT算法的算法的分解思路完全相同，故根据分解思路完全相同，故根据DIT-FFT蝶形运蝶形运算式，可得算式，可

38、得第第4 4章章 快速傅里叶变换快速傅里叶变换63/864.6.2 离散哈特曼变换1 1、离散哈特曼变换的定义、离散哈特曼变换的定义设设x x( (n n) )为一为一N N点的实序列，则其离散哈特曼变换点的实序列，则其离散哈特曼变换（DHTDHT）为）为1H02( )DHT ( )( )cas()01NnnkXkx nx nkNN逆变换为逆变换为1HH012( )IDHT( )( )cas()01Nknkx nXkXknNNN式中式中222cas()cos()sin()nknknkNNN第第4 4章章 快速傅里叶变换快速傅里叶变换64/862. DHT与与DFT的关系的关系DFT和和IDFT

39、用欧拉公式展开可表示成用欧拉公式展开可表示成12/010( )DFT ( )( )22( )cos()sin()01Njnk NnNnX kx nx n enknkx njkNNN12/010( )IDFT( )( )22( )cos()sin()01Njnk NkNkx nX kX k enknkX kjnNNN第第4 4章章 快速傅里叶变换快速傅里叶变换65/86可以看出，可以看出，DHT的核函数的核函数是是DFT核函数核函数的实部和虚部之和。的实部和虚部之和。将将 分解为奇对称分量和偶对称分量之和分解为奇对称分量和偶对称分量之和222cas()cos()sin()nknknkNNN2/2

40、2cos()sin()jnk NnknkejNNH( )XkHHH( )( )( )eoXkXkXk第第4 4章章 快速傅里叶变换快速傅里叶变换66/86其中其中由由DHT的定义可得的定义可得HHHHHH1( )( )()21( )( )()2eoXkXkXNkXkXkXNk1H01H02( )( )cos()2( )( )sin()NenNonnkXkx nNnkXkx nN第第4 4章章 快速傅里叶变换快速傅里叶变换67/86所以所以因此因此 如果不考虑因子如果不考虑因子1/2，只要增加，只要增加2N次实数加法次实数加法运算就能由运算就能由DHT求出求出DFT。HHH( )( )( )(

41、)Re( )Im( )eoX kXkjXkXkX kX kHHHH1( )( )()( )()22jX kXkXNkXkXNk第第4 4章章 快速傅里叶变换快速傅里叶变换68/863. DHT的快速算法（的快速算法（FHT）n仿照快速仿照快速DFT的分解方法，可通过时域的抽取的分解方法，可通过时域的抽取或频域抽取方式实现快速或频域抽取方式实现快速DHT。将。将N=2M点的点的实序列实序列 进行奇偶抽取进行奇偶抽取n则则12( )(2 )0,1,/2 1( )(21)x rxrrNx rxr/2 1/2 1H0022( )DHT ( )(2 )cas(2)(21)cas(21) )NNrrXkx

42、 nxrrkxrrkNN第第4 4章章 快速傅里叶变换快速傅里叶变换69/86根据三角公式根据三角公式令令 ， ，根据，根据DHT的周期性，在的周期性，在 时时cas()cascoscas()sin/2 1H10/2 1/2 122002( )( )cas()/22222cos()( )cas()sin()( )cas()/2/2NrNNrrXkx rrkNkx rrkkx rrkNNNN1H1( )DHT( )Xkx n2H2( )DHT( )Xkx n0k H1H2222( )( )cos()( )sin()()2HHNXkXkk Xkk XkNN第第4 4章章 快速傅里叶变换快速傅里叶变

43、换70/86类似于类似于DIT基基-2FFT分解的同址运算思想，可用分解的同址运算思想，可用 、 、 和和 四个四个点同址运算得出点同址运算得出 、 、和和 。这种运算构成了一个运算蝶形，。这种运算构成了一个运算蝶形，称为称为“哈特曼蝶形哈特曼蝶形”。设设将上式中将上式中k分别取分别取k、N/2+k、N/2-k和和N-k四个值四个值，并考虑，并考虑 和和 以以N/2为周期，当为周期，当时，得到时，得到 1H( )Xk2H( )Xk1H(/2)XNk2H(/2)XNkH( )XkH(/2)XNkH(/2)XNkH()XNk2cos()kckN2sin()kskN1H( )Xk2H( )Xk8N

44、第第4 4章章 快速傅里叶变换快速傅里叶变换71/86H1H22H1H22H1H22H1H22( )( )( )(/2)(/2)( )( )(/2)04(/2)(/2)( )(/2)()(/2)( )(/2)kHkHkHkHkHkHkHkHXkXkc Xks XNkXNkXkc Xks XNkNkXNkXNks Xkc XNkXNkXNks Xkc XNkH1H2H1H2(0)(0)(0)0(/2)(0)(0)HHXXXkXNXXH1H2H1H2(/4)(/4)(/4)(3/4)(/4)(/4)4HHXNXNXNNkXNXNXN第第4 4章章 快速傅里叶变换快速傅里叶变换72/86上述的运算可

45、用哈德曼蝶形表示如图上述的运算可用哈德曼蝶形表示如图4.5.1所示所示第第4 4章章 快速傅里叶变换快速傅里叶变换73/86四点的四点的FHT的蝶形图如图的蝶形图如图4.5.2所示。所示。第第4 4章章 快速傅里叶变换快速傅里叶变换74/86图图4.5.3显示了显示了8点的点的DIT基基-2FHT流图流图。第第4 4章章 快速傅里叶变换快速傅里叶变换75/86运算量运算量乘法次数乘法次数加法次数加法次数211(2)34MLHLmNNMN33222HaNMN2MN 第第4 4章章 快速傅里叶变换快速傅里叶变换76/864.7 线性卷积和线性相关的线性卷积和线性相关的FFT算法算法10( )( )

46、* ( )() ()Mmy nx nh nh m x nm 4.7.1 有限长序列线性卷积的有限长序列线性卷积的FFT算法算法序列：序列：x(n)（n=0L-1） h(n)（n=0M-1）线性卷积：线性卷积：y(n)的长度：的长度：1MLdmLM 计算量（乘法次数）：计算量（乘法次数）：第第4 4章章 快速傅里叶变换快速傅里叶变换77/86用用FFT算法也就是用算法也就是用循环循环卷积来代替线性卷积卷积来代替线性卷积时，为了不产生混叠，其必要条件是使时，为了不产生混叠，其必要条件是使x(n)，h(n)都补零值，补到至少都补零值，补到至少N=M+L-1，即，即( ) 01( )0 1h nnMh

47、 nMnN ( ) 01( )0 1x nnLx nLnN 这时，这时，y(n)就能代表线性卷积的结果。就能代表线性卷积的结果。( )( )y nx n ( )h nN第第4 4章章 快速傅里叶变换快速傅里叶变换78/86利用利用FFT进行线性卷积的步骤：进行线性卷积的步骤：1.将序列将序列x(n)和和h(n)补零，使其长度补零，使其长度NL+M-1 2.做做x(n)和和h(n)的的N点点FFT得到得到X(k)和和H(k)，并计，并计算算Y(k)=X(k)H(k)3.求求Y(k)的的IFFT获得线性卷积的结果为获得线性卷积的结果为( )IFFT ( )01y nY knN第第4 4章章 快速傅

48、里叶变换快速傅里叶变换79/86整个过程中，共需要进行三次整个过程中，共需要进行三次FFT运算，共需运算，共需 次相乘，还有步骤（次相乘，还有步骤（2）的）的N次相乘，次相乘，因此共需相乘次数为因此共需相乘次数为23log2NN2233log(1log)22FmNNNNN22332(1log)2(1)1log (1)22dmFmMLMLKmNNMLML 第第4 4章章 快速傅里叶变换快速傅里叶变换80/86x(n)与与h(n)点数差不多，点数差不多，设设M=L，当当M=L且且M超过超过64以后，以后，M越长越长循环循环卷积的卷积的好处越明显。因而将好处越明显。因而将循环循环卷积称为卷积称为快速

49、卷积快速卷积。2-12NMM2253106log4(log)22mMMKMM 第第4 4章章 快速傅里叶变换快速傅里叶变换81/86 ，则，则 ，有，有LM 1NLML223logmMKL 4.7.2 有限长序列无限长序列线性卷积的有限长序列无限长序列线性卷积的FFT算法算法基本思路：基本思路： 把有限序列和无限序列之间的线性卷积把有限序列和无限序列之间的线性卷积变成若干个有限序列之间的线性卷积。变成若干个有限序列之间的线性卷积。具体方法：具体方法： 重叠相加法重叠相加法(Overlap-add method) (Overlap-add method) 重叠保留法重叠保留法(Overlap-s

50、ave method)(Overlap-save method)第第4 4章章 快速傅里叶变换快速傅里叶变换82/86一、重叠相加法（一、重叠相加法（overlap-add methodoverlap-add method）1.1.主要思想主要思想( ) ()ih i x ni 为无限序列，为无限序列， 为为M M点有限序列，计算点有限序列，计算( )x n( )h n10( ) () ()Mih i x nin 思路：将长序列思路：将长序列 分段分段( )x n( )( )( )y nh nx n 第第4 4章章 快速傅里叶变换快速傅里叶变换83/86a.a.将将 均匀分段，每段长度为均匀分

51、段，每段长度为N N10( )( )Piix nx n 10( )( )* ( )( )* ( )Piiy nx nh nx nh n ( ) (1)1( )0 ix niLniLx nelse 10( )Piiy n (分段过程无重叠分段过程无重叠)( )x n( )iy n 第第4 4章章 快速傅里叶变换快速傅里叶变换84/86用线性卷积或循环卷积用线性卷积或循环卷积 b. b. 计算子段卷积计算子段卷积( )( )* ( )iiy nx nh n , (1)2niLiLM01iP , (1)1 0, 1iLiLM 计算方法：计算方法：用循环卷积计算需要注意什么？用循环卷积计算需要注意什么

52、？ 思考：思考：第第4 4章章 快速傅里叶变换快速傅里叶变换85/86c. c. 相加相加10( )( )Piiy ny n ( )iy n长度为长度为L+M-1相加时，相邻相加时，相邻 重叠重叠M-1个点。个点。( )iy n( )iy n起点为起点为iL第第4 4章章 快速傅里叶变换快速傅里叶变换86/86 由于由于xi (n)的长度为的长度为L，h(n)的长度为的长度为M，所以所以yi (n)的长度为的长度为即即yi (n)的范围为的范围为将上式将上式xi (n)的范围比较，显然的范围比较，显然yi (n)的范围比的范围比xi (n)的的范围大，超出的点数为范围大，超出的点数为而而yi

53、+1(n)的范围为的范围为1NML2(1)2iLniLLMiLM (1)2(1)11iLMiLM (1)(1)2(2)2iLniLLMiLM第第4 4章章 快速傅里叶变换快速傅里叶变换87/86 可以看出，由可以看出，由(i+1) L到到(i+1) L+M-2这这M-1点上，点上， yi (n)的后部分与的后部分与yi +1(n)的前部分发生了重叠。这样的前部分发生了重叠。这样对于在此范围内的每一个对于在此范围内的每一个n值，原序列值，原序列x(n)与与h(n)的的卷积结果应该是卷积结果应该是即并不是将各段线性卷积的结果简单地拼接在一起即并不是将各段线性卷积的结果简单地拼接在一起，在某些点上需要前后两端的结果重叠相加。，在某些点上需要前后两端的结果重叠相加。1( )( )( )iiy ny nyn 第第4 4章章 快速傅里叶变换快速傅里叶变换88/86第第4 4章章 快速傅里叶变换快速傅里叶变换89/86为了得到两个序列为了得到两个序列x(n)和和h(n)最终的线性卷积结果，最终的线性卷积结果，需将相邻两段的需将相邻两段的M-1个重叠点的值相加，故称为个重叠点的值相加，故称为重叠重叠相加法相加法。在求出各。在求出各yi (n)后，后，x(n)和和h(n)的线性卷积的线性卷积y(n)可分段表示为