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

1. 传递函数严格离散化

  傅里叶变换、拉普拉斯变换，以及$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}$$

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

  对于一个连续函数$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$域$z<1$区域，显然包含不稳定区域，因此稳定的连续系统变换所得离散系统可能不稳定，这种方法很少用。

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

  同理还有后向差分，将$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$域以0.5为圆心0.5为半径的圆，也就是说，原来连续系统稳定，则离散系统稳定。

4. 离散化方法，双线性变换，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}}$$
  该方法将且s左半平面映射到z域单位圆内，这是我们很希望看到的结果。

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

5.1 问题描述

  上述方法均可实现连续系统的离散化，但因为离散化过程中采用了近似化手段，不可避免地引入了误差，使得离散化后的频域特性并不完全等价于连续系统，以下以双线性变换为例，进行说明。设计电网系统的五次谐波消抑制传递函数：
$$C_a(s)=\frac{K_{i}s}{s^2+{(2\pi \ast 250)}^2}$$
  取$K_i=1$，采用双线性变换$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}}=\frac{K_i(\frac{2}{T}\frac{z-1}{z+1})}{{\left(\frac{2}{T}\frac{z-1}{z+1}\right)}^2+{(2\pi \ast 250)}^2}=\frac{2000z^2-2000}{\left[{2000}^2+{(2\pi \ast 250)}^2\right]z^2+2\left[{(2\pi \ast 250)}^2-{2000}^2\right]z+\left[{(2\pi \ast 250)}^2+{2000}^2\right]}$$

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

a=(2*pi*250)^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 (2*pi*250)^2]);
Cd = c2d(Ca,0.001,'tustin');
bode(Ca, Cd);
legend('Ca','Cd')

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

5.2 解决方案

  对离散化所引入的误差进行分析，$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{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)$时，对于$2\pi \ast 250$ rad/s的角速度变化，传递函数取得最大值。由于离散系统的$\omega$会映射到连续系统的$\frac{2}{T}\tan \left(\frac{\omega T}{2}\right)$，利用$\frac{2}{T}\tan \left(\frac{\omega T}{2}\right)=2\pi \ast 250$求解离散系统的原始$\omega$，解得$\omega =1331.5$ rad/s，这与上述波特图结果是一致的。

  解决方案：离散系统的${\omega }_d$会映射到连续系统的${\omega }_a=\frac{2}{T}\tan \left(\frac{{\omega }_{d}T}{2}\right)$，基于${\omega }_a$完成连续传递函数的设计，采用双线性变换得到目标离散传递函数。正确的设计思路为：（1）首先，明确所需要的离散传递函数的特性，如希望在${\omega }_d=2\pi \ast 250$ rad/s时，离散传递函数取最大值；（2）计算上述离散频率所映射的连续频率，${\omega }_a=\frac{2}{T}\tan \left(\frac{{\omega }_{d}T}{2}\right)=2000$ rad/s = 318 Hz；(3) 基于映射的连续角频率${\omega }_a$设计连续传递函数，$C_a(s)=\frac{K_{i}s}{s^2+2000^2}$，并完成离散化，$C_d(z)=C_a(s){|}_{s=\frac{2}{T}\frac{z-1}{z+1}}=\frac{K_i(\frac{2}{T}\frac{z-1}{z+1})}{{\left(\frac{2}{T}\frac{z-1}{z+1}\right)}^2+2000^2}$。

  基于如下代码重绘波特图。可以看出，尽管连续传递函数在2000 rad/s (318Hz)处实现最大增益，离散传递函数在目标频率点1570 rad/s （250Hz）处取得最大增益，与预期一致。

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

6. 总结

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