线性卷积的fft算法及其matlab实现.docx

线性卷积的FFT算法线性卷积是求离散系统响应的主要方法之一,许多重要应用都建立在这一理论基础上,如卷积滤波等。以前曾讨论了用循环卷积计算线性卷积的方法归纳如下将长为N2的序列xn延长到L,补L-N2个零将长为N1的序列hn延长到L,补L-N1个零如果L≥N1N2-1,则循环卷积与线性卷积相等,此时,可有FFT计算线性卷积,方法如下a.计算XkFFT[xn]b.求HkFFT[hn]c.求YkHkYkk0～L-1d.求ynIFFT[Yk]n0～L-1可见，只要进行二次FFT,一次IFFT就可完成线性卷积计算。计算表明,L32时,上述计算线性卷积的方法比直接计算线卷积有明显的优越性,因此,也称上述循环卷积方法为快速卷积法。上述结论适用于xn,hn两序列长度比较接近或相等的情况,如果xn,hn长度相差较多,例如,hn为某滤波器的单位脉冲响应,长度有限,用来处理一个很长的输入信号xn,或者处理一个连续不断的信号，按上述方法,hn要补许多零再进行计算,计算量有很大的浪费，或者根本不能实现。为了保持快速卷积法的优越性,可将xn分为许多段后处理,每小段的长与hn接近,其处理方法有两种1重叠相加法由分段卷积的各段相加构成总的卷积输出假定xin表示图中第i段x（n）序列则输入序列可表为于是输出可分解为其中由此表明，只要将x（n）的每一段分别与h（n）卷积，然后再将这些卷积结果相加起来就可得到输出序列，这样，每一段的卷积都可用上面讨论的快速卷积来计算。先对h（n）及xi（n）补零，补到具有N点长度，NN1N2-1。一般选择N2M，然后用基2FFT算法通过正反变换计算yi（n）xi（n）*h（n）由于yi（n）的长度为N，而xi（n）的长度为N2，因此相邻两yi（n）序列必然有N-N2N1-1点发生重叠，这个重叠部分应该相加起来才能构成最后的输出序列。计算步骤a.事先准备好滤波器参数H（k）DFT[h（n）]，N点b.用N点FFT计算Xi（k）DFT[xin]c.YikXikHkd.用N点IFFT求yi（n）IDFT[Yik]e.将重叠部分相加2重叠保存法这种方法和第一种方法稍有不同，即将上面分图序列中补零的部分不是补零，而是保留原来的输入序列值，且保留在各段的前端，这时，如利用DFT实现h（n）和xi（n）的循环卷积，则每段卷积结果的前N1-1个点不等于线性卷积值需舍去。为了清楚地看出这点，研究一下x（n）中一段长为N的序列xi（n）与h（n）（长为N1）的循环卷积情况由于h（n）的长度为N1，当0≤n≤N1-2时，hn-mN将在xi（m）的尾部出现有非零值，如图n1的情况就是如此，所以0≤n≤N1-2这部分yin值中将混入xim尾部与hn-mN的卷积值，从而使yin不同于线性卷积结果，但当nN1-1～N-1时，则有hn-mNh（n-m），因此从nN1-1点开始圆周圈卷积值完全与线性卷积值一样，yi（n）的后面N2点才是正确的卷积值，而每一段卷积运算结果的前N1-1点个值需去掉。为了不造成输出信号遗漏，对x（n）分段时，需使相邻两段有N1-1个点的重叠（对于第一段，xn由于没有前一段保留信号，在其前填补N1-1点个零点）。为此将xi（n）定义为每段和h（n）的循环卷积以yi（n）表示，,由FFT算出，去掉yi（n）的前N1-1点，再把相邻各段输出顺次连接起来就构成了最终的输出序列y（n）。重叠保留法每一输入段均由N-N11N2个新点和前一段保留下来的N1-1个点所组成NN1N2-1。值得注意的是，对于有限长时间序列xn（长度为LMN2），在结束段（iM-1）做完后，我们所得到的只是L点的线性卷积，还少了N1-1点，实际上就是h-n移出xn尾部时的不完全重合点，或者说是最后一段的重叠部分（N1-1个点，即每次都要删除的长度，但由于其位于序列尾部，故在其后增加N2个0以再行计算最后一次）少做了一次卷积，为此，应再补做这一段N1-1点，在其后填补N2点个零点保证长度仍为N点，一样舍去前取N1-1点，并从N1-1点开始，保留N2-1点。重叠保留法与重叠相加法的计算量差不多，但省去了重叠相加法最后的相加运算。一般来说，用FFT作信号滤波，只用于FIR滤波器阶数h（n）大于32的情况下，且取N2（5～10）N1，这样可接近于最高效的运算。Matlab实现function[y]ovrlpsavx,h,NOverlap-Saveofblockconvolution----------------------------------------[y]ovrlpsavx,h,NyoutputsequencexsequencehimpulseresponseNblocklengthLenxlengthx;Mlengthh;M1M-1;LN-M1;没有发生混叠的部分（实际的线性卷积结果）长h[hzeros1,N-M];x[zeros1,M1,x,zeros1,L];preappendM-1zerosKfloorLenxM-1/L;ofblocksYzerosK1,N;convolutionwithsuccesiveblocksfork0Kxkxk*L1k*LN;Yk1,circconvxk,h,N;circularconvolutionendYY,MN;discardthefirstM-1samplesyY;assembleoutput其中用到的函数function[y,n]circconvx1,x2,Nimplementcircleconvolutionofx1andx2x1andx2aresequencesofn1,n2N|lengthx2Nerrorlengthofsequencesisoverrange.;endx1[x1zeros1,N-lengthx1];x2[x2zeros1,N-lengthx2];m0N-1;x2x2mod-m,N1;HzerosN,N;forn11NHn,circshiftx2,n-1,N;endyx1*H;endfunction[y,n]circshiftx,m,NcircleshiftofsequencesUsage[y,n]circshiftx,m,NxsequencetobecircleshiftedmnumbertobeshifttorightNlengthofcircleshiftyoutputsequencesnoriginalsignalindexci