Matlab-信号系统笔记

备注：MOOC上西电《信号与系统》课程中的MATLAB编程部分笔记,chapter为对应章节。

chapter1

1.连续信号的绘制

%连续信号
%f(t) = 5e.^(-0.8t)sin(πt) 0<t<5
b = 5;a = 0.8;t = 0:0.001:5;
x = b*exp(-a*t).*sin(pi*t)
plot(t,x)

2.离散信号的绘制

%离散信号
%f(k) = 2(0.8).^k
c = 2;d = 0.8;k = -5:5;
y = c*d.^k;
stem(k,y)

chapter2

1.求解系统响应

求LTI系统的零状态相应的函数为lsim，调用格式为：

y = lsim(sys,f,t)

其中，t表示计算系统相应的抽样点向量；f是系统输入信号，sys是LTI系统模型，表示微分方程。
系统模型sys要借助tf函数，调用格式为：

sys = tf(b,a)

b和a分别为微分方程的右端和左端的各项系数。
例如： ${y}''(t)+5{y}'(t)+6y(t) = {f}''(t) + 2f(t)$

%微分方程
a=[1,5,6];b=[1,0,2];
sys=tf(b,a)

例题：系统的微分方程为 ${y}''(t)+2{y}'(t)+77y(t) = f(t)$ ，在t≥0接入激励 $f(t)=10sin(2\pi t)$ ，求零状态响应。

b=[1];a=[1 2 77];
sys=tf(b,a);
t = 0:0.01:5;
f = 10*sin(2*pi*t);
y = lsim(sys,f,t);
plot(t,y);
xlabel('Time(sec)')
ylabel('y(t)')

2.求解冲激响应和阶跃响应

LTI系统的微分方程为：

$\sum_{i=1}^{n}a_{i}y^{(i)}(t) = \sum_{j=1}^{m}b_{i}f^{(i)}(t)$

冲激响应的函数为：impulse(b,a)

阶跃响应的函数为：step(b,a)

例题：求下列系统的冲激响应和阶跃响应。

$7{y}''(t)+4{y}'(t)+6y(t) = {f}'(t) + f(t)$

a=[7 4 6];
b=[1 1];
subplot(2,1,1)
impulse(b,a);
subplot(2,1,2)
step(b,a)

3.卷积的求解

MATLAB中处理离散卷积的函数为：conv(f1,f2)，对序列做卷积运算。

处理连续信号的卷积时，需要对连续信号取相同的步长，结果再乘实际步长。

例题：

已知两个连续时间信号为：

$f1(t)=\left\{\begin{matrix} 2,0<t<1\\ 0,else\end{matrix}\right.$ $f2(t)= \left\{\begin{matrix} t,0<t<2\\ 0,else\end{matrix}\right.$

利用MATLAB画出f1(t)* f2(t)的时域波形图。

t11=0;t12=1;t21=0;t22=2;
t1=t11:0.001:t12;
ft1=2*rectpuls(t1-0.5,1); %幅值为1,宽度为1,t1=0.5点对称的方波信号
t2=t21:0.001:t22
ft2=t2;
t3=t11+t21:0.001:t12+t22;
ft3=conv(ft1,ft2);
ft3=ft3*0.001;
plot(t3,ft3)
title('ft1(t)*ft2(t) ')

chapter3

1.求解单位脉冲响应

MATLAB提供了专门用于求LTI离散系统的单位脉冲响应的函数：

h=impz(b,a,k)

其中a=[a0 ,a1 ,a2 ,…, an ]，b=[b0 ,b1 ,b2 ,…, bn ]分别是差分方程左右的系数向量，k表示输出序列的取值范围，h 就是系统的单位脉冲响应。

例题：

求离散系统的单位脉冲响应。

$y(k)+3y(k-1)+2y(k-2)=f(k)$

k=0:10;
a=[1 3 2]; %构造方程左右系数向量
b=[1];
h=impz(b,a,k); %求单位脉冲响应并绘图
stem(k,h,'.')

2.卷积和的求解（离散）

MATLAB中用于计算离散序列卷积的函数为： conv(x1,x2)

例：求以下两个离散序列的卷积。

$x_{1}(k)=sin(k),0\leq k\leq 10,x_{2}(k)=0.8^{k},0\leq k\leq 15$

k1=0:10; %x1的变量取值范围
x1=sin(k1); %构建x1序列
k2=0:15; %x2的变量取值范围
x2=0.8.^k2; %构建x2序列
y=conv(x1,x2);
subplot(3,1,1)
stem(k1,x1)
subplot(3,1,2)
stem(k2,x2)
subplot(3,1,3)
stem(y)

chapter4

1.求解系统响应

例题：

已知系统的频率响应函数和输入分别为：

$H(j\omega )=\frac{1-j\omega }{1+j\omega },f(t)=sin(t)+sin(3t)$

试用Matlab求解系统的响应。

t=0:pi/100:4*pi;
b=[-1,1]; %频率响应函数的分子系数
a=[1 1]; %频率响应函数的分母系数
ft=sin(t)+sin(3*t); % 激励
yt=lsim(b,a,ft,t); %求解响应
plot(t,ft,'r')
hold on
plot(t,yt,'b')

2.Sa信号的采样和恢复

例：信号Sa(t)作为被采样信号，信号带宽B=1, 采样频率ωs=2B，此频率下的采样为Nyquist采样，对采样及恢复过程用Matlab进行仿真。


B=1; %信号带宽
wc=B; %滤波器截止频率
Ts=pi/B; %采样间隔
ws=2*pi/Ts %采样角频率
N=100; %滤波器时域采样点数
n=-N:N;
nTs=n.*Ts; %采样数据的采样时间
fs=sinc(nTs/pi); %函数的采样点
Dt=0.005; %恢复信号的采样间隔
t=-15:Dt:15; %恢复信号的范围
fa=fs*Ts*wc/pi*sinc((wc/pi)*(ones(length(nTs),1)*t-nTs'*ones(1,length(t))));
%信号重构
error=abs(fa-sinc(t/pi)); %求重构信号与原信号的归一化误差
subplot(2,1,1)
plot(t,fa)
title('由 sa(t)=sinc(t/pi)的 Nyquist抽样信号重构 sa(t)')
subplot(2,1,2)
plot(t,error)
title('Nyquist抽样信号与原信号的误差 error(t)')

从误差图形中可以看到，用采样信号恢复的信号与原始信号的误差是极小的，而且这种误差来自于计算过程的数值计算误差。

表明:对一个信号进行Nyquist采样后，能够无误差地从采样信号中恢复出原信号。

------------------------------------

以下的部分在：MOOC西电《工程信号与系统》

1.拉普拉斯变换的求解

MATLAB内置了符号函数laplace 和ilaplace函数分别对应拉普拉斯变换与反变换。

例: 计算 $f1(t)=e^{-t}\varepsilon (t),f2(t)=te^{-0.5t}\varepsilon (t)$ 的卷积。

这里先求拉普拉斯变换，然后时域卷积等于频域乘积。最后再反变换得到时域的卷积。

syms t;t=sym('t','positive'); %t定义为正的“符号”变量
fs1=laplace(exp(-t)); %f1(t)的Laplace变换
fs2=laplace(t*exp(-t/2)); %f2(t)的Laplace变换
yt=simple(ilaplace(fs1*fs2)); %利用Laplace反变换求时域解
fs1,fs2,yt

运行结果：

fs1 =1/(s + 1)fs2 =1/(s + 1/2)^2yt =4*exp(-t) - 4*exp(-t/2) + 2*t*exp(-t/2)

2.Matlab绘制零极点图、判断稳定

例题：

利用MATLAB画出系统 $H(s)=\frac{s^{2}+4s+3}{s^{4}+3s^{3}+4s^{2}+6s+4}$ 的零极点图，并判断系统的稳定性。

b=[1 4 3]; %分子系数，按降幂顺序排列
a=[1 3 4 6 4]; %分母系数，按降幂顺序排列
sys=tf(b,a)
pzmap(sys);sgrid;
azp=roots(a); %求出极点azp，在左半平面即为稳定
%根据参量wd的值判断稳定：1表示稳定，0表示不稳定
wd=1;
for k=1:length(azp)if real(azp(k))>-0.000001wd=0;  end
if wd==0title('不稳定系统');
elseif wd==1title('稳定系统');
end
end

运行结果：

sys =s^2 + 4 s + 3-----------------------------s^4 + 3 s^3 + 4 s^2 + 6 s + 4Continuous-time transfer function.

系统的零极图，如图所示，可见系统有4个极点， 2个零点，其中在虚轴上有一对共轭极点，故该系统是不稳定的。

（知识回顾：系统稳定需要H(s)极点全部都在左半开平面）

3.求频率响应函数，判断稳定

例题：

已知系统函数 $H(s)=\frac{1}{s^{3}+2s^{2}+3s+1}$ ,画出其零极点分布，求系统的单位冲激响应h(t) ，和频率响应H(jw) ，并判断系统是否稳定。

num=[1];
den=[1 2 3 1];
sys=tf(num,den);
poles=roots(den);
figure(1);
pzmap(sys);
t=0:0.02:10;
h=impulse(num,den,t);
figure(2);plot(t,h);
title('Impulse Respone')
[H,w]=freqs(num,den);
figure(3);plot(w,abs(H));

poles=-0.4302 -0.7849 +1.3071i -0.7849 -1.3071i

极点都在左半开平面，系统稳定。