离散傅里叶变换的矩阵表示及其运算量

1、 第第4 4章章 FFTFFT 4.1 4.1 引言引言 4.1.1 离散傅里叶变换的矩阵表示及其运算量离散傅里叶变换的矩阵表示及其运算量 n DFT在数字信号处理中起着非常重要的作用，在数字信号处理中起着非常重要的作用， 这是与这是与DFT存在着高效算法，存在着高效算法， 即快速傅里叶变换即快速傅里叶变换(FFT) 分不开的。快速分不开的。快速运算的关键是减少运算量。运算的关键是减少运算量。n离散傅里叶变换对为：离散傅里叶变换对为： (4.1) (4.2)n式中式中 。 下面要用矩阵来表示下面要用矩阵来表示DFT关系。关系。 101, 1 , 0,)()(:NnnkNNkWnxkXDFT1,

2、 1 , 0,)(1)(:10NnWkXNnxIDFTNknkNNjNeW2n 一般情况下，信号序列一般情况下，信号序列x(n) 及其频谱序列及其频谱序列X(k) 都是用复数都是用复数来表示的，来表示的，WN当然也是复数。因此，计算当然也是复数。因此，计算DFT的一个值的一个值X(k) 需要进行需要进行N次复数乘法（与次复数乘法（与1相乘也包括在内）和相乘也包括在内）和N-1次复数次复数加法；所以，直接计算加法；所以，直接计算N点的点的DFT需要进行需要进行N2 次复数乘法和次复数乘法和N(N-1) 复数加法。复数加法。 n显然，直接计算显然，直接计算N点的点的IDFT所需的复乘和复加的次数也

3、是这所需的复乘和复加的次数也是这么多。当么多。当N足够大时，足够大时，N2 N(N-1), 因此，因此，DFT与与IDFT的运的运算次数与算次数与N2 成正比，随着成正比，随着N的增加，运算量将急剧增加，而的增加，运算量将急剧增加，而在实际问题中，在实际问题中，N往往是较大的，因此有必要对往往是较大的，因此有必要对DFT与与IDFT的计算方法予以改进。的计算方法予以改进。 4.1.2 因子的特性因子的特性n DFT和和IDFT的快速算法的导出主要是根据的快速算法的导出主要是根据 因子的特因子的特性。性。 1周期性：周期性： n对离散变量对离散变量n有同样的周期性。有同样的周期性。nkNnNNn

4、kNNknNWWWW)(nkNWnkNW 2对称性：对称性： n或或 3. 其它其它: knNNknNnkNnkNWWWW)()(*nkNNnkNknNnkNWWWW)()(*kNNNjkNNNkNNkNWeWWWW222)2(kNkNjkNjkNWeeW 2222224.2 4.2 基基2 2时间抽选的时间抽选的FFTFFT算法算法 4.2.1 算法推导算法推导n 已经知道：已经知道： n令令DFT的长度的长度N=2M，M为正整数。为正整数。1-N,0,1,k )()()(10nkNNnWnxnxDFTkXn令：令： n于是有：于是有：12N,0,1,r ) 12()()2()(rxrprx

5、rg1-N,0,1,k )()()()() 12()2()(120rk 2/120 2/120)12(1202kPWkGWrpWWrgWrxWrxkXkNNrNkNNrrkNNrkrNNrrkNn其中，其中， n是由是由x(n)的偶数抽样点形成的的偶数抽样点形成的DFT；而；而 kr 2/120120kr 2/)2()()(NNrNrNWrxWrgkG120kr 2/kr 2/120) 12()()(NrNNNrWrxWnpkPn是由是由x(n)的奇数抽样点形成的的奇数抽样点形成的DFT。但是这两个式子并不完。但是这两个式子并不完全是全是N/2点的点的DFT，因为，因为k的范围仍然是由的范围仍

6、然是由0到到N-1，因此，还，因此，还应该进一步考虑应该进一步考虑k由由N/2到到N-1范围的情况。范围的情况。n现在令现在令 ，故对于后半段有：，故对于后半段有： n同理：同理： n又知：又知： 12, 1 , 0Nk)()()()()2(120kr 21202120)2(2/22kGWrgWWrgWrgNkGNrNNrrkNrNrNkrNNN)()2(kPNkPk 2NNkNWW 图图 4.1 N点点DFT分解为两个分解为两个N/2点的点的DFT (N=8) 图图 4.2 N/2 点的点的DFT 分解为两个分解为两个N/4点的点的DFT (N=8)n综上所述，可以得到：综上所述，可以得到：

7、 n其中其中G(k)、P(k) 分别是分别是x(n)的偶数点和奇数点的的偶数点和奇数点的N/2点点DFT。 )()()2()()()(kPWkGNkXkPWkGkX kNkN12, 1 , 0Nkn 这样，我们就将一个这样，我们就将一个N点的点的DFT分解成了两个分解成了两个N/2点的点的DFT，由于由于DFT的运算量与其点数的平方成正比，因此使运算量减的运算量与其点数的平方成正比，因此使运算量减少了。但是，还应该将每一个少了。但是，还应该将每一个N/2点的点的DFT再分解为两个再分解为两个N/4点的点的DFT，如此下去，直到分解为，如此下去，直到分解为2点的点的DFT为止，总共需为止，总共需

8、要进行要进行loglog2 2N-1=logN-1=log2 2(N/2)(N/2)次分解。次分解。图图 4.3 2 点点 DFT 信号流图（蝶形图）信号流图（蝶形图）n对于对于2点点DFT，有：，有： n所以所以2点点DFT的运算只需一次加法和一次减法，这样的运算的运算只需一次加法和一次减法，这样的运算叫做蝶形运算，这样的信号流图叫做蝶形图。叫做蝶形运算，这样的信号流图叫做蝶形图。1111111122WTn该算法每次分解都是将时域序列按奇偶分为两组，因此要求该算法每次分解都是将时域序列按奇偶分为两组，因此要求N等于等于2的正整数幂，故将这种的正整数幂，故将这种FFT算法叫做基算法叫做基2时间

9、抽选法。时间抽选法。 4.2.2 算法特点算法特点 1. 倒序重排倒序重排n这种这种FFT算法的每次分解都是将输入序列按照奇偶分为两组，算法的每次分解都是将输入序列按照奇偶分为两组，故要不断地将每组输入数据按奇偶重排，直到最后分解为故要不断地将每组输入数据按奇偶重排，直到最后分解为2点的点的DFT，输入数据才不再改变顺序。这样做的结果，使得，输入数据才不再改变顺序。这样做的结果，使得作作FFT运算时，输入序列的次序要按其序号的倒序进行重新运算时，输入序列的次序要按其序号的倒序进行重新排列。排列。 n现在将图现在将图4.4中输入序号以及重排后的序号按二进制写出如下中输入序号以及重排后的序号按二进

10、制写出如下（注：下标（注：下标“2”表示二进制数）。可以看出，将输入序号表示二进制数）。可以看出，将输入序号的二进制表示（的二进制表示（n2n1n0）位置颠倒，得到（）位置颠倒，得到（n0n1n2），就是），就是相应的倒序的二进制序号。因此，输入序列按倒序重排，实相应的倒序的二进制序号。因此，输入序列按倒序重排，实际上就是将序号为（际上就是将序号为（n2n1n0）的元素与序号为（）的元素与序号为（n0n1n2）的）的元素的位置相互交换。元素的位置相互交换。 2. 同址计算同址计算 n 从图从图4.4可以看出，整个算法流图可以分为四段，（可以看出，整个算法流图可以分为四段，（0）段）段为倒序重排

11、，后面为倒序重排，后面3段为段为3(log28=3)次迭代运算：首先计算次迭代运算：首先计算2点点DFT，然后将，然后将2点点DFT的结果组合成的结果组合成4点点DFT，最后将，最后将4点点DFT组合为组合为8点点DFT。因此，对于。因此，对于N点点FFT，只需要一列存储，只需要一列存储N个复数的存储器。个复数的存储器。 3. 运算量运算量n观察图观察图4.4可知，图可知，图4.3所示的蝶形图实际上代表了所示的蝶形图实际上代表了FFT的基的基本运算，它实际上只包含了两次复数加法运算。一个本运算，它实际上只包含了两次复数加法运算。一个N(N=2M) 点的点的FFT，需要进行，需要进行M=log2

12、N次迭代运算。每次迭次迭代运算。每次迭代运算包含了代运算包含了N/2个蝶型，因此共有个蝶型，因此共有N次复数加法；此外，除次复数加法；此外，除了第一次的了第一次的2点点DFT之外，每次迭代还包括了之外，每次迭代还包括了N/2次复数乘法次复数乘法（即乘（即乘WN的幂）。因此，一个的幂）。因此，一个N点的点的FFT共有复数乘法的次共有复数乘法的次数为：数为： n复数加法的次数为：复数加法的次数为： n因此，因此，FFT算法比直接计算算法比直接计算DFT大大减少了运算量，尤其是大大减少了运算量，尤其是当当N较大时，计算量的减少更为显著。比如，当较大时，计算量的减少更为显著。比如，当N=1024时，时

13、，采用采用FFT算法时复数乘法的次数，不超过直接计算算法时复数乘法的次数，不超过直接计算DFT时复时复乘次数的千分之五。乘次数的千分之五。 2log2) 1(log2222NNNNMcNNNNAc222loglog224.3 4.3 基基2 2频率抽选的频率抽选的FFTFFT算法算法n 时间抽选法是在时域内将输入序列时间抽选法是在时域内将输入序列x(n)逐次分解为偶数点逐次分解为偶数点子序列和奇数点子序列，通过求子序列的子序列和奇数点子序列，通过求子序列的DFT而实现整个序而实现整个序列的列的DFT。而频率抽选法则是在频域内将。而频率抽选法则是在频域内将X(k)逐次分解成偶逐次分解成偶数点子序

14、列和奇数点子序列，然后对这些分解得越来越短的数点子序列和奇数点子序列，然后对这些分解得越来越短的子序列进行子序列进行DFT运算，从而求得整个的运算，从而求得整个的DFT。当然，同样要。当然，同样要求求N为为2的正整数幂。的正整数幂。 n 设设 ，则可以分别表示出，则可以分别表示出k为偶数和奇数时为偶数和奇数时的的X(k)。 n其中，其中， 12, 1 , 0Nr120rn 2/120rn 2/120rn 2/120120)2(22102)()2()()2()()()2(NnNNnrNNNNnNNnNnrNnNnrNNnrnNWngWWNnxWnxWNnxWnxWnxrX11)2()()(2N,

15、0,n Nnxnxng120rn 2/120rn 2/120rn 2/12022120rn 2/120120)2)(12(2)() 1()2()()2()()2()() 12 (NnNnNNnnNNNnnNNNnNNrNNnNnrNNnnNNNnNnNnrNnNnrNWWnpWWNnxWWnxWWWWNnxWWnxWNnxWWnxrXn其中，其中， n 显然，显然，X(2r) 为为g(n) 的的N/2点点DFT，X(2r+1) 为为p(n)WNn 的的N/2点点DFT。这样，一个。这样，一个N点的点的DFT分解为两个分解为两个N/2点的点的DFT。将分解继续下去，直到分解为将分解继续下去，直到

16、分解为2点的点的DFT为止。当为止。当N=8时，时，基基2频率抽选的频率抽选的FFT算法的整个信号流图如图算法的整个信号流图如图4.6所示。所示。11)2()()(2N,0,n Nnxnxnpn将图将图4.6与图与图4.4比较，可知频率抽选法的计算量与时间抽选比较，可知频率抽选法的计算量与时间抽选法相同，而且都能够同址计算。时间抽选法是输入序列按奇法相同，而且都能够同址计算。时间抽选法是输入序列按奇偶分组，故偶分组，故x(n)的顺序要按倒序重排，而输出序列按前后分的顺序要按倒序重排，而输出序列按前后分半，故半，故X(k) 的顺序不需要重排；频率抽选法则是输出序列按的顺序不需要重排；频率抽选法则

17、是输出序列按奇偶分组，故奇偶分组，故X(k) 的顺序要按倒序重排，而输入序列按前后的顺序要按倒序重排，而输入序列按前后分半，故分半，故x(n) 不需要重排。不需要重排。 4.5 4.5 快速傅里叶反变换（快速傅里叶反变换（IFFTIFFT）nIFFT是是IDFT的快速算法。由于的快速算法。由于DFT的正变换和反变换的表的正变换和反变换的表达式相似，因此达式相似，因此IDFT也有相似的快速算法。可以用三种不也有相似的快速算法。可以用三种不同的方法来导出同的方法来导出IFFT算法。算法。n方法方法1n 设设 ， 则有：则有： ， n=0,1,N-1n=0,1,N-1)()(kXnxDFT 10)(

18、1)(NknkNWkXNnxn根据基根据基2 FFT的时间抽选法的第一次分解的结果：的时间抽选法的第一次分解的结果：)()()2()()()(kPWkGNkXkPWkGkX kNkN12, 1 , 0Nkn可以得到：可以得到：)2()(21)()2()(21)(NkXkXWkPNkXkXkG kN12, 1 , 0Nk 图图 4.8 由由 X(k)、X(k+N/2) 得到得到 G(k)、P(k)n再由再由N/2点的点的DFT求得求得N/4点的点的DFT，依次类推下去，就可推，依次类推下去，就可推到求出到求出x(n)的各点，如图的各点，如图4.9所示。将此流图与图所示。将此流图与图4.4比较，相

19、比较，相当于整个流向反过来，此外，因子当于整个流向反过来，此外，因子WNk成为成为WN-k-k，还增加了，还增加了因子因子1/2。但是，图。但是，图4.9的的IFFT算法不能直接利用按照图算法不能直接利用按照图4.4编编写的写的FFT算法程序，却可以利用图算法程序，却可以利用图4.6的频率抽选的频率抽选FFT算法的算法的程序，只需要将程序，只需要将X(k)作为输入序列，因子作为输入序列，因子WNk变为变为WN-k-k，并，并且将最后所得的输出序列的每个元素都除以且将最后所得的输出序列的每个元素都除以N。n方法方法2n 将将DFT的正变换式：的正变换式： n与其反变换式：与其反变换式： 10)(

20、)(NnnkNWnxkX10)(1)(NknkNWkXNnxn比较，很容易知道，可以利用比较，很容易知道，可以利用FFT算法的程序来算法的程序来计算计算IFFT，只需要将，只需要将X(k) 作为输入序列，作为输入序列，x(n) 则则是输出序列，另外将因子是输出序列，另外将因子WNk 变为变为WN-k-k， 当然，当然，最后还必须将输出序列的每个元素除以最后还必须将输出序列的每个元素除以N。n 方法方法3n 对对DFT的反变换式取共轭，有：的反变换式取共轭，有： 1, 1 , 0, )(1)(1)(10*10*NnWkXNWkXNnxNknkNNknkNn与与DFT的正变换式比较，可知完全可以利

21、用的正变换式比较，可知完全可以利用FFT的计算程序，只需要将的计算程序，只需要将X*(k) 作为输入序列，并将作为输入序列，并将最后结果取共轭，再除以最后结果取共轭，再除以N就得到就得到x(n)。4.7.1 两个长度相同的实序列两个长度相同的实序列n 可以将两个长度相同的实序列组合成一个复序列来进行可以将两个长度相同的实序列组合成一个复序列来进行FFT运算，从而一次完成这两个实序列的运算，从而一次完成这两个实序列的FFT，减少了总，减少了总的计算量。的计算量。n 设设p(n) 和和g(n) 是两个长度均为是两个长度均为N的实序列，并设：的实序列，并设： n又设又设 n ， ， n则由则由DFT

22、的线性有：的线性有： (4.36)()()(njgnpny)()(kPnpDFT )()(kGngDFT )()(kYnyDFT )()()(kjGkPkYnP(k) 和和G(k) 这两个复序列的实部都是周期性的偶序列，而这两个复序列的实部都是周期性的偶序列，而其虚部都是周期性的奇序列。其虚部都是周期性的奇序列。n 对复序列对复序列Y(k) 又有又有 (4.37)n这里下标这里下标 r、i 分别表示实部和虚部，分别表示实部和虚部，Y(k) 与其实部、虚部与其实部、虚部的长度都为的长度都为 N。现将。现将 (4.37) 式中各序列作周期延拓，有式中各序列作周期延拓，有 (4.38)(4.38)(

23、)()(kjYkYkYir)()()(kYjkYkYirn由周期性有：由周期性有： (4.39) (4.40)()(),()(kNYkY kNYkYrrrr)()(),()(kNYkY kNYkYiiiin现在将序列现在将序列 与与 作如下分解：作如下分解： (4.41) (4.42)(kYr)(kYi)()(21)()(21)(kNYkYkNYkYkYrrrrr)()(21)()(21)(kNYkYkNYkYkYiiiiin 由（由（4.394.39）式和（）式和（4.404.40）式，容易证明在）式，容易证明在(4.41)(4.41)和和(4.42)(4.42)这两个式子中，前一项都是偶序列，而后一项都是奇序列。这两个式子中，前一项都是偶序列，而后一项都是奇序列。n 将（将（4.414.41）式和（）式和（4.424.42）式代入（）式代入（4.384.38）式，并将各项）式，并将各项进行重新组合，得到：进行重新组合，得到：)( )( )()(21)()(21)()(21)()(21)(kGjkPkNYkYjkNY