Scientific Computing · Signal Processing

Zero-Phase Filtering by Forward-Backward Passes

Running a filter twice in opposite directions cancels its phase at the cost of squaring its magnitude, and the initial-condition choice decides how much edge transient survives.

Offline analysis of a recorded signal has a freedom that a real-time system does not: the whole record exists before processing starts, so the filter may look ahead. That freedom makes exact zero phase achievable with a causal, stable filter, by the simple expedient of applying it twice, once in each direction. The technique is standard enough to be a library function in most numerical environments, and its properties are often taken on trust. It is worth knowing precisely what the second pass does to the magnitude response, and why the initialization of the state variables dominates the quality of the result at the record boundaries.

The second pass cancels the first one’s phase

Let a stable filter have frequency response $H(\omega) = |H(\omega)| e^{j\phi(\omega)}$. Applying it to the record produces an output whose spectrum is $X(\omega)H(\omega)$. Reversing that output in time conjugates its spectrum, so the reversed intermediate has spectrum $X^(\omega)H^(\omega)$.

Applying the same filter to the reversed record multiplies by $H(\omega)$ again, giving $X^(\omega)H^(\omega)H(\omega) = X^*(\omega)|H(\omega)|^2$. Reversing the result a second time conjugates once more and yields

$$Y(\omega) = X(\omega),|H(\omega)|^2$$

The phase terms have cancelled exactly, leaving a real, nonnegative magnitude response. The output is phase-aligned with the input, and the effective impulse response is the autocorrelation of the filter’s impulse response, which is symmetric about zero and therefore has exactly zero group delay at every frequency.

That exactness is not an approximation or a small residual: it is algebraic. What the operation cannot do is separate the phase cancellation from the magnitude squaring. The effective filter is $|H|^2$, so its passband gain is squared, its stopband attenuation is doubled in decibels, and its transition band is narrowed. A Butterworth filter with a $-3$ dB point at $f_c$ becomes a filter with a $-6$ dB point at $f_c$ and a response falling at twice the original slope. A 60 dB stopband becomes 120 dB, which may be more attenuation than the coefficient word length can represent faithfully, so the nominal improvement can exceed the filter’s numerical accuracy.

The narrowing of the transition band is a genuine benefit for offline work: a fourth-order Butterworth applied twice has the selectivity of an eighth-order filter with none of the phase distortion, and without the pole-clustering problems of an actual eighth-order design. The cost is that the effective response is no longer a standard prototype, so any specification written against the original filter’s response must be reinterpreted.

import numpy as np
from scipy import signal

fs = 1000.0
b, a = signal.butter(4, 20.0, btype='low', fs=fs)
x = np.random.default_rng(0).standard_normal(4000)

y_fwd = signal.lfilter(b, a, x)
y_zp = signal.filtfilt(b, a, x)

# Zero phase: cross-correlation peaks at zero lag, not at the group delay.
lags = np.arange(-40, 41)
for y in (y_fwd, y_zp):
    c = np.correlate(y, x, mode='full')   # index len(x)-1 is lag zero
    c /= np.max(np.abs(c))
    window = c[len(x) - 1 + lags[0]: len(x) + lags[-1]]
    print(lags[np.argmax(window)])

The first pass peaks at a lag of about twenty samples for this filter, which is the group delay at low frequency; the second peaks at zero. That lag is set by the order and the cutoff rather than by the number of coefficients, since an IIR filter has no fixed length, and it is constant across frequency only for a linear-phase FIR filter. A single peak location is therefore only a summary of an IIR filter’s dispersion: the delay varies with frequency, and a broadband transient is spread rather than merely shifted.

Initial conditions decide the edge behaviour

A causal filter initialized to zero state assumes the signal was zero before the record started. That assumption is wrong for any record taken from a signal that was already running, and the mismatch produces a transient at the start whose duration is set by the filter’s slowest pole. The transient appears in the forward pass at the start and in the backward pass at the end, so a naively implemented two-pass scheme has artifacts at both edges.

The standard mitigation is to initialize the filter state to the steady state that the first samples imply, rather than to zero. For a first-order section this is a single multiplication; for a biquad it is a small linear solve. The state that produces a constant output equal to the first input, given the section’s coefficients, is found by solving the section’s state update at equilibrium. With that initialization the filter behaves as though the input had been constant for a long time before the record began, which removes the step at the boundary.

A complementary approach pads the record before filtering. Reflecting the first few samples about the boundary creates a continuation whose derivative is continuous, so the filter sees no step; replicating the boundary value creates a flat extension, which is simpler but leaves a discontinuity in the derivative. Odd reflection, which negates the continuation, preserves both value and first difference at the boundary and is the best of the simple choices. The padding length should exceed the filter’s settling time, which for an IIR filter is several times the inverse of the smallest pole radius margin; a length of three to ten times the filter order is a common working rule.

pad = 3 * max(len(a), len(b)) * 4
xp = np.concatenate([x[pad:0:-1], x, x[-2:-pad-2:-1]])  # odd reflection
yp = signal.lfilter(b, a, xp)
yp = signal.lfilter(b, a, yp[::-1])[::-1]
y = yp[pad:pad + x.size]

The reflection above is the mechanism behind the padding that library implementations perform, and doing it explicitly makes the discarded region and the transient length observable. Without the padding, or with too little of it, the transient at the record ends contaminates a region of the output as long as the filter’s settling time, which for a sharp low-order filter can be tens of samples and for a high-order one can be hundreds.

What squaring does to the frequency response

Because the effective response is a square, the parameters that specify the filter must be translated. A design targeting a $-3$ dB cutoff in the two-pass result needs the single-pass filter to have its $-1.5$ dB point at that frequency. For a Butterworth that means choosing the analog cutoff higher than the nominal value, because the single-pass response reaches $-1.5$ dB before it reaches its own $-3$ dB point; the correction factor is modest for a low order and shrinks toward one as the order rises, since a steeper response spends less of its transition between $-1.5$ and $-3$ dB.

Ripple specifications behave differently. A Chebyshev passband ripple of 0.5 dB in a single pass becomes 1.0 dB in the two-pass result, which may cross the specification limit. A stopband ripple specified as a maximum deviation becomes twice that deviation in decibels, so a nominal 60 dB stopband becomes 120 dB, but a specification written as “at least 60 dB” was already satisfied at 60 and the extra margin is wasted. The practical reading is that the two-pass technique is not a free upgrade: it doubles the decibel quantities that describe the magnitude response, so passband and stopband specifications must be halved before design.

There is also a structural consequence for the pole locations. The effective filter has each pole of the original together with its reciprocal, so the poles occur in pairs with one member of each pair inside the unit circle and the other outside; the zeros behave the same way. A pole outside the circle makes the filter non-causal: its impulse response is symmetric about zero and extends in both directions. That is exactly why it can be realized offline and why it cannot be realized in real time. A real-time system that needs zero phase must use a symmetric FIR filter instead, accepting the delay of half its length as the price of causality; the symmetry and length rules for that filter are the ones in FIR Design by Windowing and the Gibbs Ripple.

When the technique should not be used

Two-pass filtering is legitimate for offline analysis and actively wrong for streaming. A system that must produce output as samples arrive cannot wait for the reverse pass, so any claim of zero phase in such a system refers to a symmetric FIR design rather than to this method. The distinction matters when a processing pipeline is described as phase-free without saying which of the two constructions it uses.

The technique is also unsuitable when the signal has a genuine step at the boundary. If the record begins at the onset of a transient that is part of the signal rather than an artifact of recording, the backward pass will smear that onset in the reverse direction, producing a precursor that did not exist. A symmetric filter with a long impulse response spreads energy in both time directions, so a sharp event at the boundary acquires a tail before it. Reflecting the boundary does not fix this, because the reflection assumes the signal outside the record mirrors the inside, which is false for an isolated event.

Finally, the doubling of attenuation can hide a design error. A stopband specified at 40 dB that measures 80 dB after two passes satisfies the specification, but the margin came from the method rather than from the design, and the same filter used in a single pass would fail. Any verification of a filter’s performance has to be performed on the structure that is actually deployed, with the number of passes it actually receives, which is why the same filter can behave differently in analysis code and in an embedded implementation of the same nominal design. The related question of how a filter’s transition band interacts with subsequent resampling is taken up in Polyphase Decimation and the Cost of Changing Sample Rate.

References

  1. Discrete-Time Signal Processing — A. V. Oppenheim and R. W. Schafer, Pearson
  2. Biomechanics and Motor Control of Human Movement — D. A. Winter, Wiley
  3. SciPy signal module reference — SciPy project