数字信号处理实验.doc

《数字信号处理实验.doc》由会员分享，可在线阅读，更多相关《数字信号处理实验.doc（20页珍藏版）》请在装配图网上搜索。

数字信号处理实验 电信1班 3120503005 詹筱珊 1-11给出如下差分方程： （1）计算并画出单位脉冲响应。 （2）由此规定的系统是否稳定？ 解：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 已知某模拟信号，将它分别用不同的采样频率进行采样得到离散时间信号，试分析在以下两种采样情况下对信号频率的影响。 （1） 采样频率； （2） 采样频率 解： （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('离散信号')；hold on stem(n*Ts*1000,x);hold off subplot(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('离散信号')；hold on stem(n*Ts*1000,x);hold off subplot(2,1,2)；plot(w/pi,X);xlabel('以pi为单位的频率')；ylabel('X1(w)');title('离散信号傅里叶变换') 运行结果： 2-25 已知某系统的系统函数为 求其零、极点并绘出零、极点图。 解： MATLAB实现程序： b=[0.3 0.1 0.3 0.1 0.2];a=[1 -1.2 1.5 -0.8 0.3]; r1=roots(a) r2=roots(b) zplane(b,a) 运行结果： r1 = 0.1976 + 0.8796i 0.1976 - 0.8796i 0.4024 + 0.4552i 0.4024 - 0.4552i r2 = 0.3236 + 0.8660i 0.3236 - 0.8660i -0.4903 + 0.7345i -0.4903 - 0.7345i 2-26 已知因果系统，绘出的幅频和相频特性曲线。 解： 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);grid xlabel('');ylabel('幅度');title('幅频响应') subplot(2,1,2);plot(w/pi,phaH/pi);grid xlabel('频率( 单位：pi)');ylabel('相位（单位:pi）');title('相频响应') 运行结果： 2-27 一个线性时不变系统，描绘它的差分方程为 （1）在之间求得并画出系统的脉冲响应，从脉冲响应确定系统的稳定性； （2）画出该系统的幅频、相频特性； （3）如果此系统的输入为，在间求系统的输出； （4）讨论当输入改为时，输出波形如何变化？为什么？试根据系统的幅频特性解释。 解： 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.5714 3-4 是一4点序列： （1）计算离散时间傅里叶变换（DTFT）即,并画出它的幅度和相位。 （2）计算的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);grid xlabel('');ylabel('|X|');title('DTFT的幅度') subplot(2,1,2);plot(w/pi,phaH/pi*180);grid xlabel('以pi为单位的频率');ylabel('度');title('DTFT的相角') 运行结果： （2）MATLAB实现程序： N=4;k=0:N-1; X=dft(x,N); magX=abs(X),phaX=angle(X)*180/pi subplot(2,1,1);plot(w*N/(2*pi),magH,'--'); axis([-0.1,4.1,0,5]);hold on stem(k,magX); ylabel('|X(k)|');title('DFT的幅度：N=4')；text(4.3,-1,'k') hold off subplot(2,1,2);plot(w*N/(2*pi),phaH*180/pi,'--'); axis([-0.1,4.1,-200,200]);hold on stem(k,phaX); ylabel('度')；title('DFT相角：N=4')；text(4.3,-200,'k') 运行结果： 3-12 ，利用MATLAB程序求如下、。 （1）取的前10点数据，求N=10点的、并作图； （2）将补零至100点，求N=100点的、并作图。 解： （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<=90 zeros');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 利用脉冲响应不变法设计一个数字巴特沃斯低通滤波器，滤波器的技术要求为：通带截止频率rad/s，阻带最小衰减。研究不同采样频率对所设计数字滤波器频率响应的影响，设采样频率分别取1kHz,2kHz,4kHz。 解： （1）fs=1000Hz MATLAB实现程序： 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]);hold on fs=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]);hold off 运行结果： (2)fs=2000Hz MATLAB实现程序： 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]);hold on fs=2000;[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]);hold off 运行结果： （3）fs=4000Hz MATLAB实现程序： 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]);hold on fs=4000;[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]);hold off 运行结果： 5-9 已知利用双线性变换法设计一个ChebyshevI型数字高通滤波器。 解： MATLAB实现程序： Rp=1;Rs=20;T=0.001;fp=300;fs=200; wp=2*pi*fp*T;ws=2*pi*fs*T; wp1=(2/T)*tan(wp/2);ws1=(2/T)*tan(ws/2); [n,wn]=cheb1ord(wp1,ws1,Rp,Rs,'s');[b,a]=cheby1(n,Rp,wn,'high','s'); [bz,az]=bilinear(b,a,1/T); [db,mag,pha,grd,w]=freqz_m(bz,az);plot(w/pi,db);axis([0,1,-30,2]) 运行结果： 6-5 根据下列技术指标，设计一个FIR低通滤波器。 通带截止频率rad，通带允许波动； 阻带截止频率rad，阻带衰减。 解： MATLAB实现程序： Wp=0.2*pi;Ws=0.3*pi;Rp=0.25;Rs=50; N=80;n=[0:1:N-1]; Wc=(Ws+Wp)/2; hd=ideal_lp(Wc,N); w_han=(hamming(N))'; h=hd.*w_han [db,mag,pha,grd,w]=freqz_m(h,1); subplot(2,2,1);plot(w/pi,db);xlabel('\omega^pi');ylabel('dB'); subplot(2,2,2);plot(w/pi,pha);xlabel('\omega^pi');ylabel('pha'); subplot(2,2,3);stem(n,h,'.');xlabel('n');ylabel('h(n)'); subplot(2,2,4);stem(n,w_han,'.');xlabel('n');ylabel('Hamming window'); 运行结果： 6-6 利用频率采样法去设计一个线性相位FIR低通滤波器，已知 （1）采样点数； （2）采样点数；设置一个过渡点。 （3）采样点数；设置两个过渡点，。 解： （1）MATLAB实现程序： N=33;alpha=(N-1)/2;k=0:N-1;wk=(2*pi/N)*k; Hk=[ones(1,9),zeros(1,16),ones(1,8)]; angH=-alpha*(2*pi)/N*k;H=Hk.*exp(i*angH); h=real(ifft(H,N)); [db,mag,pha,grd,w]=freqz_m(h,1); [Hr,ww,a,L]=Hr_Type1(h); subplot(3,2,1); plot(ww/pi,Hr,wk(1:33)/pi,Hk(1:33),'o');title('幅度特性'); axis([0,1,-0.5,1.2]); subplot(3,2,2);plot(w/pi,db);axis([0,1,-100,10]);title('幅度特性'); 运行结果： （2）MATLAB实现程序： N=33;alpha=(N-1)/2;k=0:N-1;wk=(2*pi/N)*k; Hk=[ones(1,9),0.39,zeros(1,14),0.39,ones(1,8)]; angH=-alpha*(2*pi)/N*k;H=Hk.*exp(i*angH); h=real(ifft(H,N)); [db,mag,pha,grd,w]=freqz_m(h,1); [Hr,ww,a,L]=Hr_Type1(h); subplot(3,2,1); plot(ww/pi,Hr,wk(1:33)/pi,Hk(1:33),'o');title('幅度特性'); axis([0,1,-0.5,1.2]); subplot(3,2,2);plot(w/pi,db);axis([0,1,-100,10]);title('幅度特性'); 运行结果： （3）MATLAB实现程序： N=34;alpha=(N-1)/2;k=0:N-1;wk=(2*pi/N)*k; Hk=[ones(1,9),0.5925,0.1099,zeros(1,13),-0.1099,-0.5925,-ones(1,8)]; angH=-alpha*(2*pi)/N*k;H=Hk.*exp(i*angH); h=real(ifft(H,N)); [db,mag,pha,grd,w]=freqz_m(h,1); [Hr,ww,a,L]=Hr_Type2(h); subplot(3,2,1); plot(ww/pi,Hr,wk(1:34)/pi,Hk(1:34),'o');title('幅度特性'); axis([0,1,-0.5,1.2]); subplot(3,2,2);plot(w/pi,db);axis([0,1,-100,10]);title('幅度特性'); 运行结果：