bilinear transform

The bilinear transform is a first-order Padé approximant of the natural logarithm function that is an exact mapping of the z-plane to the s-plane. When the Laplace transform is performed on a discrete-time signal (with each element of the discrete-time sequence attached to a correspondingly delayed unit impulse), the result is precisely the Z transform of the discrete-time sequence with the substitution of

:$$

\begin{align}

z &= e^{sT} \\

&= \frac{e^{sT/2}}{e^{-sT/2}} \\

&\approx \frac{1 + s T / 2}{1 - s T / 2}

\end{align}

where $T$ is the numerical integration step size of the trapezoidal rule used in the bilinear transform derivation;{{cite book |title=Discrete Time Signal Processing Third Edition |last=Oppenheim |first=Alan |year=2010 |publisher=Pearson Higher Education, Inc. |location=Upper Saddle River, NJ |isbn=978-0-13-198842-2 |page=504}} or, in other words, the sampling period. The above bilinear approximation can be solved for $s$ or a similar approximation for $s\; =\; (1/T)\; \backslash ln(z)$ can be performed.

The inverse of this mapping (and its first-order bilinear approximation) is

:$$

\begin{align}

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{align}

The bilinear transform essentially uses this first order approximation and substitutes into the continuous-time transfer function, $H\_a(s)$

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

That is

:$H\_d(z)\; =\; H\_a(s)\; \backslash bigg|\_\{s\; =\; \backslash frac\{2\}\{T\}\; \backslash frac\{z\; -\; 1\}\{z\; +\; 1\}\}=\; H\_a\; \backslash left(\; \backslash frac\{2\}\{T\}\; \backslash frac\{z-1\}\{z+1\}\; \backslash right).\; \backslash$

A continuous-time causal filter is stable if the poles of its transfer function fall in the left half of the complex s-plane. A discrete-time causal filter is stable if the poles of its transfer function fall inside the unit circle in the complex z-plane. The bilinear transform maps the left half of the complex s-plane to the interior of the unit circle in the z-plane. Thus, filters designed in the continuous-time domain that are stable are converted to filters in the discrete-time domain that preserve that stability.

Likewise, a continuous-time filter is minimum-phase if the zeros of its transfer function fall in the left half of the complex s-plane. A discrete-time filter is minimum-phase if the zeros of its transfer function fall inside the unit circle in the complex z-plane. Then the same mapping property assures that continuous-time filters that are minimum-phase are converted to discrete-time filters that preserve that property of being minimum-phase.

A general LTI system has the transfer function

H_a(s) = \frac{b_0 + b_1s + b_2s^2 + \cdots + b_Qs^Q}{a_0 + a_1s + a_2s^2 + \cdots + a_Ps^P}

The order of the transfer function {{math|N}} is the greater of {{math|P}} and {{math|Q}} (in practice this is most likely {{math|P}} as the transfer function must be proper for the system to be stable). Applying the bilinear transform

s = K\frac{z - 1}{z + 1}

where {{math|K}} is defined as either {{math|2/T}} or otherwise if using frequency warping, gives

H_d(z) = \frac{b_0 + b_1\left(K\frac{z - 1}{z + 1}\right) + b_2\left(K\frac{z - 1}{z + 1}\right)^2 + \cdots + b_Q\left(K\frac{z - 1}{z + 1}\right)^Q}

{a_0 + a_1\left(K\frac{z - 1}{z + 1}\right) + a_2\left(K\frac{z - 1}{z + 1}\right)^2 + \cdots + b_P\left(K\frac{z - 1}{z + 1}\right)^P}

Multiplying the numerator and denominator by the largest power of {{math|(z + 1)⁻¹}} present, {{math|(z + 1)^−N}}, gives

H_d(z) = \frac{b_0(z+1)^N + b_1K(z-1)(z+1)^{N-1} + b_2K^2(z-1)^2(z+1)^{N-2} + \cdots + b_QK^Q(z-1)^Q(z+1)^{N-Q}}

{a_0(z+1)^N + a_1K(z-1)(z+1)^{N-1} + a_2K^2(z-1)^2(z+1)^{N-2} + \cdots + a_PK^P(z-1)^P(z+1)^{N-P}}

It can be seen here that after the transformation, the degree of the numerator and denominator are both {{math|N}}.

Consider then the pole-zero form of the continuous-time transfer function

H_a(s) = \frac{(s - \xi_1)(s - \xi_2) \cdots (s - \xi_Q)}{(s - p_1)(s - p_2) \cdots (s - p_P)}

The roots of the numerator and denominator polynomials, {{math|ξ_i}} and {{math|p_i}}, are the zeros and poles of the system. The bilinear transform is a one-to-one mapping, hence these can be transformed to the z-domain using

z = \frac{K + s}{K - s}

yielding some of the discretized transfer function's zeros and poles {{math|ξ'_i}} and {{math|p'_i}}

\begin{aligned}

\xi'_i &= \frac{K + \xi_i}{K - \xi_i} \quad 1 \leq i \leq Q \\

p'_i &= \frac{K + p_i}{K - p_i} \quad 1 \leq i \leq P

\end{aligned}

As described above, the degree of the numerator and denominator are now both {{math|N}}, in other words there is now an equal number of zeros and poles. The multiplication by {{math|(z + 1)^−N}} means the additional zeros or poles are

{{cite web|url=http://www.ee.ic.ac.uk/hp/staff/dmb/courses/DSPDF/00800_TransIIR.pdf

|last=Bhandari |first=Ayush

|title=DSP and Digital Filters Lecture Notes

|access-date=16 August 2022

|archive-url=https://web.archive.org/web/20220303144755/http://www.ee.ic.ac.uk/hp/staff/dmb/courses/DSPDF/00800_TransIIR.pdf

|archive-date=3 March 2022}}

\begin{aligned}

\xi'_i &= -1 \quad Q < i \leq N \\

p'_i &= -1 \quad P < i \leq N

\end{aligned}

Given the full set of zeros and poles, the z-domain transfer function is then

H_d(z) = \frac{(z - \xi'_1)(z - \xi'_2) \cdots (z - \xi'_N)}

{(z - p'_1)(z - p'_2) \cdots (z - p'_N)}

As an example take a simple low-pass RC filter. This continuous-time filter has a transfer function

:$\backslash begin\{align\}$

H_a(s) &= \frac{1/sC}{R+1/sC} \\

&= \frac{1}{1 + RC s}.

\end{align}

If we wish to implement this filter as a digital filter, we can apply the bilinear transform by substituting for $s$ the formula above; after some reworking, we get the following filter representation:

:

The coefficients of the denominator are the 'feed-backward' coefficients and the coefficients of the numerator are the 'feed-forward' coefficients used for implementing a real-time digital filter.

It is possible to relate the coefficients of a continuous-time, analog filter with those of a similar discrete-time digital filter created through the bilinear transform process. Transforming a general, first-order continuous-time filter with the given transfer function

:$H\_a(s)\; =\; \backslash frac\{b\_0\; s\; +\; b\_1\}\{a\_0\; s\; +\; a\_1\}\; =\; \backslash frac\{b\_0\; +\; b\_1\; s^\{-1\}\}\{a\_0\; +\; a\_1\; s^\{-1\}\}$

using the bilinear transform (without prewarping any frequency specification) requires the substitution of

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

where

:$K\; \backslash triangleq\; \backslash frac\{2\}\{T\}$.

However, if the frequency warping compensation as described below is used in the bilinear transform, so that both analog and digital filter gain and phase agree at frequency $\backslash omega\_0$, then

:$K\; \backslash triangleq\; \backslash frac\{\backslash omega\_0\}\{\backslash tan\backslash left(\backslash frac\{\backslash omega\_0\; T\}\{2\}\backslash right)\}$.

This results in a discrete-time digital filter with coefficients expressed in terms of the coefficients of the original continuous time filter:

:$H\_d(z)=\backslash frac\{(b\_0\; K\; +\; b\_1)\; +\; (-b\_0\; K\; +\; b\_1)z^\{-1\}\}\{(a\_0\; K\; +\; a\_1)\; +\; (-a\_0\; K\; +\; a\_1)z^\{-1\}\}$

Normally the constant term in the denominator must be normalized to 1 before deriving the corresponding difference equation. This results in

:$H\_d(z)=\backslash frac\{\backslash frac\{b\_0\; K\; +\; b\_1\}\{a\_0\; K\; +\; a\_1\}\; +\; \backslash frac\{-b\_0\; K\; +\; b\_1\}\{a\_0\; K\; +\; a\_1\}z^\{-1\}\}\{1\; +\; \backslash frac\{-a\_0\; K\; +\; a\_1\}\{a\_0\; K\; +\; a\_1\}z^\{-1\}\}.$

The difference equation (using the Direct form I) is

:$$

y[n] = \frac{b_0 K + b_1}{a_0 K + a_1} \cdot x[n] + \frac{-b_0 K + b_1}{a_0 K + a_1} \cdot x[n-1] - \frac{-a_0 K + a_1}{a_0 K + a_1} \cdot y[n-1] \ .

A similar process can be used for a general second-order filter with the given transfer function

:$H\_a(s)\; =\; \backslash frac\{b\_0\; s^2\; +\; b\_1\; s\; +\; b\_2\}\{a\_0\; s^2\; +\; a\_1\; s\; +\; a\_2\}\; =\; \backslash frac\{b\_0\; +\; b\_1\; s^\{-1\}\; +\; b\_2\; s^\{-2\}\}\{a\_0\; +\; a\_1\; s^\{-1\}\; +\; a\_2\; s^\{-2\}\}\; \backslash \; .$

This results in a discrete-time digital biquad filter with coefficients expressed in terms of the coefficients of the original continuous time filter:

:$H\_d(z)=\backslash 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\}\}$

Again, the constant term in the denominator is generally normalized to 1 before deriving the corresponding difference equation. This results in

:$H\_d(z)=\backslash frac\{\backslash frac\{b\_0\; K^2\; +\; b\_1\; K\; +\; b\_2\}\{a\_0\; K^2\; +\; a\_1\; K\; +\; a\_2\}\; +\; \backslash frac\{2b\_2\; -\; 2b\_0\; K^2\}\{a\_0\; K^2\; +\; a\_1\; K\; +\; a\_2\}z^\{-1\}\; +\; \backslash frac\{b\_0\; K^2\; -\; b\_1\; K\; +\; b\_2\}\{a\_0\; K^2\; +\; a\_1\; K\; +\; a\_2\}z^\{-2\}\}\{1\; +\; \backslash frac\{2a\_2\; -\; 2a\_0\; K^2\}\{a\_0\; K^2\; +\; a\_1\; K\; +\; a\_2\}z^\{-1\}\; +\; \backslash frac\{a\_0\; K^2\; -\; a\_1\; K\; +\; a\_2\}\{a\_0\; K^2\; +\; a\_1\; K\; +\; a\_2\}z^\{-2\}\}.$

The difference equation (using the Direct form I) is

:$$

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] + \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] - \frac{a_0 K^2 - a_1 K + a_2}{a_0 K^2 + a_1 K + a_2} \cdot y[n-2] \ .

To determine the frequency response of a continuous-time filter, the transfer function $H\_a(s)$ is evaluated at $s\; =\; j\; \backslash omega\_a$ which is on the $j\; \backslash omega$ axis. Likewise, to determine the frequency response of a discrete-time filter, the transfer function $H\_d(z)$ is evaluated at $z\; =\; e^\{\; j\; \backslash omega\_d\; T\}$ which is on the unit circle, $|z|\; =\; 1$. The bilinear transform maps the $j\; \backslash omega$ axis of the s-plane (which is the domain of $H\_a(s)$) to the unit circle of the z-plane, $|z|\; =\; 1$ (which is the domain of $H\_d(z)$), but it is not the same mapping $z\; =\; e^\{sT\}$ which also maps the $j\; \backslash omega$ axis to the unit circle. When the actual frequency of $\backslash omega\_d$ is input to the discrete-time filter designed by use of the bilinear transform, then it is desired to know at what frequency, $\backslash omega\_a$, for the continuous-time filter that this $\backslash omega\_d$ is mapped to.

:$H\_d(z)\; =\; H\_a\; \backslash left(\; \backslash frac\{2\}\{T\}\; \backslash frac\{z-1\}\{z+1\}\backslash right)$

:

This shows that every point on the unit circle in the discrete-time filter z-plane, $z\; =\; e^\{\; j\; \backslash omega\_d\; T\}$ is mapped to a point on the $j\; \backslash omega$ axis on the continuous-time filter s-plane, $s\; =\; j\; \backslash omega\_a$. That is, the discrete-time to continuous-time frequency mapping of the bilinear transform is

:$\backslash omega\_a\; =\; \backslash frac\{2\}\{T\}\; \backslash tan\; \backslash left(\; \backslash omega\_d\; \backslash frac\{T\}\{2\}\; \backslash right)$

and the inverse mapping is

:$\backslash omega\_d\; =\; \backslash frac\{2\}\{T\}\; \backslash arctan\; \backslash left(\; \backslash omega\_a\; \backslash frac\{T\}\{2\}\; \backslash right).$

The discrete-time filter behaves at frequency $\backslash omega\_d$ the same way that the continuous-time filter behaves at frequency $(2/T)\; \backslash tan(\backslash omega\_d\; T/2)$. Specifically, the gain and phase shift that the discrete-time filter has at frequency $\backslash omega\_d$ is the same gain and phase shift that the continuous-time filter has at frequency $(2/T)\; \backslash tan(\backslash omega\_d\; T/2)$. This means that every feature, every "bump" that is visible in the frequency response of the continuous-time filter is also visible in the discrete-time filter, but at a different frequency. For low frequencies (that is, when $\backslash omega\_d\; \backslash ll\; 2/T$ or $\backslash omega\_a\; \backslash ll\; 2/T$), then the features are mapped to a slightly different frequency; $\backslash omega\_d\; \backslash approx\; \backslash omega\_a$.

One can see that the entire continuous frequency range

:

is mapped onto the fundamental frequency interval

:

The continuous-time filter frequency $\backslash omega\_a\; =\; 0$ corresponds to the discrete-time filter frequency $\backslash omega\_d\; =\; 0$ and the continuous-time filter frequency $\backslash omega\_a\; =\; \backslash pm\; \backslash infty$ correspond to the discrete-time filter frequency $\backslash omega\_d\; =\; \backslash pm\; \backslash pi\; /\; T.$

One can also see that there is a nonlinear relationship between $\backslash omega\_a$ and $\backslash omega\_d.$ This effect of the bilinear transform is called frequency warping. The continuous-time filter can be designed to compensate for this frequency warping by setting $\backslash omega\_a\; =\; \backslash frac\{2\}\{T\}\; \backslash tan\; \backslash left(\; \backslash omega\_d\; \backslash frac\{T\}\{2\}\; \backslash right)$ for every frequency specification that the designer has control over (such as corner frequency or center frequency). This is called pre-warping the filter design.

It is possible, however, to compensate for the frequency warping by pre-warping a frequency specification $\backslash omega\_0$ (usually a resonant frequency or the frequency of the most significant feature of the frequency response) of the continuous-time system. These pre-warped specifications may then be used in the bilinear transform to obtain the desired discrete-time system. When designing a digital filter as an approximation of a continuous time filter, the frequency response (both amplitude and phase) of the digital filter can be made to match the frequency response of the continuous filter at a specified frequency $\backslash omega\_0$, as well as matching at DC, if the following transform is substituted into the continuous filter transfer function.{{cite book |last=Astrom |first=Karl J. |date=1990 |title=Computer Controlled Systems, Theory and Design |edition=Second |publisher=Prentice-Hall |page=212 |isbn=0-13-168600-3}} This is a modified version of Tustin's transform shown above.

:$s\; \backslash leftarrow\; \backslash frac\{\backslash omega\_0\}\{\backslash tan\backslash left(\backslash frac\{\backslash omega\_0\; T\}\{2\}\backslash right)\}\; \backslash frac\{z\; -\; 1\}\{z\; +\; 1\}.$

However, note that this transform becomes the original transform

:$s\; \backslash leftarrow\; \backslash frac\{2\}\{T\}\; \backslash frac\{z\; -\; 1\}\{z\; +\; 1\}$

as $\backslash omega\_0\; \backslash to\; 0$.

The main advantage of the warping phenomenon is the absence of aliasing distortion of the frequency response characteristic, such as observed with Impulse invariance.

• Impulse invariance
• Matched Z-transform method

{{reflist}}

• [http://ocw.mit.edu/courses/mechanical-engineering/2-161-signal-processing-continuous-and-discrete-fall-2008/lecture-notes/lecture_19.pdf MIT OpenCourseWare Signal Processing: Continuous to Discrete Filter Design]
• [http://web.cecs.pdx.edu/~tymerski/ece452/6.pdf Lecture Notes on Discrete Equivalents]
• [https://www.native-instruments.com/fileadmin/ni_media/downloads/pdf/VAFilterDesign_2.1.0.pdf#page=69 The Art of VA Filter Design]

{{DSP}}

{{DEFAULTSORT:Bilinear Transform}}

Category:Digital signal processing

Category:Transforms

Category:Control theory