双线性变换

在数字信号处理和离散时间的控制理论中，双线性变换（也称为塔斯廷[Tustin]变换法，以阿诺德·塔斯廷命名）用于在连续时间系统和离散时间系统的表示之间进行转换。

双线性变换是共形映射（具体来说是莫比乌斯变换）的一种特例，常用于将连续时域中的线性时不变滤波器（通常被称为模拟滤波器）的传递函数$H_{a}(s)$转换为离散时域中的线性移不变滤波器（通常被称为数字滤波器，但有些使用开关电容构建的模拟滤波器实际上也是离散时间滤波器）的传递函数$H_{d}(z)$。双线性变换将s平面中虚轴上的点$j\omega$（$\mathrm {Re} [s]=0$）映射到z平面中的单位圆 $|z|=1$。此外，双线性变换可用于对离散时间线性系统的频率响应进行扭曲（例如用于逼近人类听觉系统的非线性频率分辨率），并可在离散域中通过替换系统的单位延迟$z^{-1}$来实现一阶全通滤波器。

双线性变换能够保持系统的稳定性，并将连续时间滤波器的频率响应$H_{a}(j\omega _{a})$的每一点映射到离散时间滤波器的频率响应 $H_{d}(e^{j\omega _{d}T})$中的对应点，尽管映射后频率略有不同（具体参见下文的频率扭曲部分）。这意味着，在模拟滤波器的频率响应中看到的每个特性，在数字滤波器的频率响应中都有一个对应的特性，其增益和相移完全相同，但频率可能略有不同。频率的变化在低频时几乎不可察觉，但在接近奈奎斯特频率时非常明显。

离散时间估计

双线性变换是自然对数函数的一阶估计法，也就是将z平面映射到s平面，当拉普拉斯变换被用在离散时间信号上（将离散时间串行中的每个元素附在对应的延迟狄拉克δ函数），其结果确实为将离散时间串行的Z变换替代成

${\begin{aligned}z&=e^{sT}\\&={\frac {e^{sT/2}}{e^{-sT/2}}}\\&\approx {\frac {1+sT/2}{1-sT/2}}\end{aligned}}$

其中$T\$是用在推导双线性变换的梯形公式中数值积分每阶的大小，换句话说就是采样间距。上述的双线性估计可以透过$s\$来解或是产生一个近似估计$s=(1/T)\ln(z)\ \$。

逆映射则为

${\begin{aligned}s&={\frac {1}{T}}\ln(z)\\&={\frac {2}{T}}\left[{\frac {z-1}{z+1}}+{\frac {1}{3}}\left({\frac {z-1}{z+1}}\right)^{3}+{\frac {1}{5}}\left({\frac {z-1}{z+1}}\right)^{5}+{\frac {1}{7}}\left({\frac {z-1}{z+1}}\right)^{7}+\cdots \right]\\&\approx {\frac {2}{T}}{\frac {z-1}{z+1}}\\&={\frac {2}{T}}{\frac {1-z^{-1}}{1+z^{-1}}}\end{aligned}}$

双线性变换的本质是使用这种一阶估计法且将连续时间传递函数$H_{a}(s)\$中的$s\$替换成

$s\leftarrow {\frac {2}{T}}{\frac {z-1}{z+1}}.$

也就是说

$H_{d}(z)=H_{a}(s){\bigg |}_{s={\frac {2}{T}}{\frac {z-1}{z+1}}}=H_{a}\left({\frac {2}{T}}{\frac {z-1}{z+1}}\right).\$

保留稳定性及最小相位性质

如果有一个连续时间且有因果性的滤波器，其传递函数的极点落在复数S平面的左半边，此滤波器则为稳定的。如果有一个离散时间且有因果性的滤波器，其传递函数的极点落在复数Z平面的单位圆内，此滤波器则为稳定的。双线性变换将复数S平面的左半边映射到复数Z平面的单位圆内，因此稳定的连续时间滤波器被转变成离散时间滤波器后也保有稳定性。

同样地，如果有一个连续时间的滤波器，其传递函数的零点落在复数S平面的左半边，此滤波器则有最小相位性质。如果有一个离散时间且有因果性的滤波器，其传递函数的零点落在复数Z平面的单位圆内，此滤波器则有最小相位性质。透过相同的映射性质，可以保证有最小相位性质的连续时间滤波器被变换成离散时间滤波器后也保有最小性质。

例子

以一个简单的低通RC电路当例子，这种连续时间滤波器的传递函数为

${\begin{aligned}H_{a}(s)&={\frac {1/sC}{R+1/sC}}\\&={\frac {1}{1+RCs}}.\end{aligned}}$

如果我们想将这种滤波器应用成数字滤波器，我们可以将上式中的$s$做替换，因此可以得到下列表示式

$H_{d}(z)\$	$=H_{a}\left({\frac {2}{T}}{\frac {z-1}{z+1}}\right)\$
	$={\frac {1}{1+RC\left({\frac {2}{T}}{\frac {z-1}{z+1}}\right)}}\$
	$={\frac {1+z}{(1-2RC/T)+(1+2RC/T)z}}\$
	$={\frac {1+z^{-1}}{(1+2RC/T)+(1-2RC/T)z^{-1}}}.\$

在应用在即时数字滤波器时，分母的系数为’反馈系数’而分子的系数为’前馈系数’。

一般的双二阶变换

将连续时间的模拟滤波器的系数对应到由双线性变换展成的相似的离散时间数字滤波器是有可能的，假设有一个传递函数为下式的一般二阶连续时间滤波器

$H_{a}(s)={\frac {b_{0}s^{2}+b_{1}s+b_{2}}{a_{0}s^{2}+a_{1}s+a_{2}}}={\frac {b_{0}+b_{1}s^{-1}+b_{2}s^{-2}}{a_{0}+a_{1}s^{-1}+a_{2}s^{-2}}}$

利用下列替换方法做双线性变换

$s\leftarrow K{\frac {1-z^{-1}}{1+z^{-1}}}$

其中$K\triangleq {\frac {2}{T}}$.

其结果为一个离散时间的数字双二阶滤波器，且由原本连续时间滤波器的系数所组成的表达式如

$H_{d}(z)={\frac {(b_{0}K^{2}+b_{1}K+b_{2})+(2b_{2}-2b_{0}K^{2})z^{-1}+(b_{0}K^{2}-b_{1}K+b_{2})z^{-2}}{(a_{0}K^{2}+a_{1}K+a_{2})+(2a_{2}-2a_{0}K^{2})z^{-1}+(a_{0}K^{2}-a_{1}K+a_{2})z^{-2}}}$

一般而言，在推导对应的差分方程序前，分母的常数项会被标准化为1

$H_{d}(z)={\frac {{\frac {b_{0}K^{2}+b_{1}K+b_{2}}{a_{0}K^{2}+a_{1}K+a_{2}}}+{\frac {2b_{2}-2b_{0}K^{2}}{a_{0}K^{2}+a_{1}K+a_{2}}}z^{-1}+{\frac {b_{0}K^{2}-b_{1}K+b_{2}}{a_{0}K^{2}+a_{1}K+a_{2}}}z^{-2}}{1+{\frac {2a_{2}-2a_{0}K^{2}}{a_{0}K^{2}+a_{1}K+a_{2}}}z^{-1}+{\frac {a_{0}K^{2}-a_{1}K+a_{2}}{a_{0}K^{2}+a_{1}K+a_{2}}}z^{-2}}}.$

差分方程序则为

${\begin{aligned}y[n]={}&{\frac {b_{0}K^{2}+b_{1}K+b_{2}}{a_{0}K^{2}+a_{1}K+a_{2}}}\cdot x[n]+{\frac {2b_{2}-2b_{0}K^{2}}{a_{0}K^{2}+a_{1}K+a_{2}}}\cdot x[n-1]\\[8pt]&{}+{\frac {b_{0}K^{2}-b_{1}K+b_{2}}{a_{0}K^{2}+a_{1}K+a_{2}}}\cdot x[n-2]-{\frac {2a_{2}-2a_{0}K^{2}}{a_{0}K^{2}+a_{1}K+a_{2}}}\cdot y[n-1]\\[8pt]&{}-{\frac {a_{0}K^{2}-a_{1}K+a_{2}}{a_{0}K^{2}+a_{1}K+a_{2}}}\cdot y[n-2].\end{aligned}}$

频率扭曲

为了要计算连续时间滤波器的频率响应，会去计算传递函数$H_{a}(s)\$在$j\omega \$轴上的值（$s=j\omega \$）。相同地，为了要计算离散时间滤波器的频率响应，会去计算在单位圆上传递函数$H_{d}(z)\$的值（$z=e^{j\omega T}\$，$|z|=1\$）。当真正的频率$\omega \$被代入由双线性变换产生的离散时间滤波器，可以透过下列式子得到连续时间滤波器的频率$\omega _{a}\$

$H_{d}(z)=H_{a}\left({\frac {2}{T}}{\frac {z-1}{z+1}}\right)\$

$H_{d}(e^{j\omega T})\$	$=H_{a}\left({\frac {2}{T}}{\frac {e^{j\omega T}-1}{e^{j\omega T}+1}}\right)\$
	$=H_{a}\left({\frac {2}{T}}\cdot {\frac {e^{j\omega T/2}\left(e^{j\omega T/2}-e^{-j\omega T/2}\right)}{e^{j\omega T/2}\left(e^{j\omega T/2}+e^{-j\omega T/2}\right)}}\right)\$
	$=H_{a}\left({\frac {2}{T}}\cdot {\frac {\left(e^{j\omega T/2}-e^{-j\omega T/2}\right)}{\left(e^{j\omega T/2}+e^{-j\omega T/2}\right)}}\right)\$
	$=H_{a}\left(j{\frac {2}{T}}\cdot {\frac {\left(e^{j\omega T/2}-e^{-j\omega T/2}\right)/(2j)}{\left(e^{j\omega T/2}+e^{-j\omega T/2}\right)/2}}\right)\$
	$=H_{a}\left(j{\frac {2}{T}}\cdot {\frac {\sin(\omega T/2)}{\cos(\omega T/2)}}\right)\$
	$=H_{a}\left(j{\frac {2}{T}}\cdot \tan \left(\omega T/2\right)\right)\$

由此可知，离散时间滤波器在z平面单位圆中的每一点（$z=e^{j\omega T}\$）可以被映射到连续时间滤波器在s平面$j\omega \$轴上的一点（$s=j\omega _{a}\$）。也就是说，双线性变换将离散时间滤波器的频率映射到连续时间滤波器的方法为下式

$\omega _{a}={\frac {2}{T}}\tan \left(\omega {\frac {T}{2}}\right)$

反之则为

$\omega ={\frac {2}{T}}\arctan \left(\omega _{a}{\frac {T}{2}}\right).$

离散时间滤波器在频率为$\omega \$的表现和连续时间滤波器在频率为$(2/T)\tan(\omega T/2)\$的表现相同，具体来说，离散时间滤波器在频率为$\omega \$的增益和相位平移与连续时间滤波器在频率为$(2/T)\tan(\omega T/2)\$的增益和相位平移相同。也就是说，在连续时间滤波器的频率响应所看到的每一个特征，都可以在离散时间滤波器得频率响应中看到，但频率位置可能会不同。对于低频而言（也就是当$\omega \ll 2/T$或$\omega _{a}\ll 2/T$），$\omega \approx \omega _{a}\$。

连续时间滤波器的频率范围是

$-\infty <\omega _{a}<+\infty \$

对应到在离散时间滤波器的频率区间是

$-{\frac {\pi }{T}}<\omega <+{\frac {\pi }{T}}.\$

当连续时间滤波器的频率$\omega _{a}=0\$，对应到离散时间滤波器的频率$\omega =0\$；当连续时间滤波器的频率<$\omega _{a}=\pm \infty \$，对应到离散时间滤波器的频率$\omega =\pm \pi /T.\$

可以看到$\omega _{a}\$和$\omega .\$之间是非线性的关系，这个由双线性变换产生的影响称为频率扭曲。设计连续时间滤波器时可以透过设定$\omega _{a}={\frac {2}{T}}\tan \left(\omega {\frac {T}{2}}\right)\$来补偿频率扭曲，这在滤波器设计中称作为预先扭曲。

当设计一个数字滤波器去估计连续时间滤波器时，如果将下列变换式代入连续时间滤波器的传递函数中，这个数字滤波器的频率响应（包含幅度跟相位）可以被做成符合连续时间滤波器在$\omega _{0}$的频率响应，这是一种修改过的Tustin变换。然而，当$\omega _{0}\to 0$时，这种变换方式就会变成上面所说的Tustin变换。也就是说，上面的变换使得数字滤波器的响应在直流分量时会对应到模拟滤波器响应

$s\leftarrow {\frac {\omega _{0}}{\tan({\frac {\omega _{0}T}{2}})}}{\frac {z-1}{z+1}}.$

这种扭曲现象的主要优点是去除频率响应的混叠失真。然而，还需要透过预先扭曲给定的连续时间系统频率能补偿所造成的频率扭曲，这些被预先扭曲的频率用在双线性变换上可以得到想要的离散时间系统。

参见

• Z变换
• S变换