Bilinear Transform

The bilinear transform is a transformation from continuous-time systems (in the Laplace domain) to discrete-time systems (in the Z-domain). It uses the trapezoidal rule for numerical integration.

Derivation

Consider a continuous-time integrator $y(t)={\int }_0^{t}u(\tau )d\tau$. (Where $y(t)$ is the output of the integrator for input $u(t)$).
In the Laplace domain, we can represent these signals as follows: $$U(s)\triangleq \mathcal{L}\{u(t)\}$$ $$Y(s)\triangleq \mathcal{L}\{y(t)\}=\mathcal{L}\{{\int }_0^{t}u(\tau )d\tau \}=\frac{1}{s}U(s)$$ Therefore, the transfer function of an integrator is: $$\begin{array}{}\text{(1)}& H_a(s)\triangleq \frac{Y(s)}{U(s)}=\frac{1}{s}\frac{U(s)}{U(s)}=\frac{1}{s}\end{array}$$

If we evaluate $y(t)$ at a specific time $t=k{T}_s$, $y(t)$ can be written as follows: $$y[k]\triangleq y(k{T}_s)={\int }_0^{k{T}_s}u(\tau )d\tau$$ By splitting up the integral: $$y[k]={\int }_0^{(k-1)T_s}u(\tau )d\tau +{\int }_{(k-1)T_s}^{k{T}_s}u(\tau )d\tau$$ This results in a recursive formula for $y[k]$: $$\begin{array}{}\text{(2)}& y[k]=y[k-1]+{\int }_{(k-1)T_s}^{k{T}_s}u(\tau )d\tau \end{array}$$ The second term expresses the area under the curve of $u(t)$ between $t=(k-1)T_s$ and $t=k{T}_s$. It can be approximated in several different ways, using numerical integration techniques. In the case of the bilinear transform, the trapezoidal rule is used, with a step size (sample time) of $T_s$. This is demonstrated in the following figure.

The light red area on the left is the first term in Equation $\text{(2)}$, and the darker red trapezoid on the right approximates the second term of the equation.
The area of this trapezoid is given by: $$A=T_s\frac{u((k-1)T_s)+u(k{T}_s)}{2}=\frac{T_s}{2}(u[k-1]+u[k])$$ $$u[k]\triangleq u(k{T}_s)$$ Therefore, we can approximate Equation $\text{(2)}$: $$\begin{array}{}\text{(3)}& y[k]\approx y[k-1]+\frac{T_s}{2}(u[k-1]+u[k])\end{array}$$ We now have an approximation of the recurrence relation for the integrator.

In Equation $\text{(1)}$, we found that the transfer function of an integrator in the Laplace domain was $H_a(s)=\frac{1}{s}$. We can now apply the Z-transform to the recurrence relation in Equation $\text{(3)}$ to relate the Laplace domain to the Z-domain. $$\begin{array}{rl}& \mathcal{Z}\{y[k]\approx y[k-1]+\frac{T_s}{2}(u[k-1]+u[k])\}\\ \iff & Y(z)\approx z^{-1}Y(z)+\frac{T_s}{2}(z^{-1}U(z)+U(z))\\ \iff & (1-z^{-1})Y(z)\approx \frac{T_s}{2}(1+z^{-1})U(z)\end{array}$$ $$\begin{array}{}\text{(4)}& H_d(z)\triangleq \frac{Y(z)}{U(z)}\approx \frac{T_s}{2}\frac{z+1}{z-1}\end{array}$$ If we compare the continuous-time transfer function $H_a(s)$ from Equation $\text{(1)}$ to the approximated discrete-time transfer function $H_d(z)$ from Equation $\text{(4)}$, we find an approximation of $s$ in function of $z$, that can be used to discretize a continuous-time transfer function: $$s\leftrightarrow \frac{2}{T_s}\frac{z-1}{z+1}$$ We can also express $z$ in function of $s$: $$\begin{array}{r}s\leftrightarrow \frac{2}{T_s}\frac{z-1}{z+1}\\ sz+s\leftrightarrow \frac{2}{T_s}z-\frac{2}{T_s}\\ z(s-\frac{2}{T_s})\leftrightarrow -s-\frac{2}{T_s}\\ z\leftrightarrow \frac{-s-\frac{2}{T_s}}{s-\frac{2}{T_s}}\\ z\leftrightarrow \frac{1+\frac{T_s}{2}s}{1-\frac{T_s}{2}s}\end{array}$$

Alternatively, the identity $z=e^{s{T}_s}$ can be used, together with the Taylor expansion of $e^x=1+x+\frac{x^2}{2!}+\frac{x^3}{3!}+\dots$: $$\begin{array}{rl}z& =e^{s{T}_s}\\ & =\frac{e^{s{T}_s/2}}{e^{-s{T}_s/2}}\\ & \approx \frac{1+\frac{T_s}{2}s}{1-\frac{T_s}{2}s}\end{array}$$

Stability

In the Laplace domain, poles that lie to the right of the imaginary axis are unstable. In the Z-domain, poles that lie outside of the unit circle are unstable.
We wish to determine what happens to the stability of poles when applying the bilinear transform.
Let's look at the image of the imaginary axis $(s=j\omega )$: $$\begin{array}{rl}z=& \frac{1+\frac{T_s}{2}j\omega }{1-\frac{T_s}{2}j\omega }\\ =& \frac{1+\frac{T_s}{2}j\omega }{1-\frac{T_s}{2}j\omega }\cdot \frac{1+\frac{T_s}{2}j\omega }{1+\frac{T_s}{2}j\omega }\\ =& \frac{{(1+\frac{T_s}{2}j\omega )}^2}{1+{\left(\frac{T_s}{2}\omega \right)}^2}\\ =& \frac{1+2\frac{T_s}{2}j\omega -{\left(\frac{T_s}{2}\omega \right)}^2}{1+{\left(\frac{T_s}{2}\omega \right)}^2}\end{array}$$ $$\begin{array}{}\text{(5)}& z=\frac{1-{\left(\frac{T_s}{2}\omega \right)}^2}{1+{\left(\frac{T_s}{2}\omega \right)}^2}+j\frac{2\frac{T_s}{2}\omega }{1+{\left(\frac{T_s}{2}\omega \right)}^2}\end{array}$$ We'll now calculate ${\left|z\right|}^2$, you'll see why in a moment. $$\begin{array}{rl}{\left|z\right|}^2=& \mathrm{\Re }(z{)}^2+\mathrm{\Im }(z{)}^2\\ =& \frac{{(1-{\left(\frac{T_s}{2}\omega \right)}^2)}^2}{{(1+{\left(\frac{T_s}{2}\omega \right)}^2)}^2}+\frac{{\left(2\frac{T_s}{2}\omega \right)}^2}{{(1+{\left(\frac{T_s}{2}\omega \right)}^2)}^2}\\ =& \frac{1-2{\left(\frac{T_s}{2}\omega \right)}^2+{\left(\frac{T_s}{2}\omega \right)}^4+4{\left(\frac{T_s}{2}\omega \right)}^2}{{(1+{\left(\frac{T_s}{2}\omega \right)}^2)}^2}\\ =& \frac{1+2{\left(\frac{T_s}{2}\omega \right)}^2+{\left(\frac{T_s}{2}\omega \right)}^4}{{(1+{\left(\frac{T_s}{2}\omega \right)}^2)}^2}\\ =& \frac{{(1+{\left(\frac{T_s}{2}\omega \right)}^2)}^2}{{(1+{\left(\frac{T_s}{2}\omega \right)}^2)}^2}\\ {\left|z\right|}^2=& 1\end{array}$$ In other words, the imaginary axis in the s-plane (Laplace domain) maps to the unit circle in the z-plane (Z-domain).

Let's now look at the image of the real axis $(s=\sigma )$: $$z=\frac{1+\frac{T_s}{2}\sigma }{1-\frac{T_s}{2}\sigma }$$ Substituting $x\triangleq \frac{T_s}{2}\sigma$ results in: $$z=\frac{1+x}{1-x}$$ Note that $T_s$ (the sample time) is positive by definition, so $x$ is positive when $\sigma$ is positive, and $x$ is negative when $\sigma$ is negative. As you can see, the right half of the real axis $(\sigma <0)$ is mapped to the interval $z\in (-1,1)$, so inside the unit circle; the left half of the real axis $(\sigma >0)$ is mapped to the interval $z\in (-\mathrm{\infty },-1)\cup (1,+\mathrm{\infty })$.

We can conclude that the left half of the s-plane is mapped to the part of the z-plane inside of the unit circle, and the right half of the s-plane is mapped to the part of the z-plane outside of the unit circle.
This has a very important consequence: stable poles in the Laplace domain map to stable poles in the Z-domain, and unstable poles in the Laplace domain map to unstable poles in the Z-domain.

Frequency response

In the Laplace domain, we determine the frequency response of a system by evaluating the transfer function at $s=j{\omega }_a$. In the Z-domain, on the other hand, we evaluate the transfer function at $z=e^{j{\omega }_d}$.

When designing a filter in the Laplace domain with a certain corner-frequency, we want the corner-frequency to be the same after discretization.
However, we'll find that there is no linear mapping from the frequency response of the continuous-time system to the frequency response of the discrete-time system.

In the previous section, we already found that the imaginary axis of the s-plane maps to the unit circle of the z-plane.
The phase of any $z$ on the unit circle in the Z-domain is the frequency ${\omega }_d$ $\left(\left[\text{rad}/\text{sample}\right]\right)$ of a specific point in the frequency response of the discrete-time system.
The imaginary part of any $s$ on the imaginary axis in the Laplace domain is the frequency ${\omega }_a$ $\left(\left[\text{rad}/s\right]\right)$ of a specific point in the frequency response of the continuous-time system.
$${\omega }_d\triangleq \mathrm{\angle }z=\arctan \left(\frac{\mathrm{\Im }(z)}{\mathrm{\Re }(z)}\right)$$ We'll be using Equation $\text{(5)}$, substituting $\alpha \triangleq \frac{T_s}{2}$ to make it easier to read. $$\begin{array}{rl}& {\omega }_d=\arctan \left(\frac{2\alpha {\omega }_a}{1-{\alpha }^2{\omega }_a^2}\right)\\ \iff & \tan \left({\omega }_d\right)=\frac{2\alpha {\omega }_a}{1-{\alpha }^2{\omega }_a^2}\\ \iff & (1-{\alpha }^2{\omega }_a^2)\tan \left({\omega }_d\right)=2\alpha {\omega }_a\\ \iff & {\alpha }^2\tan \left({\omega }_d\right){\omega }_a^2+2\alpha {\omega }_a-\tan \left({\omega }_d\right)=0\end{array}$$ We can solve this quadratic polynomial in ${\omega }_a$ using the discriminant method: $$\begin{array}{rl}\mathrm{\Delta }& =4{\alpha }^2+4{\alpha }^2{\tan }^2{\omega }_d\\ & =4{\alpha }^2(1+{\tan }^2{\omega }_d)\\ & =4{\alpha }^2{\sec }^2{\omega }_d\end{array}$$ $$\begin{array}{rl}{\omega }_a& =\frac{-2\alpha \pm \sqrt{\mathrm{\Delta }}}{2{\alpha }^2\tan {\omega }_d}\\ & =\frac{-2\alpha \pm 2\alpha \sec {\omega }_d}{2{\alpha }^2\tan {\omega }_d}\\ & =\frac{-1}{\alpha }\cdot \frac{1\mp \sec {\omega }_d}{\tan {\omega }_d}\\ & =\frac{-1}{\alpha }\cdot \frac{\cos {\omega }_d\mp 1}{\sin {\omega }_d}\\ & =\frac{-1}{\alpha }\cdot \frac{2{\cos }^2\frac{{\omega }_d}{2}-1\mp 1}{2\sin \frac{{\omega }_d}{2}\cos \frac{{\omega }_d}{2}}\\ {\omega }_{a,1}& =\frac{-1}{\alpha }\cdot \frac{2{\cos }^2\frac{{\omega }_d}{2}-2}{2\sin \frac{{\omega }_d}{2}\cos \frac{{\omega }_d}{2}}\\ & =\frac{-1}{\alpha }\cdot \frac{2({\cos }^2\frac{{\omega }_d}{2}-1)}{2\sin \frac{{\omega }_d}{2}\cos \frac{{\omega }_d}{2}}\\ & =\frac{1}{\alpha }\cdot \frac{{\sin }^2\frac{{\omega }_d}{2}}{\sin \frac{{\omega }_d}{2}\cos \frac{{\omega }_d}{2}}\\ & =\frac{1}{\alpha }\cdot \frac{\sin \frac{{\omega }_d}{2}}{\cos \frac{{\omega }_d}{2}}\\ & =\frac{1}{\alpha }\tan \frac{{\omega }_d}{2}\\ & =\frac{2}{T_s}\tan \frac{{\omega }_d}{2}\\ {\omega }_{a,2}& =\frac{-1}{\alpha }\cdot \frac{2{\cos }^2\frac{{\omega }_d}{2}}{2\sin \frac{{\omega }_d}{2}\cos \frac{{\omega }_d}{2}}\\ & =\frac{-1}{\alpha }\cdot \frac{\cos \frac{{\omega }_d}{2}}{\sin \frac{{\omega }_d}{2}}\\ & =\frac{-1}{\alpha }\cot \frac{{\omega }_d}{2}\\ & =\frac{-2}{T_s}\cot \frac{{\omega }_d}{2}\end{array}$$ I used some trigonometric double-angle identities to simplify the result.

We get two possible solutions for ${\omega }_a$, however, only one is correct.
First, let's look at the case where ${\omega }_a\ge 0$: in that case, $\mathrm{\Im }(z)$ is always positive (see Equation $\text{(5)}$), so ${\omega }_d\triangleq \mathrm{\angle }z\in [0,\pi ]$.
The first solution ${\omega }_{a,1}=\frac{2}{T_s}\tan \frac{{\omega }_d}{2}$ is positive, because the tangent function is positive for the argument $\frac{{\omega }_d}{2}\in [0,\pi /2]$.
The second solution ${\omega }_{a,2}=\frac{-2}{T_s}\cot \frac{{\omega }_d}{2}$ is negative, because the cotangent function is positive for the argument $\frac{{\omega }_d}{2}\in [0,\pi /2]$. This is in contradiction to our initial assumption that ${\omega }_a\ge 0$.
Now consider the case where ${\omega }_a\le 0$: in that case, $\mathrm{\Im }(z)$ is always negative, so ${\omega }_d\triangleq \mathrm{\angle }z\in [0,-\pi ]$.
The first solution ${\omega }_{a,1}=\frac{2}{T_s}\tan \frac{{\omega }_d}{2}$ is negative, because the tangent function is negative for the argument $\frac{{\omega }_d}{2}\in [0,-\pi /2]$.
The second solution ${\omega }_{a,2}=\frac{-2}{T_s}\cot \frac{{\omega }_d}{2}$ is positive, because the cotangent function is negative for the argument $\frac{{\omega }_d}{2}\in [0,-\pi /2]$. This is in contradiction to our initial assumption that ${\omega }_a\le 0$.

In conclusion, if we want our digital system to have a specific characteristic at a given frequency ${\omega }_d$, we have to use a different frequency ${\omega }_a$ when designing the system in the analog domain. This action is called frequency prewarping. The prewarped frequency is given by: $$\begin{array}{}\text{(6)}& {\omega }_a=\frac{2}{T_s}\tan \frac{{\omega }_d}{2}\end{array}$$ Where ${\omega }_d$ is the normalized angular frequency in radians per sample. The actual frequency $f_d$ is given by: $$f_d=\frac{{\omega }_d}{2\pi T_s}$$ Note the units of the quantities involved: ${\omega }_d$ is given in $\text{rad}/\text{sample}$, dividing by $2\pi$ cancels out the radians. The unit of $T_s$ is $s/\text{sample}$, the samples cancel out, and the result $f_d$ is given in $Hz$ or $s^{-1}$.