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

信号与信号处理实验参考答案
实验⼀熟悉MATLAB 环境
2、（2）粗略描绘下列各函数的波形说明：MA TLAB 中有函数t
t
t c ππsin )(sin = ④ f(t)=sint/t
t=-3*pi:0.01*pi:3*pi; t1=t/pi; y=sinc(t1); plot(t,y); hold on
plot(t,0)
⑤在⾃⼰的⼯作⽬录work 下创建Heaviside 函数的M ⽂件，该⽂件如下：function f=Heaviside(t)
f=(t>0) %t>0时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 on plot(t,0)
4、分别⽤for 和while 循环结构编写程序，求出
s=
∑=63
2
k k
=1+2+22+23+…+262+263
并考虑⼀种避免循环的简洁⽅法来进⾏求和。

程序如下： s=1; for k=1:63
s=s+2^k;
s
运⾏结果是：
s =
1.8447e+019
（2）s=1;
k=1;
while k<=63
s=s+2^k;
k=k+1;
end
s
运⾏结果：
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. 利⽤MA TLAB绘制下列信号的卷积积分f1(t)*f2(t)的时域波形。

（1）f1(t)=2[u(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-2t u(t), f2(t)=e-t u(t)
(1) 先编写实现连续信号卷积的通⽤函数sconv(),程序如下：function[f,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的卷积f
f=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;
k2=0:p:10;
f2=ones(1,length(k2)); [f,k]=sconv(f1,f2,k1,k2,p)
第（2）题图上
实验⼆信号的卷积与系统的响应1．
n=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));
f2=u(k2+2)-u(k2-2);
[f,k]=sconv(f1,f2,k1,k2,p)
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)
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 =
-15.0000
-1.3292
0.6302
0.3941
0.1265
6. 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分析模拟信号频谱1．T=0.1;
n=0:20;
xn=exp(-2*n.*T);
k=0:31
XM=T.*Xk
MagXk=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)
2．fs=400;
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)
3．t=-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*Yn
subplot(1,2,2);
stem((1/T)/N*n,Ym)
实验五离散信号与离散系统分析基础参考答案
function[z,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; %确定横坐标范围 clc hold on
axis([-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,'jIm[z]') text(y,0.1,'Re[z]')
plot(real(p),imag(p),'x') %画极点 plot(real(z),imag(z),'o') %画零点 title('离散系统零极点图') %标注标题 hold off
2. 离散时间系统的系统函数如下所⽰，绘制系统的零极点图，并判断系统的稳定性。

(1)1235223)(2323++++++=z z z z z z z H (2)1
422
)(232+-++=z z z z z H
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)
3．B=[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 on
subplot(1,2,2)
plot(w/pi,angle(H))
xlabel('频率/pi')
ylabel('频率响应相位')
grid on
5.
B=[1 1 1];
A=[1 1 -2 0];
[r p k]=residue(B,A) r =
0.5000
1.0000
-0.5000
p =
-2
1
k =
[]
f nδ5
u
n
2
n
0-
=
-
u
+
5
.
(
)
(
)
)
(
(n
)
(.
)
n
实验六离散傅⾥叶变换参考程序1．dft1.m: ⽤for循环计算DFT function [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(X[exp(jw)]/radians');
2．dft2.m: ⽤MA TLAB矩阵运算计算DFT function [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库函数计算DFT function [Am,pha]=dft3(x)
Xk=fft(x);
Am=abs(Xk);
pha=angle(Xk)
t1 =
0.3280
2．采样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) end
subplot(2,2,3)
x=0:N-1;
stem(x,yn)
xlabel('n');ylabel('y(n)=x(n)*h(n)')
实验⼋IIR数字滤波器的设计
1．wp=0.2*pi; %digital Passband freq in rad
ws=0.3*pi; %digital Stopband freq in rad
Fs=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));
grid
xlabel('频率/弧度')
ylabel('频率响应幅度')
[H,W]=freqz(bz,az,50)
figure
stem(W,abs(H))
2．
3. 参考程序
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 =
2
wn =
0.6630
b =
0.1578 -0.3155 0.1578
a =
1.0000 0.6062 0.2373
4．参考程序
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)。