The Gold Standard of Z-Transform | Prewarping a Digital Cutoff Before the Bilinear Transform

An analog angular frequency Ω is measured in radians per second. A digital angular frequency ω is measured in radians per sample. A digital frequency f in hertz is related to ω by ω = 2πf/fₛ, where fₛ is the sampling rate in samples per second.

At fₛ = 1000 samples per second, a digital frequency of 200 Hz corresponds to ω = 0.4π radians per sample. Its ordinary physical angular frequency is 400π radians per second. The analog prototype frequency used before a bilinear transform will generally be a different number because the mapping is nonlinear.

That last distinction is the reason for prewarping. The desired digital frequency describes where a feature should appear in the completed filter. The analog prototype frequency is an intermediate design coordinate chosen so the transformation places that feature correctly. Calling both simply “frequency” can make an incorrect substitution look reasonable.

We use capital Ω for the analog coordinate and lower-case ω for the digital coordinate throughout. The sampling interval is T = 1/fₛ. For positive-frequency examples, 0 < f < fₛ/2, or 0 < ω < π. The endpoint fₛ/2 is the Nyquist frequency. Keeping the units and the allowed interval beside the symbols prevents a missing factor of 2π from becoming an unexplained design error.

Contents · Next chapter

The standard bilinear substitution is s = 2fₛ(z−1)/(z+1), equivalently s = 2fₛ(1−z⁻¹)/(1+z⁻¹). Apply it to a rational analog transfer function Hₐ(s) to obtain a rational digital function H_d(z).

The factor 2fₛ has units of inverse seconds, matching the units of s. The ratio involving z is dimensionless. If the sample interval is used instead, the same factor is 2/T. Writing 2fₛ in one place and 2T in another would describe different substitutions and fail the unit check.

SciPy’s bilinear documentation gives this substitution and explicitly states that the function performs no prewarping. A caller who requires a particular digital cutoff must therefore supply an appropriately designed analog prototype or use a different design interface whose documented behaviour includes that step. The name of a transformation routine does not imply that it knows the intended cutoff.

The substitution changes the independent variable of the transfer function. It is not the same operation as sampling an analog impulse response or applying a zero-order hold model. Those methods answer different conversion questions. We will compare a simple pole mapping later, but first derive the frequency relation directly so that the source of the cutoff shift is visible.

Contents · Previous chapter · Next chapter

Digital sinusoidal frequency response is evaluated at z = eⁱω. Substitute that into (z−1)/(z+1). Multiplying numerator and denominator by e⁻ⁱω/² gives (eⁱω/²−e⁻ⁱω/²)/(eⁱω/²+e⁻ⁱω/²), which is i tan(ω/2).

Consequently s = i2fₛ tan(ω/2). Along the analog imaginary axis s = iΩ, so Ω = 2fₛ tan(ω/2). In hertz, replace ω by 2πf/fₛ to obtain Ω = 2fₛ tan(πf/fₛ). This is the prewarping relation when the desired digital f is known.

The inverse relation is ω = 2 arctan[Ω/(2fₛ)], or f = (fₛ/π) arctan[Ω/(2fₛ)]. It tells us where an existing analog feature lands after the standard transform. Julius O. Smith’s frequency-warping discussion and the University of Amsterdam’s signal-processing notes develop this same nonlinear mapping.

For 0 < ω < π, the tangent is positive and increasing. It maps the finite digital interval from zero to Nyquist onto the unbounded positive analog axis. Equal digital frequency increments therefore cannot correspond to equal analog increments everywhere. The geometry already tells us why one linear rescaling cannot preserve the entire frequency axis.

Contents · Previous chapter · Next chapter

Suppose the desired digital cutoff is 200 Hz at fₛ = 1000 Hz. A tempting analog choice is Ω₀ = 2π × 200 ≈ 1256.6371 radians per second. Under the standard bilinear transform, that analog frequency maps to f = (1000/π) arctan(1256.6371/2000).

The result is approximately 178.5662 Hz, not 200 Hz. The digital cutoff is about 21.4338 Hz lower than intended, a relative difference of approximately 10.72% of the target. The transform has not lost track of the frequency; it has applied the nonlinear map exactly.

The direction of the error can be predicted before calculating. For positive u, arctan(u) is less than u. Thus setting Ω₀ = 2πf꜀ gives a mapped digital frequency below f꜀, except in the limiting low-frequency approximation where the difference becomes small.

This is a useful diagnostic because a filter can be stable and have the expected qualitative low-pass shape while still missing its cutoff. A pole check alone will not detect that placement error. Keep separate questions for stability, gain and frequency location. Each has its own calculation, and a successful answer to one does not automatically answer the others.

Contents · Previous chapter · Next chapter

Use the inverse design relation instead: Ω꜀ = 2000 tan(0.2π). This gives approximately 1453.0851 radians per second. It is larger than 1256.6371 because the subsequent bilinear map compresses the analog frequency axis into the finite digital interval.

Substitute the new Ω꜀ back into the forward map. The argument of arctan is tan(0.2π), and 0.2π lies in the interval where arctan reverses tangent. The result is f = (1000/π) × 0.2π = 200 Hz exactly in the mathematical design.

The number Ω꜀/(2π), approximately 231.2657 Hz, is the analog prototype’s frequency expressed in cycles per second. It is not a revised digital target. The target remains 200 Hz; the intermediate prototype has been changed to make the transformed feature land there.

Keep both values in a design record: target digital frequency and prewarped analog angular frequency. Recording only “cutoff 231.27” would be ambiguous and could lead a later user to apply another conversion incorrectly. The whole calculation has a direction: digital requirement, analog prototype coordinate, transformed digital filter, then verification at the original digital requirement.

Contents · Previous chapter · Next chapter

Use the first-order analog low-pass prototype Hₐ(s) = Ω꜀/(s+Ω꜀). Its DC gain is one. On the imaginary axis, its magnitude is Ω꜀/√(Ω²+Ω꜀²), so at Ω = Ω꜀ the magnitude is 1/√2.

That point corresponds to a power ratio of one half and an amplitude level of about −3.0103 dB relative to the DC gain. Here “cutoff” means exactly that half-power frequency for this prototype. Other filter families can use different edge conventions, so the word cutoff should not replace the actual specification.

For example, a passband edge defined by one decibel of attenuation is not automatically the natural-frequency parameter of every second-order prototype. Prewarping accurately maps whichever analog feature has been identified, but it cannot correct a mismatch between a parameter and the response feature the designer intended.

The first-order choice keeps that issue transparent. We can derive the complete digital coefficients and verify the target magnitude without a numerical optimiser. A more elaborate design should preserve the same discipline: define the desired response feature, identify its analog counterpart, prewarp the relevant frequency and check the resulting digital response. The transformation is one step in a specification-driven calculation, not a substitute for defining the specification.

Contents · Previous chapter · Next chapter

Substitute s = 2fₛ(1−q)/(1+q), with q = z⁻¹, into Hₐ. Multiplying numerator and denominator by 1+q gives H_d = Ω꜀(1+q)/[(2fₛ+Ω꜀)+(Ω꜀−2fₛ)q]. Divide through by 2fₛ+Ω꜀ to normalise the denominator constant to one.

Write H_d = b(1+q)/(1−aq), where a = (2fₛ−Ω꜀)/(2fₛ+Ω꜀) and b = Ω꜀/(2fₛ+Ω꜀). For the prewarped 200 Hz design at 1000 Hz sampling, a ≈ 0.1583844 and b ≈ 0.4208078.

The corresponding zero-state difference equation is y[n] = ay[n−1] + bx[n] + bx[n−1]. In a denominator coefficient list, the feedback coefficient is −a, because the transfer denominator is 1−aq. In the rearranged update, the coefficient multiplying the previous output is +a. Mixing those two conventions reverses the feedback sign.

An exact identity provides a useful check: 2b = 1−a. It follows immediately from the formulas. If rounded coefficients violate that relation slightly, the DC gain may shift slightly too. Retain adequate precision for the intended calculation, then verify the implemented coefficients rather than only the ideal symbolic expressions from which they were derived.

Contents · Previous chapter · Next chapter

The digital pole is z = a ≈ 0.1583844, strictly inside the unit circle. The numerator has a zero at z = −1, the Nyquist point. At DC, q = 1, so H_d(1) = 2b/(1−a) = 1. At Nyquist, q = −1, the numerator vanishes and the response is zero.

These checks catch several common mistakes. Unequal numerator coefficients can spoil the exact Nyquist zero. A wrong feedback sign can change both DC gain and pole location. A missing normalisation can scale the entire response even if the shape looks approximately correct.

For any positive finite Ω꜀ and fₛ, |2fₛ−Ω꜀| is smaller than 2fₛ+Ω꜀. Therefore |a| < 1. This proves stability of the declared causal first-order digital filter, not merely of our numerical example. The sign of a can become negative for a sufficiently high cutoff, but its magnitude remains below one.

Stability here concerns the exact linear recurrence and its causal impulse response. Quantised state updates, coefficient errors and overflow require additional implementation analysis. The existing Z-Transform specialists on rounding and coefficient uncertainty develop those questions. An exact pole inside the unit circle is an important check, but it should not be described as a certificate for every possible stored-arithmetic implementation.

Contents · Previous chapter · Next chapter

At digital frequency f, let ω = 2πf/fₛ. The transformed response equals the analog prototype response at Ω = 2fₛ tan(ω/2). For our target f = 200 Hz, this Ω is exactly the prewarped Ω꜀, so the magnitude is 1/√2.

A second check can use the digital coefficients directly. The squared numerator magnitude is b²|1+e⁻ⁱω|² = b²(2+2cosω). The squared denominator magnitude is |1−ae⁻ⁱω|² = 1+a²−2acosω. Their ratio at ω = 0.4π should be one half.

This direct formula tests both the coefficient algebra and the frequency conversion. It does not merely repeat a software function’s reported cutoff. Using the unrounded a and b gives the exact half-power result within arithmetic precision. Using rounded coefficients gives a nearby value whose discrepancy can be measured explicitly.

Also evaluate at frequencies on either side, such as 100 Hz and 300 Hz. The response should be larger at the lower frequency and smaller at the higher one for this prototype. Those extra points check that the correctly placed feature belongs to the intended low-pass response. A single matching magnitude could otherwise occur in a wrongly interpreted or differently scaled filter.

Contents · Previous chapter · Next chapter

Apply a unit impulse from zero state. With x[−1] = 0 and y[−1] = 0, the first output is h[0] = b. At index one, h[1] = ab+b = b(1+a). At later indices, the input terms vanish and h[n] = a h[n−1].

Thus h[n] = b(1+a)aⁿ⁻¹ for n ≥ 1, together with h[0] = b. For our positive a, the tail decays without alternating signs. The absolute sum is finite because |a| < 1, consistent with the pole check.

The sum of the impulse response is b + b(1+a)/(1−a) = 2b/(1−a) = 1, matching the DC gain. This provides a time-domain verification of the endpoint calculation. For a unit step from zero state, the first output is b and later outputs approach one; the direct current-input term means the first sample is not zero.

If a program returns a delayed first response, check whether it implements a different timing convention or discretisation. The bilinear filter derived here has the explicit numerator b(1+q), which uses both current and previous input. A valid implementation must match that equation, including its initial histories. The unilateral-shift lesson explains how nonzero histories alter startup without changing the zero-state transfer function.

Contents · Previous chapter · Next chapter

For a small argument θ, tanθ is approximately θ. With θ = πf/fₛ, the prewarped Ω is therefore approximately 2πf when the target frequency is small relative to the sampling rate. That is why an unprewarped design can appear nearly correct in a low-frequency example.

The next term is tanθ = θ + θ³/3 + higher powers. Hence Ω/(2πf) is approximately 1 + θ²/3. The relative correction grows with the square of the small frequency ratio at leading order. The approximation should be judged against a stated tolerance, not by calling the ratio “small” without quantifying it.

At f/fₛ = 0.01, θ is about 0.0314 and the leading relative correction is about 0.000329, or 0.0329%. At f/fₛ = 0.2, θ is about 0.6283 and the exact correction is much larger: tanθ/θ is approximately 1.1563.

These are arithmetic comparisons, not universal recommendations for when to omit prewarping. A task with a tight frequency-placement requirement may care about a small difference. Since the exact tangent formula is simple, retaining it also avoids having to justify an approximation merely for convenience. The low-frequency expansion is most useful for understanding the error’s direction and scale.

Contents · Previous chapter · Next chapter

Keep the desired digital cutoff at 200 Hz but increase fₛ to 2000 Hz. The digital angle becomes ω = 0.2π, and the prewarped analog frequency becomes Ω꜀ = 4000tan(0.1π), approximately 1299.6788 radians per second.

The new coefficient a is approximately 0.5095254 and b is approximately 0.2452373. The pole is numerically closer to one than in the 1000 Hz design, yet each sample now represents half as much time. Comparing pole values without sample intervals would confuse discrete index behaviour with physical timescale.

The prewarped analog frequency is now closer to the ordinary 2π × 200 value because the target occupies a smaller fraction of Nyquist. If we had omitted prewarping at this higher sampling rate, the cutoff shift would be smaller, though still nonzero.

For a changed-context exercise, calculate the direct magnitude at f = 200 Hz using the new ω and coefficients. It should still equal 1/√2. Do not reuse the old 0.4π angle simply because the physical cutoff stayed at 200 Hz. The conversion from hertz to radians per sample depends on the sampling rate, and that dependency is exactly what this exercise is designed to expose.

Contents · Previous chapter · Next chapter

At fₛ = 1000 Hz, choose a new cutoff of 100 Hz. The prewarped Ω꜀ is 2000tan(0.1π), approximately 649.8394 radians per second. The resulting coefficients are again a ≈ 0.5095254 and b ≈ 0.2452373.

They match the previous chapter because both designs have the same normalised ratio f꜀/fₛ = 0.1. The digital transfer function describes behaviour per sample. The physical interpretation differs: one design runs at 2000 Hz with a 200 Hz cutoff, while the other runs at 1000 Hz with a 100 Hz cutoff.

This equivalence is useful when checking coefficient generation. For this first-order design, substitute Ω꜀ = 2fₛt with t = tan(πf꜀/fₛ). Then a = (1−t)/(1+t) and b = t/(1+t). The explicit sampling-rate factors cancel, leaving only the normalised cutoff ratio.

Do not conclude that sample rate is irrelevant. It remains essential when translating the coefficients into physical frequency and time. A coefficient file detached from its sampling rate cannot by itself specify a cutoff in hertz. Preserve the rate as part of the design metadata, and recalculate the physical meaning if the same coefficients are used in a different processing environment.

Contents · Previous chapter · Next chapter

As f꜀ approaches fₛ/2 from below, tan(πf꜀/fₛ) grows without bound. The required analog prototype frequency therefore becomes very large. A target exactly at Nyquist does not correspond to a finite Ω꜀ under this mapping.

For a high but valid target of 400 Hz at fₛ = 1000 Hz, t = tan(0.4π) ≈ 3.0776835. The coefficients become a ≈ −0.5095254 and b ≈ 0.7547627. The pole is negative but still strictly inside the circle. Its transient can alternate rather than decay with one sign.

The response remains the mapped first-order low-pass prototype, with its half-power point at 400 Hz and a zero at Nyquist. A negative feedback pole in the rearranged coefficient convention does not automatically make the filter high-pass or unstable. The full transfer function and frequency response determine those properties.

Near the endpoint, frequency mapping is steep. The derivative dΩ/df = 2πsec²(πf/fₛ) grows strongly as the target approaches Nyquist. Small changes in the requested digital edge can therefore require large changes in the analog design coordinate. This sensitivity is another reason to retain units, precision and an explicit target range rather than extrapolating a low-frequency approximation to the edge of the digital band.

Contents · Previous chapter · Next chapter

For a continuous pole p with negative real part, the bilinear map gives z = (2fₛ+p)/(2fₛ−p). The squared magnitude of the denominator minus that of the numerator is −8fₛ Re(p), which is positive. Therefore the mapped pole lies inside the unit circle.

That proves a stability relationship for the substitution. It does not say that every sample of the transformed digital response equals a sample of the original analog response. A sampled autonomous exponential eᵖᵗ advances by eᵖᵀ each sample, whereas the bilinear pole uses the rational expression above.

For the original analog pole p = −2π × 200 at fₛ = 1000 Hz, exact sampling of that exponential gives e⁻¹·²⁵⁶⁶…, approximately 0.28461. The unprewarped bilinear pole is approximately 0.22826. They differ because the conversion methods preserve different relationships.

Prewarping adds another deliberate design choice: it changes the analog prototype frequency before applying the substitution so a selected digital frequency feature lands correctly. It should not be described as restoring the complete physical sampled dynamics of the original analog system. If the task is to model a plant with a specified hold and sampling mechanism, choose and justify the appropriate discretisation rather than assuming a frequency-design transformation answers that modelling question.

Contents · Previous chapter · Next chapter

Suppose a design has two desired digital edges, f₁ = 100 Hz and f₂ = 300 Hz, at fₛ = 1000 Hz. Map each edge separately: Ω₁ = 2000tan(0.1π) ≈ 649.8394 and Ω₂ = 2000tan(0.3π) ≈ 2752.7638 radians per second.

Their ratio is approximately 4.2361, not the digital hertz ratio of three. A linear scaling of both analog frequencies by one common factor cannot generally turn an analog ratio of three into these two required coordinates. The nonlinear mapping changes relative spacing as well as absolute placement.

An analog prototype designed using both prewarped edge specifications can then have those two features map to the desired digital edges. This does not contradict the statement that one simple rescaling cannot preserve every frequency. The prototype design is being given transformed specifications before the substitution.

For a bandpass construction that uses an analog centre and bandwidth, those quantities should be derived from the prewarped edges according to that construction’s definitions. Do not insert an arithmetic digital midpoint where an analog geometric centre is required, or subtract unwarped digital edges and call the result an analog angular bandwidth. This chapter identifies the bookkeeping issue; a complete higher-order bandpass design also needs attenuation requirements, prototype choice and its own verified response.

Contents · Previous chapter · Next chapter

For the exact substitution, H_d(eⁱω) equals Hₐ(iΩ(ω)). This is an equality of complex responses at corresponding mapped frequencies. Magnitude features and phase values are carried along that correspondence. They are not generally preserved at equal numerical values of analog and digital physical frequency.

The distinction also affects group delay, which involves a derivative of phase. If φ_d(ω) = φₐ(Ω(ω)), the chain rule gives −dφ_d/dω = [−dφₐ/dΩ]dΩ/dω. Since dΩ/dω = fₛsec²(ω/2), the conversion includes a frequency-dependent factor.

The analog group delay has units of seconds, while the derivative with respect to digital radians per sample gives delay in samples. Multiplying by the mapping derivative reconciles the units. Simply copying an analog group-delay curve onto a digital frequency axis would omit that factor.

This is a useful limit on the phrase “the shape is preserved.” Under a nonlinear reparameterisation, response values correspond, but slopes and widths measured along the frequency coordinate change. State which property is being compared. The original first-order example is deliberately simple enough to verify magnitude at the target directly, while the derivative calculation shows why matching one cutoff is not a promise to preserve every timing descriptor unchanged.

Contents · Previous chapter · Next chapter

Begin with the target f꜀ and fₛ in the same hertz convention. Compute Ω꜀ using a tangent whose argument is in radians. Build the analog numerator [Ω꜀] and denominator [1, Ω꜀] for the first-order prototype, then apply the standard bilinear substitution at the stated sampling rate.

For a routine using SciPy’s documented bilinear convention, the returned digital numerator should contain two equal coefficients and the denominator should match [1, −a] after normalisation. Check the actual returned values against the closed-form expressions. Do not prewarp again if a different high-level design routine already interprets the supplied edge as a digital specification and performs the necessary conversion internally.

Evaluate the implemented frequency response at DC, the intended cutoff and Nyquist. Check the pole magnitude, then compute several impulse-response samples from the actual recurrence. These tests address different failure modes: gain normalisation, frequency units, numerator symmetry, feedback sign and timing.

Keep an exact or high-precision reference separate from a deliberately rounded implementation. If coefficients are shortened for storage, repeat the response checks using those shortened values. A plot of the ideal filter does not validate a different stored filter. For higher orders, numerical representation and conditioning deserve additional attention; the clean first-order algebra here should not be stretched into a universal implementation guarantee.

Contents · Previous chapter · Next chapter

Suppose a calculation passes the number 200 directly as the analog angular cutoff, while intending a digital cutoff of 200 Hz at 1000 Hz sampling. The units are already inconsistent: 200 radians per second is not 200 cycles per second. The mapped digital frequency becomes (1000/π)arctan(200/2000), approximately 31.73 Hz.

That error is much larger than the 178.57 Hz result caused by omitting prewarping after an otherwise correct conversion to 2π × 200. The two failures have different causes. One omitted the conversion from cycles to radians; the other used a linear-frequency analogy where the bilinear map requires a tangent.

A third mistake is to supply a tangent argument in degrees to a function expecting radians. The correct dimensionless argument for the 200 Hz design is 0.2π, about 0.6283 radians. It is also 36 degrees, but the numerical value 36 is not interchangeable with 0.6283 unless the function’s angle convention changes with it.

Diagnose these errors in order. First identify the physical units of every input. Then check the normalised frequency and angle convention. Finally check whether prewarping is required by the selected interface. Repeatedly adjusting coefficients until a graph looks plausible can conceal a unit error and make the same mistake reappear when the sample rate changes.

Contents · Previous chapter · Next chapter

The standard method in this guide keeps the substitution factor at 2fₛ and changes the analog prototype cutoff to Ω꜀ = 2fₛtan(ω꜀/2). Another convention keeps a chosen analog cutoff Ω₀ and uses a substitution s = c(z−1)/(z+1), selecting c = Ω₀/tan(ω꜀/2).

For the first-order prototype, the resulting coefficients depend on the ratio Ω₀/c. With that choice, Ω₀/c = tan(ω꜀/2), exactly the same ratio produced by our standard prewarped prototype. The digital coefficients therefore agree. This explains why two sources can present different-looking prewarping formulas without contradicting one another.

Do not combine both adjustments blindly. If the prototype has already been moved to the prewarped frequency and the substitution factor is then independently altered to compensate again, the intended correspondence may be lost. State which quantity is held fixed and which is adjusted before comparing formulas.

The equivalence here is demonstrated for matching the specified feature through the ratio entering this first-order design. It is not a licence to mix arbitrary prototype and substitution conventions in a higher-order calculation. Preserve the complete mapping, including sample-time interpretation, and verify the final digital response at the originally requested frequency. An apparently familiar formula is only useful when its symbols describe the same design choices.

Contents · Previous chapter · Next chapter

Design the same first-order half-power low-pass at f꜀ = 150 Hz with fₛ = 1200 Hz. The normalised ratio is 0.125, so πf꜀/fₛ = π/8. Use tan(π/8) = √2−1 to obtain Ω꜀ = 2400(√2−1), approximately 994.1125 radians per second.

The coefficients simplify to a = √2−1 ≈ 0.4142136 and b = 1−1/√2 ≈ 0.2928932. Verify 2b = 1−a. The difference equation is y[n] = ay[n−1] + b(x[n]+x[n−1]), with a pole inside the circle and a zero at minus one.

At the target, ω = π/4. Substituting into the direct magnitude formula gives a power ratio of one half. If the prototype had instead used the unprewarped Ω = 2π × 150, the mapped cutoff would be (1200/π) arctan(π/8), approximately 142.93 Hz. The exact value can be recomputed from that expression before deciding how many decimal places to report.

This exercise changes both the rate and target, uses an exact trigonometric value and requires independent coefficient and response checks. Explain the route from the original digital requirement to the intermediate analog coordinate and back. A learner who can reproduce that route can diagnose a new design without relying on memorised coefficients from the 200 Hz example.

Contents · Previous chapter · Continue to questions