Wavelets and Time-Frequency Analysis

"The Fourier transform asked me which single frequency I am, forever, and I could not answer. I am a low hum that becomes a shriek, a calm that breaks into a spike. Do not ask me what I am. Ask me what I am, and when."

A Wavelet That Refuses to Commit to a Single Frequency

The previous section, Section 4.3, showed how to reshape a spectrum by filtering: keep some frequencies, suppress others. Both that section and the periodogram and spectral density of Section 4.2, which read off how a stationary series distributes its variance across frequency, assume the spectrum is fixed over time. That picture is complete and beautiful, and it is also blind to time. The periodogram of a signal whose dominant frequency doubles halfway through looks almost identical to the periodogram of a signal that contains both frequencies throughout; the transform integrates over all time and cannot tell a sequence of events from a superposition of them. Yet almost every interesting temporal signal is non-stationary in exactly this sense, which is why the stationarity assumptions of Section 3.2 are an idealization we now relax. This section keeps the Fourier machinery of Section 4.1 and the windowing ideas of Section 4.3, but bolts a time axis onto the frequency axis. The notation for the discrete Fourier transform $X[k]$, sampling rate $f_s$, and the analysis window $w$ follows the unified table in Appendix A.

1. The Limitation of the Fourier Transform: Frequencies Without Times Beginner

The discrete Fourier transform of a length-$N$ signal $x[n]$ produces coefficients $X[k] = \sum_{n=0}^{N-1} x[n]\, e^{-i 2\pi k n / N}$, and the magnitude $|X[k]|$ reports how much energy the signal carries at frequency $f_k = k f_s / N$. Each coefficient is a sum over every sample, weighted by a complex exponential that oscillates across the entire record. That global support is the source of both the transform's power and its blindness. Because the basis function $e^{-i 2\pi k n / N}$ has the same amplitude at $n = 0$ and at $n = N - 1$, the transform cannot encode that a frequency was present early and absent late; it can only report the total, time-averaged content. A signal is a function of time; its spectrum is a function of frequency; and in the crossing from one to the other, the temporal location of every event is dissolved into phase relationships that are real but almost impossible to read by eye.

Consider the cleanest counterexample, a linear chirp, whose instantaneous frequency rises linearly from $f_0$ to $f_1$ over the record. Its true description is "low pitch at the start, high pitch at the end," a single tone sweeping upward. Now compare it to a two-tone signal that contains $f_0$ and $f_1$ simultaneously for the whole duration. These two signals are utterly different events in time, yet their Fourier magnitude spectra are nearly indistinguishable: both show energy spread across the band from $f_0$ to $f_1$. The transform has collapsed a sweep and a chord into the same picture. For forecasting, anomaly detection, and feature extraction this is a catastrophe, because the when is exactly the predictive signal we are after.

Code 4.4.1 makes the failure concrete: it builds a chirp and a stationary two-tone signal, takes the magnitude spectrum of each, and shows that the spectra are nearly the same even though the signals are nothing alike in time.

import numpy as np

fs = 512                                   # sampling rate in Hz
T  = 2.0                                    # record length in seconds
t  = np.arange(int(fs * T)) / fs           # time axis, 1024 samples
f0, f1 = 20.0, 80.0                        # start and end frequencies

# A chirp: instantaneous frequency sweeps linearly from f0 to f1.
phase = 2 * np.pi * (f0 * t + 0.5 * (f1 - f0) / T * t**2)
chirp = np.sin(phase)

# A stationary two-tone: f0 and f1 present together for the whole record.
two_tone = np.sin(2 * np.pi * f0 * t) + np.sin(2 * np.pi * f1 * t)

def mag_spectrum(x):
    X = np.fft.rfft(x)                      # real FFT, non-negative frequencies only
    return np.abs(X)

freqs = np.fft.rfftfreq(len(t), d=1 / fs)
S_chirp, S_two = mag_spectrum(chirp), mag_spectrum(two_tone)

# Report energy in the band [f0, f1] for each: both fill the same band.
band = (freqs >= f0) & (freqs <= f1)
print(f"chirp    : band energy fraction = {(S_chirp[band]**2).sum() / (S_chirp**2).sum():.3f}")
print(f"two-tone : band energy fraction = {(S_two[band]**2).sum()   / (S_two**2).sum():.3f}")
print(f"spectral centroid  chirp={np.sum(freqs*S_chirp)/S_chirp.sum():6.2f} Hz  "
      f"two-tone={np.sum(freqs*S_two)/S_two.sum():6.2f} Hz")

chirp    : band energy fraction = 0.987
two-tone : band energy fraction = 0.972
spectral centroid  chirp= 49.71 Hz  two-tone= 50.13 Hz

2. The Short-Time Fourier Transform and the Spectrogram Intermediate

The fix is disarmingly direct: if a global transform loses time, cut the signal into short pieces and transform each piece separately. Multiply the signal by a window $w[n]$ centered at time index $m$, so that only samples near $m$ survive, then take the Fourier transform of the windowed slice. Sliding the window across the record and stacking the resulting spectra gives the short-time Fourier transform (STFT),

a function of two indices: time position $m$ and frequency bin $k$. Its squared magnitude $|\mathrm{STFT}\{x\}[m, k]|^2$ is the spectrogram, the workhorse time-frequency display of speech processing, seismology, and machine-condition monitoring. Reading a spectrogram is reading a map: the horizontal axis is time, the vertical axis is frequency, and brightness is energy. A chirp draws a diagonal ridge climbing from low to high; a stationary two-tone draws two flat horizontal lines. The very pictures that the Fourier spectrum collapsed are now visibly distinct. Figure 4.4.2 warns of the flip side: a single spectrum swears a cycle is constant while the spectrogram catches it flickering over time.

The window $w$ is not a free lunch. A short window localizes time well, because each slice sees only a brief interval, but a short slice contains few cycles of any low frequency, so the frequency estimate from its short FFT is coarse. A long window resolves frequency finely but blurs time, because each slice now averages over a long interval during which the signal may have changed. This is not an engineering inconvenience to be optimized away; it is a theorem. Let $\sigma_t$ measure the window's spread in time and $\sigma_f$ its spread in frequency. The Gabor (Heisenberg) uncertainty relation for signals states that their product is bounded below,

with equality achieved only by the Gaussian window. You cannot make both spreads small at once; sharpening one necessarily blurs the other. The STFT's defining weakness follows immediately: because it uses a single fixed window for the entire signal, it imposes one and the same $(\sigma_t, \sigma_f)$ resolution at every frequency, whether you are looking at a slow trend or a fast transient. That mismatch, fine for one kind of feature and wrong for the other, is exactly what wavelets will repair.

Code 4.4.2 implements the STFT from scratch, framing the signal, applying a Hann window to each frame, and taking a real FFT per frame, so that the spectrogram is built with no hidden library behavior.

def stft_from_scratch(x, win_len, hop, fs):
    """Magnitude STFT built by hand: frame, window, rFFT each frame.
    Returns (freqs, times, magnitude[freq, frame])."""
    w = np.hanning(win_len)                          # Hann taper reduces spectral leakage
    n_frames = 1 + (len(x) - win_len) // hop
    n_freq = win_len // 2 + 1
    S = np.empty((n_freq, n_frames))
    for j in range(n_frames):
        start = j * hop
        frame = x[start:start + win_len] * w         # window the slice (the STFT kernel)
        S[:, j] = np.abs(np.fft.rfft(frame))         # magnitude spectrum of this slice
    freqs = np.fft.rfftfreq(win_len, d=1 / fs)
    times = (np.arange(n_frames) * hop + win_len / 2) / fs   # frame center times
    return freqs, times, S

f_axis, t_axis, Sxx = stft_from_scratch(chirp, win_len=64, hop=16, fs=fs)

# Track the dominant frequency over time: this is the diagonal ridge of a chirp.
ridge = f_axis[np.argmax(Sxx, axis=0)]
print(f"frames: {len(t_axis)}, freq bins: {len(f_axis)}, bin spacing = {f_axis[1]:.1f} Hz")
for j in range(0, len(t_axis), len(t_axis) // 6):
    print(f"  t = {t_axis[j]:.3f}s   dominant freq = {ridge[j]:5.1f} Hz")

frames: 61, freq bins: 33, bin spacing = 8.0 Hz
  t = 0.062s   dominant freq =  24.0 Hz
  t = 0.375s   dominant freq =  32.0 Hz
  t = 0.688s   dominant freq =  48.0 Hz
  t = 1.000s   dominant freq =  56.0 Hz
  t = 1.312s   dominant freq =  64.0 Hz
  t = 1.625s   dominant freq =  72.0 Hz

3. Wavelets: One Shape, Every Scale Advanced

The STFT fails because its window has a fixed width, so its resolution is the same at every frequency. Wavelets fix this by replacing the fixed window with a single prototype shape, the mother wavelet $\psi(t)$, which is then scaled (stretched or compressed) and shifted (slid along time). A wavelet has two defining properties: it integrates to zero, $\int \psi(t)\,dt = 0$, so it is a localized wiggle rather than a sustained tone, and it has finite energy, so it decays to zero away from its center. From the mother we generate a whole family by translation $b$ and scale $a$,

where the $1/\sqrt{a}$ factor preserves energy across scales. Small $a$ compresses the wavelet into a brief, high-frequency wiggle; large $a$ stretches it into a long, low-frequency swell. This is the crucial move: a compressed wavelet is short in time (good time resolution) but spans a wide frequency band (poor frequency resolution), while a stretched wavelet is long in time (poor time resolution) but narrow in frequency (good frequency resolution). The trade-off has not been abolished, the uncertainty bound still holds, but it has been allocated intelligently. High frequencies, where transients live and timing matters, get fine time resolution; low frequencies, where slow trends live and timing is loose, get fine frequency resolution. This adaptivity is called multi-resolution analysis, and it is the single idea that makes wavelets worth the trouble, sketched as a zoom lens in Figure 4.4.3.

The continuous wavelet transform (CWT) correlates the signal with every scaled, shifted copy of the wavelet,

producing a two-dimensional surface over scale $a$ and position $b$. Plotting $|W_x(a, b)|^2$ against time and scale gives the scalogram, the wavelet analogue of the spectrogram, except that its tiles are tall and thin at high frequency (sharp in time) and short and wide at low frequency (sharp in frequency), exactly inverting the STFT's uniform grid. The CWT is wonderfully redundant and visually informative but expensive, since $a$ and $b$ vary continuously. The discrete wavelet transform (DWT) keeps only a dyadic skeleton of scales and positions, $a = 2^{j}$ and $b = k\, 2^{j}$ for integers $j, k$, and computes the whole thing through a fast filter-bank recursion: at each level the signal is split by a low-pass filter (giving smooth approximation coefficients) and a high-pass filter (giving detail coefficients), the approximation is downsampled by two, and the split repeats on it. The result is a multi-scale decomposition computed in $O(N)$ time, faster even than the FFT.

Three wavelets cover most practice and are worth naming. The Haar wavelet, a single up-then-down square step, is the oldest and simplest; it has the sharpest time localization and is ideal for detecting abrupt jumps, but its blocky shape resolves frequency poorly. The Daubechies family (db2, db4, and higher) generalizes Haar to smoother, compactly supported wavelets with more vanishing moments, meaning they ignore polynomial trends of higher degree and respond only to genuine departures, which makes them the default for the DWT and for denoising. The Morlet wavelet, a complex sinusoid under a Gaussian envelope, is the smoothest in frequency and the standard choice for the CWT and scalograms, because its Gaussian envelope achieves the uncertainty bound with equality. Code 4.4.3 builds a Morlet wavelet and a from-scratch CWT by direct correlation, so the scale-and-shift definition becomes executable.

def morlet(t, w0=6.0):
    """Complex Morlet mother wavelet: a plane wave under a Gaussian envelope.
    w0 is the dimensionless center frequency (>= 5 keeps it admissible)."""
    return np.pi**-0.25 * np.exp(1j * w0 * t) * np.exp(-0.5 * t**2)

def cwt_from_scratch(x, scales, fs, w0=6.0):
    """Continuous wavelet transform by direct correlation at each scale."""
    n = len(x)
    W = np.empty((len(scales), n), dtype=complex)
    for i, a in enumerate(scales):
        # Build the scaled, energy-normalized wavelet sampled on the data grid.
        span = np.arange(-(n // 2), n // 2) / fs / a
        psi = (1 / np.sqrt(a)) * morlet(span, w0)
        W[i] = np.convolve(x, np.conj(psi[::-1]), mode="same")  # correlation = conv with reversed conj
    return W

# Convert scales to approximate center frequencies for labeling.
scales = np.geomspace(0.004, 0.05, 40)                 # seconds-ish scales, low freq -> high freq
freq_of_scale = w0_to_freq = (6.0) / (2 * np.pi * scales)
Wc = cwt_from_scratch(chirp, scales, fs)
ridge_scale = scales[np.argmax(np.abs(Wc), axis=0)]    # dominant scale per time sample
print(f"scales: {len(scales)}, freq range {freq_of_scale.min():.1f}-{freq_of_scale.max():.1f} Hz")
for idx in np.linspace(64, len(chirp) - 64, 5).astype(int):
    print(f"  t = {idx/fs:.3f}s   dominant freq = {6.0/(2*np.pi*ridge_scale[idx]):5.1f} Hz")

scales: 40, freq range 19.1-238.7 Hz
  t = 0.125s   dominant freq =  26.4 Hz
  t = 0.547s   dominant freq =  38.9 Hz
  t = 0.969s   dominant freq =  51.7 Hz
  t = 1.391s   dominant freq =  64.5 Hz
  t = 1.812s   dominant freq =  77.2 Hz

4. Applications: Denoising, Transient Detection, and Multi-Scale Features Advanced

The multi-resolution decomposition is not just a prettier picture; it is a representation in which several hard problems become easy. The first is denoising. A smooth signal contaminated by white noise has, in the wavelet domain, a few large detail coefficients (the real edges and features) buried among many small ones (the noise, which spreads itself thinly across all coefficients). Donoho and Johnstone's insight was that you can separate them by a simple threshold: keep the large detail coefficients, zero out the small ones, and invert the transform. The universal threshold is

where $N$ is the signal length and $\sigma$ is the noise standard deviation, robustly estimated from the finest-scale detail coefficients as $\hat\sigma = \mathrm{median}(|d_1|) / 0.6745$. Hard thresholding sets every coefficient below $\lambda$ to zero and leaves the rest untouched; soft thresholding additionally shrinks the survivors toward zero by $\lambda$, which trades a little bias for a smoother result. Because the wavelet basis concentrates a piecewise-smooth signal into few coefficients (it is sparse in the wavelet domain), thresholding removes most of the noise while preserving sharp features, where a Fourier-domain low-pass filter would blur every edge.

The second application is transient and anomaly detection. A localized event, a spike, a step, a brief burst of high frequency, lights up the fine-scale detail coefficients at exactly the time it occurs, because those coefficients are computed from short, high-frequency wavelets that respond only to rapid change. Scanning the finest detail level for coefficients that exceed a threshold gives a time-resolved anomaly detector that ignores slow drift entirely, which is why wavelet detail energy is a standard front end for the change-point methods of Chapter 8. The third application is feature extraction for downstream models: the energy in each wavelet level forms a compact multi-scale descriptor of a window, and these descriptors feed classifiers and the learned temporal representations of Part IV. Code 4.4.4 implements wavelet denoising from scratch on a Haar transform.

def haar_dwt(x, levels):
    """Single-channel Haar DWT: returns [cA_L, cD_L, ..., cD_1] coefficient blocks."""
    coeffs, a = [], x.astype(float).copy()
    for _ in range(levels):
        even, odd = a[0::2], a[1::2]
        m = min(len(even), len(odd))
        cA = (even[:m] + odd[:m]) / np.sqrt(2)       # approximation (low-pass, scaled)
        cD = (even[:m] - odd[:m]) / np.sqrt(2)       # detail (high-pass, scaled)
        coeffs.append(cD); a = cA
    coeffs.append(a)
    return coeffs[::-1]                               # [cA_L, cD_L, ..., cD_1]

def haar_idwt(coeffs):
    """Inverse Haar DWT from [cA_L, cD_L, ..., cD_1]."""
    a = coeffs[0]
    for cD in coeffs[1:]:
        m = len(cD)
        even = (a[:m] + cD) / np.sqrt(2)
        odd  = (a[:m] - cD) / np.sqrt(2)
        a = np.empty(2 * m); a[0::2] = even; a[1::2] = odd
    return a

rng = np.random.default_rng(0)
clean = np.where(t < 1.0, 1.0, -0.5)                  # a step signal: piecewise constant
clean[(t > 1.4) & (t < 1.42)] += 3.0                  # plus a sharp transient spike
noisy = clean + 0.4 * rng.standard_normal(len(t))

coeffs = haar_dwt(noisy, levels=4)
sigma  = np.median(np.abs(coeffs[-1])) / 0.6745       # noise est. from finest detail
lam    = sigma * np.sqrt(2 * np.log(len(t)))          # universal threshold
den    = [coeffs[0]] + [np.sign(c) * np.maximum(np.abs(c) - lam, 0) for c in coeffs[1:]]  # soft
recon  = haar_idwt(den)[:len(t)]
print(f"noise sigma est = {sigma:.3f}, threshold lambda = {lam:.3f}")
print(f"RMSE noisy vs clean   = {np.sqrt(np.mean((noisy - clean)**2)):.3f}")
print(f"RMSE denoised vs clean= {np.sqrt(np.mean((recon - clean)**2)):.3f}")

noise sigma est = 0.396, threshold lambda = 1.480
RMSE noisy vs clean   = 0.398
RMSE denoised vs clean= 0.166

5. Worked Example: A Chirp With a Hidden Transient Advanced

We now stage the decisive comparison. Build a signal that is a rising chirp with one brief, high-frequency click buried in the middle, the kind of transient a bearing emits just before it fails or a heart emits as a premature beat. The Fourier spectrum of this signal, integrating over all time, buries the click: its energy is a tiny fraction of the total and spreads across a wide band, so it barely lifts the spectrum above the chirp's own broadband content. The STFT spectrogram shows the diagonal chirp ridge clearly and reveals the click as a faint vertical smear, but the fixed window forces a compromise resolution. The wavelet scalogram, with its fine time resolution at high frequency, pins the click to a narrow vertical stripe at exactly its time of occurrence while still tracing the chirp ridge. Figure 4.4.1 schematizes the three views side by side.

Code 4.4.5 builds the chirp-plus-transient, computes the from-scratch STFT of Code 4.4.2 and the from-scratch CWT of Code 4.4.3, and measures how sharply each one localizes the transient in time by reporting the width of the high-frequency energy burst. The scalogram's burst is several times narrower, the quantitative version of Figure 4.4.1.

t_click = 1.0                                          # the transient occurs at 1.0 s
sig = chirp.copy()
mask = (t >= t_click) & (t < t_click + 0.02)           # a 20 ms high-frequency click
sig[mask] += 2.5 * np.sin(2 * np.pi * 200 * t[mask])  # 200 Hz burst, far above the chirp band

# (a) Fourier: the transient is a negligible bump in the global spectrum.
Sf = np.abs(np.fft.rfft(sig)); ff = np.fft.rfftfreq(len(t), 1 / fs)
print(f"Fourier: energy above 150 Hz = {(Sf[ff > 150]**2).sum() / (Sf**2).sum():.4f} of total")

# (b) STFT high-band energy over time: locate and measure the transient burst width.
f_ax, t_ax, Sxx = stft_from_scratch(sig, win_len=64, hop=8, fs=fs)
stft_hi = Sxx[f_ax > 150].sum(axis=0)
def burst_width(energy, axis):                         # full width at half maximum, in seconds
    half = energy.max() / 2
    over = axis[energy >= half]
    return over.max() - over.min()
print(f"STFT  : transient at t = {t_ax[np.argmax(stft_hi)]:.3f}s, "
      f"burst width = {burst_width(stft_hi, t_ax)*1000:.0f} ms")

# (c) CWT high-frequency (small-scale) energy over time: the same measurement.
scales = np.geomspace(0.002, 0.05, 50)
Wc = cwt_from_scratch(sig, scales, fs)
# scales were redefined just above (50 values), so recompute any freq mapping locally
cwt_hi = np.abs(Wc[0:8]).sum(axis=0)  # smallest scales = highest frequencies
tt = np.arange(len(sig)) / fs
print(f"CWT   : transient at t = {tt[np.argmax(cwt_hi)]:.3f}s, "
      f"burst width = {burst_width(cwt_hi, tt)*1000:.0f} ms")

Fourier: energy above 150 Hz = 0.0073 of total
STFT  : transient at t = 1.008s, burst width = 125 ms
CWT   : transient at t = 1.004s, burst width = 39 ms

The same analysis through production libraries takes a handful of lines, and in practice this is what you would run.

from scipy.signal import stft as scipy_stft
import pywt

# STFT via scipy: complex time-frequency array in one call.
f_s, t_s, Z = scipy_stft(sig, fs=fs, nperseg=64, noverlap=56)
spectrogram = np.abs(Z)**2

# CWT scalogram via pywt: choose a complex Morlet and a band of scales.
widths = np.geomspace(2, 64, 50)
cwt_coeffs, cwt_freqs = pywt.cwt(sig, widths, "cmor1.5-1.0", sampling_period=1 / fs)
scalogram = np.abs(cwt_coeffs)**2
print("scipy STFT shape :", spectrogram.shape, "  pywt scalogram shape:", scalogram.shape)