When we analyse signals in physics — gravitational waves, neutrino detectors, particle collisions — we almost never have an infinite, perfectly known time series. We have a finite chunk of noisy, random data. The question becomes: how do we reliably extract frequency information from it? The answer lies in spectral estimation.

This post covers the key ideas from the Laboratory of Data Analysis course at Sapienza (Prof. Paola Leaci), building up from the basic periodogram to the Welch method — the workhorse of modern PSD estimation. Each section includes an interactive plot and the equivalent Python code.

1 — Why spectral estimation, not calculation?

The power spectrum \(S(\Omega)\) is a theoretical quantity defined over an infinitely long, continuous signal. What we compute from a finite dataset is its estimate — the periodogram \(\hat{S}(\Omega)\).

Because our signals are random, not deterministic, we need a statistical framework. In the time domain we characterise random signals via the autocorrelation function; its Fourier transform gives the PSD:

Two broad families of methods exist:

• Non-parametric (Periodogram, Bartlett, Welch): minimal assumptions, used when signal structure is unknown.
• Parametric (AR, MA, ARMA): assume a generative model — more accurate when the model matches, e.g. acoustic signals.

2 — The Periodogram

Given \(N\) samples \(\{x_1,\ldots,x_N\}\), compute the Discrete Fourier Transform:

The periodogram estimate is:

3 — Spectral Resolution

The minimum resolvable frequency separation depends only on the length of data, not on zero-padding:

Zero-padding increases DFT length \(L\) and interpolates between bins — it improves appearance but not true resolution. To actually resolve two frequencies, we need more data.

import numpy as np
import matplotlib.pyplot as plt

# ── parameters ────────────────────────────────────────────────────
N  = 64       # try 16, 64, 256 ...
df = 0.05     # frequency separation between the two components

# ── signal ────────────────────────────────────────────────────────
n  = np.arange(N)
x  = np.sin(2*np.pi*0.12*n) + np.cos(2*np.pi*(0.12 + df)*n)

# ── periodogram (zero-padded for display) ─────────────────────────
L    = N * 4
xdft = np.fft.fft(x, L)
psd  = np.abs(xdft[:L//2])**2 / N
freq = np.arange(L//2) / L

delta_f = 1.0 / N   # true resolution limit

# ── plot ──────────────────────────────────────────────────────────
fig, ax = plt.subplots(figsize=(9, 4))
ax.plot(freq, psd / psd.max(), linewidth=1.5, label='Periodogram')
ax.axvline(delta_f, color='green', linestyle='--', label=f'Δf = 1/N = {delta_f:.3f}')
ax.set_xlabel('Normalised frequency (cycles/sample)')
ax.set_ylabel('Normalised PSD')
ax.set_title(f'Spectral Resolution — N={N}, df={df}')
ax.legend()
ax.grid(True, linestyle='--', linewidth=0.4, alpha=0.5)
plt.tight_layout()
plt.show()

# When df < delta_f, the two peaks are indistinguishable.
# Increase N (more data) to improve actual resolution.

4 — Windowing and Spectral Leakage

Any finite measurement multiplies the signal by a rectangular window. Multiplication in time = convolution in frequency, so the true spectrum is convolved with the window's Fourier transform — causing spectral leakage: energy bleeds into neighbouring bins.

The windowed periodogram becomes:

Three common windows:

• Rectangular: \(w_i = 1\) — narrowest main lobe (best resolution), largest side lobes (worst leakage).
• Von Hann (Hanning): \(w_i = \frac{1+\cos(2\pi(i-N/2)/N)}{2}\) — excellent side-lobe suppression, slightly wider main lobe.
• Bartlett (Triangular): linear taper — slightly narrower main lobe than Hanning but wider side lobes.

import numpy as np
import matplotlib.pyplot as plt
from scipy.fft import fft, ifft
from scipy.signal import windows

N   = 100
fs  = 8
A   = 1
fr  = 0.1
ts  = 1.0 / fs

t     = np.arange(N) * ts
x1    = A * np.sin(2 * np.pi * fr * t)

nfft  = 1024
xdft  = fft(x1, nfft)

def to_db(x):
    return 20 * np.log10(np.abs(x) + 1e-12)

def cconv_custom(a, b, n):
    return ifft(fft(a, n) * fft(b, n), n)

# ── windows ───────────────────────────────────────────────────────
wr = np.ones(N)
wh = windows.hann(N)
wt = windows.triang(N)

# ── circular convolution of signal DFT with window DFT ───────────
conv_rect   = to_db(1/nfft * cconv_custom(xdft, fft(wr, nfft), nfft))
conv_hann   = to_db(1/nfft * cconv_custom(xdft, fft(wh, nfft), nfft))
conv_triang = to_db(1/nfft * cconv_custom(xdft, fft(wt, nfft), nfft))

# ── one-sided frequency axis [0, fs/2] ───────────────────────────
dfr  = 1.0 / (nfft * ts)
frs  = np.arange(nfft) * dfr
sel  = slice(0, nfft // 2 + 1)   # even nfft
frsR = frs[sel]

# ── plot ──────────────────────────────────────────────────────────
fig, ax = plt.subplots(figsize=(10, 6))
ax.plot(frsR, conv_rect[sel],   label='Rectangular', linewidth=1.5)
ax.plot(frsR, conv_hann[sel],   label='Hanning',     linewidth=1.5)
ax.plot(frsR, conv_triang[sel], label='Triangular',  linewidth=1.5)
ax.set_xlabel('Frequency (Hz)')
ax.set_ylabel('Magnitude (dB)')
ax.set_title('Magnitude of the Power Spectrum (Circular Convolution)')
ax.legend()
ax.grid(True, linestyle='--', alpha=0.7)
plt.tight_layout()
plt.show()

5 — Spectral Estimate Uncertainty

The raw periodogram has a fundamental problem: its variance does not decrease as \(N\) increases. Adding more data improves frequency resolution, but each bin remains exponentially distributed with standard deviation equal to its mean — regardless of \(N\).

For white noise with \(S(\Omega)=1\), the periodogram values follow an exponential distribution with mean 1. Increasing \(N\) just stretches the Fourier-transformed vector — it does not give us independent repeated measurements.

6 — Bartlett and Welch: Averaging to Reduce Variance

The solution: average \(M\) independent periodograms. The standard deviation reduces by \(\sqrt{M}\) and the distribution becomes a normalised \(\chi^2\) with \(2M\) degrees of freedom.

Bartlett's method: split \(N\) samples into \(K\) non-overlapping rectangular windows of \(M=N/K\) samples each, compute and average their periodograms. Resolution worsens by \(K\), variance reduces by \(K\).

Welch's method: same as Bartlett but with 50% overlap and a Hanning window. Averages \(\approx 2K-1\) periodograms — further variance reduction. Slightly wider main lobe but much better leakage control. This is the standard in gravitational-wave data analysis.

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import periodogram, welch

# ── signal (from pwelchMethod.m) ──────────────────────────────────
Fs  = 1000
t   = np.arange(0, 1, 1/Fs)
x   = np.cos(2 * np.pi * 100 * t) + np.random.randn(len(t))
N   = len(x)

# ── periodogram (rectangular window, no averaging) ────────────────
freq, per = periodogram(x, fs=Fs, window='boxcar', nfft=N, scaling='density')
# 'boxcar' = rectangular, equivalent to rectwin(N) in MATLAB

# ── Welch's method ────────────────────────────────────────────────
nfft     = 500
win_len  = 500          # Hamming window by default in scipy (same as MATLAB pwelch)
noverlap = win_len // 2 # 50% overlap — set None for scipy default

f, pxx = welch(x, fs=Fs, window='hamming', nperseg=win_len,
               noverlap=noverlap, nfft=nfft, scaling='density')

# Number of averaged segments (with 50% overlap):
# n_segments ≈ 2 * (N / win_len) - 1 = 2*2-1 = 3

# ── comparison plot (dB) ──────────────────────────────────────────
fig, ax = plt.subplots(figsize=(9, 5))
ax.plot(freq, 10*np.log10(per), linewidth=1.5, label='Periodogram (no avg)',  color='#378ADD')
ax.plot(f,    10*np.log10(pxx), linewidth=2,   label='Welch (50% overlap)',   color='#832434')
ax.set_xlabel('Frequency (Hz)')
ax.set_ylabel('PSD (dB/Hz)')
ax.set_title('Periodogram vs Welch — variance comparison')
ax.legend()
ax.grid(True, linestyle='--', linewidth=0.4, alpha=0.5)
plt.tight_layout()
plt.show()

# ── variance comparison ───────────────────────────────────────────
print(f"var(periodogram) = {np.var(per):.6f}")
print(f"var(welch)       = {np.var(pxx):.6f}")
print(f"variance ratio   = {np.var(per)/np.var(pxx):.1f}x")
# Expected: var(welch) << var(periodogram)

7 — Whitening

In many physics applications the noise is coloured — it has a strongly frequency-dependent PSD \(S_N(f)\). Before applying a matched filter or searching for weak signals, we apply a whitening filter \(G(f) = S_N(f)^{-1/2}\):

After whitening the noise floor is flat and standard periodogram methods apply cleanly. This pipeline — estimate PSD → build whitening filter → apply to data → matched filter — underpins LIGO/Virgo gravitational-wave searches. In Python this is done with scipy.signal.lfilter or GWpy's built-in TimeSeries.whiten().

Summary

Key takeaways: more averaging always reduces variance, but shorter windows reduce frequency resolution. Hanning windows suppress leakage but widen the main lobe. Welch with Hanning and 50% overlap is the standard choice in physics data analysis.