信号与信号处理实验参考答案_新0_20070611.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);%输出一个波形，以t为纵坐标，y函数为横坐标。 hold on%不用的话 前面的图会消失的，要在第一个波形后面用。plot(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;%离散的所以是用矩阵赋值，连续的用HEAVISIDE。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)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)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）题图上，第（3）题图下，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)plot(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)另一个：fs=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);Xk=fftshift(Xk);k=0:255;plot(-fs/2+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 radians);ylabel(|X(exp(jw)|)subplot(2,1,2),plot(w,pha1,r);grid;gtext(Phase part);xlabel(frequency in radians);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. N=16;f=0:8/8;a=1,1,1,1,0.389,0,0,0,0;h=fir2(N,f,a);freqz(h,1)实验十 信号滤波：参考程序1 IIR滤波器实现clcclearFs=1000;T=1/10;t=0:1/Fs:1;x=square(2*pi*10*t);%产生幅度为-1,1的方波plot(t,x)title(矩形波)axis(0,0.2,-1.5,1.5)N=T*Fs; %方波的频谱t1=0:1/Fs:(0.1-1/Fs);xn=square(2*pi*10*t1);fx=fftshift(fft(x,N)/N;figuref=(-N/2:(N/2-1)*Fs/N;stem(f,abs(fx)title(方波信号的频谱)n,Wn=buttord(40/500,50/500,3,15); %确定butterworth滤波器阶次b,a=butter(n,Wn); %确定滤波器传递函数figurefreqz(b,a); %滤波器频率特性title(低通滤波器频率特性)y=filter(b,a,x); %滤波figureplot(t,y)title(矩形波通过低通滤波器后的波形)axis(0,0.2,-1.5,1.5)N=512; %滤波后的频谱fx=fftshift(fft(y,N)*(1/Fs);figuref=(-N/2:(N/2-1);f=f*Fs/N;plot(f,abs(fx)title(滤波后信号的频谱)s=randn(1,1001); %信号叠加噪声xs=x+s;figureplot(t,xs)title(信号叠加噪声的波形)%axis(0,0.5,-2,2)n,Wn=buttord(80/500,100/500,1,15); %确定butterworth滤波器阶次b,a=butter(n,Wn); %确定滤波器传递函数figurefreqz(b,a); %滤波器频率特性title(低通滤波器频率特性)y1=filter(b,a,xs); %采用filter函数滤波figureplot(t,y1) %滤波后的波形title(滤波后的波形)%axis(0,0.5,-2,2)参考程序2 FIR滤波器实现%采用FIR滤波器clcclearFs=1000;T=1/10;t=0:1/Fs:1;x=square(2*pi*10*t);%产生幅度为-1,1的方波figureplot(t,x)title(矩形波)axis(0,0.5,-1.5,1.5)N=T*Fs; %方波的频谱t1=0:1/Fs:(0.1-1/Fs);xn=square(2*pi*10*t1);fx=fftshift(fft(x,N)/N;figuref=(-N/2:(N/2-1)*Fs/N;stem(f,abs(fx)title(方波信号的频谱)n=48win=chebwin(n+1,40);%窗函数法设计FIR滤波器(切比雪夫窗），也可以采用其他窗bt=fir1(n,45/500,win); %截止频率45hzfigurefreqz(bt,1); %滤波器频率特性title(低通滤波器频率特性)y=filter(bt,1,x); %滤波figureplot(t,y)title(矩形波通过低通滤波器后的波形)axis(0,0.5,-1.5,1.5)N=512; %滤波后的频谱fx=fftshift(fft(y,N)*(1/Fs);figuref=(-N/2:(N/2-1);f=f*Fs/N;plot(f,abs(fx)title(滤波后信号的频谱)s=randn(1,1001); %信号叠加噪声xs=x+s;figureplot(t,xs)title(信号叠加噪声的波形)axis(0,0.5,min(xs),max(xs)n=48win=chebwin(n+1,40);%窗函数法设计FIR滤波器(切比雪夫窗），也可以采用其他窗bt=fir1(n,80/500,win); %截止频率80hz（可自己调整）figurefreqz(bt,1); %滤波器频率特性title(低通滤波器频率特性)y1=filter(bt,1,xs); %采用filter函数滤波figureplot(t,y1) %滤波后的波形title(滤波后的波形)axis(0,0.5,min(y1),max(y1)fxs=fftshift(fft(y1,N)*(1/Fs); %滤波后波形的频谱figuref=(-N/2:(N/2-1);f=f*Fs/N;plot(f,abs(fxs) title(滤波后波形的频谱)二参考程序clcclearmusic=wavread(d:matlab6p5toolboxdspblksdspdemosdspafsx_mono.wav);%可根据自己电脑情况调整music=music;a=size(music)%利用media player观察音乐持续时间为3s，信号的采样点数为60000，可知信号的采样频率为20KHzmusic=decimate(music,2);%音乐的采样频率降低为10KHzt=0:0.0001:(3-0.0001);x=sin(2*pi*1000*t)+sin(2*pi*3000*t);mx=music+x;xa=0.1;%显示x轴上限，可调整figuresubplot(3,1,1),