数字信号处理-第4章

1第四章离散傅里叶变换快速算法FastFourierTransform(FFT)2快速傅里叶变换（FFT）是计算DFT的一种快速有效方法，并不是新的傅立叶变换式。从前面的讨论中看到，有限长序列在数字技术中占有很重要的地位。有限长序列的一个重要特点是其频域也可以离散化，即离散傅里叶变换（DFT）。3DFT是信号分析与处理中的一种重要变换。但直接计算DFT的计算量与变换区间长度N的平方成正比，当N较大时，计算量太大，直接用DFT算法进行谱分析和信号的实时处理是不切实际的。1965年库里和图基发现了DFT的一种快速算法，使DFT的运算效率提高1－2个数量级，为数字信号处理技术应用于各种信号的实时处理创造了条件，推动了数字处理技术的发展。1984年，提出了分裂基快速算法，使运算效率进一步提高；4一、直接计算DFT的问题及改进的方法DFT的定义两者形式类似，差别只在于的指数符号不同，及常数因子运算量是相同的10,)(1)(10,)()(1010NnWkXNnxNkWnxkXNknkNNnnkN51、运算量运算量复数乘法复数加法一个X(k)NN–1N个X(k)(N点DFT)N2N(N–1)实数乘法实数加法一次复乘42一次复加2一个X(k)4N2N+2(N–1)=2(2N–1)N个X(k)(N点DFT)4N22N(2N–1)(a+jb)(c+jd)=(ac-bd)+j(ad+bc)∑10)(NnnkNWnx672、减少运算量的途径之一nkNW具有如下特性：使DFT运算中的有些项可以合并。对称性周期性kNnNjkNnNnkNeWW)(2)(可约性mnkmNnkNnmkmNnkNnkNnkNWW)(knNNkNnNWW)()(kNNkNWW)2/(NNW2/,18减少运算量的途径之二可将长序列的DFT分解为短序列的DFT。一个长序列DFT乘加运算比多个短序列DFT乘加运算多长序列DFT可由多个短序列DFT组合而成9二、按时间抽取(DIT)的基-2FFT算法设MN2M为自然数将长度为N的序列x(n)按n的奇偶分成两组)12/(,1,0),12()()12/(,1,0),2()(21NrrxrxNrrxrx……1、算法原理10则x(n)的DFT为12/02212/02112/0)12(12/02)()()12()2()()()(NrkrNkNNrkrNNrrkNNrkrNnnknNknNWrxWWrxWrxWrxWnxWnxkX偶数奇数)()(21kXWkXkN=12/0Nk∑∑=12/2/212/2/1)()(NkrNkNNkrNWrxWWrxr=0r=01、算法原理11其中12/02/12/02/11)2()()(NrkrNNrkrNWrxWrxkX12/02/12/02/22)12()()(NrkrNNrkrNWrxWrxkX是x(2r)与x(2r+1)的N/2点DFT。1、算法原理12由于)()2()()()2(22112/0)2(2/11kXNkXkXWrxNkXNrNkrN得)()(22221221kXWkXNkXWNkXNkXkNNkN12,,1,0Nk…1、算法原理13信号流图将X1(k)和X2(k)合成X(k)运算可归结为：bWabWa蝶形运算的简化Wa+bWa-bW-W-1a+bWa-bWWabab(a)(b)X1(k)X2(k)X1(k)+WNkX2(k)X1(k)WNkX2(k)WNk图：蝶形运算符号144点基2时间抽取FFT算法流图x[0]x[2]x[1]x[3]X1[0]X1[1]X2[0]X2[1]2点DFT2点DFT04W14W02W02WX[0]X[1]X[2]X[3]1,0],[][][241mmXWmXmXm1,0],[][]2[241mmXWmXmXm15最后剩下的是2点DFT，它可以用一个蝶形结表示：1)1()0()1(1)1()0()0(12120202WxWxXWxWxXx(0)x(1)X(0)X(1)02W168点基2FFT算法N/2点DFTN/2点DFTx(0)x(2)x(4)x(6)x(1)x(3)x(5)x(7)X1(0)X1(1)X1(2)X1(3)X2(0)X2(1)X2(2)X2(3)WN0WN1WN2WN3X(0)X(1)X(2)X(3)X(4)X(5)X(6)X(7)1212()()()0,1,12()()()0,1,122kNkNNXkXkWXkkNNXkXkWXkk1212()()()0,1,12()()()0,1,122kNkNNXkXkWXkkNNXkXkWXkk｛N=8点的DIT-2FFT(时域抽取图)一次分解图17(3)第二次分解：将x1(r)按r取奇、偶可分解成2个长度为N/4的子序列x3(l)=x1(2l)、x4(l)=x1(2l+1)，根据上面推导可得：X1(k)=X3(k)+WN/2kX4(k)，k=0,1,…,N/2-1将x2(r)按r取奇、偶可分解成2个长N/4的子序列x5(l)=x2(2l)，l=0,1,…，N/4x6(l)=x2(2l+1)；同理得l=0,1,…，N/4-1；X1(k+N/2)=X3(k)WN/2kX4(k)，k=0,1,…,N/4-1；X1(k)=X3(k)+WN/2kX4(k)，k=0,1,…,N/4-1；X2(k)=X5(k)+WN/2kX6(k)，k=0,1,…,N/4-1；X2(k+N/2)=X5(k)WN/2kX6(k),k=0,1,…,N/4-1；18N/4点DFTx(0)x(4)x(2)x(6)x(1)x(5)x(3)x(7)X1(0)X1(1)X1(2)X1(3)X2(0)X2(1)X2(2)X2(3)WN0WN1WN2WN3X(0)X(1)X(2)X(3)X(4)X(5)X(6)X(7)N/4点DFTN/4点DFTN/4点DFTX3(0)X3(1)X4(0)X4(1)X5(0)X5(1)X6(0)X6(1)WN/20WN/21WN/21WN/20N=8点DFT的二次时域抽取分解图X1(k+N/2)=X3(k)WN/2kX4(k)X1(k)=X3(k)+WN/2kX4(k)X2(k)=X5(k)+WN/2kX6(k)X2(k+N/2)=X5(k)WN/2kX6(k)k=0,1,…,N/4-1；19再次分解，对N=8点，可分解三次。X(5)N=8点DIT-FFT运算流图x(0)x(4)x(2)x(6)x(1)x(5)x(3)X1(1)X1(2)X1(3)X2(0)X2(1)X2(2)X2(3)WN0WN1WN2WN3X(0)X(1)X(2)X(3)X(4)X(6)X3(1)X4(0)X4(1)X5(0)X5(1)X6(0)X6(1)WN/20WN/21WN/20WN/40WN/40WN/40x(7)WN/21WN/40L=1级L=2L=3X(7)X3(0)X1(0)202点DFT的蝶形结表示：)1()0()1()0()1()1()0()1()0()0(1202xWxxWxXxWxxWxXoNoNx(0)x(1)X(0)X(1)0NW218点基2DIT-FFT算法N=8点DIT-FFT运算流图x(0)x(4)x(2)x(6)x(1)x(5)x(3)WN0WN1WN2WN3X(0)X(1)X(2)X(3)X(4)X(5)X(6)WN0WN2WN0WN0WN0WN0x(7)WN2WN0L=1级L=2L=3X(7)22基2DIT-FFT算法由于这种方法每一步分解都是按输入时间序列是属于偶数还是奇数来抽取的，所以称为“按时间抽取法”或“时间抽取法”。23DIT―FFT算法与直接计算DFT运算量的比较1、直接DFT运算N点运算：复数乘次数：N×N复数加次数：N×(N-1)2、用DIT-FFT作N点运算：复数乘次数：M×N/2=N/2×log2N；复加次数：2×N/2×M=N×log2N。整个运算流图中有M级蝶形，每一级运算流图中有N/2个蝶形，每个蝶形需一次复乘和两次复数加运算。2425时间抽取法FFT的运算特点（1）蝶形运算（2）原位计算（3）序数重排（4）蝶形类型随迭代次数成倍增加26（1）蝶形运算对于N=2M，总是可以通过M次分解最后成为2点的DFT运算。这样构成从x(n)到X(k)的M级运算过程。从上面的流图可看到，每一级运算都由N/2个蝶形运算构成。因此每一级运算都需要N/2次复乘和N次复加，这样经过时间抽取后M级运算总共需要的运算：复乘复加而直接运算时则与N2成正比。NNMN2log22NNMN2log27（2）原位计算当数据输入到存储器中以后，每一级运算的结果仍然储存在同一组存储器中，直到最后输出，中间无需其它存储器，这叫原位计算。每一级运算均可在原位进行，这种原位运算结构可节省存储单元，降低设备成本，还可节省寻址的时间。28（3）序数重排对按时间抽取FFT的原位运算结构，当运算完毕时，正好顺序存放着x(0)，x(1)，x(2)，…x(7)，因此可直接按顺序输出，但这种原位运算的输入x（n）却不能按这种自然顺序存入存储单元中，而是按x(0)，x(4)，x(2)，x(6)，…，x(7)的顺序存入存储单元，这种顺序看起来相当杂乱，然而它也是有规律的。当用二进制表示这个顺序时，它正好是“码位倒置”的顺序。例如，原来的自然顺序应是x(1)的地方，现在放着x(4)，用二进制码表示这一规律时，则是在x（001）处放着x（100），x（011）处放着x（110）。29自然顺序二进码表示码位倒置码位倒序0000000010011004201001023011110641000011510110156110010371111117在实际运算中，一般直接将输入数据x（n）按码位倒置的顺序排好输入很不方便，总是先按自然顺序输入存储单元，然后再通过变址运算将自然顺序的存储转换成码位倒置顺序的存储，然后进行FFT的原位计算。目前有许多通用DSP芯片支持这种码位倒置的寻址功能。30x(000)x(100)x(010)x(110)010101x(001)x(101)x(011)x(111)01010101x(n2n1n0)n2n1n0造成倒位序的原因是输入序列x(n)，按标号(序号)n的奇偶不断地分组造成的。倒位序的形成31（4）蝶形类型随迭代次数成倍增加观察8点FFT的三次迭代运算:第一级迭代，有一种类型的蝶形运算系数W80，两个数据点间隔为1第二级迭代，有二种类型的蝶形运算系数W80、W82，参加运算的两个数据点间隔为2。第三级迭代，有四类蝶形运算系数W80、W81、W82、W83，参加运算的两个数据点间隔为4。结论：每迭代一次，蝶形类型增加一倍，数据点间隔也增大一倍。每一级的取数间隔和蝶形类型种类均为2i-1，i=1,2,…M。32N点DIT—FFT运算流图中，每级都有N/2个蝶形。每个蝶形都要乘以因子，称其为旋转因子，p称为旋转因子的指数。但各级的旋转因子和循环方式都有所不同，为了编写计算程序，应先找出旋转因子与运算级数的关系。用L表示从左到右的运算级数（L=1，2，、、，M）第L级共有个不同的旋转因子。N=8时各级旋转因子表示如下：3210,310,20,1222/24/，，，时，，时，时，JL33对的一般情况，第L级的旋转因子为MN212,,2,1,0,222212,,2,1,0,12J212MLJNNpNML