数字信号处理实验

数字信号处理实验电信1班3120503005詹筱珊1-11给出如下差分方程：nnxnnyny);()2(9.0)1()(（1）计算并画出单位脉冲响应12020);(，，nnh。（2）由此)(nh规定的系统是否稳定？解：MATLAB实现程序：a=[1,-1,0.9]；b=[1]；x=impseq(0,-20,120)；%输入x为单位脉冲序列n=[-20:120]；h=filter(b,a,x)；%系统输出为单位脉冲响应stem(n,h,'.')；axis（[-20,120,-1.1,1.1]）title('脉冲响应')；text(125,-1.1,'n')；ylabel('h(n)')sum(abs(h))%对单位脉冲响应的模值求和程序运行结果：ans=14.8785系统是稳定的。1-12已知某模拟信号taetx1000)(，将它分别用不同的采样频率进行采样得到离散时间信号，试分析在以下两种采样情况下对信号频率的影响。（1）采样频率Hzfs5000；（2）采样频率Hzfs1000解：（1）MATLAB实现程序：dt=0.00005;t=-0.005:dt:0.005;xa=exp(-1000*abs(t));fs=5000,Ts=1/fs;n=-25:1:25;x=exp(-1000*abs(n*Ts));N=500;k=0:1:N;w=pi*k/N;X=x*exp(-j*n'*w);X=real(X);W=[-fliplr(w),w(2:N+1)];X=[fliplr(X),X(2:N+1)];subplot(2,1,1)；plot(t*1000,xa)；xlabel('时间：ms');ylabel('x1(n)')title('离散信号')；holdonstem(n*Ts*1000,x);holdoffsubplot(2,1,2)；plot(w/pi,X);xlabel('以pi为单位的频率')；ylabel('X1(w)');title('离散信号傅里叶变换')运行结果：（2）MATLAB实现程序：dt=0.00005;t=-0.005:dt:0.005;xa=exp(-1000*abs(t));fs=1000,Ts=1/fs;n=-25:1:25;x=exp(-1000*abs(n*Ts));N=500;k=0:1:N;w=pi*k/N;X=x*exp(-j*n'*w);X=real(X);W=[-fliplr(w),w(2:N+1)];X=[fliplr(X),X(2:N+1)];subplot(2,1,1)；plot(t*1000,xa)；xlabel('时间：ms');ylabel('x1(n)')title('离散信号')；holdonstem(n*Ts*1000,x);holdoffsubplot(2,1,2)；plot(w/pi,X);xlabel('以pi为单位的频率')；ylabel('X1(w)');title('离散信号傅里叶变换')运行结果：2-25已知某系统的系统函数为432143213.08.05.12.112.01.03.01.03.0)(zzzzzzzzzH求其零、极点并绘出零、极点图。解：MATLAB实现程序：b=[0.30.10.30.10.2];a=[1-1.21.5-0.80.3];r1=roots(a)r2=roots(b)zplane(b,a)运行结果：r1=0.1976+0.8796i0.1976-0.8796i0.4024+0.4552i0.4024-0.4552ir2=0.3236+0.8660i0.3236-0.8660i-0.4903+0.7345i-0.4903-0.7345i2-26已知因果系统)()1(9.0)(nxnyny，绘出)(jeH的幅频和相频特性曲线。解：MATLAB实现程序：b=[1,0];a=[1,-0.9][H,w]=freqz(b,a,100,'whole');magH=abs(H);phaH=angle(H);subplot(2,1,1),plot(w/pi,magH);gridxlabel('');ylabel('幅度');title('幅频响应')subplot(2,1,2);plot(w/pi,phaH/pi);gridxlabel('频率(单位：pi)');ylabel('相位（单位:pi）');title('相频响应')运行结果：2-27一个线性时不变系统，描绘它的差分方程为)2()1(2)()2(25.0)1(5.0)(nxnxnxnynyny（1）在1000n之间求得并画出系统的脉冲响应，从脉冲响应确定系统的稳定性；（2）画出该系统的幅频、相频特性；（3）如果此系统的输入为)()]3.0sin(4)2.0cos(35[)(1nunnnx，在2000n间求系统的输出)(ny；（4）讨论当输入改为)()]8.0sin(4)2.0cos(35[)(2nunnnx时，输出波形如何变化？为什么？试根据系统的幅频特性解释。解：MATLAB实现程序：b=[1,2,1];a=[1,-0.5,0.25];x=impseq(0,0,100);n=0:100;h=filter(b,a,x);stem(n,h,'.');title('脉冲响应波形');sum(abs(h))figure;[H,w]=freqz(b,a,200);magH=abs(H);phaH=angle(H);subplot(2,1,1),plot(w/pi,magH);grid;xlabel('');ylabel('幅度')；title('幅频响应')subplot(2,1,2);plot(w/pi,phaH/pi);grid;xlabel('频率(单位：pi)');ylabel('相位（单位：pi）');title('相频响应');figure;n1=0:200;x1=5+3*cos(0.2*pi*n1)+4*sin(0.3*pi*n1);y1=filter(b,a,x1);subplot(2,1,1);plot(n1,x1);title('输入x1波形');subplot(2,1,2);plot(n1,y1);title('输出y1波形')；figure;n1=0:200;x2=5+3*cos(0.2*pi*n1)+4*sin(0.8*pi*n1);y2=filter(b,a,x2);subplot(2,1,1);plot(n1,x2);title('输入x2波形');subplot(2,1,2);plot(n1,y2);title('输出y2波形')运行结果：ans=6.57143-4)(nx是一4点序列：,030,1)(nnx（1）计算离散时间傅里叶变换（DTFT）即jeX,并画出它的幅度和相位。（2）计算)(nx的4点DFT。解：（1）MATLAB实现程序：x=[1,1,1,1];w=[0:1:500]*2*pi/500;[H]=freqz(x,1,w);magH=abs(H);phaH=angle(H);subplot(2,1,1);plot(w/pi,magH);gridxlabel('');ylabel('|X|');title('DTFT的幅度')subplot(2,1,2);plot(w/pi,phaH/pi*180);gridxlabel('以pi为单位的频率');ylabel('度');title('DTFT的相角')运行结果：（2）MATLAB实现程序：N=4;k=0:N-1;X=dft(x,N);magX=abs(X),phaX=angle(X)*180/pisubplot(2,1,1);plot(w*N/(2*pi),magH,'--');axis([-0.1,4.1,0,5]);holdonstem(k,magX);ylabel('|X(k)|');title('DFT的幅度：N=4')；text(4.3,-1,'k')holdoffsubplot(2,1,2);plot(w*N/(2*pi),phaH*180/pi,'--');axis([-0.1,4.1,-200,200]);holdonstem(k,phaX);ylabel('度')；title('DFT相角：N=4')；text(4.3,-200,'k')运行结果：3-12)]52.0cos()48.0[cos()(nnnx，利用MATLAB程序求如下jeX、kX。（1）取)(nx的前10点数据，求N=10点的jeX、kX并作图；（2）将)(nx补零至100点，求N=100点的jeX、kX并作图。解：（1）MATLAB实现程序：n=[0:1:9];x=cos(0.48*pi*n)+cos(0.52*pi*n);w=[0:1:500]*2*pi/500;X=x*exp(-j*n'*w);magX=abs(X);x1=fft(x);magX1=abs(x1(1:1:10));k1=0:1:9;w1=2*pi/10*k1;subplot(3,1,1);stem(n,x);title('x(n),0=n=9');xlabel('n')axis([0,10,-2.5,2.5]);line([0,10],[0,0]);subplot(3,1,2);plot(w/pi,magX);title('DFT幅度');xlabel('频率（单位：pi）');axis([0,1,0,10])subplot(3,1,3);stem(w1/pi,magX1);title('DFT幅度')；xlabel('频率（单位：pi）');axis([0,1,0,10])运行结果：（2）MATLAB实现程序：n=[0:1:9];y=cos(0.48*pi*n)+cos(0.52*pi*n);n1=[0:1:99];x=[y(1:1:10)zeros(1,90)];w=[0:1:500]*2*pi/500;X=x*exp(-j*n1'*w);magx=abs(X);x1=fft(x);magx1=abs(x1(1:1:50));k1=0:1:49;w1=2*pi/100*k1;subplot(3,1,1);stem(n,y);title('x(n),0=n=90zeros');xlabel('n')axis([0,100,-2.5,2.5]);line([0,100],[0,0]);subplot(3,1,2);plot(w/pi,magx);title('DTFT幅度')；xlabel('频率（单位：pi）');axis([0,1,0,10])subplot(3,1,3);stem(w1/pi,magx1);title('DFT幅度');xlabel('频率（单位：pi）');axis([0,1,0,10])运行结果：5-7利用脉冲响应不变法设计一个数字巴特沃斯低通滤波器，滤波器的技术要求为：通带截止频率600srad/s，阻带最小衰减dBs12。研究不同采样频率对所设计数字滤波器频率响应的影响，设采样频率sf分别取1kHz,2kHz,4kHz。解：（1）fs=1000HzMATLAB实现程序：wp=200*pi;ws=600*pi;Rp=3;Rs=12;[n,wn]=buttord(wp,ws,Rp,Rs,'s')[b,a]=butter(n,wn,'s')[db,mag,pha,w]=freqs_m(b,a,500*2*pi);plot(w/(2*pi),db);axis([0,500,-20,1]);holdonfs=1000;[bz,az]=impinvar(b,a,fs);[db,mag,pha,grd,w]=freqz_m(bz,az);plot(0.5*fs*w/pi,db);axis([0,500,-20,1]);h