Octave(或Matlab)分析RLC二阶电路系统（万字全面解析：包含时域、频域、S域）

Octave介绍

\hspace{2em} 学习经典控制理论或者信号与系统时，会涉及到很多微分方程的计算分析，有时是对机械系统建模分析、有时是对电路系统建模分析、甚至是机电一体的系统进行建模分析。 当系统建立微分方程的阶数过高时，手工计算十分困难，因此掌握一款数值分析工具是必不可少的。相信只要是理工科类毕业的工程师或者在校学生，一定听过大名鼎鼎的Matlab，其代码编写方式简单，分析理工科类相关的工程问题时，分析工具齐全，深受广大用户欢迎。但是，其价格昂贵，体积臃肿，让许多初学者不能快速上手练习。因此对于初学者，作者向各位推荐学习Octave，这是GNU项目下的一款开源软件、内存不到1GB、支持跨平台运行，其编码语法格式几乎完全适配Matlab，许多代码只需稍加修改甚至无需修改，即可在Matlab下运行。所以，Octave也是一款优秀的数值计算和科学工程分析工具，完全可以胜任一些数值计算的代码编写。那是不是Octave就可以取代Matlab的位置了呢？了解Matlab的读者可能知道，Matlab强大的Simulink仿真建模工具，才是难以替代的功能，Octave目前也只是应用于数值分析的脚本代码编写。但是，其开源和内存小的优点，绝对是初学者为后续学习Matlab或者替代Matlab数值分析的优秀工具。关于Octave的下载安装方式非常的简单，这里贴出网址，各位自行下载安装：(下载只需要注意平台的选择即可，Octave支持Linux、MacOS、Windows平台运行)
Octave下载
\hspace{2em} 安装完毕后，点击左上角文件–>新建–>新建脚本，创建M文件以后，读者可以先将以下代码复制粘贴进代码编辑区，点击运行观察能否成功得出正弦信号波形，运行成功说明安装完成。

clc
clear
x=0:0.01:100;%设置变量范围-->起点:步进:终点
y=sin(x);
plot(x,y);%绘图

RLC电路的分析方法

时域分析

\hspace{2em} 如图1所示，为RLC二阶电路系统的示意图。研究RLC电路对实际开发工作有什么重要指示？这样典型的RLC电路可能在实际电路原理图中不容易看到，但是这样典型的二阶系统却隐藏在PCB的寄生参数中、RC低通滤波器、通讯线缆传输过程等等一系列场合，都能找到这样简化的二阶系统，这些场合经常会观察到RLC组合产生的过冲、振铃现象。有读者可能认为过冲、振铃现象是阻抗不匹配，信号反射导致的，需要端接电阻，用于阻抗匹配解决过冲现象。然后，可能就会有工程师记住这样的结论，遇到过冲就开始调节电阻，发现有些场合下，把电阻加大就没有过冲现象了，持续加大也没有过冲，也就是反射现象消失了？可是根据反射的定义，阻抗不匹配就会带来反射，那么电阻持续加大超过了匹配值，也不会过冲振荡了，这该如何解释呢？说明在该场合下反射的影响是很小的，其过冲现象是由RLC的阻尼状态占主导。其实界定两者的主导情况，通常需要判断电路系统的电尺寸和波长之间的关系，对于电大尺寸下高频情况的电路是分布式参数模型，电路理论不再适用，反射就会占主导。但是在电小尺寸的低频情况的电路是集总式参数模型，也就是像RLC这样的电路模型占主导。在《信号完整性与电源完整性分析》（作者Eric Bogatin）一书中，该书作者提出对于传输线反射问题，也可以构建n阶集总式参数模型进行分析，但是分析难度也会变高，通常要借助仿真工具分析，但是仿真中无法直接观察到表达式之间的关系。因此需要一个具有代表性的阶数，其分析精确度可接受，分析难度较低的模型，便于我们理解信号在经过这些无源器件时发生了什么，RLC二阶电路系统就是这样的一个典型。（学习过微波技术的读者可能知道，简单传输线理论模型就是一个典型的二阶系统，只是比RLC电路多了一个电导，用于表征位移电流带来的介质损耗）

图1 RLC二阶电路系统

\hspace{2em} 接下来，我们进入正题，分析RLC电路的时域响应。为了观察RLC的动态性能，假设系统初始条件为0，我们对RLC系统施加一个单位阶跃信号$u(t)$，然后求解其阶跃响应$y(t)$，根据KVL列写微分方程如下：
$$i(t)R+L\frac{di(t)}{dt}+\frac{1}{C}\int i(t)dt=u(t)$$
又由于$i(t)=C\frac{dy(t)}{dt}$，初始条件为0，微分方程可以写成如下形式:
$$\frac{dy(t)}{dt}RC+LC\frac{d{y}^2(t)}{d{t}^2}+y(t)=u(t)$$
这是一个非齐次二阶微分方程，要求通解，则需要求解齐次解和特解。在信号与系统的学习中，我们知道通解就是系统的全响应，齐次解为瞬态响应，特解为稳态响应，两者相加就是通解。当$t>0$时，输入$u(t)=1$为常数，而多项式左边的输出应保持相同的阶数，可以设特解为常数$y_p(t)=A$，带回微分方程可解得$A=1$。接下来重点是求齐次解，由齐次二阶微分方程可以列写特征方程如下：
$${\lambda }^{2}LC+\lambda RC+1=0$$
根据判别式$\Delta =\sqrt{(RC{)}^2-4LC}$确定解的形式：
$$\Delta >0，存在不同实数根{\lambda }_1、{\lambda }_2:y_i(t)=C_1{e}^{{\lambda }_{1}t}+C_2{e}^{{\lambda }_{2}t}$$
$$\Delta =0，存在两个相同实数根，二重根\lambda :y_i(t)=(C_1+C_{2}t)e^{\lambda t}$$
$$\Delta <0，存在一对共轭复数根，\sigma \pm j{\omega }_d:y_i(t)=e^{\sigma t}(C_{1}cos{\omega }_{d}t+C_{2}sin{\omega }_{d}t)$$
所以，阶跃响应$y(t)=y_p(t)+y_i(t)$，根据阶跃响应的结果可以判断，RLC过冲的秘密就藏在共轭复数中，因为解的形式中存在三角函数这样周期波动的信号。为了得到最终的$y(t)$表达式，还需要求出特征方程的根表达式。无论实数还是复数，求根公式依旧适用，整理特征根的表达式如下：
$${\lambda }_1、{\lambda }_2=-\frac{R}{2L}\pm \sqrt{(\frac{R}{2L}{)}^2-\frac{1}{LC}}$$
从根的表达式中可以看到，还可以从根号的定义域来判别不同实根、二重根、复数根。由此我们做出这样的定义：（无阻尼频率在后文频域再详细分析来源）
$$阻尼系数：\alpha =\frac{R}{2L}无阻尼角频率：{\omega }_n=\frac{1}{\sqrt{LC}}阻尼比：\xi =\frac{\alpha }{{\omega }_n}=\frac{R}{2}\sqrt{\frac{C}{L}}$$
如果是复数根，无阻尼频率与复数根的实部、虚部还存在如下关系：(后文在S域再详细分析来源)
$$实部：\sigma =-\frac{R}{2L}虚部（阻尼角频率）：{\omega }_d=\sqrt{{\omega }_n^2-{\sigma }^2}$$

所以，可以使用阻尼比这个参数对系统的瞬态响应特性做判别：
$$\xi >1：过阻尼，系统存在不同实数根，无振荡、缓慢衰减$$
$$\xi =1：临界阻尼，系统存在二重根，无振荡、快速衰减$$
$$\xi <1：欠阻尼，系统存在共轭复根，振荡、衰减$$
\hspace{2em} 如果想将图1中的器件参数用于计算阶跃响应，可以先快速验证解的形式。通过计算阻尼比发现$\xi =0.158<1$，存在共轭复根，到这里我们就能初步预估这个系统的瞬态响应会出现过冲振荡衰减的现象。不过，要求出具体表达式并绘制波形，可以预估其计算量并不小，这时就能体现数值分析工具的优势了，我们将利用Octave编写程序求解这组RLC参数的阶跃响应并绘制波形。不过，在编写程序之前，我们还得搞懂齐次解形式里面的常数如何计算，因为程序可以帮助我们解决数值计算的重复繁琐问题，但是算法还是需要工程师提供。当然也有数值解的计算方法，数值解的好处就是在于无需理论推导，只需构建微分方程，直接可以得到结果，但是缺点就是没有表达式，无法定量分析洞察系统的性能，只能通过大量数值定性分析推演系统性能（等同于仿真软件得出的结果）。因此，我将会编写两版程序，一版使用理论解析编写程序，一版采用ode45微分方程求解工具，求出数值解对比。
\hspace{2em} 求解常数的基本思路，就是依据初始条件进行求解，由于是二阶微分方程，解的形式存在两个未知常数，那就需要构建两个方程进行求解。根据初始条件为0，可构建如下方程：
$$y(0^{+})=0y^{\prime }(0^{+})=0$$
求导对于不同实根、二重根的计算可能不算繁琐，但是对于复数根解的形式进行求导，其步骤比较繁琐，有指数和三角函数，看上去就不是那么容易处理。虽然一个训练有素的工程师对于微积分的运算是信手拈来的，但是我们希望更多时间花费在设计和分析上，计算的工作应该丢给数值分析工具。所以，可以使用Octave提供的符号工具箱进行求导计算（Matlab同样适用），帮助我们把常数的约束方程确定下来，以下代码以复数根为例子，展示求导并得到约束方程的程序设计:

clc
clear
pkg load symbolic;%导入符号工具箱
syms t C1 C2 a w;%定义符号变量
y_p=1;%特解
y_t = exp(a*t)*(C1*cos(w*t)+C2*sin(w*t))+y_p;%定义函数
dy_dt=diff(y_t,t);%y(t)对t求导
fprintf("dy_dt=:\r\n");
disp(dy_dt);%显示结果
%---------带入初始条件，求得未知常数的约束方程-------
y_t_0=subs(y_t,t,0)==0;%初始条件y(0)=0带入y(t)方程
dy_dt_0=subs(dy_dt,t,0)==0;%初始条件y'(0)=0带入y'(t)方程
fprintf("y(0)=0:\r\n");
disp(y_t_0);
fprintf("y'(0)=0:\r\n")
disp(dy_dt_0);

图2 程序运行结果

通过程序计算，我们可以得到不同解形式下的未知常数约束方程组：
$$不同实根：\{\begin{array}{l}{C}_1+C_2+1=0\\ \\ C_1{\lambda }_1+C_2{\lambda }_2=0\end{array}$$
$$二重根：\{\begin{array}{l}{C}_1+1=0\\ \\ C_1\lambda +C_2=0\end{array}$$
$$复数根：\{\begin{array}{l}{C}_1+1=0\\ \\ C_1\sigma +C_2{\omega }_d=0\end{array}$$
所有的方程已经推导完毕，剩下的计算就交给计算机处理，RLC求解阶跃响应程序如下：

clc
clear
R=10;L=1e-3;C=1e-6;%定义器件参数
t=0:1e-6:2*1e-3;%定义计算范围
y_p=1;%特解
figure
hold on %开启图像叠加
plot(t,t-t+1);%绘制t>0的单位阶跃信号，用于对比
for i = 1:1:5  %用于观察电阻变化，阶跃响应的变化情况
Alpha=R./(2*L);%阻尼系数
Omega_0=1./(sqrt(L*C));%无阻尼频率
xi=Alpha./Omega_0;%阻尼比
x1=-Alpha+sqrt(Alpha.^2-Omega_0.^2);x2=-Alpha-sqrt(Alpha.^2-Omega_0.^2);%实数根
Sigma=-Alpha;Omega_d=sqrt(Omega_0.^2-Sigma.^2);%复数根的实部和虚部
if xi > 1 %不同实根
  A=[1,1;x1,x2];B=[-1;0];%根据未知实数的约束方程组，可以提取系数矩阵和常数矩阵
  Cn=A\B;%矩阵左除运算求解方程组
  C1=Cn(1);C2=Cn(2);%提取结果
  y_t=C1*exp(x1*t)+C2*exp(x2*t)+y_p; %得到阶跃响应
  elseif xi == 0 %二重根
    A=[1,0;x1,1];B=[-1;0];%根据未知实数的约束方程组，可以提取系数矩阵和常数矩阵
    Cn=A\B;%矩阵左除运算求解方程组
    C1=Cn(1);C2=Cn(2);%提取结果
    y_t=(C1+C2*t).*exp(x1*t)+y_p; %得到阶跃响应
  else    %复数根
    A=[1,0;Sigma,Omega_d];B=[-1;0];%根据未知实数的约束方程组，可以提取系数矩阵和常数矩阵
    Cn=A\B;%矩阵左除运算求解方程组
    C1=Cn(1);C2=Cn(2);%提取结果
    y_t=exp(Sigma*t).*(C1*cos(Omega_d*t)+C2*sin(Omega_d*t))+y_p; %得到阶跃响应
endif
plot(t,y_t);%绘制阶跃响应曲线
R=R+20;%每次循环电阻值增大20，重新绘制图像
fprintf("...
         Alpha=%.4f\r\n...
         Omega_0=%.4f\n...
         xi=%.4f\r\n...
         Sigma=%.4f\r\n...
         Omega_d=%.4f\r\n...
         C1=%.4f\r\n...
         C2=%.4f\r\n",...
         Alpha,Omega_0,xi,Sigma,Omega_d,C1,C2);%打印计算过程中的变量
         disp('x1=');disp(x1);disp('x2=');disp(x2);%打印特征根
endfor

图3 不同阻尼状态下的RLC阶跃响应程序运行结果（解析解）

\hspace{2em} 接下来，我们再简单讨论一下数值解的计算方法。这里主要是想让读者掌握ode45微分方程数值计算工具的使用方法。ode45可以直接用于求解一阶微分方程，但是RLC系统是二阶微分方程，因此需要进行降阶处理，将非齐次二阶微分方程改写成两个一阶微分方程的组合即可：（换元法）
$$\{\begin{array}{l}\frac{dy(t)}{dt}=Y(t)\\ \\ \frac{dY(t)}{dt}=\frac{u(t)}{LC}-\frac{R}{L}Y(t)-\frac{y(t)}{LC}\end{array}$$
因此，代码里将上述两个方程写成匿名函数，在ode45求解器中调用即可：

clc
clear
R=10;L=1e-3;C=1e-6;%定义器件参数
u_t=1;%t>0时的单位阶跃信号
%y(1)表示y(t),y(2)表示Y(t)
RLC_ODE=@(t,y)[y(2);u_t/(L*C)-R*y(2)/L - y(1)/(L*C)];%列写微分方程匿名函数
y0=[0;0];%给定初始条件，y(0)=0和Y(0)=0
t_span=[0,2*1e-3];%定义计算范围
ode_opset=odeset('RelTol',1e-6,'AbsTol',1e-6);%设置微分方程求解精度
[t,y]=ode45(RLC_ODE,t_span,y0,ode_opset);
y_t=y(:,1);%求得通解
figure
hold
plot(t,t-t+1);%绘制单位阶跃信号用于对比
plot(t,y_t);%绘制阶跃响应曲线

图4 RLC阶跃响应程序运行结果（数值解）

其运行结果与理论分析的解析解一致，说明理论分析的过程和代码编写都是正确的。数值解代码编写简单，其优势就是只需要列出系统的微分方程，就可以直接算出结果，对于大规模复杂的系统，解析解求解困难时，数值解分析复杂系统灵敏度、稳定性等方面也大有用处。不过，如果可以求得解析解，很显然解析解对于洞察系统的各种行为和性能，会更有优势。
\hspace{2em} 图1的电路图是使用Psim仿真软件绘制的，我们不妨顺便运行一下仿真结果。如图5所示，其运行结果与解析解、数值解的结果一致，说明上述分析过程正确。

图5 RLC阶跃响应Psim仿真结果（R=10Ω）

频域分析

\hspace{2em} 有时我们还关心周期性的信号经过RLC会发生什么，而这种周期性的信号是非常常见的。例如在开关电源系统或者从市电中取电的设备中，都不可避免的存在开关信号谐波或者市电50Hz工频信号的谐波。如果将丰富的谐波都放在时域的时间轴上观察，势必非常的混乱，难以观察。因此，需要将时域变换到频域进行观察，以频率作为变量，可以直观的看到不同频率的信号经过系统将会发生什么。伟大的傅里叶先生就提供了这个转换的方法，称为傅里叶变换，用于研究系统的稳态响应。傅里叶级数（周期傅里叶变换）向世人揭示了，任意一个周期信号都可以分解为多个不同频率的正弦信号叠加（因此频域中唯一存在的信号就只有正弦信号），而稳态响应的定义就是信号激励对系统有着长期稳定的影响，正弦信号就具备这样的特点，所以傅里叶变换适用于求解周期信号的稳态响应。
\hspace{2em} 周期信号虽然可以分解为多个正弦信号，但经过系统后，又该如何计算？那就得保证系统是线性时不变的，线性时不变系统满足齐次叠加性，可以先单独计算分解后的正弦信号经过系统，然后再进行求和叠加，得到正弦信号响应或称为频率响应。接下来我们就进一步围绕RLC电路进行分析，求解RLC电路的频率响应，由于篇幅有限这里不加证明的给定RLC是线性时不变系统。基于上述分析，常规的分析方法就是依据时域的微分方程进行傅里叶变换转换到频域分析，进而求得激励和响应之间的频率关系，可见在频域中计算还会涉及到数学变换的计算，计算变得难以处理了。但是对于工程计算上，我们更加常用的是阻抗分析法，这里运用傅里叶变换公式的时域微分性质，简单推导一下电阻、电感、电容的阻抗表达式。
$$傅里叶变换时域微分公式：\frac{df(t)}{dt}\stackrel{傅里叶变换}{\to }j\omega F(\omega )$$
$$电阻：Z_R(t)=\frac{U(t)}{I(t)}=R\stackrel{傅里叶变换}{\to }{Z}_R(\omega )=\frac{U(\omega )}{I(\omega )}=R$$
$$电感：Z_L(t)=\frac{U(t)}{I(t)}=\frac{L\frac{dI(t)}{dt}}{I(t)}\stackrel{傅里叶变换}{\to }{Z}_L(\omega )=\frac{j\omega LI(\omega )}{I(\omega )}=j\omega L$$
$$电容：Z_C(t)=\frac{U(t)}{I(t)}=\frac{U(t)}{C\frac{dU(t)}{dt}}\stackrel{傅里叶变换}{\to }{Z}_C(\omega )=\frac{U(\omega )}{j\omega CU(\omega )}=\frac{1}{j\omega C}$$
有了阻抗表达式，那么要计算激励和响应之间的关系就变得非常简单了，只需使用阻抗的串联分压关系就可以计算频率响应：
$$U_o(\omega )=\frac{U_i(\omega )}{1+j\omega RC+(j\omega {)}^{2}LC}$$
为了更直观的分析频率响应的特性，1930年荷兰科学家伯德发明了Bode图。他将输出和输入作为一个比值，也就是频域下的传递函数$H(\omega )$，取其模值和相位放置在同一张图，考虑到频率范围变化会引起模值大幅度变化，便于绘制图像还需对模值取对数。(后文会解释$20lo{g}_{10}$的来源)
$$\frac{U_o(\omega )}{U_i(\omega )}=H(\omega )=\frac{1}{1+j\omega RC+(j\omega {)}^{2}LC}$$
$$对数模值:20lo{g}_{10}|H(\omega )|相位：\angle H(\omega )$$
接下来，使用Octave编写代码观察Bode图，Bode图里面蕴含了很多丰富的信息，这里可能无法面面俱到，只是尽量让读者能理解和看懂Bode图：

clc
clear
R=10;L=1e-3;C=1e-6;%定义器件参数
f=logspace(log10(1),log10(10^7),1000);%定义扫频范围--->起点：终点：点数
w=2*pi*f;%角频率
H_w=1./(1+j*w*R*C+(j*w).^2*L*C);%传递函数
H_abs=abs(H_w);%取模值
H_pha=angle(H_w)*180/pi;%取相位并转换为度
figure
ax=plotyy(f,20*log10(H_abs),f,H_pha);
set(ax(1),'XScale','log');%坐标轴取对数显示
ylabel(ax(1),'dB');%标注单位
set(ax(2),'XScale','log');
ylabel(ax(2),'°');
xlabel('Hz');
title('Bode');

图6 RLC电路Bode图程序运行结果

图7 Psim仿真Bode图运行结果

仿真结果与计算结果一致，说明分析过程和程序设计正确。纵坐标的左侧为传递函数的幅值，右侧为相位，纵坐标的颜色分别对应曲线的颜色。首先观察幅频响应曲线，整体看上去，在某频率点之前幅值基本不变，而随着输入信号频率的增大，传递函数的幅值在下降，这表明输出信号的幅值相比于输入信号越来越小。这点就体现了RLC有低通的特性，也就是可以用作低通滤波器，即低频信号可以通过，高频信号会被抑制。观察相频曲线，发现在某频率点之前都是0°，接近该点相位开始变化，之后快速翻转到-180°，说明信号超过一定频率后，经过RLC系统会发生相位翻转的现象。幅频曲线和相频曲线放在同一张图中，很大的原因是两者之间存在某些关系，例如方便观察穿越频率、相位裕度等指标，用于判断系统的一些性能，不过由于篇幅有限这些内容我计划在后续的文章中再进行展开讲解。
\hspace{2em} 本文关于RLC电路在频域的分析，重点关注峰值频率、谐振频率、截止频率，接下来我们先讨论一下峰值频率。Bode图的幅频响应曲线存在一个尖峰点，这个点对应的频率就是RLC的输出电压峰值频率点，如果在频域中存在这个尖峰，就意味着在时域表现为输出电压在这个频率下会高于输入电压，如果输入信号是阶跃信号，阶跃响应会出现明显过冲。可能有不少读者认为峰值角频率就是前文提到的无阻尼角频率${\omega }_n=1/\sqrt{LC}$，换算成频率为$f_n=1/(2\pi \sqrt{LC})$，可事实真是如此吗？我们接下来进行简单的数学推导，尝试找到这个最大值点。要使得RLC传递函数对应最大值，则意味着传递函数的分母最小，将分母单独命为一个函数并取模，令模函数的导数为0 ，找到极值点，这里不加证明的给出极值点是极小值点且为最小值点。推导过程如下：
$$D(\omega )=-{\omega }^{2}LC+j\omega RC+1\underset{等式两边取平方}{\overset{取模}{\to }}|D(\omega ){|}^2=(1-{\omega }^{2}LC{)}^2+(\omega RC{)}^2$$
$$上式\underset{令导数｜D(\omega ){｜}^{\prime }=0}{\overset{等式两边同时求导}{\to }}整理得:(RC{)}^2=2LC(1-{\omega }^{2}LC)$$
$$化简求得：\omega =\sqrt{\frac{1}{LC}-\frac{R^2}{2L^2}}={\omega }_0(定义为峰值角频率)$$
$$RLC输出电压峰值频率公式：f_0=\frac{\sqrt{\frac{1}{LC}-\frac{R^2}{2L^2}}}{2\pi }$$
$$假设为无阻尼情况R=0，则有f_0=f_n，因此f_n为输出电压无阻尼频率$$
$$假设阻尼过大，使f_0表达式分子根号内出现负值，表明电压无峰值频率$$
\hspace{2em} 因此，用$f_n$表示输出电压的峰值频率并不完整。我们验证一下推导的结论是否正确，验证方法依旧是使用Octave的数值解验证解析解，同时使用Psim仿真标定幅频特性曲线的最大值点对应的频率值，是否和推导的一致：

clc
clear
R=10;L=1e-3;C=1e-6;%定义器件参数
f=logspace(log10(1),log10(10^7),1e6);%点数选1e6提高运算精度
H=@(w)1./(1+j*w*R*C+(j*w).^2*L*C);%匿名函数定义传递函数方便调用
w=2*pi*f;%角频率
H_abs=abs(H(w));%取模值
H_pha=angle(H(w))*180/pi;%取相位并转换为度
w_n=1./(sqrt(L*C));%无阻尼角频率
w_0=sqrt(w_n^2-R^2/(2*L^2));%峰值角频率
f_n=w_n/(2*pi);
f_0=w_0/(2*pi);
[H_max,H_mx]=max(20*log10(H_abs));%使用max函数提取最大值
fprintf("f_x=%.1f  H_max=%.3f\r\n",f(H_mx),H_max);%最大值数值解
fprintf("f_0=%.1f  H_0=%.3f\r\n",f_0,20*log10(abs(H(w_0))));%峰值点解析解
fprintf("f_n=%.1f  H_n=%.3f\r\n",f_n,20*log10(abs(H(w_n))));%无阻尼频率

图8 峰值频率和无阻尼频率计算程序运行结果

图9 Psim仿真标定最大值运行结果
程序运行结果中的 $f_x$表示数值解，其计算结果与理论解析解 $f_0$完全相等，无阻尼频率 $f_n$与其存在几十 $Hz$的误差，同时Psim仿真的结果与理论计算的一致，说明上述分析过程正确。
\hspace{2em} 前文我为什么会认为读者可能会将 $f_n$等同于峰值频率，因为 $f_n$其实就是谐振频率，而RLC电路发生谐振时，有一个特点就是其回路电流达到最大，可能会想当然的认为电流最大，使得输出电压最大，经过峰值频率的计算可以得知并非如此。无论是RLC串联谐振还是并联谐振，谐振的本质就是系统中存在一个固有的频率，使得LC储能器件的能量达到一个周期性的相互交换，能够形成稳定的振荡。本文描述的是串联RLC，回路电流是相等，根据 $P(\omega )=I^2(\omega )Z(\omega )$，只要感抗和容抗相等，在同一时间尺度下，就可以表征两者交换的能量相等。可以推导谐振频率公式：
$$|I(\omega )|=\frac{|U_i(\omega )|}{|Z(\omega )|}=\frac{|U_i(\omega )|}{\sqrt{R^2+(\omega L-\frac{1}{\omega C}{)}^2}}$$
感抗和容抗相等时，阻抗最小，呈纯电阻，使得电流最大：
$$\omega L=\frac{1}{\omega C}\stackrel{整理得}{\to }\omega =\frac{1}{\sqrt{LC}}={\omega }_n(谐振角频率)$$
$$RLC谐振频率公式：f_n=\frac{1}{2\pi \sqrt{LC}}$$
\hspace{2em} 在设计滤波器时，最重要的指标之一就是截止频率，接下来我们讨论RLC截止频率的计算方法。前文提到由于频率变化范围广，需要使用对数进行压缩处理，方便计算。那么应该只需要将其取 $lo{g}_{10}|H(\omega )|$就可以达到压缩的目的了，为什么演变成 $20lo{g}_{10}|H(\omega )|$？这实际涉及到单位的统一，1928年分贝( $dB$)单位被正式定义用于表示和计算电话线路中信号的增益和损耗，以 $10lo{g}_{10}(|\frac{P_1}{P_2}|)$作为换算公式，而功率和电压的平方存在正比关系 $P\propto U^2$,因此统一到以 $dB$为单位表示电压，使用换算公式表示为 $20lo{g}_{10}(|\frac{U_1}{U_2}|)$，这就是Bode图中纵坐标用分贝单位表示的由来。
\hspace{2em} 截止频率的定义来源于半功率点，从广义的角度分析，半功率点是指输出功率相对于输入功率衰减一半的点（输入幅度保持不变），此时有 $P_i=2P_o$，注意我在功率和电压平方之间使用了正比关系，说明功率换算到电压后，两者可能并不是严格相等的，那么在半功率点处电压是否会存在 $U_o=\sqrt{2}{U}_i$关系式？什么意思呢？我们做一下严格的数学推导，看看转换过程发生了什么？
$$10lo{g}_{10}(\frac{P_o}{P_i})=10lo{g}_{10}(\frac{\frac{U_o^2}{Z_o}}{\frac{U_i^2}{Z_i}})=20lo{g}_{10}(\frac{U_o}{U_i})+10lo{g}_{10}(\frac{Z_i}{Z_o})$$
$$若Z_i=Z_o，则有等式成立，10lo{g}_{10}(\frac{P_o}{P_i})=20lo{g}_{10}(\frac{U_o}{U_i})$$
通过转换计算，可见除非输入阻抗和输出阻抗相等，才能使得两者严格相等。由此我们得到截止频率的严格定义：输入幅度保持不变，输入功率和输出功率作用在同一负载下，截止频率点处输出功率衰减为输入功率的一半。根据这个定义，我们才能理所当然的写出如下等式：
$$\frac{1}{2}=|\frac{P_o}{P_i}|\stackrel{设相同负载为R}{\to }\frac{1}{2}=|\frac{\frac{U_o^2}{R}}{\frac{U_i^2}{R}}|=|H({\omega }_c){|}^2$$
$$|H({\omega }_c)|=\frac{1}{\sqrt{2}}\approx 0.707\stackrel{等式两边同时取对数}{\to }20lo{g}_{10}|H({\omega }_c)|\approx -3.01(dB)$$
由此，我们解释了截止频率的定义，以及 $-3dB$的来源，唯一需要特别注意的就是这里的截止频率处半功率点的取法，前提是在相同阻抗下取得的，而非广义上系统实际的半功率点。接下来，就依据这个关系式，推导RLC的截止频率 $f_c$：
$$由截止频率定义得：(1-{\omega }_c^{2}LC{)}^2+({\omega }_{c}RC{)}^2=2$$
$$展开整理得：{\omega }_c^4{L}^2{C}^2+{\omega }^2(R^2{C}^2-2LC)=1$$
$$做变量代换，令x={\omega }_c^2：x^2(LC{)}^2+x(R^2{C}^2-2LC)-1=0$$
$$求根公式，保留正根：x=\frac{-(R^2{C}^2-2LC)+\sqrt{(R^2{C}^2-2LC{)}^2+4L^2{C}^2}}{2L^2{C}^2}$$
到这一步，换元后还元即可得到，RLC电路截止角频率的公式 ${\omega }_c=\sqrt{x}$，但是这个公式似乎有些臃肿，不够优雅。我们再做一些化简，思路就是使用分子有理化，最终整理得如下公式：
$$分子有理化，提公因式：x=\frac{{\omega }_n^2}{\sqrt{(2{\xi }^2-1{)}^2+1}+(2{\xi }^2-1)}$$
$$RLC电路截止频率计算公式：f_c=\frac{\sqrt{x}}{2\pi }$$
为了保证这个公式的正确性和直观性，使用Mathcad可视化数学分析工具，计算出截止频率结果，然后将这个结果输入到Psim仿真测量工具中。计算过程和Psim仿真分别如图10、11所示：

图10 Mathcad可视化计算截止频率

图11 Psim仿真标定截止频率

由此证明上述推导的截止频率计算公式正确。

S域分析

\hspace{2em} 经过前文的时域和频域分析，相信读者已经对RLC电路系统已经有了较为深刻的理解，对比一下时域和频域不难发现，频域中引入阻抗分析法求解传递函数后，大大简化了计算量，而且如果输入是绝对可积的连续时间信号，还可以使用傅里叶逆变换将其变换回时域，可以说傅里叶变换引入频域为分析信号和系统带来极大的方便。但是，美中不足的就是傅里叶变换需要满足绝对可积的条件（即狄利克雷条件）才能直接使用。例如在我们在做时域分析时，为了观察系统的动态性能，输入的阶跃信号是一种奇异信号，并不满足绝对可积的条件，数学上通过广义函数理论可求解阶跃信号的傅里叶变换，因此想要使用傅里叶变换研究系统的阶跃响应是存在困难的。为了保证傅里叶变换是绝对可积的，法国数学家拉普拉斯通过引入一个指数衰减因子$e^{-\sigma t}$修正了傅里叶变换，使得频域过渡到复频域，定义了复频域的变量为$s=\sigma +j\omega$，所以复频域也称为S域。这样经过修正后使得不绝对可积的信号依旧可以适用变换的处理方法研究系统，这就是拉普拉斯变换。因此，拉普拉斯变换可以处理阶跃信号这种瞬态变换的信号，也可以处理稳态信号，适用于分析系统的瞬态响应和稳态响应。在工程领域大多数工程师都倾向于使用拉普拉斯变换分析线性时不变系统，不过时域、频域、S域都有其各自的特点，作者建议最好是都能够掌握不同域的分析方法。
\hspace{2em} 接下来，我们就使用拉普拉斯变换研究RLC电路初始条件为零的阶跃响应求解问题。同频域分析类似，在S域可以使用复阻抗分析的方法求解传递函数，只是注意在S域下传递函数的变量是s，依据频域分析的方法，我们可以直接写出S域下的传递函数：
$$H(s)=\frac{1}{1+sRC+s^{2}LC}$$
传递函数在时域也被称为系统函数$h(t)$，在时域中输入经过系统再到输出，可以用卷积运算表示这个过程，输出信号等于输入信号和系统函数的卷积，时域卷积变换到频域和S域中有一个很重要结论，就是在频域和S域中输出信号等于输入信号和传递函数的乘积：
$$y(t)=f(t)\ast h(t)\stackrel{拉普拉斯变换}{\to }Y(s)=F(s)H(s)$$
因此在S域中要求解RLC电路的阶跃响应非常的简单，只需要将时域阶跃信号变换到S域再乘上传递函数即可，单位阶跃信号的拉普拉斯变换等于$\frac{1}{s}$，所以写出RLC电路阶跃响应的S域表达式如下:
$$Y(s)=\frac{H(s)}{s}$$
如果系统是稳定的，拉普拉斯变换还有一个很有用的结论用于求解系统的稳态响应，叫做终值定理：
$$终值定理：\underset{t\to +\infty }{\lim }y(t)=\underset{s\to 0}{\lim }sY(s)$$
从时域表达式中可以看出RLC电路是稳定系统，因此使用终值定理的公式很快就能算出，单位阶跃信号作用在RLC电路，其输出最终等于1Ｖ，可见终值定理用于求解系统稳态响应是非常方便的。这里插一句题外话，关于系统在S域下的稳定性判断，在后文会从极点分布的角度解释。稳定性判断是很大的一个话题，本文篇幅有限，无法在此展开。不过预计将会在后续的文章中分析S域下的各种重要的分析方法，结合实例分享例如零极点分布与配置、稳定性判据、根轨迹法等经典控制理论相关知识。所以，这里介绍S域也只是浮于表面，目的是让读者先具备一个简单的认知，更重要的是这里的知识深度，也足够分析RLC二阶电路系统了。
\hspace{2em} 回到正题上，求解时域阶跃响应，需要计算拉普拉斯逆变换，手工计算一般思路就是采用部分分式展开再利用公式求解。这种繁琐的工作，我们依旧采用数值分析工具进行计算，依据S域的分析思路，编写程序：

clc
clear
pkg load control%导入控制工具箱
R=10;L=1e-3;C=1e-6;%定义器件参数
t=0:1e-6:2*1e-3;%定义计算范围
num=[1];%分子多项式的系数
den=[L*C,R*C,1];%分母多项式系数:LC*s^2+RC*s+1
H_s=tf(num,den);%构建传递函数
figure
step(H_s,t);%根据传递函数绘制阶跃响应
figure
pzmap(H_s);%绘制零极点分布图
xlim([-8000,4000]);

图12 传递函数S域求解阶跃响应

可见求解结果与前文时域分析是一致的，说明分析过程正确。同时代码中还调用了绘制零极点分布图的函数，这里简单介绍一下，零极点的求解方法，顺便解释一下前文提到的无阻尼频率和特征方程复数根实部与虚部之间的关系，以及S平面的阻尼比。
零极点的求解方法，也很简单（其来源就复杂许多，零极点有许多特性，计划在今后更新的文章中再剖析），将分子视为零点多项式，分母视为极点多项式。令分子多项式等于零求出来的根就是零点（RLC的例子中分子是常数，表明没有零点），令分母多项式等于零求出来的根就是极点，在RLC电路中极点方程其实就等同于时域中的二阶特征方程，表明有两个根。因此，RLC电路是一个没有零点，存在两个极点的系统，其极点分布如图13所示：

图13 RLC电路（欠阻尼）零极点分布图

在复频域中用S平面进行绘制零极点分布，横坐标是实部轴，纵坐标是虚部轴，平面中用o表示零点，用×表示极点。可见图中并无零点，存在两个极点，当系统极点完全分布在S平面的左半平面时，就可以判断系统是稳定的，因为实部对应时域表达式的指数项，负指数信号最终会衰减到0 。有读者可能会疑惑，RLC阶跃响应刚开始不是会振荡吗？怎么是稳定的呢？这里一定要搞定系统稳定的概念，系统稳定是有界输入得到有界输出，也就是在某个时间节点之前系统在振荡，经过某时间节点后稳定，系统最终也是稳定的。真正不稳定的振荡系统，即使时间推向无穷，输出也不是有界的，在S平面中体现为极点落在虚轴以及右半平面。例如RLC电路无阻尼的情况极点为纯虚数，就可以构成LC振荡器形成稳定振荡。当极点在右半平面，实际体现在时域表达式中指数项为正指数信号，数学上会将信号放大至无穷，不过实际系统中可能会放大到非线性区域，最终受到限制，显然这是不稳定的。
\hspace{2em} 在S平面我们还会使用极坐标的方式表示零极点，也就是使用模和相角进行表示。正如图中标注三角形斜边的箭头和角度，显然根据勾股定理，可以计算模值的大小。还有一个重要的信息，就是模值是由实部和虚部两个分量构成的，而且这个模值正如前文所说是等于谐振角频率的，这里有什么实际的物理意义吗？假设以模值为半径，可以绘制一个左半圆，圆半径是不变的，而虚部和实部的分量是此消彼长的，这就对应着系统中存在一个固定的频率，只由系统的储能元件决定，储能元件能在这个固定的频率下恰好完成能量的交换，这个频率就是系统的固有频率，也称为谐振频率，也可以是前文中写的无阻尼频率。当系统中存在其他的耗能元件，比如电阻器，能量势必会重新分布，一部分会耗散在实部代表的耗能元件上，一部分能量还在虚部代表的储能器件之间交换，因此虚部称为阻尼角频率，这两个分量的和是保持不变的，也就是圆的半径是不变的，当然随着电阻的持续加大，系统会进入过阻尼状态，这种状态下就不存在虚部了，自然也就没有这个关系了。为了表征系统由于耗能器件影响系统阻尼频率与固有频率的偏差程度，可以使用相角进行描述，取相角的正弦值进行观察，也就是阻尼比。在S平面的第二象限中，模和虚轴之间的夹角（相角）越大，说明阻尼频率和固有频率之间偏差得越多。可以将电阻变化的过程在S平面绘制成一个轨迹，也就是根轨迹图。绘制过程中不断的计算模值的大小，计算结果可以发现$R<63$，模值转换为频率值一直等于前文计算的$f_n=5032.9$，当$R>63$不再出现虚部，这里的电阻值就是通过前文的临界阻尼条件公式进行计算得来，计算结果为$R=63.246$。绘制根轨迹的程序以及运行结果如下：

clc
clear
pkg load control%导入控制工具箱
R=0;L=1e-3;C=1e-6;%定义器件参数
figure
for i=0:0.5:70
hold on
num=[1];%分子多项式的系数
den=[L*C,R*C,1];%分母多项式系数:LC*s^2+RC*s+1
poles=roots(den);%求极点
H_s=tf(num,den);%构建传递函数
pzmap(H_s);%也可以用根轨迹函数进行绘制rlocus(H_s);
xlim([-40000,40000]);
ylim([-40000,40000]);
R
f_n=abs(poles)./(2*pi)
R=R+0.5;
endfor

图14 RLC电路电阻变化时的根轨迹图

图15 RLC电路电阻变化时模值转换频率运算结果

\hspace{2em} 关于RLC二阶电路系统的分析就到此结束了，如果读者很仔细的从头看到尾，说明您大约阅读了两万字！本文我历时差不多半个月完成，力求文章知识点的系统完整性和连贯性，希望对各位读者能够有所帮助和提升。由于作者水平有限，所阅读的资料有限，文章中难免有疏漏和不妥之处，望各位读者批评指正！