Power Spectral Density

Power Spectral Density (cross-cutting/numerics/power-spectral-density)

Use when the task is estimating the one-sided power spectral density (units^2/Hz) of a stochastic sampled signal by Welch's averaged modified periodogram: Hann-windowed, overlapping segments, density-scaled, with the equivalent noise bandwidth and the integrated total power. This is the measurement-to-spectrum step for random-vibration and noise surveys, complementing the single-record spectrum leaf that transforms one record without windowing, segment averaging, or Hz scaling. Implemented in pure stdlib Python, offline and deterministic.

Domain quick reference

• Welch method: split the record into segments of M samples, multiply each by the Hann window w[n] = 0.5 - 0.5 cos(2 pi n / M), transform each windowed segment with the leaf's iterative radix-2 Cooley-Tukey transform, scale each to a one-sided density periodogram, and average over segments. Segments start every hop = M * (1 - overlap) samples; the default overlap is 50 percent, giving 63 segments for an 8192-sample record at M = 256: (8192 - 256) / 128 + 1.
• One-sided density scaling for a real segment, X the transform of the windowed segment and w the window: P[k] = 2 |X[k]|^2 / (fs * sum(w^2)) on interior bins k = 1..M/2 - 1. The DC bin (k = 0) and the Nyquist bin (k = M/2) are NOT doubled: P[0] and P[M/2] keep the single-sided numerator |X|^2 / (fs * sum(w^2)).
• Frequency axis: freqs[k] = k * fs / M for k = 0..M/2, so the bin width df = fs / M. A 60 Hz tone at fs = 1024 Hz with M = 256 lands on bin 15 exactly (60 = 15 * 4).
• Equivalent noise bandwidth: ENBW = fs * sum(w^2) / (sum(w))^2, the Hz width of an ideal rectangular filter passing the same white-noise power as the window. Hann at M = 256: 1024 * 96 / 128^2 = 6.0 Hz.
• Total power = sum_k P[k] * df. With this scaling the integral equals the variance for a zero-mean signal and the mean square otherwise, the energy-conservation check that catches scaling bugs.
• Peak density of a pure sine of amplitude A on an interior bin: P_peak = A^2 / (2 * ENBW); the integrated one-sided power is A^2 / 2 regardless of window, because sum(w^2) cancels in the density scale.
• Units: fs in Hz, amplitudes in the signal unit, PSD in (signal unit)^2 / Hz, e.g. g^2/Hz for an acceleration survey.
• All functions are deterministic and stdlib-only (math); no network, no third-party numerical libraries. NACA Report 824 anchors the pack's public-domain reference set; Welch's method is generic signal-processing methodology, summary-only.

Workflow

1. Collect the sampled record x (real values), the sample rate fs in Hz, and choose the segment length M, a power of two.
2. Form the window with hann_window(M); check the noise-width cost of the window with equivalent_noise_bw(window, fs).
3. Estimate the density: welch_psd(x, fs, M, overlap=0.5) returns the frequency axis and the averaged PSD. Reduce overlap to lower the segment count, or toward 1.0 for more averaging at the cost of correlated segments (the hop never drops below one sample).
4. Identify the tone or response peak: psd_summary returns the peak density and peak frequency together with enbw_hz, df_hz, total power, and the full freqs/psd arrays.
5. Integrate the level: psd_total_power(psd, df) gives the band or total power; compare with the variance of a zero-mean record as an energy check.
6. For a random-vibration survey, report the peak g^2/Hz level, the frequency of the peak, and the total rms from sqrt(total power), the quantities a downstream response analysis consumes.
7. Confirm the deterministic checks with the contract test scripts/test_power_spectral_density.py.

Worked example

fs = 1024 Hz, M = 256, Hann window, 50% overlap, 8192-sample record (63 segments), 60 Hz sine at amplitude A = 1.0 (bin 15 exactly). Real outputs of this module:

• hann_window(256) sums to 128.000000 (M/2) with sum of squares 96.000000 (3M/8).
• equivalent_noise_bw = 6.000000 Hz = 1024 * 96 / 128^2.
• welch_psd returns 129 bins, df = 4.000000 Hz; the peak sits at bin 15, 60.000000 Hz, with density 0.08333333333 = A^2 / (2 * ENBW) = 1/12 (ratio to theory 1.0 to 1e-15).
• Integrated total power = 0.50000000000 = A^2 / 2 (energy check).
• Amplitude A = 0.5: peak density 0.02083333333, exactly one quarter of the A = 1 peak (6.02 dB down), integrated power 0.125000000.
• DC bin not doubled: a constant record of 1.0 integrates to 1.000000000 (= c^2, its mean square) with P[0] = 0.166666667 = (sum w)^2 / (fs * sum w^2), not the doubled 0.333333333.
• Nyquist bin not doubled: a (-1)^n record integrates to 1.000000000 with P[128] = 0.166666667, again undoubled.
• Tone A = 1 plus seeded Gaussian noise of variance 0.2 (sigma 0.447), fixed seed 21: peak density 0.084140 at 60.0 Hz against the clean tone 0.083333, off-peak noise floor 5.66e-4 (max bin beyond three bins of the peak), signal-to-noise ratio about 149, and the run is bit-identical on repetition.

Pitfalls

• Doubling the DC or Nyquist bins: P[0] and P[M/2] are NOT doubled, so a constant record of 1.0 integrates to 1.0 with P[0] = 0.166666667, not the doubled 0.333333333.
• Using a non-power-of-two segment length or other non-physical inputs: seg_len not a power of two, fs <= 0, overlap outside [0, 1), x shorter than seg_len, empty or zero-sum windows, and df <= 0 all raise ValueError.
• Comparing peak densities across different windows: the peak of a pure sine sits at A^2 / (2 * ENBW), so the density depends on the window's noise bandwidth, while the integrated one-sided power A^2 / 2 does not.
• Reading the peak as the total power: the energy-conservation check integrates sum_k P[k] * df to the variance (or mean square); the peak height alone undercounts broadband content.
• Confusing this Welch density with the single-record spectrum: the Hann windowing, 50%-overlap segment averaging, and Hz density scaling are what define this estimate; the FFT sibling transforms one record without them.
• Expecting the noisy PSD to repeat without its seed: the tone-plus- noise worked case is deterministic only with the fixed seed, and the noiseless anchors are the repeatable checks.

Verification

• Confirm the window identities: sum(w) = 128.0 and sum(w^2) = 96.0 for M = 256, and ENBW = fs * sum(w^2) / (sum(w))^2 = 6.0 Hz.
• Confirm the pure-sine anchors to 1e-6 relative: peak density equals A^2 / (2 * ENBW) and the integrated power equals A^2 / 2; halving the amplitude quarters the peak density.
• Confirm the DC and Nyquist bins are not doubled: a constant record and a (-1)^n record each integrate to their mean square with P[0] and P[M/2] at the undoubled value.
• Confirm determinism: two calls on the same record return identical PSDs, and a fixed seed regenerates an identical noisy PSD.
• Confirm psd_summary returns exactly the keys freqs, psd, enbw_hz, df_hz, total_power, peak_density, peak_freq_hz.
• Confirm every non-physical input raises ValueError: seg_len not a power of two, fs <= 0, overlap outside [0, 1), x shorter than seg_len, empty or zero-sum windows, and df <= 0.
• Run the contract test offline: python3 scripts/test_power_spectral_density.py (35 tests, deterministic).

Related leaves

• cross-cutting/numerics/fast-fourier-transform: the single-record spectrum sibling; this leaf adds windowing, segment averaging and density scaling on top of the radix-2 transform.
• cross-cutting/numerics/fir-filter-design: Hamming-window truncation for FIR coefficient design, a different use of windows.
• structures/loads/random-vibration-analysis: a PSD consumer that turns an input density spectrum into a response; this leaf estimates the input PSD from a measured time history.

Behavior contract (gate 3)

Run the deterministic contract test (stdlib unittest, offline):

python3 scripts/test_power_spectral_density.py

The test covers the Hann window identities (sum M/2, sum of squares 3M/8, symmetry), the equivalent noise bandwidth (6.0 Hz anchor and the rectangular-window limit fs/M), the internal radix-2 transform against the transform definition, the one-sided periodogram (axis, peak density A^2/(2 ENBW), integrated power A^2/2, the quarter-amplitude law, undoubled DC and Nyquist bins), the Welch average (63 segments at default overlap, mean-of-periodograms identity, frequency axis, full 8192-sample anchors, zero-overlap path, determinism, seeded-noise reproducibility and the peak/floor/SNR magnitude bounds), the total power integral, the summary dict keys and values, and ValueError rejection of every non-physical input.

Compliance

• Standards referenced, not reproduced: NACA TR-824 anchors the pack's public-domain reference set; Welch's averaged periodogram method is standard signal-processing methodology, summary-only per standards-map.yaml.
• compliance: STANDARDS-REF, gated: false.