Bilinear transform
The bilinear transform, also known as Tustin's method, is a mathematical mapping technique in digital signal processing (DSP) and discrete-time control theory that converts continuous-time systems represented in the Laplace domain (s-plane) to discrete-time systems in the z-domain.[1][2] Introduced by Arnold Tustin in 1947 for analyzing linear systems via time series, it employs the substitution s = \frac{2}{T} \frac{1 - z^{-1}}{1 + z^{-1}}, where T is the sampling period, to derive discrete equivalents while mapping the imaginary axis in the s-plane (continuous frequencies) to the unit circle in the z-plane (discrete frequencies).[1][3] This transform is particularly valued in IIR filter design, where it enables the direct conversion of analog prototypes—such as Butterworth, Chebyshev, or elliptic filters—into stable digital filters without aliasing, as the entire frequency axis is compressed into the range from 0 to the Nyquist frequency.[2][4] Key properties include stability preservation, since the left-half s-plane maps to the interior of the unit circle in the z-plane, and a one-to-one correspondence between poles and zeros, ensuring no information loss from partial fraction expansions.[3][4] However, it induces nonlinear frequency warping, described by \Omega = \frac{2}{T} \tan\left(\frac{\omega T}{2}\right), where \Omega is the analog frequency and \omega is the digital frequency, which distorts the response unless mitigated by prewarping specific band-edge frequencies before design.[2][4] Beyond filter design, the bilinear transform finds applications in control systems for discretizing controllers and in audio signal processing for frequency-domain mappings, such as Bark-scale warping, due to its conformal nature that maps lines and circles to lines and circles while fixing DC (s=0 to z=1) and infinite frequency (to z=-1).[3][5] Its advantages over methods like impulse invariance include avoiding aliasing and simplifying implementation in tools like MATLAB, though the warping limitation makes it most suitable for bandlimited signals or lowpass/highpass filters with narrow transition bands.[2][4]Introduction
Definition and Purpose
The bilinear transform is a mathematical mapping technique employed in digital signal processing to convert continuous-time systems represented in the s-plane of the Laplace domain to equivalent discrete-time systems in the z-plane of the z-transform domain.[3] This substitution enables the discretization of linear time-invariant (LTI) systems while approximating their dynamic behavior.[2] The primary purpose of the bilinear transform is to design infinite impulse response (IIR) digital filters by starting from well-established analog prototypes, such as Butterworth or Chebyshev filters, and transforming them into digital implementations that retain core performance attributes.[6] By avoiding the need to derive digital filters from scratch, it streamlines the development process in applications requiring precise frequency-selective processing.[3] Beyond filter design, the bilinear transform is widely applied in discrete-time control systems to adapt continuous-time controllers for implementation on digital hardware, ensuring compatibility with sampled-data environments.[7] The transformation explicitly incorporates the sampling period T as a scaling parameter, which relates the continuous-time frequency response to the discrete-time domain.[2] This method is fundamentally based on the trapezoidal rule for numerical integration.[7]Historical Development
The bilinear transform, also known as Tustin's method, was first introduced by British engineer Arnold Tustin in the mid-1940s as a technique for analyzing linear systems using time series data in control applications. Tustin developed it during World War II to model the behavior of continuous-time systems through discrete approximations, particularly for servo mechanisms and feedback control in time-domain simulations. His seminal work appeared in the 1947 paper "A Method of Analysing the Behaviour of Linear Systems in Terms of Time Series," published in the Journal of the Institution of Electrical Engineers, where he derived the transformation from the trapezoidal rule of numerical integration to approximate differential equations with difference equations.[8] Parallel to Tustin's contributions, German engineers Rudolf Oldenbourg and Hans Sartorius independently advanced similar ideas using difference equations to model dynamic systems during the same wartime period. Their work focused on sampled-data control and the discretization of continuous processes for automatic regulation, laying groundwork for digital approximations in industrial control. This was detailed in their 1948 book The Dynamics of Automatic Controls, which emphasized practical applications in mechanical and electrical engineering contexts. Following the war, the method gained traction in control theory during the 1950s for numerical simulations and early digital control systems, such as those implemented in industrial processes. With the proliferation of digital computers in the 1960s, it transitioned into digital signal processing (DSP), where it became a cornerstone for designing infinite impulse response (IIR) filters by mapping analog prototypes to discrete-time equivalents. By the 1970s, its role in early numerical integration evolved into a standard tool in DSP, as evidenced by its prominent inclusion in foundational texts like Oppenheim and Schafer's Digital Signal Processing (1975), which popularized it for filter design amid the field's rapid growth.[9]Mathematical Derivation
Trapezoidal Rule Approximation
The bilinear transform arises from the discretization of continuous-time systems, particularly through the approximation of the integrator, which is fundamental to linear time-invariant (LTI) systems described by differential equations.[10] In continuous time, an integrator computes the output y(t) as the integral of the input u(t), i.e., y(t) = \int u(t) \, dt. To discretize this for digital implementation with sampling period T, numerical integration methods approximate the integral over each interval [nT, (n+1)T]. The trapezoidal rule provides a second-order accurate approximation by assuming a linear variation of the integrand between samples, using the average of the input values at the endpoints of the interval.[4] The step-by-step derivation begins with the integral over one sampling interval:y((n+1)T) = y(nT) + \int_{nT}^{(n+1)T} u(\tau) \, d\tau.
Applying the trapezoidal rule, the integral is approximated as the area of a trapezoid:
\int_{nT}^{(n+1)T} u(\tau) \, d\tau \approx \frac{T}{2} \left( u(nT) + u((n+1)T) \right).
Thus, the difference equation for the discrete integrator becomes
y(n+1) = y(n) + \frac{T}{2} \left( u(n) + u(n+1) \right),
where y(n) = y(nT) and u(n) = u(nT). This implicit relation solves for the output at the next step using both current and future input values, making it suitable for offline computation or when combined with other equations in a system.[10][11] To connect this to the s-domain, consider the continuous-time integrator in the Laplace domain: Y(s) = \frac{1}{s} U(s), where s represents differentiation. The trapezoidal approximation implies a mapping from the continuous derivative to a discrete difference. Taking the z-transform of the difference equation yields
Y(z) = \frac{T}{2} \frac{1 + z^{-1}}{1 - z^{-1}} U(z),
which rearranges to approximate the continuous relation as
s Y(s) \approx \frac{2}{T} \frac{1 - z^{-1}}{1 + z^{-1}} Y(z).
This substitution, known as the bilinear mapping, effectively replaces s with \frac{2}{T} \frac{1 - z^{-1}}{1 + z^{-1}} in the transfer function of the continuous system.[4][11] The trapezoidal rule is preferred over simpler methods like backward Euler (s \approx \frac{1 - z^{-1}}{T}) or forward Euler (s \approx \frac{1 - z^{-1}}{T z^{-1}}) because it provides a more accurate frequency response, particularly at higher frequencies. While Euler methods introduce significant distortion and potential instability for certain mappings, the trapezoidal approximation symmetrically averages the endpoints, yielding a mapping that better preserves the phase and gain characteristics up to the Nyquist frequency, albeit with nonlinear frequency warping.[10][11]