数字信号处理实验吴镇扬版matlab程序.docx

1、数字信号处理实验吴镇扬版matlab程序（1）数组的加、减、乘、除和乘方运算。输入A=1 2 3 4，B=3 4 5 6，求C=A+B，D=A-B，E=A.*B，F=A./B，G=A.B并用stem语句画出A、B、C、D、E、F、G。clear all;a=1 2 3 4; b=3 4 5 6;c=a+b;d=a-b;e=a.*b;f=a./b;g=a.b;n=1:4;subplot(4,2,1);stem(n,a);xlabel(n);xlim(0 5);ylabel(A);subplot(4,2,2);stem(n,b);xlabel(n);xlim(0 5);ylabel(B);subp

2、lot(4,2,3);stem(n,c);xlabel(n);xlim(0 5);ylabel(C);subplot(4,2,4);stem(n,d);xlabel(n);xlim(0 5);ylabel(D);subplot(4,2,5);stem(n,e);xlabel(n);xlim(0 5);ylabel(E);subplot(4,2,6);stem(n,f);xlabel(n);xlim(0 5);ylabel(F);subplot(4,2,7);stem(n,g);xlabel(n);xlim(0 5);ylabel(G);（2）用MATLAB实现下列序列：a) x(n)=0.8n

3、 0n15b) x(n)=e(0.2+3j)n 0n15c) x(n)=3cos(0.125n+0.2)+2sin(0.25n+0.1) 0n15d) 将c)中的x(n)扩展为以16为周期的函数x16(n)=x(n+16),绘出四个周期。e) 将c)中的x(n)扩展为以10为周期的函数x10(n)=x(n+10),绘出四个周期。clear all;N=0:15;xa=0.8.N;figure;subplot(2,1,1);stem(N,xa); xlabel(n);xlim(0 16);ylabel(xa);xb=exp(0.2+3*j)*N);subplot(2,1,2);stem(N,xb

4、);xlabel(n);xlim(0 16);ylabel(xb);figure;xc=3*cos(0.125*pi*N+0.2*pi)+2*sin(0.25*pi*N+0.1*pi);subplot(3,1,1);stem(N,xc);xlabel(n);xlim(0 16);ylabel(xc);k=0:3;m=0;for i=1:4for j=1:16m=m+1;n(m)=N(j)+16*k(i);x16(m)=3*cos(0.125*pi*n(m)+0.2*pi)+2*sin(0.25*pi*n(m)+0.1*pi);endend subplot(3,1,2);stem(n,x16);

5、xlabel(n);ylabel(x16);for j=1:10 x10(j)=x16(j);endfor i=1:3for m=1:10 x10(i*10+m)=x10(m);end endn=1:40;subplot(3,1,3);stem(n,x10); xlabel(n);ylabel(x10);（3）x(n)=1,-1,3,5,产生并绘出下列序列的样本:a) x1(n)=2x(n+2)-x(n-1)-2x(n)b) clear alln=1:4;T=4;x=1 -1 3 5;x(5:8)=x(1:4);subplot(2,1,1);stem(1:8,x);grid;for i=1:4

6、if i-10x1(i)=2*x(i+2)-x(i-1)-2*x(i);elsex1(i)=2*x(i+2)-x(i-1+T)-2*x(i);endendx1(5:8)=x1(1:4);subplot(2,1,2);stem(1:8,x1);grid; （4）绘出下列时间函数的图形，对x轴、y轴以及图形上方均须加上适当的标注：a) x(t)=sin(2t) 0t10sb) x(t)=cos(100t)sin(t) 0t4sta=0:0.05:10;xa=sin(2*pi*ta);subplot(2,1,1);plot(ta,xa);xlabel(t);ylabel(幅度);tb=0:0.01:

7、4;xb=cos(100*pi*tb).*sin(pi*tb);subplot(2,1,2);plot(tb,xb);xlabel(t);ylabel(幅度); （5）编写函数stepshift(n0,n1,n2)实现u(n-n0),n1n0=0;stem(n,x);（6）给一定因果系统求出并绘制H(z)的幅频响应与相频响应。clear all;b=1,sqrt(2),1;a=1,-0.67,0.9;h,w=freqz(b,a);am=20*log10(abs(h);subplot(2,1,1);plot(w,am);ph=angle(h);subplot(2,1,2);plot(w,ph);

8、（7）计算序列8 -2 -1 2 3和序列2 3 -1 -3的离散卷积，并作图表示卷积结果。clear all;a=8 -2 -1 2 3;b=2 3 -1 -3;c=conv(a,b); %计算卷积M=length(c)-1;n=0:1:M;stem(n,c);xlabel(n);ylabel(幅度); （8）求以下差分方程所描述系统的单位脉冲响应h(n),0n50y(n)+0.1y(n-1)-0.06y(n-2)=x(n)-2x(n-1)clear all;N=50;a=1 -2;b=1 0.1 -0.06;x=1 zeros(1,N-1);k=0:1:N-1;y=filter(a,b,x

9、);stem(k,y);xlabel(n);ylabel(幅度 ); (1)、观察高斯序列的时域和幅频特性，固定信号xa(n)中参数p=8，改变q的值，使q分别等于2，4，8，观察它们的时域和幅频特性，了解当q取不同值时，对信号序列的时域幅频特性的影响；固定q=8，改变p，使p分别等于8，13，14，观察参数p变化对信号序列的时域及幅频特性的影响，观察p等于多少时，会发生明显的泄漏现象，混叠是否也随之出现？记录实验中观察到的现象，绘出相应的时域序列和幅频特性曲线。解：程序见附录程序一：n=0:1:15;%p=8不变，q变化（2,4,8）;p=8;q=2; %p=8;q=2; xa1=exp(-

10、(n-p).2)/q);subplot(5,2,1);plot(n,xa1,-*);xlabel(t/T);ylabel(xa(n);title(p=8 q=2)xk1=abs(fft(xa1);subplot(5,2,2);stem(n,xk1)xlabel(k);ylabel(Xa(k);title(p=8 q=2)p=8;q=4; %p=8;q=4; xa1=exp(-(n-p).2)/q);subplot(5,2,3);plot(n,xa1,-*);xlabel(t/T);ylabel(xa(n);title(p=8 q=4)xk1=abs(fft(xa1);subplot(5,2,4

11、);stem(n,xk1)xlabel(k);ylabel(Xa(k);title(p=8 q=4)p=8;q=8; %p=8;q=8; xa1=exp(-(n-p).2)/q);subplot(5,2,5);plot(n,xa1,-*);xlabel(t/T);ylabel(xa(n);xk1=abs(fft(xa1);title(p=8 q=8)subplot(5,2,6);stem(n,xk1)xlabel(k);ylabel(Xa(k);title(p=8 q=8)%q=8不变,p变化(8,13,14);p=8;q=8; %p=8;q=8; xa1=exp(-(n-p).2)/q);s

12、ubplot(5,2,5);plot(n,xa1,-*);xlabel(t/T);ylabel(xa(n);xk1=abs(fft(xa1);title(p=8 q=8)subplot(5,2,6);stem(n,xk1)xlabel(k);ylabel(Xa(k);title(p=8 q=8)p=13;q=8; %p=13;q=8; xa1=exp(-(n-p).2)/q);subplot(5,2,7);plot(n,xa1,-*);xlabel(t/T);ylabel(xa(n);xk1=abs(fft(xa1);title(p=13 q=8)subplot(5,2,8);stem(n,x

13、k1)xlabel(k);ylabel(Xa(k);title(p=13 q=8)p=14;q=8; %p=14;q=8; xa1=exp(-(n-p).2)/q);subplot(5,2,9);plot(n,xa1,-*);xlabel(t/T);ylabel(xa(n);title(p=14 q=8)xk1=abs(fft(xa1);subplot(5,2,10);stem(n,xk1)xlabel(k);ylabel(Xa(k);title(p=14 q=8)分析：由高斯序列表达式知n=p为期对称轴；当p取固定值时，时域图都关于n=8对称截取长度为周期的整数倍，没有发生明显的泄漏现象；但

14、存在混叠，当q由2增加至8过程中，时域图形变化越来越平缓，中间包络越来越大，可能函数周期开始增加，频率降低，渐渐小于fs/2,混叠减弱； 当q值固定不变，p变化时，时域对称中轴右移，截取的时域长度渐渐地不再是周期的整数倍，开始无法代表一个周期，泄漏现象也来越明显，因而图形越来越偏离真实值，p=14时的泄漏现象最为明显，混叠可能也随之出现； (2)、观察衰减正弦序列xb(n)的时域和幅频特性，a=0.1，f=0.0625，检查谱峰出现位置是否正确，注意频谱的形状，绘出幅频特性曲线，改变f，使f分别等于0.4375和0.5625，观察这两种情况下，频谱的形状和谱峰出现位置，有无混叠和泄漏现象？说明

15、产生现象的原因。n1=0:1:15;xb1=exp(-0.1*n1).*sin(2*pi*0.0625*n1);subplot(3,2,1);plot(n1,xb1,-*);xlabel(n);ylabel(x(n);title(f=0.0625);xk1=abs(fft(xb1);subplot(3,2,2);stem(n1,xk1)xlabel(k);ylabel(X(k);title(f=0.0625);n2=0:1:15;xb2=exp(-0.1*n2).*sin(2*pi*0.4375*n2);subplot(3,2,3);plot(n2,xb2,-*);xlabel(n);ylab

16、el(x(n);title(f=0.4375);xk2=abs(fft(xb2);subplot(3,2,4);stem(n2,xk2)xlabel(k);ylabel(X(k);title(f=0.4375);n3=0:1:15;xb3=exp(-0.1*n3).*sin(2*pi*0.5625*n3);subplot(3,2,5);plot(n3,xb3,-*);xlabel(n);ylabel(x(n);title(f=0.5625);xk3=abs(fft(xb3);subplot(3,2,6);stem(n3,xk3)xlabel(k);ylabel(X(k);title(f=0.5

17、625);分析：当f=f1=0.0625时，谱峰位置出现正确，存在在混叠现象，时域采样为一周期，不满足采样定理。当f=0.4375和0.5625时，时域图像关于Y轴对称，频域完全相同。这是因为频域图是取绝对值的结果，所以完全相同。另外由于时域采样为6个半周期，满足采样定理，无混叠；但由于截取长度不是周期整数倍，出现泄漏。(3)、观察三角波和反三角波序列的时域和幅频特性，用N=8点FFT分析信号序列xc(n)和xd(n)的幅频特性，观察两者的序列形状和频谱曲线有什么异同？绘出两序列及其幅频特性曲线。在xc(n)和xd(n)末尾补零，用N=16点FFT分析这两个信号的幅频特性，观察幅频特性发生了什

18、么变化？两情况的FFT频谱还有相同之处吗？这些变化说明了什么？三角波序列:反三角波序列：clear all;n=0:3;k=1:8;%定义三角波序列Xc(n+1) = n;Xc(n+5) =4-n;%定义反三角波序列Xd(n+1) = 4-n;Xd(n+5) =n;% 三角波特性 %subplot(2,2,1);plot(k-1,Xc);xlabel(n);ylabel(时域特性);text(1,3,三角波);subplot(2,2,2);plot(k-1,abs(fft(Xc);xlabel(k);ylabel(幅频特性);text(4,10,三角波);% 反三角波特性 %subplot(2

19、,2,3);plot(k-1,Xd);xlabel(n);ylabel(时域特性);text(3,3,反三角波);subplot(2,2,4);plot(k-1,abs(fft(Xd);xlabel(k);ylabel(幅频特性);text(4,10,反三角波);%末尾补0，计算32点FFTXc(9:32)=0;Xd(9:32)=0;k=1:32;figure;% 三角波特性 %subplot(2,2,1);plot(k-1,Xc);xlabel(n);ylabel(时域特性);text(1,3,三角波);subplot(2,2,2);plot(k-1,abs(fft(Xc);xlabel(k

20、);ylabel(幅频特性);text(4,10,三角波);% fan三角波特性 %subplot(2,2,3);plot(k-1,Xd);xlabel(n);ylabel(时域特性);text(3,3,反三角波);subplot(2,2,4);plot(k-1,abs(fft(Xd);xlabel(k);ylabel(幅频特性);text(4,10,反三角波);N=8时域和幅度频谱图：分析： 由图知，三角波序列和反三角波序列的时域图像成镜像关系，但频域图像完全一样，只是因为幅频图是对x（k）的值取绝对值。N=32时域和幅度频谱图：分析：由实验所得的图形知，N=32点时和的幅频特性都更加密集，

21、更多离散点的幅值显示，“栅栏效应”减小，分辨率提高，而对于来说变化更加明显。在原序列的末端填补零值，变动了DFT的点数，人为的改变了对真实频谱采样的点数和位置，相当于搬动了“尖桩栅栏”的位置，从而使得频谱的峰点和谷点暴露出来。N=32时，和的频谱差别较大，但总体趋势仍然都是中间最小，两侧呈对称。 (4)、一个连续信号含两个频率分量，经采样得x(n)=sin2*0.125n+cos2*(0.125+f)n n=0,1,N-1已知N=16，f分别为1/16和1/64，观察其频谱；当N=128时，f不变，其结果有何不同，为什么？clear all;% N = 16 %N=16;detf=1/16;n

22、=0:N-1;x1(n+1)=sin(2*pi*0.125.*n)+cos(2*pi*(0.125+detf).*n);detf = 1/64;x2(n+1)=sin(2*pi*0.125.*n)+cos(2*pi*(0.125+detf).*n);% N = 16,detf = 1/16 %subplot(2,2,1);stem(n,x1);hold;plot(n,x1);xlabel(n);ylabel(时域特性);text(6,1,N=16,detf=1/16);subplot(2,2,2);stem(n,abs(fft(x1);xlabel(n);ylabel(幅值特性);text(6

23、,4,N=16,detf=1/16);% N = 16,detf = 1/64 %subplot(2,2,3);stem(n,x2);xlabel(n);ylabel(时域特性);text(6,1,N=16,detf=1/64);subplot(2,2,4);stem(n,abs(fft(x2);xlabel(n);ylabel(幅值特性);text(6,4,N=16,detf=1/64);% N = 128 %N=128;detf=1/16;n=0:N-1;x3(n+1)=sin(2*pi*0.125.*n)+cos(2*pi*(0.125+detf).*n);detf = 1/64;x4(

24、n+1)=sin(2*pi*0.125.*n)+cos(2*pi*(0.125+detf).*n);% N = 128,detf = 1/16 %figure;subplot(2,2,1);stem(n,x3);xlabel(n);ylabel(时域特性);axis(0 128 -2 2);text(6,1.5,N=128,detf=1/16);subplot(2,2,2);stem(n,abs(fft(x3);xlabel(n);ylabel(幅值特性);axis(0 128 -10 70);text(40,60,N=128,detf=1/16);% N = 128,detf = 1/64

25、%subplot(2,2,3);stem(n,x3);xlabel(n);ylabel(时域特性);axis(0 128 -2 2);text(6,1.5,N=128,detf=1/16);subplot(2,2,4);stem(n,abs(fft(x4);xlabel(n);ylabel(幅值特性);axis(0 128 -10 70);text(40,60,N=128,detf=1/16); (5)、用FFT分别实现xa(n)（p8，q2）和 xb(n)（a0.1，f0.0625）的16点圆周卷积和线性卷积。clear all;N=16;n=0:N-1;p=8;q=2;Xa(n+1)=ex

26、p(-(n-p).2./q);a=0.1;f=0.0625;Xb(n+1)=exp(-a.*n).*sin(2*pi*f.*n);%16点循环卷积Fa=fft(Xa); Fb=fft(Xb);Fx=Fa.*Fb;X51=ifft(Fx);stem(n,X51);%16点线性卷积Xa(N+1:2*N-1)=0;Xb(N+1:2*N-1)=0;Fa=fft(Xa); Fb=fft(Xb);Fc=Fa.*Fb;X52=ifft(Fc);figure;stem(1:2*N-1,X52);（7）用FFT分别计算xa(n)(p=8,q=2)和xb(n)(a=0.1,f=0.0625)的16点循环相关和线性

27、相关，问一共有多少种结果，他们之间有何异同点。clear all;N=16;n=0:N-1;p=8;q=2;Xa(n+1)=exp(-(n-p).2./q);a=0.1;f=0.0625;Xb(n+1)=exp(-a.*n).*sin(2*pi*f.*n); N=length(Xa); %16点循环相关Fa=fft(Xa,2*N); Fb=fft(Xb,2*N);Fx=conj(Fa).*Fb;X71=real(ifft(Fx);X71=X71(N+2:2*N) X71(1:N);n=(-N+1):(N-1); stem(n,X71);%16点线性相关Xa(N+1:2*N-1)=0;Xb(N+

28、1:2*N-1)=0;Fa=fft(Xa); Fb=fft(Xb);Fc=conj(Fa).*Fb;X72=real(ifft(Fc);figure;stem(1:2*N-1,X72);（8）用FFT分别计算xa(n)(p=8,q=2)和xb(n)(a=0.1,f=0.0625)的自相关函数。clear all;N=16;n=0:N-1;p=8;q=2;Xa(n+1)=exp(-(n-p).2./q);a=0.1;f=0.0625;Xb(n+1)=exp(-a.*n).*sin(2*pi*f.*n); %自然对数的底:e=:2.71828 18284 59045 23536 N=length(Xa); % Xa(n) 16点自相关Fa=fft(Xa,2*N); Fb=fft(X