Editing Bilinear transform

{{Short description|Signal processing operation}}
{{more citations needed|date=June 2009}}
The '''bilinear transform''' (also known as '''Tustin's method''', after [[Arnold Tustin]]) is used in [[digital signal processing]] and discrete-time [[control theory]] to transform continuous-time system representations to discrete-time and vice versa.

The bilinear transform is a special case of a [[conformal map]]ping (namely, a [[Möbius transformation]]), often used for converting a [[transfer function]] <math> H_a(s) </math> of a [[linear]], [[time-invariant]] ([[LTI system theory|LTI]]) filter in the [[continuous function|continuous]]-time domain (often named an [[analog filter]]) to a transfer function <math>H_d(z)</math> of a linear, shift-invariant filter in the [[discrete signal|discrete]]-time domain (often named a [[digital filter]] although there are analog filters constructed with [[switched capacitor]]s that are discrete-time filters). It maps positions on the <math> j \omega </math> axis, <math> \mathrm{Re}[s]=0 </math>, in the [[s-plane]] to the [[unit circle]], <math> |z| = 1 </math>, in the [[complex plane|z-plane]].  Other bilinear transforms can be used for warping the [[frequency response]] of 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 <math> \left( z^{-1} \right) </math> with first order [[all-pass filter]]s.

The transform preserves [[BIBO stability|stability]] and maps every point of the [[frequency response]] of the continuous-time filter, <math> H_a(j \omega_a) </math> to a corresponding point in the frequency response of the discrete-time filter, <math> H_d(e^{j \omega_d T}) </math> although to a somewhat different frequency, as shown in the [[#Frequency warping|Frequency warping]] section 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 the [[Nyquist frequency]].

== Discrete-time approximation ==
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 [[Dirac delta function|unit impulse]]), the result is precisely the [[Z transform]] of the discrete-time sequence with the substitution of

:<math>
\begin{align}
z &= e^{sT}   \\
  &= \frac{e^{sT/2}}{e^{-sT/2}} \\
  &\approx \frac{1 + s T / 2}{1 - s T / 2}
\end{align}
</math>

where <math> T </math> is the [[numerical integration]] step size of the [[trapezoidal rule]] used in the bilinear transform derivation;<ref>{{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}}</ref> or, in other words, the sampling period. The above bilinear approximation can be solved for <math> s </math> or a similar approximation for <math> s = (1/T) \ln(z) </math> can be performed.

The inverse of this mapping (and its first-order bilinear [[Logarithm#Power series|approximation]]) is

:<math>
\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}
</math>

The bilinear transform essentially uses this first order approximation and substitutes into the continuous-time transfer function, <math> H_a(s) </math>

:<math>s \leftarrow \frac{2}{T} \frac{z - 1}{z + 1}.</math>

That is

:<math>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). \ </math>

== Stability and minimum-phase property preserved ==
A continuous-time causal filter is [[BIBO stability|stable]] if the [[Pole (complex analysis)|poles]] of its transfer function fall in the left half of the [[complex number|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 plane|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 [[Zero (complex analysis)|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.

== Transformation of a General LTI System ==
A general [[LTI system]] has the transfer function
<math display=block>
    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}
</math>
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 transfer function|proper]] for the system to be stable). Applying the bilinear transform
<math display=block>
    s = K\frac{z - 1}{z + 1}
</math>
where {{math|''K''}} is defined as either {{math|2/''T''}} or otherwise if using [[frequency warping]], gives
<math display=block>
    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}
</math>
Multiplying the numerator and denominator by the largest power of {{math|(''z'' + 1)<sup>−1</sup>}} present, {{math|(''z'' + 1)<sup>−''N''</sup>}}, gives
<math display=block>
   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}}
</math>
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
<math display=block>
    H_a(s) = \frac{(s - \xi_1)(s - \xi_2) \cdots (s - \xi_Q)}{(s - p_1)(s - p_2) \cdots (s - p_P)}
</math>
The roots of the numerator and denominator polynomials, {{math|''ξ<sub>i</sub>''}} and {{math|''p<sub>i</sub>''}}, 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
<math display=block>
    z = \frac{K + s}{K - s}
</math>
yielding some of the discretized transfer function's zeros and poles {{math|''ξ'<sub>i</sub>''}} and {{math|''p'<sub>i</sub>''}}
<math display=block>
    \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}
</math>
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)<sup>−''N''</sup>}} means the additional zeros or poles are
<ref>
{{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}}
</ref>
<math display=block>
    \begin{aligned}
        \xi'_i &= -1 \quad Q < i \leq N \\
          p'_i &= -1 \quad P < i \leq N
    \end{aligned}
</math>
Given the full set of zeros and poles, the z-domain transfer function is then
<math display=block>
    H_d(z) = \frac{(z - \xi'_1)(z - \xi'_2) \cdots (z - \xi'_N)}
                  {(z - p'_1)(z - p'_2) \cdots (z - p'_N)}
</math>

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

:<math>\begin{align}
H_a(s) &= \frac{1/sC}{R+1/sC} \\
&= \frac{1}{1 + RC s}.
\end{align}</math>

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

:{|
|-
|<math>H_d(z) \ </math>
|<math> =H_a \left( \frac{2}{T} \frac{z-1}{z+1}\right) \ </math>
|-
|
|<math>= \frac{1}{1 + RC \left( \frac{2}{T} \frac{z-1}{z+1}\right)} \ </math>
|-
|
|<math>= \frac{1 + z}{(1 - 2 RC / T) + (1 + 2RC / T) z} \ </math>
|-
|
|<math>= \frac{1 + z^{-1}}{(1 + 2RC / T) + (1 - 2RC / T) z^{-1}}. \ </math>
|}

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]].

== 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

:<math>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}}</math>

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

:<math>s \leftarrow K \frac{1 - z^{-1}}{1 + z^{-1}}</math>

where

:<math>K \triangleq \frac{2}{T} </math>.

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 <math>\omega_0</math>, then

:<math>K \triangleq \frac{\omega_0}{\tan\left(\frac{\omega_0 T}{2}\right)} </math>.

This results in a discrete-time digital filter with coefficients expressed in terms of the coefficients of the original continuous time filter:
:<math>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}}</math>

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

:<math>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}}. </math>

The difference equation (using the [[Digital filter#Direct form I|Direct form I]]) is

:<math>
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] \ .
</math>

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

:<math>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}} \ . </math>

This results in a discrete-time [[digital biquad filter]] with coefficients expressed in terms of the coefficients of the original continuous time filter:
:<math>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}}</math>

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

:<math>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}}. </math>

The difference equation (using the [[Digital filter#Direct form I|Direct form I]]) is

:<math>
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] \ .
</math>

== Frequency warping ==
To determine the frequency response of a continuous-time filter, the [[transfer function]] <math> H_a(s) </math> is evaluated at <math>s = j \omega_a </math> which is on the <math> j \omega </math> axis.  Likewise, to determine the frequency response of a discrete-time filter, the transfer function <math> H_d(z) </math> is evaluated at <math>z = e^{ j \omega_d T} </math> which is on the unit circle, <math> |z| = 1 </math>.  The bilinear transform maps the <math> j \omega </math> axis of the ''s''-plane (which is the domain of <math> H_a(s) </math>) to the unit circle of the ''z''-plane, <math> |z| = 1 </math> (which is the domain of <math> H_d(z) </math>), but it is '''not''' the same mapping <math> z = e^{sT} </math> which also maps the <math> j \omega </math> axis to the unit circle.  When the actual frequency of <math> \omega_d </math> is input to the discrete-time filter designed by use of the bilinear transform, then it is desired to know at what frequency, <math> \omega_a </math>, for the continuous-time filter that this <math> \omega_d </math> is mapped to.

:<math>H_d(z) = H_a \left( \frac{2}{T} \frac{z-1}{z+1}\right) </math>

:{|
|-
|<math>H_d(e^{ j \omega_d T}) </math>
|<math>= H_a \left( \frac{2}{T} \frac{e^{ j \omega_d T} - 1}{e^{ j \omega_d T} + 1}\right) </math>
|-
|
|<math>= 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) </math>
|-
|
|<math>= 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) </math>
|-
|
|<math>= 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) </math>
|-
|
|<math>= H_a \left(j \frac{2}{T} \cdot \frac{ \sin(\omega_d T/2) }{ \cos(\omega_d T/2) }\right) </math>
|-
|
|<math>= H_a \left(j \frac{2}{T} \cdot \tan \left( \omega_d T/2 \right) \right) </math>
|}

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

:<math> \omega_a = \frac{2}{T} \tan \left( \omega_d \frac{T}{2} \right) </math>

and the inverse mapping is

:<math> \omega_d = \frac{2}{T} \arctan \left( \omega_a \frac{T}{2} \right). </math>

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

One can see that the entire continuous frequency range

: <math> -\infty < \omega_a < +\infty </math>

is mapped onto the fundamental frequency interval

: <math> -\frac{\pi}{T} < \omega_d < +\frac{\pi}{T}. </math>

The continuous-time filter frequency <math> \omega_a = 0 </math> corresponds to the discrete-time filter frequency <math> \omega_d = 0 </math> and the continuous-time filter frequency <math> \omega_a = \pm \infty </math> correspond to the discrete-time filter frequency <math> \omega_d = \pm \pi / T. </math>

One can also see that there is a nonlinear relationship between <math> \omega_a </math> and <math> \omega_d.</math>  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 <math> \omega_a = \frac{2}{T} \tan \left( \omega_d \frac{T}{2} \right) </math> 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 <math> \omega_0 </math> (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 <math> \omega_0 </math>, as well as matching at DC, if the following transform is substituted into the continuous filter transfer function.<ref>{{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}}</ref> This is a modified version of Tustin's transform shown above.

:<math>s \leftarrow \frac{\omega_0}{\tan\left(\frac{\omega_0 T}{2}\right)} \frac{z - 1}{z + 1}.</math>

However, note that this transform becomes the original transform

:<math>s \leftarrow \frac{2}{T} \frac{z - 1}{z + 1}</math>

as <math> \omega_0 \to 0 </math>.

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

==See also==
* [[Impulse invariance]]
* [[Matched Z-transform method]]

==References==
{{reflist}}

==External links==
* [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]]