一阶与二阶连续时间系统

\(\newcommand{\dif}{\mathop{}\!\mathrm{d}} \newcommand{\xleftrightarrow}[1]{\stackrel{#1}{\longleftrightarrow}} \newcommand{\F}{\mathcal{F}} \newcommand{\ft}{\xleftrightarrow{\F}}\)

一阶连续时间系统

Bode图

一阶系统的微分方程往往表示为:

\[\tau \frac{\dif y(t)}{\dif t}+y(t)=x(t)\]

(注:尽管更一般的方程应该是 $a \frac{\dif y(t)}{\dif t}+ b y(t)=c x(t)$,但由于 $b,c$ 只是线性增大了输入/输出的幅度,对探讨系统性质没影响,所以我们用上面的形式来讨论。)

其频率响应为:

\[H(j\omega)=\frac{1}{1+j\omega \tau}\]

我们先考虑幅频特性。对其取对数,并在不同区间取近似:

\[\begin{align} 20\lg |H(j\omega)|&=-10\lg [(\omega \tau)^2+1]\\ &\approx \begin{cases} -10 \lg 1=0 & \omega \ll 1/\tau\\ -10\lg 2=-3\;\mathrm{dB} & \omega=1/\tau\\ -10 \lg [(\omega \tau)^2]=-20 \lg (\omega \tau) & \omega \gg 1/\tau \end{cases} \end{align}\]

类似地,对相频特性也进行近似处理(不是很清楚第二条近似的数学推导,知道的同学在评论区写写):

\[\begin{align} \angle H(j\omega) &= -\arctan(\omega\tau)\\ &\approx \begin{cases} 0 & \omega \ll 0.1/\tau\\ -(\pi/4)[\lg(\omega\tau)+1] & 0.1/\tau \leq \omega\leq 10/\tau\\ -\pi/2 & \omega \geq 10/\tau \end{cases} \end{align}\]

作出 Bode图:

我们能从图中读出以下几点:

  • 幅频曲线:
    • 频率较低时,其增益为 0
    • 频率较高时,按照每 10 倍频 20 dB 衰减
    • 在 $1/\tau$ 处,两条近似直线相交,实际衰减幅度为 $-10\lg 2=-3\;\mathrm{dB}$。故该点又称为:转折频率(break frequency)3dB点
  • 相频曲线:
    • 在转折频率上下十倍频范围内,$\angle H(j\omega)\sim\lg(\omega\tau)$ 图像近似线性下降
MATLAB代码
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
%matlab代码
tau=0.1;
omega=0:0.01:1000;
H=1./(1+i*omega*tau);

%幅频特性
subplot(2,1,1)
plot(log10(omega), 20*log10(abs(H)),'r');
hold on;
xticks(log10([0.01/tau, 0.1/tau 1/tau 10/tau 100/tau]))
xticklabels({'0.01/\tau','0.1/\tau','1/\tau','10/\tau','100/\tau'})
xlabel({'$\omega$'},'Interpreter','latex');
ylabel({'$20 \log_{10} (|H(j\omega)|)$'},'Interpreter','latex');

plot(log10(omega), 0.*omega,'b--');
plot(log10(omega), -20*log10(omega)-20*log10(tau),'b--');

plot([log10(1/tau),log10(1/tau)],[-40, 0],'b:')

axis([log10(omega(1)), log10(max(omega)),-40,5])

hold off

%相频特性
subplot(2,1,2)
plot(log10(omega), angle(H),'r');
hold on
xticks(log10([0.01/tau, 0.1/tau 1/tau 10/tau 100/tau]))
xticklabels({'0.01/\tau','0.1/\tau','1/\tau','10/\tau','100/\tau'})
yticks([-3*pi/4 -pi/2 -pi/4 0 pi/4 pi/2 3*pi/4])
yticklabels({'-3\pi/4','-\pi/2','-\pi/4','0','\pi/4','\pi/2','3\pi/4'})
xlabel({'$\omega$'},'Interpreter','latex');
ylabel({'$\angle H(j\omega)$'},'Interpreter','latex');

plot(log10(omega), 0.*omega,'b--');
plot(log10(omega), -(pi/4).*(log10(omega.*tau)+1),'b--');
plot(log10(omega), 0.*omega-pi/2,'b--')

plot([log10(0.1/tau),log10(0.1/tau)],[-3*pi/4, 0],'b:')
plot([log10(1/tau),log10(1/tau)],[-3*pi/4, -pi/4],'b:')
plot([log10(10/tau),log10(10/tau)],[-3*pi/4, -pi/2],'b:')

axis([log10(omega(1)), log10(max(omega)),-3*pi/4,pi/4])

hold off

冲激响应和阶跃响应

对应的单位冲激响应为:

\[H(j\omega)=\frac{1}{1+j\omega \tau}\\ h(t)=\frac{1}{\tau}e^{-t/\tau} u(t)\]

对应的阶跃响应为:

\[s(t)=h(t)*u(t)=[1-e^{-t/\tau}]u(t)\]

可以看出,$\tau$ 决定着冲激衰减与阶跃上升的速度,所以 $\tau$ 又称为时间常数。

最后,我们来看看不同频率的正弦信号经过一阶系统后的图像,加深对一阶系统的理解。从下面三个图可以看出,频率越高,信号输出的衰减幅度越大,同时相位差也越大。

输入 $\sin(t)$ $\sin(2t)$ $\sin(4t)$
图像 1 2 3

二阶连续时间系统

冲激响应与阶跃响应

二阶系统的微分方程可表示为:

\[\frac{\dif^2 y(t)}{\dif t^2}+2 \zeta\omega_n \frac{\dif y(t)}{\dif t}+\omega_n^2 y(t)=\omega_n^2 x(t)\]

频率响应为:

\[H(j\omega)=\frac{\omega_n^2}{(j\omega)^2+2\zeta \omega_n (j\omega)+\omega_n^2}\]

我们对分母进行多项式分解:

\[H(j\omega)=\frac{\omega_n^2}{(j\omega-c_1)(j\omega-c_2)}\\ \begin{cases} c_1=\zeta\omega_n+\omega_n\sqrt{\zeta^2-1}\\ c_2=\zeta\omega_n-\omega_n\sqrt{\zeta^2-1} \end{cases}\]

对 $\sqrt{\zeta-1}$ 分情况讨论:

  $\zeta \neq 1$ $\zeta=1$
$c_1,c_2$ $c_1 \neq c_2$ $c_1=c_2=-\omega$
频率响应 \(H(j\omega)=\frac{M}{j\omega-c_1}-\frac{M}{j\omega-c_2}\\M=\frac{\omega_n}{2\sqrt{\zeta^2-1}}\) \(H(j\omega)=\frac{\omega_n^2}{(j\omega+\omega_n)^2}\)
冲激响应 $h(t)=M[e^{c_1t}-e^{c_2t}]u(t)$ $h(t)=\omega_n^2 t e^{-\omega_n t}u(t)$

进一步将 $\zeta \neq 1$ 划分。如果 $0<\zeta<1$,那么冲激响应为:

\[\begin{align} h(t) &= \frac{\omega_n e^{-\zeta\omega_n t}}{2 j \sqrt{1-\zeta^2}}\left\{ \exp\left[ j(\omega_n\sqrt{1-\zeta^2})t \right]-\exp\left[ -j(\omega_n\sqrt{1-\zeta^2})t \right]\right\}u(t)\\ &= \frac{\omega_n e^{-\zeta\omega_n t}}{j \sqrt{1-\zeta^2}} \sin(\omega_n t\sqrt{1-\zeta^2})u(t) \end{align}\]

从上式可以看出,这种情况下的冲激响应会有振荡(式中 $\sin$ 部分),当其幅度会衰减(式中 $e^{-\zeta\omega_n t}$ 部分)

如果 $\zeta>1$,这时是两个指数信号相减,不会有振荡。

综上,我们作出不同 $\zeta$ 值的图像(进行了适当缩放,y 轴是 $h(t)/\omega_n$):

图中清楚地说明了 $\zeta$ 的影响,所以称 $\zeta$ 为 阻尼系数(damping ratio)。根据 $\zeta$ 的取值,将系统分为:

  • $0<\zeta<1$,欠阻尼(underdamped)
  • $\zeta=1$,临界阻尼(critically damped)
  • $\zeta>1$,过阻尼(overdamped)

$\omega_n$ 称为 无阻尼自然频率,$\omega_n$ 与 $\zeta$ 一同控制着时间尺度与幅度。$\omega_n$ 越大,则时间上越压缩,幅度上越大。在 $\zeta=0$ 时,振荡频率为 $\omega_n$。

最后补充一下,$\zeta<0$ 的情况与上面类似,也是分成三类振荡情况,只不过其幅度是增加而不是减小。

MATLAB代码
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
function secondOrderSystem(zeta,omegan)
    t=0:0.1/omegan:10/omegan;
    if zeta~=1
        c1=-zeta*omegan+omegan*sqrt(zeta^2-1);
        c2=-zeta*omegan-omegan*sqrt(zeta^2-1);
        h=omegan/(2*sqrt(zeta^2-1))*(exp(c1.*t)-exp(c2.*t));
    else
        h=omegan^2.*t.*exp(-omegan.*t);
    end
    plot(t.*omegan,h./omegan,'linewidth',2);
    %text=['\zeta=', num2str(zeta)];
    %legend(text);
end

hold on;
omegan=1;
secondOrderSystem(0.1,omegan);
secondOrderSystem(0.5,omegan);
secondOrderSystem(1,omegan);
secondOrderSystem(1.5,omegan);

xticks([1:10]/omegan)
xticklabels({'1/\omega_n','2/\omega_n'})
xlabel({'$t$'},'Interpreter','latex');
ylabel({'$h(t)/\omega_n$'},'Interpreter','latex');
xyplot %网上找的function,用于将 x 轴移到中间
legend('\zeta=0.1','\zeta=0.5','\zeta=1','\zeta=1.5')

对于阶跃响应,其特点与上面类似,也是分成:

$\zeta\neq1$ $$ s(t)=\left\{ 1+M\left[ \frac{e^{c_1t}}{c_1}-\frac{e^{c_2t}}{c_1} \right]\right\}u(t) $$
$\zeta=1$ $$ s(t)=\left[ 1-e^{-\omega_n t}-\omega_n t e^{-\omega_n t} \right] u(t) $$

其图像为(懒得再用matlab了):

可以看到,欠阻尼情况下,阶跃响应有 超量(overshoot),即阶跃响应超过它的终止。

Bode图

频率响应为:

\[\begin{align} H(j\omega)&=\frac{\omega_n^2}{(j\omega)^2+2\zeta \omega_n (j\omega)+\omega_n^2}\\ &=\frac{1}{(j\omega/\omega_n)^2+2\zeta (j\omega/\omega_n)+1}\\ &=\frac{1}{1-(\omega/\omega_n)^2+2\zeta (j\omega/\omega_n)} \end{align}\]

幅频特性:

\[20\lg |H(j\omega)|=-10\lg \left\{ \left[ 1- \left( \frac{\omega}{\omega_n}\right)^2 \right]^2+4\zeta^2 \left( \frac{\omega}{\omega_n} \right)^2 \right\}\]
  • 若 $\omega\ll\omega_n$,则 $\omega/\omega_n\approx0$,$20\lg H(j\omega) =0$
  • 若 $\omega\gg\omega_n$,则 $20\lg H(j\omega) \approx -10 \lg \left( \dfrac{\omega}{\omega_n} \right)^4=-40\lg \omega+40\lg \omega_n$

作出图像:

从图像中可以看出:

  • 低频时,$20\lg\vert H(j\omega) \vert\approx 0 \rm{dB}$
  • 高频时,以每十倍频 -40dB 速度下降
  • 低频与高频交于 $\omega=\omega_n$,故 $\omega_n$ 称为二阶系统的转折频率
  • 特殊地,注意到当 $\zeta$ 较小时,$\omega=\omega_n$ 处会有尖峰,我们用 品质因数 Q 来衡量峰值的尖锐程度。对于该二阶系统,$Q=1/(2\zeta)$
MATLAB代码
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
function secondOrderSystem(zeta,omegan)
    range=[-2:0.1:3]
    t=omegan.*10.^range;
    H=omegan^2./((i*omega).^2+2*zeta.*omegan.*i.*omega+omegan.^2)
    
    plot(range,20.*log10(abs(H)),'linewidth',1);
    %text=['\zeta=', num2str(zeta)];
    %legend(text);
end

hold on;
omegan=1;
secondOrderSystem(0.1,omegan);
secondOrderSystem(0.5,omegan);
secondOrderSystem(1,omegan);
secondOrderSystem(1.5,omegan);

xticks([-2:2])
xticklabels({'0.01\omega_n','0.1\omega_n','1\omega_n','10\omega_n','100\omega_n'})
xlabel({'$\omega$'},'Interpreter','latex');
ylabel({'$20\lg |H(j\omega|$'},'Interpreter','latex');
legend('\zeta=0.1','\zeta=0.5','\zeta=1','\zeta=1.5')