Filtering in the Frequency Domain

"Slow trends may pass; I wave them through without a glance. Fast noise knocks all day and I never answer. The high frequencies call me cold. I prefer the word principled."

A Low-Pass Filter With Strong Boundaries

The previous section, Section 4.2, gave us the power spectral density and consistent estimators of it (building on the periodogram of Section 4.1): the tools to see a signal's frequency content. Seeing is half the job. Once the spectrum reveals that the trend lives below some frequency, the seasonal cycle sits at a known peak, and the measurement noise sprawls across the high end, the natural next move is to act: keep what you want, remove what you do not. That act is filtering, and its theoretical engine is the Fourier transform of Section 4.1. We assume the discrete-time spectral notation of this chapter, with $X(f)$ the Fourier transform of a series $\{x_t\}$ and the frequency $f$ measured in cycles per sample (so the Nyquist limit sits at $f = 1/2$); the full symbol table lives in the unified notation of Appendix A. Filtering will return, in learned form, when frequency-domain forecasters appear in Chapter 14; the present section is where the classical idea is forged.

1. Filtering as Multiplication in Frequency Beginner

A linear time-invariant filter acts on an input series $\{x_t\}$ by convolution with a fixed kernel $\{h_k\}$, producing the output

The kernel $\{h_k\}$ is the filter's impulse response: it is literally what comes out when you feed in a single unit spike. Convolution in the time domain is a sliding, weighted sum, and it is cumbersome both to reason about and, for long kernels, to compute. The convolution theorem dissolves the difficulty. Writing $X(f)$, $H(f)$, and $Y(f)$ for the Fourier transforms of the input, the kernel, and the output, the theorem states

Convolution becomes pointwise multiplication. The function $H(f)$, the Fourier transform of the kernel, is the frequency response of the filter, and it is the single object that fully characterizes what the filter does. At each frequency $f$ the response is a complex number, $H(f) = |H(f)|\, e^{i\angle H(f)}$, whose magnitude $|H(f)|$ scales the amplitude of that frequency component and whose phase $\angle H(f)$ shifts its timing. A filter is nothing more than a recipe for the gain curve $|H(f)|$: where the gain is near one the frequency passes, where it is near zero the frequency is blocked, and the band between is the transition.

The cleanest illustration is the simple moving average, the most-used filter in all of forecasting. An $L$-term average has the kernel $h_k = 1/L$ for $k = 0, \dots, L-1$ and zero elsewhere. Summing the geometric series of its Fourier transform gives a closed form, the Dirichlet kernel,

whose magnitude is $|H(f)| = \big|\sin(\pi f L) / (L \sin(\pi f))\big|$. Inspect this gain: at $f = 0$ it equals one (the average passes the constant, the trend), and as $f$ grows it falls toward zero, with the first null exactly at $f = 1/L$. A moving average is therefore a low-pass filter, and now we know precisely how low: lengthening the window pushes the first null down to lower frequency, smoothing more aggressively. This is why a 12-month moving average annihilates a 12-sample seasonal cycle (its period sits right on a null) while a 3-month average barely touches it. The numeric example pins the numbers down.

Code 4.3.1 makes the gain curve concrete by computing it two ways that must agree: the closed-form Dirichlet expression and the discrete Fourier transform of the actual kernel. Their agreement is the convolution theorem made executable.

import numpy as np

def moving_average_response(L, freqs):
    """Magnitude gain |H(f)| of an L-term moving average, two ways."""
    freqs = np.asarray(freqs, dtype=float)
    # 1. Closed-form Dirichlet kernel magnitude (guard the f=0 limit -> 1).
    num = np.sin(np.pi * freqs * L)
    den = L * np.sin(np.pi * freqs)
    closed = np.where(np.isclose(freqs, 0.0), 1.0, np.abs(num / den))
    # 2. Brute force: DFT of the actual kernel evaluated at the same freqs.
    h = np.full(L, 1.0 / L)
    k = np.arange(L)
    dft = np.array([np.abs(np.sum(h * np.exp(-1j * 2 * np.pi * f * k)))
                    for f in freqs])
    return closed, dft

freqs = np.array([0.0, 0.01, 1 / 12, 0.4])
closed, dft = moving_average_response(12, freqs)
for f, c, d in zip(freqs, closed, dft):
    print(f"f = {f:6.4f}:  closed |H| = {c:6.4f}   DFT |H| = {d:6.4f}")

f = 0.0000:  closed |H| = 1.0000   DFT |H| = 1.0000
f = 0.0100:  closed |H| = 0.9770   DFT |H| = 0.9770
f = 0.0833:  closed |H| = 0.0000   DFT |H| = 0.0000
f = 0.4000:  closed |H| = 0.0544   DFT |H| = 0.0544

2. Filter Families and the Cost of Sharp Cutoffs Intermediate

Filters are named for the band they let pass. A low-pass filter keeps frequencies below a cutoff $f_c$ and blocks those above; it extracts trends and smooths noise. A high-pass filter does the reverse, keeping fast variation and removing the slow drift, which is how you detrend a series in one stroke ($\text{high-pass} = \text{identity} - \text{low-pass}$). A band-pass filter keeps a contiguous middle band $[f_1, f_2]$ and rejects everything outside it, isolating a business cycle or a known oscillation. A band-stop (or notch) filter does the opposite, carving out one narrow band while passing the rest, which is exactly how engineers remove the 50 or 60 hertz power-line hum from a sensor trace. The four families are duals built from the same parts; Figure 4.3.1 sketches their ideal gain curves.

The ideal gain curves above are rectangles, perfectly flat in the passband and exactly zero in the stopband with a vertical wall between. They are also impossible. The inverse Fourier transform of a rectangular gain is the sinc function, $h_k \propto \sin(2\pi f_c k)/(\pi k)$, which extends infinitely in both directions and is non-causal (it needs future samples). Any realizable filter must truncate or approximate that infinite kernel, and truncation has a price. Cutting a sinc to finite length is multiplying it by a rectangular window in time, which convolves the ideal brick-wall gain with the window's spectrum, smearing the sharp wall into a sloped transition band and wrapping ripples around it. Those ripples are the Gibbs phenomenon: overshoot near the cutoff that does not shrink as you lengthen the kernel, it only narrows. The deep trade-off of all filter design lives here.

A second axis classifies filters by their structure. A finite impulse response (FIR) filter has a kernel of finite length, like the moving average: its output is a weighted sum of a bounded window of inputs, it is always stable, and it can be made exactly linear-phase (a constant delay at every frequency, so no waveform distortion). An infinite impulse response (IIR) filter feeds its own past outputs back in, $y_t = \sum_k b_k x_{t-k} - \sum_j a_j y_{t-j}$, so a finite kernel of coefficients produces an impulse response that rings on forever. IIR filters reach a given sharpness with far fewer coefficients than FIR (the Butterworth and Chebyshev designs are IIR), but they can go unstable and they impose a nonlinear phase that distorts waveform shape unless corrected. The corrective trick, zero-phase filtering, is the subject of subsection three. The trade-offs are collected below.

The comparison table is the lookup card for choosing a filter type, and it explains why the worked example of subsection five builds a windowed-sinc FIR by hand (for its transparent linear phase) and reaches for a Butterworth IIR as the library shortcut (for its cheap steep roll-off). Both will be made zero-phase before we trust their output.

3. Designing, Applying, and Not Distorting Intermediate

Two design recipes cover most time-series needs. The windowed-sinc FIR builds a low-pass filter by taking the ideal sinc kernel and multiplying it by a smooth window that tapers its ends to zero, trading a slightly wider transition for vastly less Gibbs ringing. For a cutoff $f_c$ (in cycles per sample) and a kernel of odd length $N = 2M + 1$ centered at zero, the windowed kernel is

where $\operatorname{sinc}(x) = \sin(\pi x)/(\pi x)$ and $w_k$ is a window such as the Hamming taper $w_k = 0.54 + 0.46 \cos(\pi k / M)$. The Butterworth IIR takes a different route: rather than approximate a brick wall, it asks for the gain that is maximally flat in the passband, giving the magnitude

where $n$ is the filter order. At $f = f_c$ the gain is exactly $1/\sqrt{2}$ (the half-power point) regardless of order, and raising $n$ steepens the roll-off beyond the cutoff while keeping the passband flat. This compact form is why Butterworth filters dominate practical denoising.

Applying a filter raises a hazard that ruins more analyses than any other: phase distortion. Any causal filter delays the signal, and an IIR filter delays different frequencies by different amounts, so a single forward pass shifts peaks, troughs, and crossings in time. For a time series this is catastrophic: a filtered turning point that has moved three samples to the right is a forecast that is wrong by three samples. The cure is zero-phase filtering, implemented as filtfilt: filter the series forward, then filter the result backward. The backward pass undoes the forward pass's phase exactly (the two delays cancel), leaving a net phase of zero at every frequency while squaring the magnitude response, $|H_{\text{eff}}(f)| = |H(f)|^2$. The output is perfectly time-aligned with the input. The price is that the filter is non-causal: it uses future samples, so it is a tool for offline analysis, not for real-time forecasting where the future is unavailable.

The third hazard is edge effects. A convolution near the start or end of a finite series reaches past the data's edge, and whatever you assume lives there (zeros, the first value repeated, a mirror reflection) contaminates the first and last $M$ outputs. On a long series these few bad samples are a rounding error; on a short series they can be most of the output. The standard mitigations are reflecting the signal at the boundary (the default padding filtfilt uses) or simply discarding the unreliable edge region. For the short economic and clinical series that dominate this part of the book, edge effects are not a footnote; they are the difference between a usable trend estimate and a fabricated one at exactly the most interesting point, the recent end.

4. Time-Series Applications: Trend, Cycle, and Clean Sensors Advanced

Filtering is the workhorse behind decomposition. In macroeconomics the canonical task is to split a series like real GDP into a slow trend and a business cycle around it, and the two most-cited tools are filters in disguise. The Hodrick-Prescott (HP) filter finds the trend $\{\tau_t\}$ that minimizes a penalized least-squares objective balancing fit against smoothness,

where the smoothing parameter $\lambda$ controls how rigid the trend is: $\lambda = 0$ returns the data untouched, and $\lambda \to \infty$ forces a straight line. The conventional values are $\lambda = 1600$ for quarterly data and $\lambda = 129{,}600$ for monthly. Although it is written as an optimization, the HP filter is exactly a linear low-pass filter; its cyclical component $x_t - \tau_t$ is the matching high-pass. The Baxter-King (BK) filter is more explicit: it is a band-pass FIR designed to isolate fluctuations with periods between 6 and 32 quarters (the standard definition of a business cycle), built as the difference of two ideal low-pass filters and then windowed to finite length, which forces it to discard $K$ samples at each end. Both let an economist answer "is the economy above or below trend right now," and both inherit every caveat of subsection three, especially the unreliable endpoint.

• HP filter: a penalized-smoothness low-pass, where $\lambda$ sets how rigid the trend is; equivalent to a fixed linear filter despite its optimization dress.
• BK filter: a band-pass FIR isolating the 6 to 32 quarter business-cycle band, built from two low-pass filters and windowed, which discards $K$ samples at each end.

Seasonal filtering is the same idea aimed at a known period. A series with a 12-month cycle can have that cycle removed by a band-stop notch centered at $f = 1/12$ and its harmonics, or extracted by the matching band-pass; the official statistical agencies' X-13ARIMA-SEATS procedure is, at its core, an elaborate cascade of such moving-average filters. Sensor denoising, the running theme of Part II's later chapters, is the most direct application: an accelerometer or temperature probe produces a slow physical signal buried in high-frequency electronic noise, and a low-pass filter recovers the signal, provided the cutoff sits above the physics and below the noise. When the noise overlaps the signal band, filtering alone fails and you escalate to the model-based anomaly detection of Chapter 8; filtering is the cheap first line of defense that handles the easy ninety percent.

5. Worked Example: A Low-Pass Filter From Scratch and From a Library Advanced

We now build a low-pass filter end to end. The target is a noisy signal: a slow component at $f = 0.02$ cycles per sample (the trend we want to keep) plus a fast contaminant at $f = 0.30$ (the noise we want to remove), with a cutoff placed between them at $f_c = 0.08$. We first design a windowed-sinc FIR kernel by hand, apply it by FFT-based convolution (using the convolution theorem of subsection one so the cost is one multiply in frequency rather than a full sliding sum), and read off its frequency response. Code 4.3.2 does the design and the spectral application.

import numpy as np

def windowed_sinc_lowpass(fc, num_taps):
    """Design a low-pass FIR kernel: ideal sinc times a Hamming window.
    fc is the cutoff in cycles/sample; num_taps must be odd for symmetry."""
    if num_taps % 2 == 0:
        num_taps += 1
    M = num_taps // 2
    k = np.arange(-M, M + 1)
    ideal = 2 * fc * np.sinc(2 * fc * k)          # np.sinc is the normalized sinc
    window = 0.54 + 0.46 * np.cos(np.pi * k / M)  # Hamming taper kills Gibbs ringing
    h = ideal * window
    return h / h.sum()                            # normalize to unit DC gain

def fft_convolve(x, h):
    """Apply filter h to signal x via the convolution theorem (multiply in freq)."""
    n = len(x) + len(h) - 1
    nfft = 1 << (n - 1).bit_length()              # next power of two for a fast FFT
    Y = np.fft.rfft(x, nfft) * np.fft.rfft(h, nfft)   # Y(f) = X(f) H(f)
    y = np.fft.irfft(Y, nfft)[:len(x)]            # back to time, trim to input length
    return y

# Build the test signal: slow trend (keep) + fast noise (remove).
rng = np.random.default_rng(0)
t = np.arange(400)
trend = np.sin(2 * np.pi * 0.02 * t)             # f = 0.02, the wanted component
noise = 0.6 * np.sin(2 * np.pi * 0.30 * t)       # f = 0.30, the contaminant
signal = trend + noise + 0.1 * rng.standard_normal(400)

h = windowed_sinc_lowpass(fc=0.08, num_taps=61)
clean = fft_convolve(signal, h)

# Frequency response: gain at the two component frequencies.
H = np.abs(np.fft.rfft(h, 4096))
fr = np.fft.rfftfreq(4096)
for f in (0.02, 0.08, 0.30):
    g = np.interp(f, fr, H)
    print(f"|H({f:.2f})| = {g:.4f}")
print(f"kernel length = {len(h)} taps")

|H(0.02)| = 0.9988
|H(0.08)| = 0.5126
|H(0.30)| = 0.0007

The 61-tap kernel above is a faithful low-pass, but it is a forward convolution and therefore carries a delay of $M = 30$ samples; its output peaks lag the input. To demonstrate phase distortion and its cure, Code 4.3.3 applies a Butterworth IIR filter two ways: a single forward pass (lfilter, which delays) and a zero-phase forward-backward pass (filtfilt, which does not). This is also the library shortcut: the entire hand-built design-and-convolve pipeline of Code 4.3.2 collapses to two scipy.signal calls.

from scipy.signal import butter, lfilter, filtfilt

# Design a 4th-order Butterworth low-pass. Wn is the cutoff as a fraction of
# Nyquist (=0.5 cycles/sample), so 0.08 cycles/sample -> Wn = 0.08 / 0.5 = 0.16.
b, a = butter(N=4, Wn=0.16, btype="low")

one_pass  = lfilter(b, a, signal)    # causal: correct for real-time, but delayed
zero_phase = filtfilt(b, a, signal)  # forward+backward: zero delay, offline only

# Measure the timing error by locating the first trend peak in each output.
true_peak = np.argmax(trend[:80])
lag_one  = np.argmax(one_pass[:80])  - true_peak
lag_zero = np.argmax(zero_phase[:80]) - true_peak
print(f"true trend peak at t = {true_peak}")
print(f"one-pass  lfilter  peak shifted by {lag_one:+d} samples (phase distortion)")
print(f"zero-phase filtfilt peak shifted by {lag_zero:+d} samples (aligned)")

true trend peak at t = 12
one-pass  lfilter  peak shifted by +6 samples (phase distortion)
zero-phase filtfilt peak shifted by +0 samples (aligned)

Figure 4.3.2 sketches what the three time-domain outputs look like side by side, making the phase shift visible: the noisy input, the delayed one-pass result, and the aligned zero-phase result tracking the true trend.

The line-count contrast is stark. The from-scratch path of Code 4.3.2 spends roughly 22 lines designing the windowed-sinc kernel and implementing FFT convolution before it filters a single sample. The library path of Code 4.3.3 reaches the same low-pass result in two lines, butter to design and filtfilt to apply, an order-of-magnitude reduction, and it throws in the zero-phase correction, the stability checks, and the boundary padding for free.