信号与系统实验.doc

信号与系统实验
实验目的
1、 通过基于matlab的实验，加深对理论课程的理解
2、 初步掌握用matlab绘图、微积分、解微分差分方程的方法
3、 初步掌握matlab在信号系统中的应用

实验内容
一、信号的描述、运算、绘图
1、 用MATLAB生成下列函数，连续信号用plot，离散信号用stem绘图
(1)
t=linspace(0,5,11);y=stepfun(t,0);plot(t,y)

t=linspace(0,10,51);y=stepfun(t,0);z=sin(t).*y;plot(t,z)

t=linspace(0,10,101);y=stepfun(t,0);z=exp(-t).*y.*cos(t);plot(t,z)

t=linspace(-2,2,1001);y=rectpuls(t,2);plot(t,y)

t=-3*pi:0.001:3*pi;x=sinc(t)
;plot(t,x);grid on

(2)

k=linspace(-5,5,101);
x=cos(1/2*pi*k);
stem(k,x)

k=linspace(-5,5,51);
x=cos(1/8*pi*k);
stem(k,x)

k=linspace(-5,5,51);
x=cos(1/4*pi*k);
stem(k,x)

k=linspace(-5,5,101);
x=cos(pi*k);
stem(k,x)

k=linspace(-5,5,101);
x=cos(3/2*pi*k);
stem(k,x)

k=linspace(-5,5,51);
x=cos(7/4*pi*k);
stem(k,x)

k=linspace(-5,5,51);
x=cos(15/8*pi*k);
stem(k,x)

k=linspace(-5,5,51);
x=cos(2*pi*k);
stem(k,x)

k=linspace(-5,5,51);
x=cos(7/4*pi*k+1/3*pi);
stem(k,x)

k=linspace(-5,5,51);
y=cos(1/3*pi*k)+j*sin(1/3*pi*k);
stem(k,y)

k=linspace(-5,5,51);
x=cos(1/3*pi*k)+j*sin(1/3*pi*k);
y=0.5.^k;
z=x.*y;
stem(k,z)

k=linspace(-5,5,51);
y=exp(j*1/3*pi*k);
stem(k,y)

-5
-4
-3
-2
-1
0
1
2
3
4
5
-1
-0.8
-0.6
-0.4
-0.2
0
0.2
0.4
0.6
0.8
1

2、 熟悉ezplot、polar、mesh等指令
(1) 用ezplot绘的图。
(2) 用mesh绘的三维曲面。的取值范围是

ezplot('exp(-1*t).*cos(t)',[0,10])

x=-8:0.5:8;
y=x';
X=ones(size(y))*x;
Y=y*ones(size(x));
R=sqrt(X.^2+Y.^2)+eps;
Z=sin(R)./R;
mesh(X,Y,Z);
colormap(hot);
xlabel('x');
ylabel('y');
zlabel('z');

3、 计算数值积分：
(1)
t=linspace(0,pi,1001);
x=exp(-t).*sin(t);
y1=sum(x)*t(2),
y2=trapz(x)*t(2),
f=inline('exp(-t).*sin(t)'),
y3=quad(f,0,pi),

f=inline('sinc(2*t)');
y=quad(f,-10*pi,10*pi)

y =
0.5014
(2) ，，求
dt=0.01;
k1=0:dt:4;
f1=cos(3*k1);
k2=k1;
f2=stepfun(k2,0)-stepfun(k2,4);
f=dt*conv(f1,f2);
k0=k1(1)+k2(1);
k3=length(f1)+length(f2)-2;
k=k0:dt:k0+k3*dt;
subplot(2,2,1);plot(k1,f1);
subplot(2,2,2);plot(k2,f2);
subplot(2,2,3);plot(k,f);

二、 系统的描述、时域与变换域分析、求解
1、 解方程：
(1) ,y(-1)=1,y(-2)=2, f(k)=ε(k),求系统的h(k)、g(k)、y(k)
num=[1 1 0];den=[1 2 2];N=20;
k=0:19; y1=dimpulse(num,den,20);y2=dstep(num,den,20);
subplot(3,1,1);
stem(k,y1,'filled');
xlabel('k');ylabel('h(k)');title(' dimpulse ');
subplot(3,1,2);
stem(k,y2,'filled');
xlabel('k');ylabel('g(k)');title('dstep')
n=1:20+3;
F=stepfun(n,3);
F(1)=0;F(2)=0;
Y(1)=2;Y(2)=1;
for k=3:23
Y(k)=F(k)-F(k-1)-2*Y(k-1)-2*Y(k-2);
end
n1=n-3;
subplot(3,1,3)
stem(n1,Y,'filled');xlabel('k');ylabel('f(k)');

(2) (a)求解微分方程，初始条件为：y(0+)=y'(0+)=1，。理论计算全响应，用ezplot绘结果图。范围为区间[0,4]。
(b)用数值法求解需先化为差分方程。化出对应的差分方程。步长取T=0.1。
(c)确定此差分方程的初始条件。
(d)迭代法计算机求解。请编写程序，运行调试。比较数值解的准确程度。
Ts=0.1;
n=1:4/Ts+5;
F(n)=1;F(1)=0;F(2)=0;
Y(3)=1;Y(4)=1.1;
for k=5:length(n)
Y(k)=(13*F(k)-10*F(k-1)-100*Y(k-2)+240*Y(k-1))/143;
end
n1=(n-5)*Ts;subplot(2,1,1);plot(n1,Y);
xlabel('t');ylabel('y(t)');title('全响应(数值计算)');grid on;
subplot(2,1,2);ezplot('-0.5*exp(-3*t)+0.5*exp(-t)+1',[0,4]);
xlabel('t');ylabel('y(t)');title('理论全响应');grid on;

(3) 电压源激励的RLC串联回路微分方程为：。求当R=1Ω，C=1F，L=1H时二阶电路的冲激响应和阶跃响应并绘波形图。
a=[1 1 1];
b=[1];
subplot(2,1,1);impulse(b,a);title('冲激响应');
subplot(2,1,2);step(b,a);title('阶跃响应');

(4) 求解微分方程y''+3y'=3f(t)，初始条件y'(0)=1.5，y(0)=0，f(t)=ε(t)
y=dsolve('(D2y)+Dy*3=3*stepfun(n,0)','y(0)=0','Dy(0)=1.5')

y =
(stepfun(n, 0)/3 - 1/2)/exp(3*t) + t*stepfun(n, 0) - stepfun(n, 0)/3 + 1/2

(5) 求微分方程的冲激响应。
a=[1 4 3 ];b=[1 2];
impulse(b,a);
title('冲激响应');

(6)对连续系统y''+3y'+2y=3(ε(t)-ε(t-2))建立SIMULINK模型，并观察系统的零状态响应。

(7)求下面流图的传输函数。

syms z;
Q=vpa(zeros(7,7));
Q(2,1)=2;Q(2,3)=-1;
Q(2,5)=-2;Q(3,2)=1/z;Q(4,1)=3;Q(4,3)=1;Q(5,4)=1/z;Q(6,1)=1;Q(6,5)=2;
Q(7,6)=1;
P=[1;0;0;0;0;0;0];
I=eye(size(Q));
H=(I-Q)\P;
HS=simple(H(7))

HS =
(6*z + 10)/(z^2 + z + 2) + 1

2、 对如下系统函数绘制零极点图、判断稳定性、对稳定的绘制频率特性曲线。
(1)、、、。

b=[1 -2 -1];a=[2 0 0 -1];zplane(b,a);

w=0:0.01:10;freqz(b,a,w);

b=[1 1];a=[1 0 0 -1]; zplane(b,a);

w=0:0.01:10;freqz(b,a,w);

b=[1 0 2];a=[1 2 -4 1];zplane(b,a)

b=[1 0 0 0];a=[1 0.2 0.2 0.4];zplane(b,a)

w=0:0.01:10;freqz(b,a,w);

(2)、、。

b=[1 3];a=[1 3 6 4];pzmap(b,a);

w=0:0.01:10;freqz(b,a,w);

b=[1,2];a=[1 9 27 27];pzmap(b,a);

w=0:0.01:10;freqz(b,a,w);

b=[1 1];a=[1 2 1 1];pzmap(b,a)

w=0:0.01:10;freqz(b,a,w);

3、 用系统转换的指令得到上述系统函数的零极点增益形式、状态方程形式和多项式形式。
(1)、、、

b=[1 -2 -1];a=[2 0 0 -1];
[Z,P,K]=tf2ss(b,a)
[A,B,C,D]=tf2ss(b,a)
[b,a]=ss2tf(A,B,C,D)

Z =
Column 1
0
1.0000
0
Column 2
0
0
1.0000
Column 3
0.5000
0
0
P =
1
0
0
K =
Column 1
0.5000
Column 2
-1.0000
Column 3
-0.5000
A =
Column 1
0
1.0000
0
Column 2
0
0
1.0000
Column 3
0.5000
0
0
B =
1
0
0
C =
Column 1
0.5000
Column 2
-1.0000
Column 3
-0.5000
D =
0
b =
Column 1
0
Column 2
0.5000
Column 3
-1.0000
Column 4
-0.5000
a =
Column 1
1.0000
Column 2
-0.0000
Column 3
-0.0000
Column 4
-0.5000

b=[1 1];a=[1 0 0 -1]; [Z,P,K]=tf2ss(b,a)
[A,B,C,D]=tf2ss(b,a)
[b,a]=ss2tf(A,B,C,D)

Z =
Column 1
0
1
0
Column 2
0
0
1
Column 3
1
0
0
P =
1
0
0
K =
Column 1
0
Column 2
1
Column 3
1
A =
Column 1
0
1
0
Column 2
0
0
1
Column 3
1
0
0
B =
1
0
0
C =
Column 1
0
Column 2
1
Column 3
1
D =
0
b =
Column 1
0
Column 2
0
Column 3
1.0000
Column 4
1.0000
a =
Column 1
1.0000
Column 2
0
Column 3
0
Column 4
-1.0000

b=[1 0 2];a=[1 2 -4 1]; [Z,P,K]=tf2ss(b,a)
[A,B,C,D]=tf2ss(b,a)
[b,a]=ss2tf(A,B,C,D)

Z =
Column 1
-2
1
0
Column 2
4
0
1
Column 3
-1
0
0
P =
1
0
0
K =
Column 1
1
Column 2
0
Column 3
2
A =
Column 1
-2
1
0
Column 2
4
0
1
Column 3
-1
0
0
B =
1
0
0
C =
Column 1
1
Column 2
0
Column 3
2
D =
0
b =
Column 1
0
Column 2
1.0000
Column 3
-0.0000
Column 4
2.0000
a =
Column 1
1.0000
Column 2
2.0000
Column 3
-4.0000
Column 4
1.0000

b=[1 0 0 0];a=[1 0.2 0.2 0.4]; [Z,P,K]=tf2ss(b,a)
[A,B,C,D]=tf2ss(b,a)
[b,a]=ss2tf(A,B,C,D)

Z =
Column 1
-0.2000
1.0000
0
Column 2
-0.2000
0
1.0000
Column 3
-0.4000
0
0
P =
1
0
0
K =
Column 1
-0.2000
Column 2
-0.2000
Column 3
-0.4000
A =
Column 1
-0.2000
1.0000
0
Column 2
-0.2000
0
1.0000
Column 3
-0.4000
0
0
B =
1
0
0
C =
Column 1
-0.2000
Column 2
-0.2000
Column 3
-0.4000
D =
1
b =
Column 1
1
Column 2
0
Column 3
0
Column 4
0
a =
Column 1
1.0000
Column 2
0.2000
Column 3
0.2000
Column 4
0.4000

(2)、、。

b=[1 3];a=[1 3 6 4]; [Z,P,K]=tf2ss(b,a)
[A,B,C,D]=tf2ss(b,a)
[b,a]=ss2tf(A,B,C,D)

Z =
Column 1
-3
1
0
Column 2
-6
0
1
Column 3
-4
0
0
P =
1
0
0
K =
Column 1
0
Column 2
1
Column 3
3
A =
Column 1
-3
1
0
Column 2
-6
0
1
Column 3
-4
0
0
B =
1
0
0
C =
Column 1
0
Column 2
1
Column 3
3
D =
0
b =
Column 1
0
Column 2
0.0000
Column 3
1.0000
Column 4
3.0000
a =
Column 1
1.0000
Column 2
3.0000
Column 3
6.0000
Column 4
4.0000

b=[1,2];a=[1 9 27 27];
[Z,P,K]=tf2ss(b,a)
[A,B,C,D]=tf2ss(b,a)
[b,a]=ss2tf(A,B,C,D)

Z =
Column 1
-9
1
0
Column 2
-27
0
1
Column 3
-27
0
0
P =
1
0
0
K =
Column 1
0
Column 2
1
Column 3
2
A =
Column 1
-9
1
0
Column 2
-27
0
1
Column 3
-27
0
0
B =
1
0
0
C =
Column 1
0
Column 2
1
Column 3
2
D =
0
b =
Column 1
0
Column 2
0.0000
Column 3
1.0000
Column 4
2.0000
a =
Column 1
1.0000
Column 2
9.0000
Column 3
27.0000
Column 4
27.0000

b=[1 1];a=[1 2 1 1];
[Z,P,K]=tf2ss(b,a)
[A,B,C,D]=tf2ss(b,a)
[b,a]=ss2tf(A,B,C,D)

Z =
Column 1
-2
1
0
Column 2
-1
0
1
Column 3
-1
0
0
P =
1
0
0
K =
Column 1
0
Column 2
1
Column 3
1
A =
Column 1
-2
1
0
Column 2
-1
0
1
Column 3
-1
0
0
B =
1
0
0
C =
Column 1
0
Column 2
1
Column 3
1
D =
0
b =
Column 1
0
Column 2
-0.0000
Column 3
1.0000
Column 4
1.0000
a =
Column 1
1.0000
Column 2
2.0000
Column 3
1.0000
Column 4
1.0000

4、 用Simulink搭建上述系统
(1)、、、。

(2)、、。

三、 Fourier变换
1、 计算Fourier级数
(1)求下图周期三角波y(t)的傅立叶级数系数，绘制振幅谱和相位谱图。

t=linspace(-pi,pi,201);
n=0:9;n1=linspace(1,10,10);
gt= (n1'*((-1/pi)*abs(t)+1)).*exp(-i*n'*t);
Fn=1/(2*pi)*trapz(gt,2)*0.02;
sig=abs(Fn)>0.00001;
Fn=Fn.*sig;
figure;
subplot(2,1,1);
stem(n,abs(Fn));
xlabel('n'),ylabel('Fn');title('振幅谱');grid;
subplot(2,1,2);
stem(n,angle(Fn));
xlabel('n'),ylabel('Theta');title('相位谱');grid;

(3) 计算按照T=10周期延拓的周期信号的傅立叶系数，并绘制幅频、相频曲线。

t=linspace(-3,3,101);
n=0:9;n1=[1 1 1 1 1 1 1 1 1 1];
gt=(n1'*(cos(20*t).*rectpuls(t,6))).*exp(-i*n'*t);
Fn=1/10*trapz(gt,2)*0.02;
sig=abs(Fn)>0.00001;
Fn=Fn.*sig;
figure;
subplot(2,1,1);
stem(n,abs(Fn));
xlabel('n'),ylabel('Fn');title('振幅谱');grid;
subplot(2,1,2);
stem(n,angle(Fn));
xlabel('n'),ylabel('Theta');title('相位谱');grid;

2、 求Fourier变换与反变换（用数值方法作，计算DFT的函数fft()）：
，，，
绘制幅频和相频特性曲线。然后计算傅立叶逆变换计算时域波形并绘图，和原始函数进行对比，观察频谱混叠与频谱泄露现象

T=10;
N=512;
Ts=T/N;
PT=N/16;
t=linspace(-T/2,T/2-T/N,N);
f=0*t;
f(t>-2&t<2)=1;
subplot(3,1,1);plot(t,f);grid;
w=[0:N-1]*2*pi/(Ts*N);
Fw=fft(f)*Ts./exp(-i*0.5*T*w);
Fw=fftshift(Fw);
w=w-2*pi*(1/(Ts*2));
subplot(3,1,2);
plot(w(N/2-PT:N/2+PT),abs(Fw(N/2-PT:N/2+PT)));grid;
subplot(3,1,3);
plot(w(N/2-PT:N/2+PT),angle(Fw(N/2-PT:N/2+PT)));grid;

T=10;
N=200;
t=linspace(-T/2,T/2-T/N,N);
f=0*t;f(t>-2&t<2)=1;
W=16*pi;
K=100;
w=linspace(-W/2,W/2-W/K,K);
F=0*w;
for k=1:K
for n=1:N
F(k)=F(k)+T/N*f(n)*exp(-j*w(k)*t(n));
end
end
f1=0*t;
for n=1:N
for k=1:K
f1(n)=f1(n)+W/(2*pi*K)*F(k)*exp(j*w(k)*t(n));
end
end
subplot(2,1,1);
plot(t,f,'-k',t,f1,':k');
xlabel('t');ylabel('f(t)');legend('f(t)','反变换得f(t)');
subplot(2,1,2);
plot(w,F,'-k');
xlabel('w');ylabel('F(w)');grid;

Warning: Imaginary parts of complex X and/or Y arguments ignored
Warning: Imaginary parts of complex X and/or Y arguments ignored