雙線性轉換

在數位信號處理和離散時間的控制理論中，雙線性變換（也稱為塔斯廷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{align} z &= e^{sT} \\ <pre> &= \frac{e^{sT/2}}{e^{-sT/2}} \\ &\approx \frac{1 + s T / 2}{1 - s T / 2} </pre> \end{align}$$

其中$T \ $是用在推導雙線性變換的梯形公式中數值積分每階的大小，換句話說就是取樣間距。上述的雙線性估計可以透過$s \ $來解或是產生一個近似估計$s = (1/T) \ln(z) \ \ $。

逆映射則為

$$\begin{align} s &= \frac{1}{T} \ln(z) \\ <pre> &= \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}} </pre> \end{align}$$

雙線性變換的本質是使用這種一階估計法且將連續時間傳遞函數$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{align} H_a(s) &= \frac{1/sC}{R+1/sC} \\ &= \frac{1}{1 + RC s}. \end{align}$$

如果我們想將這種濾波器應用成數位濾波器，我們可以將上式中的$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 - 2 RC / 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{align} 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{align}$$

## 頻率扭曲

為了要計算連續時間濾波器的頻率響應，會去計算傳遞函數$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}.$$

這種扭曲現象的主要優點是去除頻率響應的混疊失真。然而，還需要透過預先扭曲給定的連續時間系統頻率能補償所造成的頻率扭曲，這些被預先扭曲的頻率用在雙線性變換上可以得到想要的離散時間系統。