Power Law Noise Whitening — Demo

Practical walkthrough of generating, characterising, and whitening \(1/f^\alpha\) noise. See index.qmd for the theory and whitening.py for the implementation.

Show the code 
import numpy as npimport matplotlib.pyplot as pltfrom whitening import (    generate_power_law_noise_psd,    kasdin_coefficients,    whiten,    acf,    estimate_lag1_correlation,)rng = np.random.default_rng(42)plt.rcParams.update({"figure.dpi": 100, "axes.grid": True})

1. Generate coloured noiseWe use generate_power_law_noise_psd() to create white (\(\alpha=0\)), pink (\(\alpha=1\)), and brown (\(\alpha=2\)) noise by shaping the spectrum of a white noise sequence. Note: this function uses np.random internally, so we set the legacy seed for reproducibility here.

Show the code 
np.random.seed(42)N = 4096alphas = [0, 1, 2]labels = ["White (α=0)", "Pink (α=1)", "Brown (α=2)"]signals = {a: generate_power_law_noise_psd(N, alpha=a) for a in alphas}fig, axes = plt.subplots(len(alphas), 2, figsize=(10, 6), constrained_layout=True)for i, a in enumerate(alphas):    x = signals[a]    # Time series    axes[i, 0].plot(x[:500], lw=0.5)    axes[i, 0].set_ylabel(labels[i])    if i == 0:        axes[i, 0].set_title("Time series (first 500 samples)")    # PSD via Welch-like periodogram    freqs = np.fft.rfftfreq(N, d=1.0)[1:]  # skip DC    psd = np.abs(np.fft.rfft(x))[1:] ** 2 / N    axes[i, 1].loglog(freqs, psd, lw=0.5, alpha=0.7)    axes[i, 1].set_ylabel("PSD")    if i == 0:        axes[i, 1].set_title("Power spectral density")    if i == len(alphas) - 1:        axes[i, 0].set_xlabel("Sample")        axes[i, 1].set_xlabel("Frequency (normalised)")plt.show()

2. Kasdin coefficientsThe AR coefficients for \(1/f^\alpha\) noise follow a recurrence relation. For white noise (\(\alpha=0\)) there is only the trivial coefficient \([1]\). For Brownian noise (\(\alpha=2\)) the coefficients are \([1, -1]\) — a simple integrator. Pink noise (\(\alpha=1\)) has a long, slowly decaying tail.

Show the code 
fig, ax = plt.subplots(figsize=(8, 3.5))for a, label in zip([0.5, 1.0, 1.5, 2.0], ["α=0.5", "α=1.0", "α=1.5", "α=2.0"]):    c = kasdin_coefficients(a, min_magnitude=1e-4)    ax.plot(c, "o-", markersize=3, label=f"{label} ({len(c)} coeffs)")ax.set_xlabel("Coefficient index")ax.set_ylabel("Value")ax.set_title("Kasdin AR coefficients for different α")ax.legend()plt.tight_layout()plt.show()# Print first few coefficients for brown noiseprint("Brown (α=2) coefficients:", kasdin_coefficients(2.0))

3. Whitening in actionTake brown noise (\(\alpha=2\)), apply the whitening filter, and compare PSD and autocorrelation before and after. The whitened signal should look approximately white.

Show the code 
np.random.seed(42)x_brown = generate_power_law_noise_psd(N, alpha=2.0)x_whitened = whiten(x_brown, alpha=2.0)fig, axes = plt.subplots(2, 2, figsize=(10, 5), constrained_layout=True)# PSDsfreqs = np.fft.rfftfreq(N, d=1.0)[1:]for j, (sig, title) in enumerate([(x_brown, "Brown noise"), (x_whitened, "After whitening")]):    psd = np.abs(np.fft.rfft(sig))[1:] ** 2 / N    axes[0, j].loglog(freqs, psd, lw=0.5, alpha=0.7)    axes[0, j].set_title(f"PSD — {title}")    axes[0, j].set_xlabel("Frequency")    axes[0, j].set_ylabel("Power")# Autocorrelationmax_lag = 50for j, (sig, title) in enumerate([(x_brown, "Brown noise"), (x_whitened, "After whitening")]):    r = acf(sig, max_lag=max_lag)    axes[1, j].stem(range(max_lag + 1), r, linefmt="C0-", markerfmt="C0o", basefmt="k-")    axes[1, j].set_title(f"ACF — {title}")    axes[1, j].set_xlabel("Lag")    axes[1, j].set_ylabel("Autocorrelation")    axes[1, j].axhline(0, color="gray", lw=0.5)plt.show()print(f"Lag-1 autocorrelation — before: {acf(x_brown, 1)[1]:.3f}, after: {acf(x_whitened, 1)[1]:.3f}")

4. Lag-1 correlation estimation (streaming)The recursive estimator from estimate_lag1_correlation() tracks the lag-1 autocorrelation sample by sample. Here we watch it converge on pink noise (\(\alpha=1\)).

Show the code 
np.random.seed(42)x_pink = generate_power_law_noise_psd(8000, alpha=1.0)r1_stream = estimate_lag1_correlation(x_pink, eta=0.99)# Batch ground truthr1_true = acf(x_pink, max_lag=1)[1]fig, ax = plt.subplots(figsize=(8, 3))ax.plot(r1_stream, lw=0.8, label="Recursive estimate")ax.axhline(r1_true, color="red", ls="--", lw=1.2, label=f"Batch ACF(1) = {r1_true:.3f}")ax.set_xlabel("Sample index")ax.set_ylabel("Lag-1 autocorrelation")ax.set_title("Streaming lag-1 estimator on pink noise")ax.legend()plt.tight_layout()plt.show()

5. Alpha estimation from lag-1 autocorrelationThe logistic relationship \(r_1 \approx 1/(1 + e^{-4\alpha+3})\) lets us estimate \(\alpha\) from a single autocorrelation coefficient. We generate noise at several known \(\alpha\) values and compare the estimate to ground truth.

Show the code 
def estimate_alpha_from_r1(r1):    """Invert the logistic fit to get alpha from lag-1 autocorrelation."""    r1 = np.clip(r1, 1e-6, 1 - 1e-6)  # avoid log(0)    return (3 - np.log(1 / r1 - 1)) / 4true_alphas = np.arange(0.2, 2.01, 0.2)estimated_alphas = []np.random.seed(42)for a in true_alphas:    x = generate_power_law_noise_psd(8192, alpha=a)    r1 = acf(x, max_lag=1)[1]    estimated_alphas.append(estimate_alpha_from_r1(r1))fig, ax = plt.subplots(figsize=(6, 4))ax.plot([0, 2.2], [0, 2.2], "k--", lw=0.8, label="Ideal")ax.plot(true_alphas, estimated_alphas, "o-", markersize=6, label="Estimated")ax.set_xlabel("True α")ax.set_ylabel("Estimated α")ax.set_title("Alpha estimation via lag-1 autocorrelation")ax.legend()ax.set_xlim(0, 2.2)ax.set_ylim(0, 2.2)ax.set_aspect("equal")plt.tight_layout()plt.show()for a_true, a_est in zip(true_alphas, estimated_alphas):    print(f"  α = {a_true:.1f}  →  estimate = {a_est:.2f}  (error = {a_est - a_true:+.2f})")

Key takeaways- Power law noise (\(1/f^\alpha\)) is common in practice. Its correlated samples break methods that assume white noise.- Kasdin’s recurrence gives AR coefficients that model the noise; inverting that AR filter gives a simple FIR whitening filter.- Whitening flattens the PSD and drives the autocorrelation toward a delta function, restoring the white-noise assumption.- The lag-1 autocorrelation provides a cheap, single-number estimate of \(\alpha\) via a logistic fit – useful for adaptive whitening in streaming applications.