Bilinear transform

Template:Short description Template:More citations needed 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 mapping (namely, a Möbius transformation), often used for converting a transfer function <math> H_a(s) </math> of a linear, time-invariant (LTI) filter in the 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-time domain (often named a digital filter although there are analog filters constructed with switched capacitors 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 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 filters.

The transform preserves 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 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 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>Template:Cite book</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 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 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.

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 Template:Math is the greater of Template:Math and Template:Math (in practice this is most likely Template:Math as the transfer function must be proper for the system to be stable). Applying the bilinear transform <math display=block>

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

</math> where Template:Math is defined as either Template:Math 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 Template:Math present, Template:Math, 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 Template:Math.

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, Template:Math and Template:Math, 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 Template:Math and Template:Math <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 Template:Math, in other words there is now an equal number of zeros and poles. The multiplication by Template:Math means the additional zeros or poles are <ref> {{#invoke:citation/CS1|citation |CitationClass=web }} </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 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 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>Template:Cite book</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

External links

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