传递函数离散化方法：差分，双线性Tustin，预矫正双线性，Matlab实现

传递函数严格离散化

  傅里叶变换、拉普拉斯变换，以及$Z$变换的定义如下所示：
$$傅里叶变换:\{\begin{array}{l}F\left(\omega \right)={\int }_{-\infty }^{\infty }f\left(t\right)e^{-j\omega t}dt\\ f\left(t\right)=\frac{1}{2\pi }{\int }_{-\infty }^{\infty }F\left(\omega \right)e^{j\omega t}d\omega \end{array}$$

$$拉普拉斯变换:\{\begin{array}{l}F\left(s\right)={\int }_{-\infty }^{\infty }f\left(t\right)e^{-st}dt\\ f\left(t\right)=\frac{1}{2\pi j}{\int }_{\sigma -j\infty }^{\sigma +j\infty }F\left(s\right)e^{st}ds\end{array}$$

$$Z变换:\{\begin{array}{l}F\left(z\right)=\sum f\left[n\right]z^{-n}\\ f\left[n\right]=\frac{1}{2\pi j}\oint F\left(z\right)z^{n-1}dz\end{array}$$

  可以利用如下映射关系，实现传递函数$F(s)$、$F(\omega )$与$F(z)$之间的等价转换，详见《信号与系统》的“离散系统的Z域分析”章节：

$$s=j\omega =\frac{1}{T}\ln z，或，z=e^{sT}=e^{j\omega T}$$

  稳定性分析： 连续型传递函数稳定的充要条件为，所有极点均位于$s$域左半平面；离散型传递函数稳定的充要条件为，所有极点均位于$z$域原点为圆心的单位圆内。对于连续传函离散化问题，$s$域左半平面所有极点，经$z=e^{sT}$映射后，均位于$z$域原点为圆心的单位圆内。因此，严格离散化前后，连续系统稳定，则离散系统稳定。

离散化方法1，前向欧拉，前向差分

  对于一个连续函数$y(t)$在$kT$时刻的微分可以做以下近似：
$$y^{\prime }(k)\approx \frac{y(k+1)-y(k)}{T}$$

  对上式的左右两端分别进行拉氏变换和Z变换，可以得到如下映射关系：
$$sY(s)\leftrightarrow \frac{z-1}{T}Y(z)$$

  因此，可以令$s=\frac{z-1}{T}$实现S域传递函数到Z域传递函数的转换。值得注意的是，这种传递函数离散化是近似的，$s=\frac{1}{T}\ln z$才是严格的，但并不方便数字实现。实际上，$s=\frac{z-1}{T}$正是$s=\frac{1}{T}\ln z$，在$s=0$即$z=1$时的近似表示。

  稳定性分析： 假设连续传递函数是稳定的，即传递函数极点位于$s$域左半平面，分析离散化后极点在$z$域平面位置，
$$s=\frac{z-1}{T}\Rightarrow z=sT+1$$
  前向欧拉离散化方法将$s$域左半平面映射到$z$域$z<1$区域，显然包含不稳定区域，稳定极点可能变为不稳定极点。因此，稳定的连续系统变换所得离散系统可能不稳定，这种方法很少用。

离散化方法2，后向欧拉，后向差分

  同理还有后向差分，将$y(t)$在$kT$时刻的微分做以下近似：
$$y^{\prime }(k)\approx \frac{y(k)-y(k-1)}{T}$$
  对上式的左右两端分别进行拉氏变换和Z变换，可以得到如下映射关系：
$$sY(s)\leftrightarrow \frac{1-z^{-1}}{T}Y(z)$$
  因此，可以令$s=\frac{1-z^{-1}}{T}$实现S域传递函数到Z域传递函数的转换。实际上，$s=\frac{1-z^{-1}}{T}$也是$s=\frac{1}{T}\ln z$，在$s=0$即$z=1$时的近似表示。

  稳定性分析： 假设连续传递函数是稳定的，即连续传函极点位于$s$域左半平面，分析离散化后极点在$z$域平面位置，
$$s=\frac{1-z^{-1}}{T}\Rightarrow z=\frac{1}{1-sT}$$
  后向欧拉离散化方法将$s$域左半平面映射到$z$域以$(0.5,0)$为圆心，以$0.5$为半径的圆内（具体分析见附录）。也就是说，原来的连续系统稳定，则离散系统稳定。

离散化方法3，双线性变换，Tustin

  由于$z=e^{sT}=\frac{e^{sT/2}}{e^{-sT/2}}\approx \frac{1+sT/2}{1-sT/2}$，分子分母进行一阶泰勒线性展开，故得名双线性变换。进一步可以得到：
$$s=\frac{2}{T}\frac{1-z^{-1}}{1+z^{-1}}$$

  稳定性分析： 假设连续传递函数是稳定的，分析离散化后极点在$z$域平面位置，
$$s=\frac{2}{T}\frac{1-z^{-1}}{1+z^{-1}}\Rightarrow z=\frac{1+sT/2}{1-sT/2}$$
  双线性变换离散化方法将$s$域左半平面映射到$z$域以$(0,0)$为圆心，以$1$为半径的单位圆内（具体分析见附录），这正是我们很希望看到的结果。因此，原来的连续传函稳定，则离散传函稳定。

离散化方法4，预矫正双线性变换，Tustin with Pre-Warping

问题描述

  上述方法均可实现连续系统的离散化，但因为离散化过程中采用了近似化手段，不可避免地引入了误差，使得离散化后的频域特性并不完全等价于连续系统，以下以双线性变换为例，进行说明。设计谐振频率${\omega }_0=1000$ $rad/s$的传递函数：
$$C_a(s)\text{ = }\frac{s}{s^2+{1000}^2}$$

  采用双线性变换$s=\frac{2}{T}\frac{z-1}{z+1}$进行离散化，其中T取0.001，得到离散化后的传递函数：
$$C_d(z)=C_a(s){|}_{s=\frac{2}{T}\frac{z-1}{z+1}}\text{ = }\frac{\left(\frac{2}{T}\frac{z-1}{z+1}\right)}{{\left(\frac{2}{T}\frac{z-1}{z+1}\right)}^2+{1000}^2}=\frac{2000z^2-2000}{\left({2000}^2+{1000}^2\right)z^2+2\left({1000}^2-{2000}^2\right)z+\left({1000}^2+{2000}^2\right)}$$

  利用如下代码，对比连续传递函数和离散传递函数的波特图：

a=1000^2;b=2000^2;
Ca = tf([1,0],[1 0 a]);
Cd = tf([2000 0 -2000],[b+a 2*a-2*b b+a],0.001);
bode(Ca, Cd);legend('Ca','Cd')

  当然，也可以利用matlab实现双线性变换，两者效果是一致的（更多函数介绍请参考matlab函数用法：传递函数tf，波特图bode，离散化c2d）：

Ca = tf([1,0],[1 0 1000^2]);
Cd = c2d(Ca,0.001,'tustin');
bode(Ca, Cd);legend('Ca','Cd')

  连续传递函数和离散传递函数的波特图对比结果如下所示。可以看出，近似的离散化过程使得离散传递函数特性发生了改变，离散化后传递函数的特性与所设计连续传递函数的特性不一致。

解决方案

  对离散化所引入的误差进行分析，$C_d(z)$代表数字控制器，$C_a(s)$代表模拟控制器，两者之间的转换关系如下：
$$C_a(s){|}_{s=\frac{2}{T}\frac{z-1}{z+1}}=C_d(z)$$

  下面探究一下$C_a(s){|}_{s=\frac{2}{T}\frac{z-1}{z+1}}$所对应的真实连续系统，注：$\sin \alpha =-j0.5\left(e^{j\alpha }-e^{-j\alpha }\right),\cos \alpha =0.5\left(e^{j\alpha }+e^{-j\alpha }\right)$:
$$\{\begin{array}{rl}{C}_a(s){|}_{s=\frac{2}{T}\frac{z-1}{z+1}}& =C_a(\frac{2}{T}\frac{e^{sT}-1}{e^{sT}+1})=C_a\left(\frac{2}{T}\frac{e^{j\omega T}-1}{e^{j\omega T}+1}\right)\\ & =C_a\left[\frac{2}{T}\frac{e^{j\frac{\omega T}{2}}\left(e^{j\frac{\omega T}{2}}-e^{-j\frac{\omega T}{2}}\right)}{e^{j\frac{\omega T}{2}}\left(e^{j\frac{\omega T}{2}}+e^{-j\frac{\omega T}{2}}\right)}\right]\\ & =C_a\left[j\ast \frac{2}{T}\tan \left(\frac{\omega T}{2}\right)\right]\end{array}$$

  $C_d(z)$并不对应$C_a(j\omega )$，而对应$C_a\left[j\frac{2}{T}\tan \left(\frac{\omega }{2}T\right)\right]$，离散系统的$\omega$会映射到连续系统的$\frac{2}{T}\tan \left(\frac{\omega T}{2}\right)$，这就是频率畸变现象。也就是说，分别以$\omega$为横轴画出$C_d(z)$的波特图，和以$\frac{2}{T}\tan \left(\frac{\omega }{2}T\right)$为横轴画出$C_a(s)$的波特图，两个波形图应该是完全对应的。当$\frac{\omega T}{2}$足够小时，$j\frac{2}{T}\tan \left(\frac{\omega }{2}T\right)$是近似等于$j\omega$的。但对于高频谐波，或者特定频率的分析，就不得不考虑这部分影响了。

  回顾上述所设计$C_a(s)$时，对于$1000$ rad/s的角速度变化，传递函数取得最大值。由于离散系统的$\omega$会映射到连续系统的$\frac{2}{T}\tan \left(\frac{\omega T}{2}\right)$，利用$\frac{2}{T}\tan \left(\frac{\omega T}{2}\right)=1000$求解离散系统的原始$\omega$，解得$\omega =927.3$ rad/s，这与上述波特图结果是一致的。

  解决方案：离散系统的$\omega$与连续系统的$\frac{2}{T}\tan \left(\frac{\omega T}{2}\right)$之间存在映射关系，考虑将连续系统进行整体比例缩放，从而实现两者在目标频率${\omega }^{\ast }$处表现一致，具体操作如下表所示。 注意，预矫正双线性变换只是通过结合比例坐标变换，实现在目标频点处的特性一致，其他频点处特性仍存在差异。

离散	变换关系	连续
$\omega$	$s=\frac{2}{T}\frac{z-1}{z+1}$	$\frac{2}{T}\tan \left(\frac{\omega T}{2}\right)$
${\omega }^{\ast }$	$s=\frac{2}{T}\frac{z-1}{z+1}$	$\frac{2}{T}\tan \left(\frac{{\omega }^{\ast }T}{2}\right)$
$\omega$	$s=\frac{\frac{{\omega }^{\ast }T}{2}}{\tan \left(\frac{{\omega }^{\ast }T}{2}\right)}\frac{2}{T}\frac{z-1}{z+1}=\frac{{\omega }^{\ast }}{\tan \left(\frac{{\omega }^{\ast }T}{2}\right)}\frac{z-1}{z+1}$	$\frac{\frac{{\omega }^{\ast }T}{2}}{\tan \left(\frac{{\omega }^{\ast }T}{2}\right)}\frac{2}{T}\tan \left(\frac{\omega T}{2}\right)=\frac{{\omega }^{\ast }}{\tan \left(\frac{{\omega }^{\ast }T}{2}\right)}\tan \left(\frac{\omega T}{2}\right)$
${\omega }^{\ast }$	$s=\frac{{\omega }^{\ast }}{\tan \left(\frac{{\omega }^{\ast }T}{2}\right)}\frac{z-1}{z+1}$	$\frac{{\omega }^{\ast }}{\tan \left(\frac{{\omega }^{\ast }T}{2}\right)}\tan \left(\frac{{\omega }^{\ast }T}{2}\right)={\omega }^{\ast }$

  取${\omega }^{\ast }=1000$ rad/s进行预矫正，使得离散传函和连续传函在${\omega }^{\ast }$处特性一致。采用预矫正双线性对$C_a(s)\text{ = }\frac{s}{s^2+{1000}^2}$进行离散化（$T=0.001$），如下所示：
$$C_d(z)=C_a(s){|}_{s=\frac{1000}{\tan \left(\frac{1000\ast 0.001}{2}\right)}\frac{z-1}{z+1}}\text{ = }\frac{\left(1830\frac{z-1}{z+1}\right)}{{\left(1830\frac{z-1}{z+1}\right)}^2+{1000}^2}=\frac{1830z^2-1830}{\left({1830}^2+{1000}^2\right)z^2+2\left({1830}^2-{1000}^2\right)z+\left({1000}^2+{1830}^2\right)}$$

  对比连续传函和离散传函波特图，如下所示：

  Matlab实现方案一：解析计算离散传递函数（如上所述），绘制波特图，

a=1000^2;b=1830^2;
Ca = tf([1,0],[1 0 a]);
Cd = tf([1830 0 -1830],[b+a 2*a-2*b b+a],0.001);
bode(Ca, Cd);legend('Ca','Cd');grid on;

  Matlab实现方案二：借助Matlab分别实现传递函数预矫正和双线性变换，

Ca = tf([1, 0], [1, 0, 1000^2]);            % 定义连续传递函数 Ca(s)

T = 0.001;                                  % 采样周期
omega_p = 1000;                             % 预矫正频率
K = (omega_p * T / 2)/tan(omega_p * T / 2); % 计算预缩放因子 K

s = tf('s');                                    % 定义 s 为传递函数变量
Ca_prewarped = (s * K) / ((s * K)^2 + 1000^2);  % s * K 取代 s，预矫正后的传递函数
Cd = c2d(Ca_prewarped, T, 'tustin');            % 使用双线性变换进行离散化

bode(Ca, Cd);legend('Ca','Cd');grid on;

  Matlab实现方案三：借助Matlab直接实现预矫正双线性变换，

Ca = tf([1, 0], [1, 0, 1000^2]);    % 定义连续传递函数 Ca(s)
opt = c2dOptions('Method','tustin','PrewarpFrequency',1000);
Cd = c2d(Ca, 0.001, opt);            % 使用双线性变换进行离散化
bode(Ca, Cd);legend('Ca','Cd');grid on;

总结

  对离散化的简单场景，优先考虑后向差分离散化方法；对于频率有较高准确度的场景，如谐振控制器，建议使用预矫正双线性变换方法。

参考文献

[1] 余成波.信号与系统(第二版)[M].清华大学出版社,2007.
[2] GeneF.Franklin,J.DavidPowell,MichaelWorkman.动态系统的数字控制(第3版)[M].清华大学出版社,2001.
[3] Matlab文档：c2d
[4] Matlab文档：c2dOption
[5] Matlab文档：Continuous-Discrete Conversion Methods

附录

附录1：后向欧拉相关证明

  对于$s$域的左半平面，其边界为$s$平面虚轴。证明：该边界在$z$域的映射为，以$(0.5,0)$为圆心，以$0.5$为半径的圆，如下所示:
$$\{\begin{array}{rl}z& =\frac{1}{1-sT}\\ & =\frac{1}{1-j\omega }\\ & =\frac{1+j\omega }{\left(1-j\omega \right)\left(1+j\omega \right)}\\ & =\frac{1}{1+{\omega }^2}+j\frac{\omega }{1+{\omega }^2}\end{array}\to \{\begin{array}{rl}{\Vert z-0.5\Vert }^2& ={\left(\frac{1}{1+{\omega }^2}-0.5\right)}^2+{\left(\frac{\omega }{1+{\omega }^2}\right)}^2\\ & =\frac{1}{4}\left[{\left(\frac{2}{1+{\omega }^2}-1\right)}^2+{\left(\frac{2\omega }{1+{\omega }^2}\right)}^2\right]\\ & =\frac{1}{4}\left[{\left(\frac{1-{\omega }^2}{1+{\omega }^2}\right)}^2+{\left(\frac{2\omega }{1+{\omega }^2}\right)}^2\right]\\ & =\frac{1}{4}，证毕\end{array}$$

  进一步证明，$s$域内任何左半平面的点，映射到上述圆内：
$$\{\begin{array}{rl}z& =\frac{1}{1-sT},Re(s)<0\\ & =\frac{1}{a-bj},a>1\\ & =\frac{a+bj}{\left(a-bj\right)\left(a+bj\right)}\\ & =\frac{a}{a^2+b^2}+j\frac{b}{a^2+b^2}\end{array}\to \{\begin{array}{rl}{\Vert z-0.5\Vert }^2& ={\left(\frac{a}{a^2+b^2}-0.5\right)}^2+{\left(\frac{b}{a^2+b^2}\right)}^2\\ & =\frac{1}{4}\left[{\left(\frac{2a}{a^2+b^2}-1\right)}^2+{\left(\frac{2b}{a^2+b^2}\right)}^2\right]\\ & <\frac{1}{4}\left[{\left(\frac{2a}{a^2+b^2}-1\right)}^2+{\left(\frac{2ab}{a^2+b^2}\right)}^2\right]\\ & <\frac{1}{4}\left[{\left(\frac{a^2-b^2}{a^2+b^2}\right)}^2+{\left(\frac{2ab}{a^2+b^2}\right)}^2\right]\\ & =\frac{1}{4}，证毕\end{array}，其中\{\begin{array}{rl}& \frac{2a}{a^2+b^2}-1<\frac{a^2-b^2}{a^2+b^2}\\ & \leftarrow 2a-a^2-b^2<a^2-b^2\\ & \leftarrow 2a<2a^2\\ & \leftarrow 1<a\end{array}$$

  Matlab绘图验证：

  Matlab代码如下：

clc; clear; close all;

% 定义参数
T = 1;

% 定义s=0（虚部从-无穷到正无穷）
omega = linspace(-1000, 1000, 20001);  % 使用大的范围代替无穷
s1 = 0.0 + 1j * omega;
s2 = -0.1 + 1j * omega;
s3 = -1.0 + 1j * omega;
s4 = -10.0 + 1j * omega;

% 映射到z平面
z1 = 1 ./ (1 - s1 * T);
z2 = 1 ./ (1 - s2 * T);
z3 = 1 ./ (1 - s3 * T);
z4 = 1 ./ (1 - s4 * T);

% 绘制z平面的值
figure;
hold on;

% 绘制s=0的映射
plot(real(z1), imag(z1), 'b', 'DisplayName', 's = 0 (Im: -\infty to +\infty)');

% 绘制s=-0.1的映射
plot(real(z2), imag(z2), 'r', 'DisplayName', 's = -0.1 + j\omega');

% 绘制s=-1.0的映射
plot(real(z3), imag(z3), 'g', 'DisplayName', 's = -1.0 + j\omega');

% 绘制s=-10.0的映射
plot(real(z4), imag(z4), 'm', 'DisplayName', 's = -10.0 + j\omega');

% 设置图形属性
xlabel('Real Part');
ylabel('Imaginary Part');
title('Mapping of s-plane lines to z-plane');
legend show;
grid on;
axis equal;
hold off;

附录2：双线性变换相关证明

  对于$s$域的左半平面，其边界为$s$平面虚轴。证明：该边界在$z$域的映射为，以$(0,0)$为圆心，以$1$为半径的圆，如下所示:
$$\{\begin{array}{rl}z& =\frac{1+sT/2}{1-sT/2}\\ & =\frac{1+j\omega }{1-j\omega }\\ & =\frac{1-{\omega }^2}{1+{\omega }^2}+j\frac{2\omega }{1+{\omega }^2}\end{array}\to \{\begin{array}{rl}{\Vert z\Vert }^2& ={\left(\frac{1-{\omega }^2}{1+{\omega }^2}\right)}^2+{\left(\frac{2\omega }{1+{\omega }^2}\right)}^2\\ & =1\end{array}$$
  对于边界内的点，类比后向欧拉方法，可证明其在$z$平面的映射在单位圆内。给出Matlab绘图结果：

  Matlab代码如下：

clc; clear; close all;
% 定义参数
T = 1;

% 定义s=0（虚部从-无穷到正无穷）
omega = linspace(-1000, 1000, 20001);  % 使用大的范围代替无穷
s1 = 0.0 + 1j * omega;
s2 = -0.1 + 1j * omega;
s3 = -1.0 + 1j * omega;
s4 = -10.0 + 1j * omega;

% 映射到z平面
z1 = (1 + s1 * T) ./ (1 - s1 * T);
z2 = (1 + s2 * T) ./ (1 - s2 * T);
z3 = (1 + s3 * T) ./ (1 - s3 * T);
z4 = (1 + s4 * T) ./ (1 - s4 * T);

% 绘制z平面的值
figure;
hold on;

% 绘制s=0的映射
plot(real(z1), imag(z1), 'b', 'DisplayName', 's = 0 (Im: -\infty to +\infty)');

% 绘制s=-0.1的映射
plot(real(z2), imag(z2), 'r', 'DisplayName', 's = -0.1 + j\omega');

% 绘制s=-1.0的映射
plot(real(z3), imag(z3), 'g', 'DisplayName', 's = -1.0 + j\omega');

% 绘制s=-10.0的映射
plot(real(z4), imag(z4), 'm', 'DisplayName', 's = -10.0 + j\omega');

% 设置图形属性
xlabel('Real Part');
ylabel('Imaginary Part');
title('Mapping of s-plane lines to z-plane');
legend show;
grid on;
axis equal;
hold off;