Matlab 连续到离散变换：任意离散化方式

Matlab 用 c2d 函数进行连续时间到离散时间的变换。c2d 提供了 zoh foh impulse tustin matched least-squares 几种离散化方法，能满足大多数应用需求。但有时你可能希望实现自己的离散化方法，这正是本文要讲的内容。

c2d 函数

常用的离散化方法是 ZOH（零阶保持）和 Tustin 近似（双线性近似），在 c2d 函数中分别称为 zoh 和 tustin。c2d 的语法是

sysd = c2d(sys,Ts)
sysd = c2d(sys,Ts,method)
sysd = c2d(sys,Ts,opts)
[sysd,G] = c2d(sys,Ts,method)
[sysd,G] = c2d(sys,Ts,opts)

ZOH（后向差分）

z=\cfrac{1}{1-sT_s}

s=\cfrac{1-z^{-1}}{T_s}

Tustin 近似

不做频率预畸变

在连续时间系统和离散化系统之间，Tustin（双线性）近似在频域上的匹配效果最好。

z=e^{sT_s}\approx\cfrac{1+sT_s/2}{1-sT_s/2}

s=\cfrac{2}{T_s}\cfrac{1-z^{-1}}{1+z^{-1}}

做频率预畸变

这种方法保证连续时间响应和离散时间响应在预畸变频率处相匹配。

s=\cfrac{\omega_{prewarp}}{\tan(\omega_{prewarp}T_s/2)}\cfrac{1-z^{-1}}{1+z^{-1}}

任意离散化

我们可以用 Matlab 的符号运算实现任意离散化。以一阶低通滤波器（LPF）为例

H_{LPF}(s) = \cfrac{\omega_c}{s+\omega_c}

假设 \omega_c=200\pi\ \mathrm{rad/s}，T_s=10^{-3}\ s，用 c2d

omega_s = 200*pi;
T_s     = 1e-3;
Hs_LPF  = tf(omega_s, [1 omega_s]);
Hz_LPF  = c2d(Hs_LPF,T_s,'tustin');

结果是

Hz_LPF =

  0.2391 z + 0.2391
  -----------------
     z - 0.5219

Sample time: 0.001 seconds
Discrete-time transfer function.

符号运算

syms s z
Ss_LPF = omega_s / (s + omega_s);
Sz_LPF = subs(Ss_LPF,s,2/T_s * (1-z^-1) / (1+z^-1));
[Snum,Sden] = numden(Sz_LPF);

num = flip(eval(coeffs(Snum)));
den = flip(eval(coeffs(Sden)));

num = num./den(1);
den = den./den(1);

Hz_LPF = tf(num,den,T_s);

得到的 Hz_LPF 与 c2d 的结果相同。

Hz_LPF =

  0.2391 z + 0.2391
  -----------------
     z - 0.5219

Sample time: 0.001 seconds
Discrete-time transfer function.

说明

subs

符号替换。subs(s,old,new) 返回 s 的一个副本，把其中所有的 old 替换为 new，然后对 s 求值。

numden

提取分子和分母。[N,D] = numden(A) 将 A 化为有理形式，其分子和分母是互质的整系数多项式。该函数返回表达式有理形式的分子和分母。

如果 A 是符号矩阵或数值矩阵，则 N 是由分子组成的符号矩阵，D 是由分母组成的符号矩阵。N 和 D 都是与 A 大小相同的矩阵。

coeffs

多项式的系数。C = coeffs(p) 返回多项式 p 关于 symvar 在 p 中确定的所有变量的系数。系数按从低次到高次的顺序排列。 例如：

syms x
c = coeffs(16*x^2 + 19*x + 11)
>> c = [ 11, 19, 16]

因此，需要用 flip 函数把结果转换成 tf(num,den) 的格式。