实验三卷积FFT

实验三卷积、FFT频谱分析方法1、实验目的（1）进一步加深对线性卷积的理解和分析能力；（2）通过编程，上机调试程序，进一步增强使用计算机解决问题的能力；（3）掌握线性卷积与循环卷积软件实现的方法，并验证二者之间的关系。(4)学习用FFT对连续信号和时域离散信号进行谱分析的方法，了解可能出现的分析误差及其原因，以便正确应用FFT。2、实验原理与方法（1）线性卷积线性时不变系统（LinearTime-InvariantSystem,orL.T.I系统）输入、输出间的关系为：当系统输入序列为)(nx，系统的单位脉冲响应为)(nh，输出序列为)(ny，则系统输出为：mnhnxmnhmxny)(*)()()()(或mnxnhmnxmhny)(*)()()()(上式称为离散卷积或线性卷积。图4.1示出线性时不变系统的输入、输出关系。)(n→L.T.I—→)(nh—→—→图4.1线性时不变系统的输入、输出关系(2)圆周卷积设两个有限长序列)(1nx和)(2nx，均为N点长)(1nx)(1kX)(2nx)(2kX如果)()()(213kXkXkX则)()(~)(~)(10213nRmnxmxnxNNm1021)()(NmNmnxmx)(nxL.T.Ih(n)mmnhmxny)()()(DFTDFT3、实验内容及要求1)卷积计算已知两个有限长序列)4(5)3(4)2(3)1(2)()(nnnnnnx编制一个计算两个序列线性卷积的通用程序，计算)(*)(nhnx。编制一个计算循环卷积的通用程序，计算两个序列)(nx与)(nh的圆周卷积。上机调试并打印或记录实验结果。2)FFT变换(1)对以下序列进行谱分析。其它nnnnnnx其它nnnnnnxnRnx,074,330,4)(,074,830,1)()()(3241选择FFT的变换区间N为8和16两种情况进行频谱分析。分别打印其幅频特性曲线。并进行对比、分析和讨论。(2)对模拟周期信号进行谱分析6()cos8cos16cos20xtttt选择采样频率HzFs64，变换区间N=16,32,64三种情况进行谱分析。分别打印其幅频特性，并进行分析和讨论。1）卷积计算（1）定义圆周卷积函数，并保存为circonv.m文件：functionyc=circonv(x1,x2,N)iflength(x1)Nerror('x1mustbelessthanlengthofN');endiflength(x2)Nerror('x2mustbelessthanlengthofN');endx1=[x1,zeros(1,N-length(x1))];x2=[x2,zeros(1,N-length(x2))];n=[0:N-1];)3(2)2()1(2)()(nnnnNHx2=x2(mod(-n,N)+1);H=zeros(N,N);forn=1:1:NH(n,:)=cirshiftd(x2,n-1,N);endyc=x1*H';（2）定义循环移位函数，并保存为cirshiftd.m文件：functiony=cirshiftd(x,m,N)iflength(x)Nerror('lengthofxmustbelessthanN');endx=[x,zeros(1,N-length(x))];n=[0:1:N-1];y=x(mod(n-m,N)+1);（3）调用相应的函数：clearall;n=[0:1:11];m=[0:1:5];N1=length(n);N2=length(m);xn=0.8.^n;hn=ones(1,N2);y1n=conv(xn,hn);ycn=circonv(xn,hn,N1);ny1=[0:1:length(y1n)-1];ny2=[0:1:length(ycn)-1];subplot(2,1,1);stem(ny1,y1n);ylabel('线性卷积');subplot(2,1,2);stem(ny2,ycn);ylabel('圆周卷积');程序运行结果如图所示：2）FFT变换(1)对序列进行谱分析%用FFT对信号作频谱分析clearallcloseallx1n=[ones(1,4)];M=8;xa=1:(M/2);xb=(M/2):-1:1;x2n=[xa,xb];x3n=[xb,xa];X1k8=fft(x1n,8);X1k16=fft(x1n,16);X2k8=fft(x2n,8);X2k16=fft(x2n,16);X3k8=fft(x3n,8);X3k16=fft(x3n,16);N=8;f=2/N*(0:N-1);subplot(2,2,1);stem(f,abs(X1k8),'.');title('(la)8点DFT[X_1(n)]');xlabel('w/pi');ylabel('幅度');axis([0,2,0,1.2*max(abs(X1k8))])N=16;f=2/N*(0:N-1);subplot(2,2,2);stem(f,abs(X1k16),'.');title('(1b)16点DFT[X_1(n)]');xlabel('w/pi');ylabel('幅度');axis([0,2,0,1.2*max(abs(X1k16))])figure(2)N=8;f=2/N*(0:N-1);subplot(2,2,1);stem(f,abs(X2k8),'.');title('(2a)8点DFT[X_2(n)]');xlabel('w/pi');ylabel('幅度');axis([0,2,0,1.2*max(abs(X2k8))])N=16;f=2/N*(0:N-1);subplot(2,2,2);stem(f,abs(X2k16),'.');title('(2b)16点DFT[X_2(n)]');xlabel('w/pi');ylabel('幅度');axis([0,2,0,1.2*max(abs(X2k16))])N=8;f=2/N*(0:N-1);subplot(2,2,3);stem(f,abs(X3k8),'.');title('(3a)8点DFT[X_3(n)]');xlabel('w/pi');ylabel('幅度');axis([0,2,0,1.2*max(abs(X3k8))])N=16;f=2/N*(0:N-1);subplot(2,2,4);stem(f,abs(X3k16),'.');title('(3b)16点DFT[X_3(n)]');xlabel('w/pi');ylabel('幅度');axis([0,2,0,1.2*max(abs(X3k16))])实验运行结果：实验结果分析：图（1a）和（1b）说明x1(n)=R4(n)的8点DFT和16点DFT分别是x1(n)的频谱函数的8点和16点采样。因为x3(n)=x2((n+3))8R8(n),所以，x3(n)和x2(n)的8点DFT的模相等，如图（2a）和（3a）所示。但是，当N=16时，x3(n)和x2(n)不满足循环移位关系，所以x3(n)和x2(n)模不同。（2）对模拟周期信号进行谱分析%模拟周期信号谱分析clearall;closeall;figure(4)Fs=64;T=1/Fs;N=16;n=0:N-1;x6nT=cos(8*pi*n*T)+cos(16*pi*n*T)+cos(20*pi*n*T);X6k16=fft(x6nT);X6k16=fftshift(X6k16);Tp=N*T;F=1/Tp;k=-N/2:N/2-1;fk=k*F;subplot(3,1,1);stem(fk,abs(X6k16),'.');boxontitle('(4a)16点|DFT[x_6(nT)]|');xlabel('f/HZ');ylabel('幅度');axis([-N*F/2-1,N*F/2-1,0,1.2*max(abs(X6k16))])N=32;n=0:N-1;x6nT=cos(8*pi*n*T)+cos(16*pi*n*T)+cos(20*pi*n*T);X6k32=fft(x6nT);X6k32=fftshift(X6k32);Tp=N*T;F=1/Tp;k=-N/2:N/2-1;fk=k*F;subplot(3,1,2);stem(fk,abs(X6k32),'.');boxontitle('(4b)32点|DFT[x_6(nT)]|');xlabel('f/HZ');ylabel('幅度');axis([-N*F/2-1,N*F/2-1,0,1.2*max(abs(X6k32))])N=64;n=0:N-1;x6nT=cos(8*pi*n*T)+cos(16*pi*n*T)+cos(20*pi*n*T);X6k64=fft(x6nT);X6k64=fftshift(X6k64);Tp=N*T;F=1/Tp;k=-N/2:N/2-1;fk=k*F;subplot(3,1,3);stem(fk,abs(X6k64),'.');boxontitle('(4c)64点|DFT[x_6(nT)]|');xlabel('f/HZ');ylabel('幅度');axis([-N*F/2-1,N*F/2-1,0,1.2*max(abs(X6k64))])实验运行结果：实验结果分析：X6(t)有三个频率成分，f1=4Hz，f2=8Hz，f3=10Hz.所以X6(t)的周期为0.5s。采样频率Fs=64Hz=16f1=8f2=6.4f3.变换区间N=16时，观察时间Tp=16T=0.25s,不是X6(t)的整数倍周期，所以所得频谱不正确，如图(4a)所示。变换区间N=32，64时，观察Tp=0.5,1s，是X6(t)的整数倍周期，所以所得频谱正确，如图（4b）和(4c)所示。图中三根谱线正好分别位于4、8、10Hz处。变换区间N=64时频谱幅度是变换区间N=32时的2倍，这种结果正好验证了用DFT对中期序列谱分析的理论。