信号与信号处理实验参考答案.doc

实验一 熟悉MATLAB环境2、（2）粗略描绘下列各函数的波形说明：MATLAB中有函数 f(t)=sint/tt=-3*pi:0.01*pi:3*pi; t1=t/pi;y=sinc(t1); plot(t,y); hold onplot(t,0)在自己的工作目录work下创建Heaviside函数的M文件，该文件如下：function f=Heaviside(t)f=(t0) %t0时f为1，否则f为0在命令窗口输入如下语句，就能绘出u(t)的波形。t=-1:0.01:3;f=Heaviside(t);plot(t,f)axis(-1 3 0.2 1.2)t=-1:0.01:2;g=Heaviside(t)-Heaviside(t-1); plot(t,g);axis(-1 2 -0.2 1.2)hold onplot(t,0)4、 分别用for和while循环结构编写程序，求出s=1+2+22+23+262+263并考虑一种避免循环的简洁方法来进行求和。 程序如下： s=1;for k=1:63 s=s+2k;ends运行结果是：s =1. 8447e+019（2）s=1;k=1;while k=63 s=s+2k; k=k+1;ends运行结果：s =1. 8447e+019(3) k=0:63;s=sum(2.k)实验二 信号的卷积与系统的响应和阶跃响应1 n=0:20;hn=0.9.n;xn=0,0 ones(1,8),0,0;yn=conv(hn,xn);stem(yn)3. 利用MATLAB绘制下列信号的卷积积分f1(t)*f2(t)的时域波形。（1） f1(t)=2u(t+1)-u(t-1), f2(t)=u(t+2)-u(t-2)（2） f1(t)=tu(t), f2(t)=u(t)（3） f1(t)=u(t)-u(t-4), f2(t)=sin(t)u(t);（4） f1(t)=e-2tu(t), f2(t)=e-tu(t)(1) 先编写实现连续信号卷积的通用函数sconv(),程序如下：functionf,k=sconv(f1,f2,k1,k2,p)%计算连续信号卷积积分f(t)=f1(t)*f2(t)%f:卷积积分f(t)对应的非零样值向量%k:f(t)的对应时间向量%f1:f1(t)非零样值向量%f2:f2(t)非零样值向量%k1:f1(t)的对应时间向量%k2:序列f2(t)的对应时间向量%p:取样时间间隔f=conv(f1,f2); %计算序列f1与f2的卷积ff=f*p;k0=k1(1)+k2(1); %计算序列f的非零样值的起点位置k3=length(f1)+length(f2)-2; %计算卷积和f的非零样值的宽度k=k0:p:(k3-(0-k0)/p)*p); %确定卷积和f非零样值的时间向量subplot(2,2,1)plot(k1,f1) %绘制f1(t)title(f1(t)xlabel(t)ylabel(f1(t)subplot(2,2,2)plot(k2,f2)title(f2(t)xlabel(t)ylabel(f2(t)subplot(2,2,3)plot(k,f);h=get(gca,position);h(3)=2.5*h(3);set(gca,position,h) %将第三个子图的横坐标范围扩为原来的2.5倍title(f(t)=f1(t)*f2(t)xlabel(t)ylabel(f(t)p=0.01;k1=-1:p:1;f1=2*ones(1,length(k1);k2=-2:p:2;f2=ones(1,length(k2);f,k=sconv(f1,f2,k1,k2,p)（2）p=0.01;k1=0:p:10;f1=k1k2=0:p:10;f2=ones(1,length(k2);f,k=sconv(f1,f2,k1,k2,p)第（2）题图上实验二 信号的卷积与系统的响应1n=0:20;hn=0.9.n;xn=stepseq(2,0,20)-stepseq(10,0,20);yn=conv(hn,xn);stem(yn)2（1）p=0.01;k1=-2:p:2;f1=2*(u(k1+1)-u(k1-1);k2=k1;f2=u(k2+2)-u(k2-2);f,k=sconv(f1,f2,k1,k2,p)（2）p=0.01;k1=-1:p:10;f1=k1.*u(k1);k2=k1;f2=u(k2);f,k=sconv(f1,f2,k1,k2,p)（3）p=0.01;k1=-4:p:10;f1=u(k1)-u(k1-4);k2=k1;f2=sin(pi*k2).*u(k2);f,k=sconv(f1,f2,k1,k2,p)（4）5已知描述某连续系统的微分方程为：y(t)+5y(t)+8y(t)=3f(t)+2f(t)绘出系统的冲激响应波形，求出t=0.5s, 1s, 1.5s, 2s系统冲激响应的数值解。6、step()函数的调用格式与impulse()函数类似。已知描述某连续系统的微分方程为：2y(t)+y(t)+8y(t)=f(t)试用MATLAB绘出该系统的冲激响应与阶跃响应的波形。5. 程序如下：a=1 5 8;b=3 0 2;impulse(b,a)y=impulse(b,a,0:0.5:2)y = -15.0000 -1.3292 0.6302 0.3941 0.12656. a=2 1 8;b=1; impulse(b,a) hold on step(b,a)legend(冲激响应,阶跃响应)7. b=1;a=1 1 0.9; n=-10:50;impz(b,a,n)实验三 利用DFT分析模拟信号频谱1T=0.1;n=0:20;xn=exp(-2*n.*T);k=0:31Xk=fft(xn,32);XM=T.*XkMagXk=abs(XM)stem(k,XM)xlabel(k);ylabel(XM)矩形脉冲的频谱T=0.1;n=-10:10;xn=Heaviside(n*T+1)-Heaviside(n*T-1)Xk=fft(xn,256);Xk=fftshift(T*abs(Xk);plot(Xk)2fs=400;n=0:160;xn=cos(2*pi*100*n/fs)+0.75*cos(2*pi*110*n/fs);Xk=fft(xn,256);Xk=(abs(Xk)/fs);k=0:255;plot(fs/256*k,Xk)3t=-0.2:0.005:0.2;x=square(2*pi*20*t,50); %产生频率为20Hz的方波subplot(1,2,1);plot(t,x);axis(-0.2 0.2 -1.5 1.5)T=0.0025;N=20;n=0:N-1;xn=square(2*pi*20.*n*T,50)yn=0.5*(xn+1)Yn=fft(yn);Yn=fftshift(abs(Yn);Ym=1/N*Ynsubplot(1,2,2);stem(1/T)/N*n,Ym)实验五 离散信号与离散系统分析基础参考答案functionz,p=ljdt(A,B)%绘离散系统零极点图的函数p=roots(A); %求系统极点z=roots(B); %求系统零点p=p; %将极点列向量转置为行向量z=z; %将零点列向量转置为行向量x=max(abs(p z 1); %确定纵坐标范围x=x+0.1;y=x; %确定横坐标范围clchold onaxis(-x x -y y) %确定坐标轴显示范围w=0:pi/300:2*pi;t=exp(i*w); %画单位圆plot(t)axis(square)plot(-x x,0 0) %画横坐标轴plot(0 0,-y y) %画纵坐标轴text(0.1,x,jImz)text(y,0.1,Rez)plot(real(p),imag(p),x) %画极点plot(real(z),imag(z),o) %画零点 title(离散系统零极点图) %标注标题hold off2. 离散时间系统的系统函数如下所示，绘制系统的零极点图，并判断系统的稳定性。(1) (2)A=1 3 2 1;B=3 2 2 5; ljdt(A,B)（2）A= 1 2 -4 1; B=1 0 2;ljdt(A,B)3B=1 2 1; A=1 -0.5 -0.005 0.3; ljdt(A,B) H,w=freqz(B,A); figure;subplot(1,2,1);plot(w/pi,abs(H) subplot(1,2,2);plot(w/pi,angle(H)4.A=1 0 0.61;B=1 -2 1;H,w=freqz(B,A);subplot(1,2,1)plot(w/pi,abs(H)xlabel(频率/pi)ylabel(频率响应幅度)grid onsubplot(1,2,2)plot(w/pi,angle(H)xlabel(频率/pi)ylabel(频率响应相位)grid on5.B=1 1 1;A=1 1 -2 0;r p k=residue(B,A)r = 0.5000 1.0000 -0.5000p = -2 1 0k =实验六 离散傅里叶变换参考程序1 dft1.m: 用for循环计算DFTfunction Am,pha=dft1(x)N=length(x);w=exp(-j*2*pi/N);for k=1:N sum=0; for n=1:N sum=sum+x(n)*w(k-1)*(n-1); end Am(k)=abs(sum); pha(k)=angle(sum);end调用dft1.m函数文件的M程序文件：x=ones(1,8),zeros(1,248);t=cputime;Am1,pha1=dft1(x);t1=cputime-t;n=0:(length(x)-1);w=(2*pi/length(x)*n;figure(1)subplot(2,1,1),plot(w,Am1,b);title(Magnitude fart);xlabel(frequency in fadians);ylabel(|X(exp(jw)|)subplot(2,1,2),plot(w,pha1,r);grid;gtext(Phase part);xlabel(frequency in fadians);ylabel(arg(Xexp(jw)/radians);2dft2.m: 用MATLAB矩阵运算计算DFTfunction Am,pha=dft2(x)N=length(x);n=0:N-1;k=0:N-1;w=exp(-j*2*pi/N);nk=n*k;wnk=w.nk;Xk=x*wnk;Am=abs(Xk);pha=angle(Xk);3. dft3.m：调用FFT库函数计算DFTfunction Am,pha=dft3(x) Xk=fft(x); Am=abs(Xk); pha=angle(Xk)t1 =0. 32802采样16点，做16点DFT,程序如下：fs=32000;n=1:16;x=cos(2*pi*6.5*1000*n/fs)+cos(2*pi*7*1000*n/fs)+cos(2*pi*9*1000*n/fs);subplot(2,2,1);stem(n,x);title(signal x(n)Xk=fft(x,16);magXk=abs(Xk);k1=0:15;w1=2*pi/16*k1;subplot(2,1,2);plot(w1,magXk);title(16点DTFT的幅值);xlabel(w(单位：) 采样16点，做256点DFT,程序如下：fs=32000;n=1:16;x=cos(2*pi*6.5*1000*n/fs)+cos(2*pi*7*1000*n/fs)+cos(2*pi*9*1000*n/fs);subplot(2,2,1);stem(n,x);title(signal x(n)Xk=fft(x,256);magXk=abs(Xk);k1=0:255;w1=2*pi/256*k1;subplot(2,1,2);plot(w1,magXk);title(16点DTFT的幅值);xlabel(w(单位：pi)结果如图二（3）采集数据长度N=256点，做256点DFT,结果如图三fs=32000;n=0:255;x=cos(2*pi*6.5*1000*n/fs)+cos(2*pi*7*1000*n/fs)+cos(2*pi*9*1000*n/fs);subplot(2,2,1);stem(n,x);title(signal x(n)Xk=fft(x,256);magXk=abs(Xk);k1=0:255;w1=2*pi/256*k1;subplot(2,1,2);plot(w1,magXk);title(256点DTFT的幅值);xlabel(w(单位：pi)3参考程序和运行结果n=0:1:15;xn=0:1:15;subplot(2,2,1)stem(n,xn);xlabel(n);ylabel(x(n);m=0:1:7;hn=ones(1,8);subplot(2,2,2)stem(m,hn)xlabel(n);ylabel(h(n);N=length(n)+length(m)-1;Xk=fft(xn,N);Hk=fft(hn,N);Yk=Xk.*Hk;yn=ifft(Yk,N);if all(imag(xn)=0)&(all(imag(hn)=0) %实序列的循环卷积仍为实序列 yn=real(yn)endsubplot(2,2,3)x=0:N-1;stem(x,yn)xlabel(n);ylabel(y(n)=x(n)*h(n)实验八 IIR数字滤波器的设计1wp=0.2*pi; %digital Passband freq in radws=0.3*pi; %digital Stopband freq in radFs=1000; wp=0.2*pi*Fs;ws=0.3*pi*Fs; Rp=1;Rs=15; N,Wn=buttord(wp,ws,Rp,Rs,s); 返回模拟滤波器的最小阶数和截止频率Z,P,K=buttap(N); 模拟低通滤波器原型Bap,Aap=zp2tf(Z,P,K);b,a=lp2lp(Bap,Aap,Wn);bz,az=bilinear(b,a,Fs);H,W=freqz(bz,az);plot(W,abs(H);gridxlabel(频率/弧度) ylabel(频率响应幅度)H,W=freqz(bz,az,50)figurestem(W,abs(H)23. 参考程序fs=1000;wp=400*2/fs;ws=200*2/fs;rp=3;rs=15;Nn=128;N,wn=buttord(wp,ws,rp,rs)b,a=butter(N,wn,high)freqz(b,a,Nn,fs)运行结果N = 2wn = 0.6630b = 0.1578 -0.3155 0.1578a = 1.0000 0.6062 0.23734参考程序fs=1000;wp=100 250*2/fs;ws=50 300*2/fs;rp=3;rs=30;Nn=128;N,wn=cheb2ord(wp,ws,rp,rs)b,a=cheby2(N,rp,wn,bandpass)freqz(b,a,Nn,fs)运行结果实验九 FIR数字滤波器的设计1参考程序clear all;Wp=0.2*pi;Ws=0.4*pi;tr_width=Ws-Wp; %过渡带宽度N=ceil(6.6*pi/tr_width)+1 %滤波器长度（选择海明窗）n=0:1:N-1;Wc=(Ws+Wp)/2; %理想低通滤波器的截止频率alpha=(N-1)/2;m=n-alpha+eps;hd=sin(Wc*m)./(pi*m); %理想低通滤波器的单位冲激响应w_ham=(hamming(N); %海明窗h=hd.*w_ham; %截取得到实际的单位脉冲响应H,w=freqz(h,1,1000,whole);H=(H(1:1:501);w=(w(1:1:501);mag=abs(H);db=20*log10(mag+eps)/max(mag); %相对幅值响应pha=angle(H);delta_w=2*pi/1000;Ap=-(min(db(1:1:Wp/delta_w+1) %实际通带文波As=-round(max(db(Ws/delta_w+1:1:501) %实际阻带文波subplot(221)stem(n,hd)title(理想单位脉冲响应hd(n)subplot(222)stem(n,w_ham)title(海明窗w(n)subplot(223)stem(n,h)title(实际单位脉冲响应h(n)subplot(224)plot(w/pi,db)title(幅度响应(dB)axis(0,1,-100,10)运行结果N = 34Ap = 0.0477As = 5291另一程序clear all;Wp=0.2*pi;Ws=0.4*pi;tr_width=Ws-Wp; %过渡带宽度N=ceil(6.6*pi/tr_width)+1 %滤波器长度（选择海明窗）n=0:1:N-1;Wc=(Ws+Wp)/2; %理想低通滤波器的截止频率h=fir1(N,Wc/pi,hamming(N+1)subplot(221)stem(h)H,w=freqz(h,1,512);mag=abs(H);db=20*log10(mag+eps)/max(mag); %相对幅值响应subplot(222)plot(w/pi,db)2window=hanning(16); b=fir1(15,0.3 0.5,window); freqz(b,1,512)window=hanning(46); b=fir1(45,0.3 0.5,window); freqz(b,1,512)4. rp=3;rs=40;Fs=2000;f=500 600;a=1 0;dev=(10(rp/20)-1)/(10(rp/20)+1) 10(-rs/20); %实际幅频响应与期望幅频响应间的最大误差n,fo,ao,w=remezord(f,a,dev,Fs);b=remez(n,fo,ao,w);h,f=freqz(b,1,1024,Fs)plot(f,20*log10(abs(h)grid实验十 参考答案Fs=400;t=0:1/Fs:1;x=square(2*pi*10*t);subplot(221)plot(t,x)title(频率为10Hz的方波信号)axis(0 1 -1.5 1.5)%求方波信号的频谱T=0.1;t1=0:1/Fs:(T-1/Fs);xn=square(2*pi*10*t1);N=Fs*T;Xk=fft(xn,N);Xk=fftshift(abs(Xk);Xk1=Xk/N;subplot(222)stem(Fs/40*0:1:39,Xk1)title(方波信号的频谱)%设计一个数字低通滤波器，截止频率为40Hzwp=40*2*pi;ws=50*2*pi;Rp=3;Rs=40;N,Wn=buttord(wp,ws,Rp,Rs,s); % 返回模拟滤波器的最小阶数和截止频率Z,P,K=buttap(N); %模拟低通滤波器原型Bap,Aap=zp2tf(Z,P,K);b,a=lp2lp(Bap,Aap,Wn);bz,az=bilinear(b,a,Fs);H,W=freqz(bz,az);subplot(223)plot(W,abs(H);gridy=filter(bz,az,x);subplot(224)plot(t,y)%滤波器输出信号的频谱yn=filter(bz,az,xn)Yk=fft(yn,40)Yk=fftshift(abs(Yk);Yk1=Yk/Nfigurestem(Fs/21*0:1:20,Yk1)title(滤波后的信号的频谱)实验十一 信号的调制与解调1Fs=1000; %被调信号x的采样频率Fc=40; %载波信号的频率N=1000;n=0:N;t=n/Fs;x=sin(2*pi*5*t); %被调信号y=modulate(x,Fc,Fs,am)subplot(2,2,1)plot(t,x)xlabel(t(s)ylabel(x)title(被调信号)axis(0 1 -1 1)w1=350;k=-N:N;w=k*w1/N;X=x*exp(-j*t*w)*(1/Fs); %求X(jw)X=real(X);Y=y*exp(-j*t*w)*(1/Fs); %求Y(jw)Y=real(Y);subplot(2,2,2);plot(t,y)xlabel(t);ylabel(y(t)=x(t)*cos(2*pi*40*t);subplot(2,2,3);plot(w,X);xlabel(w);ylabel(X(jw);subplot(2,2,4);plot(w,Y);xlabel(w);ylabel(Y(jw);2Fs=1000;Fc=100;N=1000;n=0:N;t=n/Fs;x=sin(4*pi*t);y1=modulate(x,Fc,Fs,am);y2=modulate(x,Fc,Fs,Fm);y3=modulate(x,Fc,Fs,Pm);figuresubplot(2,2,1)plot(t,x)title(sin(4*pi*t)subplot(2,2,2)plot(t,y1)title(调幅信号)subplot(2,2,3)plot(t,y2)title(调频信号)subplot(2,2,4)plot(t,y3)title(调相信号)3Fs=1000;Fc=100;N=1000;n=0:N;t=n/Fs;x=sin(4*pi*t);y1=modulate(x,Fc,Fs,am);y2=modulate(x,Fc,Fs,Fm);y3=modulate(x,Fc,Fs,Pm);x1=demod(y1,Fc,Fs,am)x2=demod(y2,Fc,Fs,Fm);x3=demod(y3,Fc,Fs,Pm);figuresubplot(2,2,1)plot(t,x1)title(调幅信号的解调信号)subplot(2,2,2)plot(t,x2)title(调频信号的解调信号)subplot(2,2,3)plot(t,x3)title(调相信号的解调信号)4t=-2:0.01:2;ft=Heaviside(t+1)-Heaviside(t-1);subplot(221)plot(t,ft)axis(-2 2 -0.5 1.5)f1t=ft.*cos(10*pi*t);subplot(222)plot(t,f1t)axis(-2 2 -1.5 1.5)%求矩形脉冲