DSP-第三章23-24

1、第三章第三章 离散傅里叶离散傅里叶变换及其快速算法变换及其快速算法3.3 快速傅里叶变换快速傅里叶变换 (FFT) 从前面的讨论中看到，有限长序列在数字技术中占从前面的讨论中看到，有限长序列在数字技术中占有很重要的地位。有限长序列的一个重要特点是其频有很重要的地位。有限长序列的一个重要特点是其频域可离散化，即离散傅里叶变换域可离散化，即离散傅里叶变换(DFT)。 虽然频谱分析和虽然频谱分析和DFT运算很重要，但在很长一段时运算很重要，但在很长一段时间里，由于间里，由于DFT运算复杂，并没有得到真正的运用，运算复杂，并没有得到真正的运用，而频谱分析仍大多采用模拟信号滤波的方法解决，直而频谱分析仍

2、大多采用模拟信号滤波的方法解决，直到到1965年首次提出年首次提出DFT运算的一种快速算法以后，情运算的一种快速算法以后，情况才发生了根本变化，人们开始认识到况才发生了根本变化，人们开始认识到DFT运算的一运算的一些内在规律，从而很快地发展和完善了一套高速有效些内在规律，从而很快地发展和完善了一套高速有效的运算方法的运算方法快速傅里叶变换快速傅里叶变换(FFT)算法。算法。 FFT的出现，使的出现，使DFT的运算大大简化，运算时间的运算大大简化，运算时间缩短缩短12个数量级，使个数量级，使DFT的运算在实际中得到广的运算在实际中得到广泛应用。泛应用。一、提出问题：设有限长序列一、提出问题：设有

3、限长序列x(n),非零值长度为非零值长度为N, 对对x(n)进行进行N点点DFT运算，共需多大的运算量运算，共需多大的运算量?10( ) ( )( )NnkNnX kDFT x nx n W k=0,1,2,N1n计算一个计算一个X(k)(一个频率成分一个频率成分)值，运算量为值，运算量为 例如，例如，k=1，则则 要进行要进行N次复数乘法次复数乘法和和(N-1)次复数加法次复数加法；n要完成整个要完成整个DFT运算，运算，需要计算需要计算N个个X(k)值值，所，所以以总总计算量为：计算量为： N*N次复数相乘和次复数相乘和N*(N-1)次复数加法次复数加法10(1)( )NnNnXx n W

4、 n在这些运算中，乘法比加法复杂，需要的运算时间在这些运算中，乘法比加法复杂，需要的运算时间多，尤其是复数相乘，每个复数相乘包括多，尤其是复数相乘，每个复数相乘包括4个实数个实数相乘和相乘和2个实数相加，例个实数相加，例 10( )( )( )( )( )NnknkeeNmmNnnknkemNmeNX kRx nRwIx nIwj Rx nIwIx nRw 又每个复数相加包括又每个复数相加包括2个实数相加，所以，每计算一个实数相加，所以，每计算一个个X(k)要进行要进行4N次实数相乘和次实数相乘和2N+2(N-1)=2(2N-1)次次实数相加。因此，整个实数相加。因此，整个DFT运算需要运算需

5、要 4N2实数相乘和实数相乘和2N(2N-1)次实数相加次实数相加 从上面的分析看到，在从上面的分析看到，在DFT计算中，不论是乘法计算中，不论是乘法和加法，运算量均与和加法，运算量均与N2成正比。因此，成正比。因此，N较大时，较大时，运算量十分可观。运算量十分可观。 反变换反变换IDFT与与DFT的运算结构相同，只是多乘一的运算结构相同，只是多乘一个常数个常数1/N，所以二者的计算量相同。，所以二者的计算量相同。例例1:当当N=1024点时，直接计算点时，直接计算DFT需要：需要： N2=220=1048576次，即一百多万次的复乘运算次，即一百多万次的复乘运算 例例2:石油勘探，石油勘探，

6、24道记录，每道波形记录长度道记录，每道波形记录长度5秒，若秒，若 每秒抽样每秒抽样500点点/秒，秒， 每道总抽样点数每道总抽样点数=500*5=2500点点 24道总抽样点数道总抽样点数=24*2500=6万点万点 DFT运算需要运算需要N2=(60000)2=36*108次次 如此之大的运算量，对计算机的处理能力和处理速度如此之大的运算量，对计算机的处理能力和处理速度都提出很高的要求；尤其在实时性很强的信号处理都提出很高的要求；尤其在实时性很强的信号处理(如雷如雷达信号处理达信号处理)中，对计算速度的要求十分苛刻；中，对计算速度的要求十分苛刻；迫切需要改进迫切需要改进DFT的计算方法，以

7、减少总的运算次数的计算方法，以减少总的运算次数二、改善二、改善DFT运算效率的基本途径运算效率的基本途径1利用利用nkNW的特点的特点2jNNWe 具有具有 1)周期性周期性 ()()nkn N kn kNNNNWWW 2)共轭性共轭性 ()()()nkNn kn NkNNNWWW 3)对称性对称性 ()2NkkNNWW 4)221NNNNWW 5)22nnNNWW ReImj8N0NW7NW6NW5NW4NW3NW2NW1NW2、FFT的思路的思路利用对称性和周期性将长序利用对称性和周期性将长序列列DFT分解为短序列分解为短序列DFTn因为因为DFT的运算量与的运算量与N2成正比的成正比的n

8、如果一个大点数如果一个大点数N的的DFT能分解为若干小点能分解为若干小点数数DFT的组合，则显然可以达到减少运算工的组合，则显然可以达到减少运算工作量的效果作量的效果n运算量分解示意图如下：运算量分解示意图如下：其运算量为：其运算量为：22)()()(NNNkXkXkX把把N点数据分成二半：点数据分成二半：2)2(N2)2(N2N再分二半：再分二半：44)()(NNkXkX44)()(NNkXkX+=22N2)4(N2)4(N2)4(N2)4(N+=24N这样一直分下去，剩下这样一直分下去，剩下两点的变换两点的变换。三、常见的三、常见的FFT算法算法n按抽取方法分：按抽取方法分： 时间抽取法时

9、间抽取法(DIT); 频率抽取法频率抽取法(DIF);n按按“基数基数”分：分： 基基-2FFT算法；算法； 基基-4FFT算法；算法； 混合基混合基FFT算法；算法；3.3.1 按时间抽取的按时间抽取的FFT1、先从、先从一一特殊情况开始，假定特殊情况开始，假定N是是2的整数次方，即的整数次方，即 N=2M，M：正整数：正整数 将序列将序列x(n)分解为两组分解为两组偶数项和奇数项偶数项和奇数项 r=0,1，N/2-1 12(2 )( )(21)( )xrx rxrx r 2101()()NNn kn kNNnnxnwxnw 偶偶奇奇 10X( )( )( )NnkNnkDFTx nx n

10、w DFT运算也相应分为两组：运算也相应分为两组：/2 1/2 12(21)00(2 )(21)NNrkrkNNrrxr wxrw /2 1/2 12200(2 )(21)NNrkkrkNNNrrxr wwxrw 因为因为 故故 其中其中 /2 102/2 102(2 )(21)NrkNrNrkNrG kxr WH kxrW 222222jnNjnnnNNNweew /2 1/2 10022(2 )(21)NNrkkrkNNNrrkNX kxr WWxrWG kW H kk=0,1,N/2-1应用系数应用系数wkN的周期性和对称性：的周期性和对称性：X(k)的后的后N/2个点（个点（ 即即N/

11、2N-1点），表示为点），表示为 ()/2222Nkr NkrkkNNNNwwWW 222222,0,1,12kNNkNG kNG kH kNH kX kNG kNWH kNNG kW H kk ,0,1,12(),0,1,122kNkNNX kG kW H kkNNX kG kW H kk依此类推，依此类推，G(k)和和H(k)可以继续分下去，这种可以继续分下去，这种按时间抽取算法是在输入序列分成越来越小的按时间抽取算法是在输入序列分成越来越小的子序列上执行子序列上执行DFT运算，最后再合成为运算，最后再合成为点的点的DFT。 一个一个N点的点的DFT被分解为两个被分解为两个N/2点的点的D

12、FT，这两，这两个个N/2点的点的DFT再合成为一个再合成为一个N点点DFT.2、蝶形信号流图、蝶形信号流图将将G(k)和和H(k) 合成合成X(k)运算可归结为：运算可归结为： abWabW Wa+bWa-bW-W1a+bWa-bWWabab蝶 形 运 算 的 简蝶 形 运 算 的 简化化(a)(a)(b)图图 (a)为实现这一运算的一般方法，它需要两次乘法、为实现这一运算的一般方法，它需要两次乘法、两次加减法。考虑到两次加减法。考虑到-bW和和bW两个乘法仅相差一负两个乘法仅相差一负号，可将图号，可将图 (a)简化成图简化成图 (b)，此时仅需一次乘法、，此时仅需一次乘法、两次加减法。两次

13、加减法。图图 (b)的运算结构像一蝴蝶通常称作蝶形运算结构简的运算结构像一蝴蝶通常称作蝶形运算结构简称蝶形结，采用这种表示法，就可以将以上所讨论称蝶形结，采用这种表示法，就可以将以上所讨论的分解过程用流图表示的分解过程用流图表示。3、N=23=8 的例子的例子。 按照这个办法，继续把按照这个办法，继续把N/2用除，由于用除，由于N=2M，仍，仍然是偶数，可以被整除，因此，可以对两个然是偶数，可以被整除，因此，可以对两个N/2点点的的DFT再分别作进一步的分解。再分别作进一步的分解。 即对即对G(k)和和H(k)的计算，又可以分别通过计算的计算，又可以分别通过计算两个长度为两个长度为N/4=2点

14、的点的DFT，进一步节省计算量，见，进一步节省计算量，见图。这样，一个点的图。这样，一个点的DFT就可以分解为四个点的就可以分解为四个点的DFT。 最后剩下的是最后剩下的是2点点DFT，它可以用一个蝶形结表示：，它可以用一个蝶形结表示： 这样，一个这样，一个8点的完整的按时间抽取运算的流图点的完整的按时间抽取运算的流图 由于这种方法每一步分解都是按输入时间序列是属由于这种方法每一步分解都是按输入时间序列是属于偶数还是奇数来抽取的，所以称为于偶数还是奇数来抽取的，所以称为“按时间抽取法按时间抽取法”或或“时间抽取法时间抽取法”。) 1 ()0() 1 ()0() 1 () 1 ()0() 1 (

15、)0()0(1202xWxxWxXxWxxWxXoNoN按时间抽取的按时间抽取的8点点FFT4、时间抽取法、时间抽取法FFT的运算特点：的运算特点：(1)蝶形运算蝶形运算(2)原位计算原位计算 (3)序数重排序数重排(4)蝶形类型随迭代次数成倍增加蝶形类型随迭代次数成倍增加(1)蝶形运算蝶形运算 对于对于N=2M，总是可以通过，总是可以通过M次分解最后成为次分解最后成为2点的点的DFT运算。这样构成从运算。这样构成从x(n)到到X(k)的的M级运算过程。级运算过程。 从上面的流图可看到，每一级运算都由从上面的流图可看到，每一级运算都由N/2个蝶形运个蝶形运算构成。因此每一级运算都需要算构成。因

16、此每一级运算都需要N/2次复乘和次复乘和N次复次复加，这样，经过时间抽取后加，这样，经过时间抽取后M级运算总共需要的运算级运算总共需要的运算： 复乘复乘 复加复加 而直接运算时则与而直接运算时则与N2 成正比。成正比。 例例N=2048，N2=4194304，(N/2)log2N=11264，N2/ (N/2)log2N=392.4。FFT显然要比直接法快得多。显然要比直接法快得多。log222NNMNlog2NMNN(2)原位计算原位计算 当数据输入到存储器中以后，每一级运算的结果当数据输入到存储器中以后，每一级运算的结果仍然储存在同一组存储器中，直到最后输出，中间仍然储存在同一组存储器中，

17、直到最后输出，中间无需其它存储器，这叫原位计算。无需其它存储器，这叫原位计算。 每一级运算均可在原位进行，这种原位运算结构每一级运算均可在原位进行，这种原位运算结构可节省存储单元，降低设备成本，还可节省寻址的可节省存储单元，降低设备成本，还可节省寻址的时间。时间。(3)序数重排序数重排 对按时间抽取对按时间抽取FFT的原位运算结构，当运算完毕时，的原位运算结构，当运算完毕时，正好顺序存放着正好顺序存放着X(0)，X(1)，X(2)，X(7)，因此可，因此可直接按顺序输出，但这种原位运算的输入直接按顺序输出，但这种原位运算的输入x(n)却不能按却不能按这种自然顺序存入存储单元中，而是按这种自然顺

18、序存入存储单元中，而是按x(0)，x(4)，x(2)，x(6)，x(7)的顺序存入存储单元，这种顺序看起来的顺序存入存储单元，这种顺序看起来相当杂乱，然而它也是有规律的。当用二进制表示这个相当杂乱，然而它也是有规律的。当用二进制表示这个顺序时，它正好是顺序时，它正好是“码位倒置码位倒置”的顺序。例如，原来的的顺序。例如，原来的自然顺序应是自然顺序应是x(1)的地方，现在放着的地方，现在放着x(4)，用二进制码表，用二进制码表示这一规律时，则是在示这一规律时，则是在 x(0 0 1)处放着处放着 x(1 0 0)， x(0 1 1)处放着处放着 x(1 1 0)。自然顺序自然顺序二进码表示二进码

19、表示码位倒置码位倒置码位倒置顺序码位倒置顺序0000000010011004201001023011110641000011510110156110010371111117 在实际运算中，一般直接将输入数据在实际运算中，一般直接将输入数据x(n)按码位倒置的顺按码位倒置的顺序排好输入很不方便，总是先按自然顺序输入存储单元，然序排好输入很不方便，总是先按自然顺序输入存储单元，然后再通过变址运算将自然顺序的存储转换成码位倒置顺序的后再通过变址运算将自然顺序的存储转换成码位倒置顺序的存储，然后进行存储，然后进行FFT的原位计算。目前有许多通用的原位计算。目前有许多通用DSP芯片芯片支持这种码位倒置的

20、寻址功能。支持这种码位倒置的寻址功能。n第一次分偶、奇，根据最低位第一次分偶、奇，根据最低位n0的的0、1状态来分，状态来分，若若n0=0，则为偶序列；，则为偶序列；n0=1则为奇序列，得到两组则为奇序列，得到两组序列：序列： 000 010 100 110 001 011 101 111n第二次对这两个偶、奇序列再分一次偶、奇序列，这第二次对这两个偶、奇序列再分一次偶、奇序列，这就要根据就要根据n1的、状态。若的、状态。若n1=0，则为偶序列；，则为偶序列；n1=1则为奇序列，得到四组序列：则为奇序列，得到四组序列： 000 100 010 110 001 101 011 111n同理，再根

21、据同理，再根据n2的、状态来分偶、奇序列，直到的、状态来分偶、奇序列，直到不能再分偶、奇时为止。对于不能再分偶、奇时为止。对于N=8, n2已是最高位，已是最高位，最后一次分得结果为最后一次分得结果为 000 100 010 110 001 101 011 111(4)蝶形类型随迭代次数成倍增加蝶形类型随迭代次数成倍增加 观察观察8点点FFT的三次迭代运算的三次迭代运算: 第一级迭代，有一种类型的蝶形运算系数第一级迭代，有一种类型的蝶形运算系数W08，两个数据点间隔为两个数据点间隔为1； 第二级迭代，有二种类型的蝶形运算系数第二级迭代，有二种类型的蝶形运算系数W08、W28 ，参加运算的两个数

22、据点间隔为参加运算的两个数据点间隔为2； 第三级迭代，有四类蝶形运算系数第三级迭代，有四类蝶形运算系数W08、W18、W28、W38，参加运算的两个数据点间隔为，参加运算的两个数据点间隔为4； 结论：每迭代一次，蝶形类型增加一倍，数据点间结论：每迭代一次，蝶形类型增加一倍，数据点间 隔也增大一倍。每一级的取数间隔和蝶形类隔也增大一倍。每一级的取数间隔和蝶形类 型种类均为型种类均为2i-1，i=1,2,M。 3.3.2按频率抽取的按频率抽取的FFT(按输出按输出X(k)在频域的在频域的顺序上属于偶数还是奇数分解为两组顺序上属于偶数还是奇数分解为两组) 对于对于N=2M情况下的另外一种普遍使用的情

23、况下的另外一种普遍使用的FFT结构是结构是频率抽取法。对于频率抽取法，输入序列不是按偶数奇频率抽取法。对于频率抽取法，输入序列不是按偶数奇数，而是按前后对半分开，这样便将数，而是按前后对半分开，这样便将N点点DFT写成前后写成前后两部分：两部分： /( )( )( )2 1102NNnknkNNnn NX kx n Wx n W/()( )()2 12 12002NNNnknkNNnnNx n Wx nW/(/ ) ( )(/ )2 1202NNknkNNnx nWx nNW又又把把X(k)进一步分解为偶数组和奇数组：进一步分解为偶数组和奇数组： X(2r+1)= 奇数为偶数1, 1) 1(,

24、 1)2/(2/kWWkkNNNN/( ) ( )()(/ )2 1012NknkNnX kx nx nNW /() ( )(/ )2 12022NnrNnXrx nx nNW/ ( )(/ )2 1202NnrNnx nx nNW/() ( )(/ )2 12102NnrNnx nx nNW/ ( )(/ )2 1202NnnrNNnx nx nNW W令令 a(n)=x(n)+x(n+N/2)，b(n)=x(n)-x(n+N/2wnN这两个序列都是这两个序列都是N/2点的序列，将其代入上两式，得点的序列，将其代入上两式，得 这正是两个这正是两个N/2点的点的DFT运算，即将一个运算，即将一

25、个N点的点的DFT分解分解为两个为两个N/2点的点的DFT，上式的运算关系可用下图表示，上式的运算关系可用下图表示. /()( )2 1202NnrNnXra n W/()( )2 12021NnrNnXrb n Wx1(n)+ x2(n)WNnx1(n)x2(n)x1(n)-x2(n) WNn-1以以N=8的频率抽取为例的频率抽取为例x(0)x(1)x(2)x(3)x(4)x(5)x(6)x(7)-1-1-1-1N/2点点DFTa(0)a(1)a(2)a(3)N/2点点DFTb(0)b(1)b(2)b(3)WN1WN2WN3X(0)X(2)X(4)X(6)X(1)X(3)X(5)X(7)按频

26、率抽取将按频率抽取将8点点DFT分解成两个分解成两个4点点DFT按频率抽取将按频率抽取将8点点DFT分解成四个分解成四个2点点DFT 与时间抽取法一样，由于与时间抽取法一样，由于N=2M，N/2仍是一个偶仍是一个偶数数,这样，一个这样，一个N=2M点的点的DFT通过通过M次分解后，最次分解后，最后只剩下全部是后只剩下全部是2点的点的DFT，2点点DFT实际上只有加实际上只有加减运算。但为了比较，也为了统一运算的结构，仍减运算。但为了比较，也为了统一运算的结构，仍用一个系数为用一个系数为W0N 的蝶形运算来表示。的蝶形运算来表示。 频率抽取法的流图是时域抽取法流图的左右翻转频率抽取法的流图是时域

27、抽取法流图的左右翻转。下图是。下图是N=8的频率抽取法的频率抽取法FFT流图。流图。 N=8的按频率抽取的按频率抽取FFT运算流图运算流图频率抽取法频率抽取法FFT的运算特点：的运算特点：(1)蝶形运算蝶形运算 对于任何一个对于任何一个2的整数幂的整数幂N=2M，总是可以通过，总是可以通过M次分解最后完全成为次分解最后完全成为2点的点的DFT运算。这样的运算。这样的M次分次分解，就构成从解，就构成从x(n)到到X(k)的的M级运算过程。将频率级运算过程。将频率抽取法的流图反转，并将输入变输出，输出变输入抽取法的流图反转，并将输入变输出，输出变输入，得到的便是时间抽取法流图，得到的便是时间抽取法

28、流图(反映了时域，频域的反映了时域，频域的对称法对称法)。频率抽取法也共有。频率抽取法也共有M级运算级运算(N=2M)，其运，其运算量与时间抽取法相同。算量与时间抽取法相同。(2)原位计算原位计算 类似于时间抽取法，当数据输入到存储器中以后类似于时间抽取法，当数据输入到存储器中以后，每一级运算的结果仍然储存在同一组存储器中，每一级运算的结果仍然储存在同一组存储器中，直到最后输出，中间无需其它存储器，所以频域抽直到最后输出，中间无需其它存储器，所以频域抽取法也可进行原位计算。取法也可进行原位计算。(3)序数重排序数重排 它的输入正好是自然顺序。但它的输出却是码位倒置它的输入正好是自然顺序。但它的

29、输出却是码位倒置的顺序。因此运算完毕后，要通过变址运算将码位倒置的顺序。因此运算完毕后，要通过变址运算将码位倒置的顺序转换为自然顺序，然后输出，变址方法同时间抽的顺序转换为自然顺序，然后输出，变址方法同时间抽取法。取法。(4)蝶形类型随迭代次数成倍减少蝶形类型随迭代次数成倍减少(与时间抽取法相反与时间抽取法相反) 第一级迭代中有第一级迭代中有N/2种蝶形运算系数，参加蝶形运算种蝶形运算系数，参加蝶形运算的两个数据相隔的两个数据相隔N/2，随后每次迭代，蝶形类形比前一，随后每次迭代，蝶形类形比前一级减少一倍，间距也减少一倍，最后一级迭代，蝶形类级减少一倍，间距也减少一倍，最后一级迭代，蝶形类形只

30、有一种形只有一种W0N，数据间隔为，数据间隔为1。由这几点规律可以看出，频率抽取法与时间抽到法是两由这几点规律可以看出，频率抽取法与时间抽到法是两种等价的种等价的FFT运算。运算。3.3.3 N为组合数的为组合数的FFT(任意基数的任意基数的FFT算法算法) 以上讨论的都是以以上讨论的都是以2为基数的为基数的FFT算法，即算法，即N=2M，这种情况实际上使用得最多。这种情况实际上使用得最多。 优点：程序简单，效率高，使用方便。优点：程序简单，效率高，使用方便。 实际应用时，有限长序列的长度实际应用时，有限长序列的长度N很大程度上由人很大程度上由人为因素确定，因此多数场合可取为因素确定，因此多数

31、场合可取N=2M，从而直接使用，从而直接使用以以2 为基数的为基数的FFT算法。算法。 如要求准确的如要求准确的N点点DFT值值,可采用任意数为基数的可采用任意数为基数的 FFT算法算法 , 其其 计算效率低于以计算效率低于以2为基数为基数FFT算法。算法。 如如N为复合数，可分解为两个整数为复合数，可分解为两个整数p与与q的乘积的乘积,像前像前面以面以2为基数时一样为基数时一样,FFT的基本思想是将的基本思想是将DFT的运算的运算尽量分小尽量分小,因此因此,在在N=pq情况下情况下,也希望将也希望将N点的点的DFT分分解为解为p个个q点点DFT或或q个个p点点DFT,以减少计算量。以减少计算

32、量。如如N不能人为确定不能人为确定 , N的数值也不是以的数值也不是以2为基数的整数次为基数的整数次方方,处理方法有两种处理方法有两种: 补零补零: 将将x(n)补零补零 , 使使N=2M. 例如例如N=30,补上补上x(30)=x(31)=0两点两点,使使N=32=25,这样可这样可直接采用以直接采用以2为基数为基数M=5的的FFT程序。程序。 3.4 FFT应用中的几个问题应用中的几个问题1)IDFT的运算方法的运算方法101( )( )( )NnkNkx nIDFT X kX k WN 10X(k)=DFTx(n) =()NnkNnx n W 以上所讨论的以上所讨论的FFT算法可用于算法

33、可用于IDFT运算运算简简称为称为IFFT。 比较比较IDFT的定义式：的定义式： IDFT： DFT：IDFT与与DFT的差别的差别 1)把把DFT中的每一个系数中的每一个系数 改为改为 ， 2)再乘以常数再乘以常数 1/N ， 则以上所讨论的时间抽取或频率抽取的则以上所讨论的时间抽取或频率抽取的FFT运算均可直接用于运算均可直接用于IDFT运算，当然，蝶形中运算，当然，蝶形中的系数的系数 应改为应改为 。nkNWnkNW nkNW nkNW b)第二种方法，完全不需要改动第二种方法，完全不需要改动FFT程序，而是直接利程序，而是直接利 用它作用它作 IFFT。 考虑到考虑到 故故 IFFT

34、计算分三步：计算分三步： 将将X(k)取共轭取共轭(虚部乘以虚部乘以-1) 对对 直接作直接作FFT 对对FFT的结果取共轭并乘以的结果取共轭并乘以1/N，得，得x(n)。101( )( )NnkNkxnXk WN 1011( )( )( )NnkNkx nXk WDFT XkNN ( )Xk 3.4 FFT应用中的几个问题应用中的几个问题2)实数序列的实数序列的FFT以上讨论的以上讨论的FFT算法都是复数运算算法都是复数运算,包括序列包括序列x(n)也认也认为是复数为是复数,但大多数场合但大多数场合,信号是实数序列信号是实数序列,任何实数都任何实数都可看成虚部为零的复数可看成虚部为零的复数,

35、例如例如,求某实信号求某实信号x(n)的复谱的复谱,可认为是将实信号加上数值为零的虚部变成复信号可认为是将实信号加上数值为零的虚部变成复信号(x(n)+j0),再用再用FFT求其离散傅里叶变换。这种作法很求其离散傅里叶变换。这种作法很不经济不经济,因为把实序列变成复序列因为把实序列变成复序列,存储器要增加一倍存储器要增加一倍,且计算机运行时且计算机运行时,即使虚部为零即使虚部为零,也要进行涉及虚部的也要进行涉及虚部的运算运算,浪费了运算量。合理的解决方法是利用复数据浪费了运算量。合理的解决方法是利用复数据FFT对实数据进行有效计算对实数据进行有效计算,下面介绍两种方法。下面介绍两种方法。 (1

36、)用用 一个一个N点点FFT同时计算两个同时计算两个N点实序列的点实序列的DFT 设设x (n)、y (n)是彼此独立的两个是彼此独立的两个N点实序列点实序列,且且 X (k)=DFTx (n)，Y (k)=DFTy(n) 则则X (k)、Y(k)可通过一次可通过一次FFT运算同时获得。运算同时获得。 首先，将首先，将x (n)、y(n)分别当作一复序列的实部及虚分别当作一复序列的实部及虚部部, 令令 g(n)=x (n)+jy(n) riririirriG kX kjY kXkjXkjYkYkXkYkjXkjYkGkjGk *DFT Re121122rriig nX kG kGNkGkGNk

37、j GkGNk 通过通过FFT运算可获得运算可获得g(n)的的DFT值值利用离散傅里叶变换的共轭对称性利用离散傅里叶变换的共轭对称性通过通过g(n)的的FFT运算结果运算结果G(k),由上式可得到由上式可得到X (k) 的值。的值。通过通过g(n)FFT运算结果运算结果G(k),由上式也可得到由上式也可得到Y (k) 的值。的值。 *DFTIm121122rriijg njY kG kGNkG kG Nkj G kG Nk 1122iirrY kG kG N kj G kG N k 做一次做一次点复序列的点复序列的FFT，再通过加、减法运算就可以，再通过加、减法运算就可以将将X(k)与与Y(k

38、)分离出来。这将使运算效率提高一倍。分离出来。这将使运算效率提高一倍。 (2)用一个用一个N点的点的FFT运算获得一个运算获得一个2N点实序列的点实序列的DFT 设设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)组成一个复序列：组成一个复序列： y(n)=x1(n)+jx2(n)通过通过N点点FFT运算可得到：运算可得到： Y(k)=X1(k)+jX2(k) ，N点点根据前面的讨论，得到根据

39、前面的讨论，得到 *1*21( )( )()2( )( )()2XkY kYNkjXkY kYNk 为求为求 2N 点点 x(n)所对应所对应 X(k)，需求出，需求出 X(k)与与 X1(k)、X2(k)的关系的关系 ： 而而 21112(21)22000( )( )(2 )(21)NNNnknknkNNnnnX kx n Wxn WxnW 11200(2 )(21)NNnkknkNNNnnxn WWxnW 1111001122100( )( )(2 )( )( )(21)NNnknkNNnnNNnknkNNnnXkx n Wxn WXkxn WxnW 所以所以1)由由x1(n)及及x2(n

40、)组成复序列，经组成复序列，经FFT运算求得运算求得 Y(k)，2)利用共轭对称性求出利用共轭对称性求出 X1(k)、X2(k)，3)最后利用上式求出最后利用上式求出 X(k), 达到用一个达到用一个N点的点的FFT计算计算一个一个2N点的实序列的点的实序列的DFT的目的。的目的。 X(k)=X1(k)+W2Nk X2(k) 3) 线性卷积的线性卷积的FFT算法算法 线性卷积是求离散系统响应的主要方法之一线性卷积是求离散系统响应的主要方法之一,许多重要许多重要应用都建立在这一理论基础上应用都建立在这一理论基础上,如卷积滤波等。如卷积滤波等。 以前曾讨论了用圆周卷积计算线性卷积方法归纳如下以前曾

41、讨论了用圆周卷积计算线性卷积方法归纳如下： 将长为将长为N2的序列的序列x(n)延长到延长到L,补补L-N2个零，个零， 将长为将长为N1的序列的序列h(n)延长到延长到L,补补L-N1个零，个零， 如果如果LN1+N2-1,则圆周卷积与线性卷积相等则圆周卷积与线性卷积相等,此时此时,可可用用FFT计算线性卷积，方法如下计算线性卷积，方法如下: a. 计算计算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就可完成线就可完

42、成线性卷积计算。性卷积计算。 计算表明计算表明, L32时时, 上述计算线性卷积的方法上述计算线性卷积的方法比直接计算线卷积有明显的优越性比直接计算线卷积有明显的优越性, 因此因此, 也称上也称上述循环卷积方法为快速卷积法。述循环卷积方法为快速卷积法。 上述结论适用于上述结论适用于 x(n)、h(n) 两序列长度比较接近或两序列长度比较接近或相等的情况相等的情况,如果如果x(n)、h(n) 长度相差较多长度相差较多,例如例如, h(n) 为某滤波器的单位脉冲响应为某滤波器的单位脉冲响应,长度有限长度有限,用来处理一个用来处理一个很长的输入信号很长的输入信号 x(n),或者处理一个连续不断的信号

43、，或者处理一个连续不断的信号，按上述方法按上述方法, h(n) 要补许多零再进行计算要补许多零再进行计算,计算量有很计算量有很大的浪费，或者根本不能实现。大的浪费，或者根本不能实现。 为了保持快速卷积法的优越性为了保持快速卷积法的优越性,可将可将 x(n) 分为许多分为许多段段,每段的长度与每段的长度与 h(n) 接近接近 , 处理方法有两种：处理方法有两种：(1) 重叠相加法重叠相加法由分段卷积的各段相由分段卷积的各段相加构成总的卷积输出加构成总的卷积输出h(n)x(n) 假定假定 xi(n) 表示表示 x(n)序列的第序列的第i段段 ： 则输入序列可表为：则输入序列可表为： 于是输出可分解

44、为：于是输出可分解为： 其中其中 22( )(1)1( )0ix niNniNx n ( )( )iix nx n ( )( )* ( )( )* ( )( )iiiiy nx nh nx nh ny n ( )( )* ( )iiy nx nh n 1)先对先对h(n)及及xi(n)补零，补到具有补零，补到具有N点长度，点长度，N=N1+N2-1。一般选。一般选 N=2M。 2)用基用基2FFT计算计算 yi(n)=xi(n)*h(n)。 3)重叠部分相加构成最后的输出序列。重叠部分相加构成最后的输出序列。 由于由于yi(n)的长度为的长度为N，而，而xi(n)的长度为的长度为N2，因此相，

45、因此相邻两段邻两段 yi(n)序列必然有序列必然有N-N2=N1-1点发生重叠。点发生重叠。 ( )( )iiy ny n 计算步骤：计算步骤：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. 将重叠部分相加将重叠部分相加(2)重叠保留重叠保留 这种方法和第一种方法稍有不同，即将上面分段序这种方法和第一种方法稍有不同，即将上面分段序列中补零的部分不是补零，而是保留原来的输入序列列中补零的部分不是补零，而是保

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

47、延差等。 两个长为两个长为N的实离散时间序列的实离散时间序列 x(n)与与y(n)的互相的互相关函数定义为关函数定义为 1100( )() ( )( ) ()NNxynnrx ny nx n y n 则可以证明，则可以证明，rxy()的离散傅里叶变换为的离散傅里叶变换为 Rxy(k)= 其中其中 X(k)=DFTx(n), Y(k)=DFTy(n), Rxy(k)=DFTrxy()， 0kN-1*X (k)Y(k)互相关函数定义为互相关函数定义为 x(n)及及y(n)的卷积公式的卷积公式 相比较，我们可以得到相关和卷积的时域关系：相比较，我们可以得到相关和卷积的时域关系： 1100()() (

48、 )( ) ()NNxynnrmx nm y nx n y nm10()() ( )()()Nnf mx mn y nx my m 1100()() ( ) () ( )()()NNxynnrmx nm y nxmn y nxmy m 211001( )( ) ()( ) ( )NNjkNxynkrx n y nXk Y k eN 因因 得得 证毕。证毕。DFT ()( )*( )NNxnRnXk 当当x(n)=y(n)时，得到时，得到x(n)的自相关函数为：的自相关函数为： 维纳维纳辛钦定理：辛钦定理： 自相关函数与信号功率谱互为傅立叶变换对。自相关函数与信号功率谱互为傅立叶变换对。 10( )( ) ()Nxxnrx n x n 21201|() |NjkNkX keN 上面的推导表明，互相关和自相关函数的计算可上面的推导表明，互相关和自相关函数的计算可利用利用FFT实现。实现。 由