写在前面

在学习过程中，接触到的都是连续系统，而计算机是离散的，因此对于连续系统进行离散处理是很有必要的。描述连续系统使用的是微分方程，而描述离散系统采用差分方程，因此将系统离散化的一般流程为

本文以一阶低通RC滤波器为例，实现从传递函数到差分方程的转化。

RC低通滤波器

滤波器的结构如图所示

滤波器的微分方程与传递函数如下

离散方法

离散化方法有很多，比如一阶向后差法、双线性变换法、零极点匹配法、保持器等价法、脉冲响应不变法、阶跃响应不变法。各个方法的特点的与优势不尽相同，下面将使用脉冲响应不变法、向后差分、双线性变换法进行系统的离散化。
离散方法的详细内容可以参考这篇文章： https://www.docin.com/p-540901626.html

脉冲响应不变法

根据脉冲响应不变法公式

可以得到系统传递函数的z变换如下：

进而得到差分方程：

一阶向后差分

一阶向后差分使用两次的差值除去采样时间T近似作为微分，其从s域到z域的变换算子如下
变换结果如下

进行z反变换求得差分方程

双线性变换法

双线性变换法变换算子：

变换结果：

解得

验证

验证方法与工具

完成了对系统的离散化并且得到了差分方程之后，需要对结果进行验证，验证的方式为使用matlab，将连续系统的时域输出与差分方程的输出进行对比。以此验证离散结果是否正确。使用 lsim 函数输出连续系统的输入与响应，使用 filter 函数来求取差分方程的结果。
lsim函数：lsim函数是针对线性时不变模型，给定任意输入，得到任意输出。lsim(sys,u,t)绘制传递函数为sys的系统对输入u的模拟时间响应。 t是仿真时间。
filter 函数：filter函数可以实现差分方程的递规求解。
设差分方程的形式为
a0y[n] + a1y[n−1] + ⋯ + amy[n−m] = b0x[n] + b1x[n−1] + ⋯ + bsx[n−s]
在零初始情况下（t<0,y[n] = 0,x[n] = 0）
使用方法为：yn = filter(B, A, xn)。其中 B=[b0,⋯,bs],A=[a0,⋯,am]；xn是输入信号；yn是零状态响应。

实施

为了方便起见。设采样时间T=1，a = 0.3；
则根据上述式子，可得：

写入matlab程序中
以xt = 1cos(0.3t)为函数，t = 0:0.01:90;为连续时间，使用lsim(sys,xt,t)打映时域响应yt（灰色是xn，蓝色是响应yn）。结果如下：

以T为采样时间，采样次数为n（90/T 次），则采样点的函数值 xn = 1cos(0.3n*T);使用filter函数计算输出值yn。结果如下：

将连续输出与离散输出相对比，可得各个方法之间的差异：

将xt 改为 1cos(0.1t)+cos(2t);则xn 变为 1cos(0.1nT)+1cos(2n*T);离散输出与时域响应符合

将采样时间T变为 T=2，查看结果，离散结果相对于T = 1时变差，符合现实情况。

将采样时间T变为 T=0.5，查看结果，除脉冲响应不变法外，离散结果相对于T = 2时变好。

matlab代码

附上整体代码：

clc;clear;close all;

%建立幕布
figure;

%1/RC 与采样时间T
a = 0.3; T = 0.1;

%生成传递函数
sys = tf(a,[1,a]);
sys
% 输入序列 为 连续
t = 0:0.01:90;
xt = 1*cos(0.1*t)+cos(2*t);
%xt = 1*cos(0.3*t);
%打映时域输入与响应
lsim(sys,xt,t);
hold on

% 输入序列 为 任意的离散
n = 0:(90/T);
xn = 1*cos(0.1*n*T)+1*cos(2*n*T);
%xn = 1*cos(0.3*n*T);
%一阶向后差分
% 线性常系数差分方程 中的 x(n) 项系数
B=(a*T)/(1+a*T);
% 线性常系数差分方程 中的 y(n) 项系数
A=[1, -1/(1+a*T)];
% 输出序列
yn=filter(B,A,xn); 
%绘制 "输出序列" 图像
plot(n*T,yn, 'cs');hold on;
%脉冲响应不变法
% 线性常系数差分方程 中的 x(n) 项系数
B=a;
% 线性常系数差分方程 中的 y(n) 项系数
A=[1, -exp(-a*T)];
% 输出序列
yn2=filter(B,A,xn); 
plot(n*T,yn2, 'b+');hold on;
%双线性不变法
% 线性常系数差分方程 中的 x(n) 项系数
B=[(a*T)/(2+a*T),(a*T)/(2+a*T)];
% 线性常系数差分方程 中的 y(n) 项系数
A=[1, (a*T-2)/(a*T+2)];
% 输出序列
yn3=filter(B,A,xn); 
plot(n*T,yn3, 'r*');hold on;


%标注
legend('tout','(y[n]-y[n-1])/T','imp','tustin')
%legend('(y[n]-y[n-1])/T','imp','tustin')
% 打开网格
grid on;