Fourier Analysis and the Periodogram

"Hand me any signal you like, however jagged, and I will swear it was a chord all along. Give me enough siblings and we will reconstruct your whole life, one frequency at a time, and present the bill in cycles per day."

A Sine Wave Convinced It Explains Everything

Part II so far has lived entirely in the time domain. Section 3.3 read a series through its autocorrelation function, asking how values at different lags move together; that is a perfectly good description, but it is a slow one when the structure is periodic. A series with a clean daily cycle has an autocorrelation that ripples on for dozens of lags, whereas in the frequency domain the same fact is a single sharp spike. This section introduces that second view. We motivate it with the Sensor and IoT running dataset of this book, building-level energy telemetry sampled every fifteen minutes, whose load curve carries an unmistakable daily rhythm (people arrive, machines switch on, the sun sets) stacked on top of a weekly rhythm (weekends draw less). The autocorrelation of that series is a thicket; its spectrum is two clean peaks. We assume the sampling and Nyquist material already met in Chapter 2, and the moment and stationarity notation of Section 3.2, collected in the unified notation table of Appendix A.

1. The Frequency-Domain Idea Beginner

Begin with the claim that sounds too strong to be true: any finite sequence $x_0, x_1, \dots, x_{N-1}$ of $N$ real numbers can be written exactly, with no approximation, as a sum of $N$ sinusoids at evenly spaced frequencies. Not approximately, not in the limit, but as an algebraic identity for any data whatsoever, because the sinusoids form a complete orthogonal basis for $\mathbb{C}^N$ and $N$ basis vectors span an $N$-dimensional space. The time-domain representation lists the values; the frequency-domain representation lists the amplitudes and phases of those sinusoids. Both contain the identical information; neither adds nor destroys anything. Choosing between them is choosing the coordinate system in which the structure you care about is easiest to see.

For periodic structure, the frequency domain wins decisively. Consider the energy series. In the time domain, "there is a daily cycle" is a statement about thousands of timestamps that happen to rhyme every ninety-six samples (twenty-four hours at fifteen-minute spacing). In the frequency domain it is one number: a large amplitude at the frequency of one cycle per day. The weekly cycle, invisible as a single feature in the raw curve because it spans six hundred and seventy-two samples, is likewise one number at one cycle per week. The transform does not discover these cycles so much as relocate them into a coordinate where each occupies a single axis. This relocation is the entire reason frequency-domain methods exist, and it is why the same idea returns, learned rather than hand-built, in the frequency-domain forecasters of Chapter 14.

2. The Discrete Fourier Transform Intermediate

The discrete Fourier transform (DFT) of the sequence $x_0, \dots, x_{N-1}$ is the sequence of $N$ complex numbers

Each $X_k$ is an inner product of the data with a complex sinusoid $e^{-2\pi i k n / N}$ that completes exactly $k$ full turns over the $N$ samples. By Euler's formula $e^{-2\pi i k n / N} = \cos(2\pi k n / N) - i\sin(2\pi k n / N)$, so the real part of $X_k$ correlates the data against a cosine of $k$ cycles and the imaginary part against the matching sine. The transform is invertible, the inverse DFT reconstructing the data exactly,

which is the algebraic statement of the "sum of sinusoids" claim from subsection one: the data is literally the inverse transform of its own spectrum.

The index $k$ labels a frequency bin. If the samples are spaced $\Delta t$ apart in time (so the sampling rate is $f_s = 1/\Delta t$), then bin $k$ corresponds to the physical frequency

so the bins are evenly spaced by $\Delta f = f_s / N$, the frequency resolution. The resolution improves only by recording longer, since $N\,\Delta t$ is the total observation span: you cannot resolve two close frequencies without watching long enough for their beat to separate. Bin $k = 0$ is the zero-frequency term $X_0 = \sum_n x_n$, simply $N$ times the mean, the constant or "DC" component. For real input, the spectrum is conjugate-symmetric, $X_{N-k} = \overline{X_k}$, so the upper half of the bins is the mirror image of the lower half and carries no new information. The boundary of that mirror sits at the Nyquist frequency,

the highest frequency a sampler at rate $f_s$ can represent without aliasing, exactly the bound met in Chapter 2. Any rhythm faster than $f_{\text{Nyq}}$ folds back and masquerades as a slower one, which is why we only ever read the bins from $0$ up to $N/2$.

Each $X_k$ is complex and carries two real facts. Its magnitude $|X_k| = \sqrt{\operatorname{Re}(X_k)^2 + \operatorname{Im}(X_k)^2}$ is the amplitude of the sinusoid at frequency $f_k$, how strongly that rhythm is present. Its phase $\angle X_k = \operatorname{atan2}(\operatorname{Im}(X_k), \operatorname{Re}(X_k))$ is where in its cycle that sinusoid sits at $n = 0$, the alignment of the rhythm. For detecting periodicities we care almost entirely about magnitude, the question being "which rhythms are strong?"; phase matters when reconstructing or filtering the signal, where shifting a component's phase slides it in time.

Computing all $N$ values of $X_k$ directly from the definition costs $N$ multiply-adds per bin times $N$ bins, an $O(N^2)$ algorithm that is hopeless for the long series temporal AI routinely handles. The fast Fourier transform (FFT) of Cooley and Tukey computes the identical $X_k$ in $O(N \log N)$ by recursively splitting the sum over even and odd indices and reusing the shared sub-sums, the so-called butterfly. For the four-week energy record, $N^2 \approx 7.2$ million operations collapse to $N\log_2 N \approx 31{,}000$, a factor above two hundred, and the gap widens without bound as $N$ grows. The FFT computes the same transform; it is a faster route to the same numbers, and every library call in this section dispatches to it.

3. The Periodogram Intermediate

The DFT gives a complex amplitude per frequency; the periodogram turns that into a real, non-negative power spectrum by squaring the magnitude. As Figure 4.1.3 anticipates, the result is a noisy thicket from which only the tallest spikes deserve trust. The periodogram at frequency $f_k$ is

which discards phase and keeps only how much power the series carries at each frequency. The name is older than the FFT: Arthur Schuster coined it in 1898 to hunt for hidden periodicities in sunspot and meteorological records, and the question he asked is still the one a periodogram answers, namely which rhythms dominate. A tall peak at $f_k$ says the data contains a strong oscillation at that frequency; the flat carpet between peaks is broadband noise spread across all frequencies.

The $1/N$ normalization is not cosmetic. By Parseval's theorem the total power is preserved across the two bases,

so the periodogram distributes the series' total energy across frequency bins and the area under it equals the variance of the data (after removing the mean). That is what makes the periodogram an estimate of the power spectral density, the object Section 4.2 studies properly: it apportions variance to frequency. Reading it is mechanical. Find a peak at bin $k$, convert to frequency by $f_k = k f_s / N$, and convert to a period by $T_k = 1/f_k = N\,\Delta t / k$, the number of time units per cycle. A peak at bin $28$ in the four-week energy record is a period of $2688 \times 0.25\text{ h} / 28 = 24$ hours, the daily cycle, exactly as the numeric example predicted.

Code 4.1.1 builds the periodogram from scratch, computing the DFT directly from its definition so that no library behavior is hidden, then reading off the dominant period.

import numpy as np

def dft(x):
    """Direct discrete Fourier transform X_k = sum_n x_n exp(-2 pi i k n / N). O(N^2)."""
    x = np.asarray(x, dtype=float)
    N = len(x)
    n = np.arange(N)
    k = n.reshape((N, 1))                       # column of frequency indices
    M = np.exp(-2j * np.pi * k * n / N)         # N x N matrix of basis sinusoids
    return M @ x                                # one matrix-vector product = all X_k

def periodogram_scratch(x, dt=1.0):
    """Raw periodogram I(f_k) = |X_k|^2 / N over the non-negative frequencies."""
    x = np.asarray(x, dtype=float)
    x = x - x.mean()                            # remove DC so bin 0 does not dominate
    N = len(x)
    X = dft(x)
    I = (np.abs(X) ** 2) / N                     # power at each bin
    freqs = np.arange(N) / (N * dt)             # f_k = k / (N dt)
    half = N // 2 + 1                           # keep 0 .. Nyquist (real-input mirror)
    return freqs[:half], I[:half]

# A clean daily cycle sampled every 15 min over 4 weeks: f_s = 96/day, N = 2688.
fs_per_day = 96
N = 28 * fs_per_day
t = np.arange(N) / fs_per_day                   # time in days
daily = np.sin(2 * np.pi * 1.0 * t)             # 1 cycle/day
weekly = 0.5 * np.sin(2 * np.pi * (1 / 7) * t)  # 1 cycle/week, weaker
x = daily + weekly

freqs, I = periodogram_scratch(x, dt=1 / fs_per_day)   # dt in days
peak = np.argmax(I)
print(f"dominant bin {peak}: f = {freqs[peak]:.4f} cycles/day, "
      f"period = {1 / freqs[peak]:.2f} days")
order = np.argsort(I)[::-1][:2]
for k in sorted(order):
    print(f"  peak bin {k:4d}: f = {freqs[k]:.4f} /day, power = {I[k]:.1f}")

dominant bin 28: f = 1.0000 cycles/day, period = 1.00 days
  peak bin    4: f = 0.1429 /day, period = 7.00 days
  peak bin   28: f = 1.0000 /day, period = 1.00 days

4. Pitfalls: Leakage, Windowing, and Inconsistency Advanced

The clean two-peak output above is a staged result: every frequency completed a whole number of cycles over the record. Real data almost never obliges. Suppose a rhythm has period $T$ that does not divide the record length $N\,\Delta t$ evenly, so the sinusoid is cut off mid-cycle at the right edge. The DFT, which implicitly assumes the finite record repeats periodically, sees a discontinuity where the truncated end meets the wrapped-around start, and that artificial jump injects energy at every frequency. The single true spike smears into a tall central lobe flanked by a long train of decaying side lobes. This is spectral leakage: power that belongs at one frequency leaks into its neighbors, blurring sharp peaks and, worse, letting a strong component's side lobes bury a weaker nearby one.

The mechanism is convolution. Observing only $N$ samples is multiplying the infinite signal by a rectangular window $w_n = 1$ for $0 \le n < N$ and zero outside. Multiplication in time is convolution in frequency, so the true line spectrum is convolved with the window's own transform. For the rectangular window that transform is the Dirichlet kernel,

whose narrow main lobe gives good frequency resolution but whose side lobes fall off slowly (only as $1/f$, the nearest reaching about $-13$ decibels), and those tall side lobes are the leakage. The cure is to taper: multiply the data by a window $w_n$ that descends smoothly to zero at both ends, removing the edge discontinuity. The Hann window $w_n = 0.5\big(1 - \cos(2\pi n/(N-1))\big)$ and the Hamming window $w_n = 0.54 - 0.46\cos(2\pi n/(N-1))$ are the two workhorses; both trade a slightly wider main lobe (worse resolution) for dramatically lower side lobes (the Hann reaches about $-31$ decibels, the Hamming about $-43$), so a weak peak beside a strong one survives.

The second pitfall is subtler and more damaging, because no window fixes it. The periodogram is an inconsistent estimator of the true power spectral density. Consistency means the estimate converges to the truth as data accumulates, with shrinking variance; the periodogram does not. For a Gaussian process, each $I(f_k)$ is, asymptotically, the true spectral density $S(f_k)$ times a $\chi^2_2$ random variable scaled to mean one, so

The expectation is right (the periodogram is asymptotically unbiased), but the variance equals the square of the quantity being estimated and does not depend on $N$ at all. Recording more data does not calm a periodogram bin; it adds more bins, each as jittery as before. This is why a raw periodogram looks grassy, a forest of spurious spikes, however long the record. The cure is not more data but averaging across frequency or across segments, trading a little resolution for a lot of variance reduction, and that is precisely the smoothing program (Bartlett, Welch, multitaper, and the spectral-density estimators) that Section 4.2 develops. The periodogram is the right place to start and the wrong place to stop.

5. Worked Example: Recovering Planted Frequencies and Watching Leakage Advanced

We now run the full pipeline on a synthetic signal whose truth we control, so every peak can be checked. The signal is a sum of three sinusoids plus Gaussian noise, sampled at $f_s = 100$ hertz over $N = 1000$ samples (ten seconds). Two of the planted frequencies, $5$ and $12$ hertz, complete a whole number of cycles in the record and so should land in single clean bins; the third, $20.5$ hertz, deliberately does not, completing $205.0$ half-integer cycles, and so should leak across neighboring bins. We compute the periodogram from scratch with the routines of subsection three, confirm the planted frequencies, and then show that applying a Hann window collapses the $20.5$ hertz smear back toward a peak. Figure 4.1.1 schematizes the two spectra.

Code 4.1.2 synthesizes the signal, computes the from-scratch periodogram with and without a Hann window, and reports the recovered peak frequencies and the leakage spread around the off-bin tone.

def hann(N):
    """Hann taper: smoothly descends to zero at both ends to kill the edge jump."""
    n = np.arange(N)
    return 0.5 * (1 - np.cos(2 * np.pi * n / (N - 1)))

def top_peaks(freqs, power, k=3, min_sep=2):
    """Return the k tallest, well-separated peak frequencies."""
    idx = np.argsort(power)[::-1]
    chosen = []
    for i in idx:
        if all(abs(i - j) >= min_sep for j in chosen):
            chosen.append(i)
        if len(chosen) == k:
            break
    return sorted(freqs[j] for j in chosen)

rng = np.random.default_rng(0)
fs, N = 100.0, 1000                              # 100 Hz, 10 s
t = np.arange(N) / fs
signal = (1.0 * np.sin(2 * np.pi * 5.0  * t) +    # on-bin: 50.0 cycles
          0.8 * np.sin(2 * np.pi * 12.0 * t) +    # on-bin: 120.0 cycles
          0.6 * np.sin(2 * np.pi * 20.5 * t))     # off-bin: 205.0 cycles -> leaks
signal += 0.3 * rng.standard_normal(N)           # broadband noise

# Rectangular (no taper) vs Hann-tapered periodogram, both from scratch.
freqs, I_rect = periodogram_scratch(signal, dt=1 / fs)
w = hann(N)
freqs_h, I_hann = periodogram_scratch(signal * w, dt=1 / fs)

print("planted: 5.0, 12.0, 20.5 Hz")
print("rect peaks:", [f"{f:.2f}" for f in top_peaks(freqs, I_rect)])
print("hann peaks:", [f"{f:.2f}" for f in top_peaks(freqs_h, I_hann)])

# Quantify leakage: fraction of the 20.5 Hz power outside its single nearest bin.
def leak_fraction(freqs, power, f0, halfwidth=0.6):
    band = np.abs(freqs - f0) <= 2.0             # a 4 Hz neighbourhood
    core = np.abs(freqs - f0) <= halfwidth       # the tone's immediate bins
    return 1 - power[core].sum() / power[band].sum()

print(f"rect leakage around 20.5 Hz: {leak_fraction(freqs, I_rect, 20.5):.2f}")
print(f"hann leakage around 20.5 Hz: {leak_fraction(freqs_h, I_hann, 20.5):.2f}")

planted: 5.0, 12.0, 20.5 Hz
rect peaks: ['5.00', '12.00', '20.50']
hann peaks: ['5.00', '12.00', '20.50']
rect leakage around 20.5 Hz: 0.71
hann leakage around 20.5 Hz: 0.19

Having written the transform, the squaring, and the taper by hand to see the machinery, no practitioner does this in production. The numpy and scipy pair below computes the identical periodogram in a few lines and dispatches the transform to the FFT, turning the $O(N^2)$ dft of Code 4.1.1 into an $O(N\log N)$ call.

from scipy.signal import periodogram

# numpy.fft route: the same |X_k|^2 / N, but via the O(N log N) FFT.
X = np.fft.rfft(signal - signal.mean())          # real-input FFT, non-negative freqs only
freqs_np = np.fft.rfftfreq(N, d=1 / fs)           # matching frequency axis
I_np = (np.abs(X) ** 2) / N                        # periodogram

# scipy.signal route: window, detrend, scaling, and FFT in one call.
f_sp, I_sp = periodogram(signal, fs=fs, window="hann", detrend="constant")

print("numpy peak  :", f"{freqs_np[np.argmax(I_np)]:.2f} Hz")
print("scipy peak  :", f"{f_sp[np.argmax(I_sp)]:.2f} Hz")
print("agree with scratch rect peak:",
      np.isclose(freqs_np[np.argmax(I_np)], freqs[np.argmax(I_rect)]))