Bilinear transform

The bilinear transform is a first-order approximation 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

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

where ${\displaystyle T}$ is the numerical integration step size of the trapezoidal rule used in the bilinear transform derivation;[1] or, in other words, the sampling period. The above bilinear approximation can be solved for ${\displaystyle s}$ or a similar approximation for ${\displaystyle s=(1/T)\ln(z)}$ can be performed.

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

${\displaystyle {\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}}}$

The bilinear transform essentially uses this first order approximation and substitutes into the continuous-time transfer function, ${\displaystyle H_{a}(s)}$

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

That is

${\displaystyle 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).\ }$

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.

Consider a continuous-time IIR filter of order ${\displaystyle N}$

${\displaystyle H_{a}(s)=k\prod _{i=1}^{N}{\frac {s-\xi _{i}}{s-p_{i}}},}$

where ${\displaystyle p_{i}}$ and ${\displaystyle \xi _{i}}$ are the transfer function poles and zeros in the s-plane. Let ${\displaystyle K=2/T}$ (or if using frequency warping as described below, let ${\displaystyle K=\omega _{0}/\tan(\omega _{0}T/2)}$).

The filter's bilinear transform is obtained by substituting ${\displaystyle s=K(z-1)/(z+1)}$:

${\displaystyle {\begin{aligned}H_{d}(z)&=H_{a}{\bigl (}K{\tfrac {z-1}{z+1}}{\bigr )}\\&=k\prod _{i=1}^{N}{\frac {K{\frac {z-1}{z+1}}-\xi _{i}}{K{\frac {z-1}{z+1}}-p_{i}}}\\&=k\prod _{i=1}^{N}{\frac {K-\xi _{i}}{K-p_{i}}}\cdot {\frac {z-{\frac {K+\xi _{i}}{K-\xi _{i}}}}{z-{\frac {K+p_{i}}{K-p_{i}}}}}\\&=H_{a}(K)\prod _{i=1}^{N}{\frac {z-\xi _{i}^{d}}{z-p_{i}^{d}}},\end{aligned}}}$

where ${\displaystyle p_{i}^{d}}$, ${\displaystyle \xi _{i}^{d}}$ are the z-plane pole and zero locations of the discretized filter,

${\displaystyle p_{i}^{d}={\frac {K+p_{i}}{K-p_{i}}},\quad \xi _{i}^{d}={\frac {K+\xi _{i}}{K-\xi _{i}}}.}$

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

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

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

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

The coefficients of the denominator are the 'feed-backward' coefficients and the coefficients of the numerator are the 'feed-forward' coefficients used to implement 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

${\displaystyle H_{a}(s)={\frac {b_{0}s+b_{1}}{a_{0}s+a_{1}}}={\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

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

where

${\displaystyle K\triangleq {\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 ${\displaystyle \omega _{0}}$, then

${\displaystyle K\triangleq {\frac {\omega _{0}}{\tan \left({\frac {\omega _{0}T}{2}}\right)}}}$.

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

${\displaystyle H_{d}(z)={\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

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

The difference equation (using the Direct Form I) is

${\displaystyle 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

${\displaystyle 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}}}\ .}$

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

${\displaystyle 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}}}}$

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

${\displaystyle 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}}}.}$

The difference equation (using the Direct form I) is

${\displaystyle 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 ${\displaystyle H_{a}(s)}$ is evaluated at ${\displaystyle s=j\omega _{a}}$ which is on the ${\displaystyle j\omega }$ axis. Likewise, to determine the frequency response of a discrete-time filter, the transfer function ${\displaystyle H_{d}(z)}$ is evaluated at ${\displaystyle z=e^{j\omega _{d}T}}$ which is on the unit circle, ${\displaystyle |z|=1}$. The bilinear transform maps the ${\displaystyle j\omega }$ axis of the s-plane (of which is the domain of ${\displaystyle H_{a}(s)}$) to the unit circle of the z-plane, ${\displaystyle |z|=1}$ (which is the domain of ${\displaystyle H_{d}(z)}$), but it is not the same mapping ${\displaystyle z=e^{sT}}$ which also maps the ${\displaystyle j\omega }$ axis to the unit circle. When the actual frequency of ${\displaystyle \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, ${\displaystyle \omega _{a}}$, for the continuous-time filter that this ${\displaystyle \omega _{d}}$ is mapped to.

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

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

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

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

and the inverse mapping is

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

The discrete-time filter behaves at frequency ${\displaystyle \omega _{d}}$ the same way that the continuous-time filter behaves at frequency ${\displaystyle (2/T)\tan(\omega _{d}T/2)}$. Specifically, the gain and phase shift that the discrete-time filter has at frequency ${\displaystyle \omega _{d}}$ is the same gain and phase shift that the continuous-time filter has at frequency ${\displaystyle (2/T)\tan(\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 ${\displaystyle \omega _{d}\ll 2/T}$ or ${\displaystyle \omega _{a}\ll 2/T}$), then the features are mapped to a slightly different frequency; ${\displaystyle \omega _{d}\approx \omega _{a}}$.

One can see that the entire continuous frequency range

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

is mapped onto the fundamental frequency interval

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

The continuous-time filter frequency ${\displaystyle \omega _{a}=0}$ corresponds to the discrete-time filter frequency ${\displaystyle \omega _{d}=0}$ and the continuous-time filter frequency ${\displaystyle \omega _{a}=\pm \infty }$ correspond to the discrete-time filter frequency ${\displaystyle \omega _{d}=\pm \pi /T.}$

One can also see that there is a nonlinear relationship between ${\displaystyle \omega _{a}}$ and ${\displaystyle \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 ${\displaystyle \omega _{a}={\frac {2}{T}}\tan \left(\omega _{d}{\frac {T}{2}}\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 ${\displaystyle \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 ${\displaystyle \omega _{0}}$, as well as matching at DC, if the following transform is substituted into the continuous filter transfer function.[2] This is a modified version of Tustin's transform shown above.

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

However, note that this transform becomes the original transform

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

as ${\displaystyle \omega _{0}\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

• MIT OpenCourseWare Signal Processing: Continuous to Discrete Filter Design
• Lecture Notes on Discrete Equivalents
• The Art of VA Filter Design