fft和卷积_fft算法

实验3利用FFT计算线性卷积 　　1、实验3利用FFT计算线性卷积、实验目的1. 掌握利用FFT计算线性卷积的原理及具体实现方法。2. 加深理解重叠相加法和重叠保留法。3. 考察利用FFT计算线性卷积各种方法的适用范围。二、实验设备与环境计算机、mat lab软件环境。三、实验基础理论1线型卷积和圆周卷积设x(n)为L点序列,h(n)为M点序列,x(n)和h(n)的线性卷积为00y/(n) = x(n) * y(n)=x(m)h(n th)m=co y】(n)的长度为L+MT。x(n)和h(n)的N点圆周卷积为00y(n) = x(n)®y(n)=x(n?)/i(n – m)N/?N(n)m =oo圆周卷积与线性卷积相等 　　2、而不产生交叠的必要条件为N > L-F M -1圆周卷积定理:根据DFT的性质，x(n)和h(n)的N点圆周卷积的DFT等于它们DFT的乘积DFTx(n)®h(n) = X(k)H(k)2. 快速卷积快速卷积算法用関周卷积实现线性卷积，根据圆周卷积定理利用 FFT算法实现圆周卷积。可以将快速卷积的步骤归纳如下：(1) 为了使线性卷积可以用圆周卷积来计算，必须选择 N>L + M-1;同时为了能使用基-2FFT完成卷积运算，要求N = 2vo 采用补零的办法使x(n)和h(n)的长度均为No(2) 计算x(n)和h(n)的N点FFTx(n)竺X(k)h(n)竺_H(k)(3 　　3、) 组成乘积Y(k) = X(k)H(k)(4) 利用IFFT计算Y(K)的IDFT,得到线性卷积y(n)Y(k)一丝_y(n)3. 分段卷积我们考察单位取样响应为h(n)的线性系统，输入为x(n),输出为 y(n),则y(n) = x(n) * h(n)当输入序列x(n)极氏时，如果要等x(n)全部集齐时再开始进行卷 积，会使输出相对输入有较大的延时，再者如果序列太长，需要大量 存贮单。为此我们把x(n)分段，分别求出每段的卷积，合在一起 得到最后总的输出。这种方法称为分段卷积分段卷积可细分为重叠保 留法和重叠相加法。重叠保留法：设x(n)的长度为N“ h(n)的长度为M。我们把序列x(n 　　4、) 分成多段N点序列焉(n),每段与前一段重叠M-1个样本。由于第一段 没有前一段保留信号，为了修正，我们在第一个输入端前面填充M-1 个零。计算每一段与h(n)的圆周卷积，则其每段卷积结果的前M-1 个样本不等于线性卷积值，不是正确的样本值。所以我们将每段卷积 结果的前M-1个样本社区，只保留后面的N-M+1个正确输出样本，把 这些输出样本合起来，得到总输出。利用FFT实现重叠保留发的步骤如下：(1) 在x(n)前面填充M-1个零，扩大以后的序列为x(n) = 0,0, .,0,x(n)M-1个(2) 将x(n)分为若干N点子段，设L二N-M+1为每一段的有效数据 长度，则第i段焉(n) ( 　　5、0 < n < N – 1)的数据为xt(n) = x(m) iL < m < iL N l,i > 0,0 < n < N 1(3) 计算每一段与h(n)的N点圆周卷积，利用FFT计算圆周卷积：兀(n)& (k)h(n)-H(k)Yi(k)=Xi(k)H(k)K(k) Z(n)(4) 舍去每一段卷积结果的前M-1个样本，连接剩下样本，得到卷 积结果y(n) o重叠相加法:设h(n)长度为M,将信号x(n)分解成长为L的子段，建议L选择与的M数量级相同，以治(n)表示没每段信号，则00x(n)=工石(町i=0Xi =x(n 4- iL), 0 　　6、<n < L 10,其它COx(n) * h(n)=Xi(n) * h(n)i=0每-段卷积为)的长度为L+M-1,所以在做求和时，相邻两段序列 有M-l个样本重叠，即前一段的最后M-1个样本喝下一段的前M-1个 序列重叠，这个重叠部分相加，再与不重叠部分共同组成输出y(n)。利用FFT实现重叠保留法的步骤如下：(1) 将x(n)分为若干L点子段片(n) o(2) 计算每一段与h(n)的卷积力)，根据快速卷积算法利用FFT 计算卷积。(3) 将各段力)相加，得到输出y(n)00y(n) = 2yI(n-/L)i=0四、实验内容假设要计算序列x ( n )二u ( n ) -u ( 　　7、nL ), 0<=n<=L和h(n) =cos (0. 2irn), 0<=n<=M的线性卷积，完成以下实验内容1. 设L二M,根据线性卷积的表达式和快速卷积的原理，分别编程实现计算两个序列线性卷积的方法，比较当序列长度分别为8,16,32,64,256,512,1024时，两种方法计算线性卷积所需时间。实验代码及结果: 程序Exp3_l代码L = input(!L=1)； M = input(* M= *)； xn = linspace(lf lf L); hn = cos(02*pi*O:M-1); ynl = zeros(lr L+M-l);ticfor k = 　　8、0:L+M-lfor n = max(0 k-M+1):min(k L-l)ynl(k+1) = ynl(k+1)+xn(n+1)*hn(k-n+1);endendtoe ticXk = fft(xnf L+M-l); Hk = fft(hnf L+M-l); Yk = Xk*Hk;yn2 = ifft(Yk); toe>> Exp3_lL=8M=8Elapsed time is 0.000136 secondsElapsed time is 0.000033 seconds>> Exp3_lL=16M=16Elapsed timeisis0.0002050.00003 　　9、0secondssecondsElapsedtime>> Exp3_1L=32M=32Elapsedtimeis0.000419secondsElapsed time is 0.000056 seconds.>> Exp3_lL=64M=64Elapsed timeElapsed timeis 0.000911is 0.000117secondsseconds>> Exp3_lL=256M=25 6Elapsed time is 0.005773 secondsElapsed time is 0.000261 seconds>> Exp3_lL= 　　10、512M=512Elapsed time is 0.017060 secondsElapsed time is 0.000150 seconds>> Exp3_lL=1024M=1024Elapsed time is 0.055688 secondsElapsed time is 0.000524 seconds由此可见，线性卷积慢得多2. 当L二2048且M二256时，比较直接计算线性卷积和快速卷积所需时间，进一步考察当L=4096且M二256时两种算法所需的时间。实验代码及结果:>> Exp3_lL=2048M=25 6Elapsed Lime is 0.03542 　　11、1 seconds.Elapsed time is 0.000448 seconds.由此可见，快速卷积快>> Exp3_lL=4096M=25 6Elapsed time is 0.069497 seconds.Elapsed time is 0.001463 seconds.由此可见，快速卷积快3. 编程实现利用重叠相加法计算两个序列的线性卷积，考察L二2048且M二256时计算线性卷积的时间,与第二题的结果进行比较。实验代码及结果: 程序Exp3_3代码L = input(1L=*)； M = input(* M= *)；xn = linspace(1,1, L);hn = 　　12、cos(02*pi*0:M-l);ticXk = fft(xnf L+M-l); Hk = fft(hnf L+M-l); Yk = Xk*Hk;yn2 = ifft(Yk); toe yn =;ticxln = linspace(0f 0r M-l) xn;Hk = fft(hn, 2*M);iend = floor(length(xln)/(M+l)+1；fori = 0: iendxin = xln (i* (M+l) +1 : min ( length (xln) i*(M+l)+2*M);Xik = fft(xinr 2*M);Yik = Xik*Hk;yin = ifft(Yik) 　　13、;yn = yn yin (M: 2*M);end yn = yn (1: L+M-l);toe>> Exp3 3L=2O48M=25 6Elapsed time is 0.000340 secondsElapsed time is 0.000513 seconds>> Exp3_3L=8152M=25 6Elapsed time is 0.003971 secondsElapsed time is 0.002105 seconds分析：当L = 2048 M=25 6时，重叠相加比较慢，原因是程序控 制流耗时较多。但足，随着序列L的增长，重叠相加变得较快4. 编程实现 　　14、利用重叠保留法计算两个序列的线性卷积，考察L二2048且M二256时计算线性卷积的时间,与第二题的结果进行比较。实验代码及结果: 程序Exp3_4代码L = input(* L=,);M = input()； xn = linspace(lf lr L);hn = cos(02*pi*O:M-1);ticXk = fft(xnf L+M-l);Hk = fft(hnr L+M-l);Yk = Xk*Hk;yn2 = ifft(Yk);toe %d±OEai»y ticHk = fft(hn, 2*M-1);iend = floor(length(xn)/M);i = 0;x 　　15、in = xn(i*M+l : min(i+1)length(xn);Xik = fft(xinr 2*M-1);Yik = Xik*Hk;yin = ifft(Yik);yn = yin;%pUO» IjEEafori = 1: iendxin = xn(i*M+l : min(i+1)length(xn);Xik = fft(xin, 2*M-1);Yik = Xik*Hk;yin = ifft(Yik);yn = yn(1:(length(yn)-M+l)yn(length(yn)-(M-2):length(yn) + yin(1:M-1) yin(M:2*M-l) ; %5：x6DdADendyn = yn(1:L+M-l);toe >> Exp3_4L=2O48M=25 6Elapsed time is 0.000441 seconds Elapsed time is 0.001055 seconds >> Exp3_4L=32768M=1024Elapsed time is 0.018384 seconds.Elapsed time is 0.016179 second