Bilinear transform

Thebilinear transform(also known asTustin's method,afterArnold Tustin) is used indigital signal processingand discrete-timecontrol theoryto transform continuous-time system representations to discrete-time and vice versa.

The bilinear transform is a special case of aconformal mapping(namely, aMöbius transformation), often used for converting atransfer function$H_{a}(s)$of alinear,time-invariant(LTI) filter in thecontinuous-time domain (often named ananalog filter) to a transfer function$H_{d}(z)$of a linear, shift-invariant filter in thediscrete-time domain (often named adigital filteralthough there are analog filters constructed withswitched capacitorsthat are discrete-time filters). It maps positions on the$j\ Omega$axis,$\mathrm {Re} [s]=0$,in thes-planeto theunit circle,$|z|=1$,in thez-plane.Other bilinear transforms can be used for warping thefrequency responseof any discrete-time linear system (for example to approximate the non-linear frequency resolution of the human auditory system) and are implementable in the discrete domain by replacing a system's unit delays$\left(z^{-1}\right)$with first orderall-pass filters.

The transform preservesstabilityand maps every point of thefrequency responseof the continuous-time filter,$H_{a}(j\ Omega _{a})$to a corresponding point in the frequency response of the discrete-time filter,$H_{d}(e^{j\ Omega _{d}T})$although to a somewhat different frequency, as shown in theFrequency warpingsection below. This means that for every feature that one sees in the frequency response of the analog filter, there is a corresponding feature, with identical gain and phase shift, in the frequency response of the digital filter but, perhaps, at a somewhat different frequency. The change in frequency is barely noticeable at low frequencies but is quite evident at frequencies close to theNyquist frequency.

Discrete-time approximation

The bilinear transform is a first-orderPadé approximantof the natural logarithm function that is an exact mapping of thez-plane to thes-plane. When theLaplace transformis performed on a discrete-time signal (with each element of the discrete-time sequence attached to a correspondingly delayedunit impulse), the result is precisely theZ transformof the discrete-time sequence with the substitution of

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

where$T$is thenumerical integrationstep size of thetrapezoidal ruleused in the bilinear transform derivation;or, in other words, the sampling period. The above bilinear approximation can be solved for$s$or a similar approximation for$s=(1/T)\ln(z)$can be performed.

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

${\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,$H_{a}(s)$

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

That is

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

Stability and minimum-phase property preserved

A continuous-time causal filter isstableif thepolesof its transfer function fall in the left half of thecomplexs-plane.A discrete-time causal filter is stable if the poles of its transfer function fall inside theunit circlein thecomplex 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 isminimum-phaseif thezerosof 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.

Transformation of a General LTI System

A generalLTI systemhas the transfer function $$H_{a}(s)={\frac {b_{0}+b_{1}s+b_{2}s^{2}+\cdots +b_{Q}s^{Q}}{a_{0}+a_{1}s+a_{2}s^{2}+\cdots +a_{P}s^{P}}}$$ The order of the transfer functionNis the greater ofPandQ(in practice this is most likelyPas the transfer function must beproperfor the system to be stable). Applying the bilinear transform $$s=K{\frac {z-1}{z+1}}$$ whereKis defined as either2/Tor otherwise if usingfrequency 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(z+ 1)⁻¹present,(z+ 1)^-N,gives $$H_{d}(z)={\frac {b_{0}(z+1)^{N}+b_{1}K(z-1)(z+1)^{N-1}+b_{2}K^{2}(z-1)^{2}(z+1)^{N-2}+\cdots +b_{Q}K^{Q}(z-1)^{Q}(z+1)^{N-Q}}{a_{0}(z+1)^{N}+a_{1}K(z-1)(z+1)^{N-1}+a_{2}K^{2}(z-1)^{2}(z+1)^{N-2}+\cdots +a_{P}K^{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 bothN.

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,ξ_iandp_i,are thezeros and polesof the system. The bilinear transform is aone-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ξ'_iandp'_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 bothN,in other words there is now an equal number of zeros and poles. The multiplication by(z+ 1)^-Nmeans the additional zeros or poles are $${\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})}}$$

Example

As an example take a simplelow-passRC filter.This continuous-time filter has a transfer function

${\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$s$the formula above; after some reworking, we get the following filter representation:

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

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-timedigital filter.

Transformation for a general first-order continuous-time 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)={\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

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

where

$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$\ Omega _{0}$,then

$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:

$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 correspondingdifference equation.This results in

$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 theDirect 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]\.$

General second-order biquad transformation

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

$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-timedigital biquad filterwith coefficients expressed in terms of the coefficients of the original continuous time filter:

$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 correspondingdifference equation.This results in

$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 theDirect 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]\.$

Frequency warping

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

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

$H_{d}(e^{j\ Omega _{d}T})$	$=H_{a}\left({\frac {2}{T}}{\frac {e^{j\ Omega _{d}T}-1}{e^{j\ Omega _{d}T}+1}}\right)$
	$=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)$
	$=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)$
	$=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)$
	$=H_{a}\left(j{\frac {2}{T}}\cdot {\frac {\sin(\ Omega _{d}T/2)}{\cos(\ Omega _{d}T/2)}}\right)$
	$=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,$z=e^{j\ Omega _{d}T}$is mapped to a point on the$j\ Omega$axis on the continuous-time filter s-plane,$s=j\ Omega _{a}$.That is, the discrete-time to continuous-time frequency mapping of the bilinear transform is

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

and the inverse mapping is

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

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

One can see that the entire continuous frequency range

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

is mapped onto the fundamental frequency interval

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

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

One can also see that there is a nonlinear relationship between$\ Omega _{a}$and$\ Omega _{d}.$This effect of the bilinear transform is calledfrequency warping.The continuous-time filter can be designed to compensate for this frequency warping by setting$\ 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 calledpre-warpingthe filter design.

It is possible, however, to compensate for the frequency warping by pre-warping a frequency specification$\ 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$\ Omega _{0}$,as well as matching at DC, if the following transform is substituted into the continuous filter transfer function.This is a modified version of Tustin's transform shown above.

$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

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

as$\ 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 withImpulse invariance.

See also

• Impulse invariance
• Matched Z-transform method

External links

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