数字信号处理实验.doc

数字信号处理实验 电信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 onstem(n*Ts*1000,x);hold offsubplot(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 onstem(n*Ts*1000,x);hold offsubplot(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.4552ir2 = 0.3236 + 0.8660i 0.3236 - 0.8660i -0.4903 + 0.7345i -0.4903 - 0.7345i2-26 已知因果系统，绘出的幅频和相频特性曲线。解：MATLAB实现程序：b=1,0;a=1,-0.9H,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 一个线性时不变系统，描绘它的差分方程为（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.57143-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);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);hold onstem(k,magX);ylabel(|X(k)|);title(DFT的幅度：N=4)；text(4.3,-1,k)hold offsubplot(2,1,2);plot(w*N/(2*pi),phaH*180/pi,-);axis(-0.1,4.1,-200,200);hold onstem(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=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);hold onfs=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=2000HzMATLAB实现程序：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 onfs=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=4000HzMATLAB实现程序：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 onfs=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_handb,mag,pha,grd,w=freqz_m(h,1);subplot(2,2,1);plot(w/pi,db);xlabel(omegapi);ylabel(dB);subplot(2,2,2);plot(w/pi,pha);xlabel(omegapi);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(幅度特性);运行结果：