Calibrated spectral analysis

A spectrum without its uncertainty is half a measurement. This page covers the Welch spectral estimators of phonometry.metrology that report, next to the spectrum itself, the statistical quality of the estimate following Bendat & Piersol, Random Data: Analysis and Measurement Procedures (4th ed., 2010): the power spectral density and cross-spectral density with the effective number of averages, the normalized random error and chi-square confidence intervals; the coherent output spectrum that splits a measured output into the part linearly explained by the input and the noise remainder, with the spectral signal-to-noise ratio; a fractional-octave smoother with a constant-power kernel; and colored-noise generators with an exact power-law slope for exercising all of the above. A Thomson multitaper estimator (Percival & Walden, 1993) completes the family for records too short to segment. Every error formula is a closed form from the sources, verified by seeded Monte Carlo in the test suite.

1. Power spectral density with its statistical error

power_spectral_density estimates the one-sided autospectral density Gxx(f) by Welch’s method: the record is split into tapered (Hann by default), 50 %-overlapped segments whose periodograms are averaged. No detrending is applied, so absolute calibration is preserved: a signal in pascals yields Pa²/Hz. Two scalings are available: 'density' (units²/Hz, integrates to the signal power) and 'spectrum' (units², reads the power of discrete tones directly).

Averaging nd independent segments gives the estimate 2·nd chi-square degrees of freedom (Eq. 8.162), from which everything else follows:

$${\epsilon }_r[{\hat{G}}_{xx}]=\frac{1}{\sqrt{n_d}},\frac{n{\hat{G}}_{xx}}{{\chi }_{n;\alpha /2}^2}\le G_{xx}\le \frac{n{\hat{G}}_{xx}}{{\chi }_{n;1-\alpha /2}^2},n=2n_d.$$

With overlapped, tapered segments the averages are correlated, so the result reports both the raw segment count (n_segments) and the effective number of independent averages (n_averages), computed with the window-correlation formula of Welch (1967) that Bendat & Piersol reference in Section 11.5.2.2; for a Hann taper at 50 % overlap it is roughly 0.95 of the raw count. The random error and the confidence interval use the effective value. At DC, and at Nyquist for an even segment length, the one-sided spectrum has a single real Fourier component, so those bins carry half the degrees of freedom and a correspondingly wider interval.

from phonometry import power_spectral_density

res = power_spectral_density(signal, fs)          # Hann, 50 % overlap, 95 % CI
print(res.n_averages, res.random_error)           # nd and 1/sqrt(nd)
print(res.ci_lower[10], res.psd[10], res.ci_upper[10])
res.plot()                                        # PSD in dB with the CI band

The resolution bias is the other half of the error budget: a finite analysis bandwidth Be (reported as resolution_bandwidth, the effective noise bandwidth of the taper) smooths sharp spectral features, always in the direction of reduced dynamic range (Eq. 8.139). For a resonance peak of half-power bandwidth Br, the first-order normalized bias is the closed form of Eq. 8.141, exposed as resolution_bias_error:

$${\epsilon }_b[{\hat{G}}_{xx}(f_r)]\approx -\frac{1}{3}{\left(\frac{B_e}{B_r}\right)}^2.$$

from phonometry import resolution_bias_error

eps_b = resolution_bias_error(res.resolution_bandwidth, 25.0)  # Br = 25 Hz peak

Narrow Be (long segments) suppresses the bias but leaves fewer averages and a larger random error; the two requirements on segment length pull in opposite directions, which is exactly the trade-off the reported numbers make visible.

Show the code for this figure

import matplotlib.pyplot as plt
import numpy as np
from phonometry import (
    fractional_octave_smoothing,
    noise_signal,
    power_spectral_density,
)

fs = 48000.0
x = noise_signal(fs, 20.0, color="pink", seed=11)
res = power_spectral_density(x, fs, nperseg=4096)
band = (res.frequencies >= 20.0) & (res.frequencies <= 20000.0)
freqs = res.frequencies[band]
smooth = fractional_octave_smoothing(res.frequencies, res.psd, 3.0)[band]

fig, ax = plt.subplots(figsize=(10, 6))
ax.fill_between(freqs, 10 * np.log10(res.ci_lower[band]),
                10 * np.log10(res.ci_upper[band]), alpha=0.3,
                label="95 % chi-square confidence interval")
ax.semilogx(freqs, 10 * np.log10(res.psd[band]), lw=1.0,
            label="Welch PSD estimate")
ax.semilogx(freqs, 10 * np.log10(smooth), lw=2.2,
            label="1/3-octave smoothed")
ax.set_xlabel("Frequency [Hz]")
ax.set_ylabel("PSD [dB re 1/Hz]")
ax.legend()
plt.show()

2. Cross-spectral density

cross_spectral_density estimates the complex Gxy(f) between two channels with the same Welch core, and reports the ordinary coherence γ²xy = |Gxy|²/(Gxx·Gyy) together with the Bendat & Piersol random errors of the magnitude and phase (Eqs. 9.33 and 9.52, with the measured coherence in place of the unknown true value, as the book recommends for measured data):

$${\epsilon }_r[|{\hat{G}}_{xy}|]=\frac{1}{|{\gamma }_{xy}|\sqrt{n_d}},s.d.[{\hat{\theta }}_{xy}]=\frac{{\left(1-{\gamma }_{xy}^2\right)}^{1/2}}{|{\gamma }_{xy}|\sqrt{2n_d}}.$$

Both shrink as the coherence approaches one: a strongly coherent pair needs far fewer averages for the same confidence. The phase is unwrapped, so its slope against frequency is the group delay τ_g = -dφ/(2π·df); for a pure delay path the phase is linear and that slope reads the propagation delay directly.

from phonometry import cross_spectral_density

res = cross_spectral_density(x, y, fs)
print(res.magnitude_random_error[100], res.phase_std[100])  # bin 100 errors
res.plot()   # magnitude, phase with ±sigma band, coherence

The cross-spectral density of a 2 ms delay path: the unwrapped phase is exactly the line −2πfτ, so its slope reads the propagation delay directly, and the ±1 s.d. band of Eq. 9.52 quantifies how far to trust it per frequency.

Show the code for this figure

import matplotlib.pyplot as plt
import numpy as np
from phonometry import cross_spectral_density, noise_signal

fs = 8000.0
tau = 0.002                                   # 2 ms = 16 samples
delay = int(tau * fs)
x = noise_signal(fs, 8.0, seed=8)
noise = noise_signal(fs, 8.0, rms=0.3, seed=9)
y = 0.9 * np.concatenate([np.zeros(delay), x[:-delay]]) + noise

res = cross_spectral_density(x, y, fs)

# One line — magnitude, phase with its ±sigma band and coherence:
res.plot()
plt.show()

# By hand, from the fields the result carries:
band = (res.frequencies >= 20) & (res.frequencies <= 3500)
freqs = res.frequencies[band]
fig, (ax_m, ax_p) = plt.subplots(2, 1, sharex=True)
ax_m.semilogx(freqs, 10 * np.log10(res.magnitude[band]),
              label="|Gxy| (Welch estimate)")
ax_m.set_ylabel("Magnitude [dB]")
ax_p.semilogx(freqs, res.phase[band], label="Unwrapped phase")
ax_p.fill_between(freqs, res.phase[band] - res.phase_std[band],
                  res.phase[band] + res.phase_std[band], alpha=0.25,
                  label="±1 s.d. (Eq. 9.52)")
ax_p.semilogx(freqs, -2 * np.pi * freqs * tau, "r--",
              label="slope -2·pi·f·tau")
ax_p.set(xlabel="Frequency [Hz]", ylabel="Phase [rad]")
for ax in (ax_m, ax_p):
    ax.legend()
plt.show()

3. Coherent output spectrum and spectral SNR

In the single-input/single-output model the measured output autospectrum splits exactly into the part linearly explained by the input and the uncorrelated remainder (Eqs. 9.55–9.57):

$$G_{vv}={\gamma }_{xy}^2{G}_{yy},G_{nn}=\left(1-{\gamma }_{xy}^2\right)G_{yy},\mathrm{SNR}(f)=\frac{{\gamma }_{xy}^2}{1-{\gamma }_{xy}^2}.$$

coherent_output_spectrum returns all three spectra, the spectral signal-to-noise ratio (linear and in dB) and the random error of the coherent output estimate (Eq. 9.73), plus the first-order propagation of the coherence error through the SNR:

$${\epsilon }_r[{\hat{G}}_{vv}]=\frac{{\left(2-{\gamma }_{xy}^2\right)}^{1/2}}{|{\gamma }_{xy}|\sqrt{n_d}},{\epsilon }_r[\widehat{\mathrm{SNR}}]=\frac{\sqrt{2}}{|{\gamma }_{xy}|\sqrt{n_d}}.$$

For additive uncorrelated output noise of known level the coherence has the closed form γ² = SNR/(1+SNR), which makes the whole chain verifiable with a synthetic signal:

import numpy as np
from phonometry import coherent_output_spectrum, noise_signal

fs = 48000.0
x = noise_signal(fs, 8.0, color="white", seed=1)
noise = noise_signal(fs, 8.0, color="white", rms=0.5, seed=2)
y = 0.8 * x + noise                      # SNR = 0.64/0.25 at every frequency

res = coherent_output_spectrum(x, y, fs)
print(np.median(res.coherence))          # -> SNR/(1+SNR) = 0.719
print(np.median(res.snr_db))             # -> 10·lg(2.56) = 4.1 dB
res.plot()                               # Gyy, Gvv, Gnn and the SNR panel

The exact split of the snippet’s model: Gvv = γ²·Gyy explains the output except for the flat noise remainder Gnn, and the spectral SNR scatters around its closed form 10·lg(0.64/0.25) = 4.1 dB at every frequency.

Show the code for this figure

import matplotlib.pyplot as plt
import numpy as np
from phonometry import coherent_output_spectrum, noise_signal

fs = 48000.0
x = noise_signal(fs, 8.0, color="white", seed=1)
noise = noise_signal(fs, 8.0, color="white", rms=0.5, seed=2)
y = 0.8 * x + noise                      # SNR = 0.64/0.25 per band

res = coherent_output_spectrum(x, y, fs, nperseg=2048)

# One line — the three spectra and the SNR panel:
res.plot()
plt.show()

# By hand, from the fields the result carries:
band = (res.frequencies >= 20) & (res.frequencies <= 20000)
freqs = res.frequencies[band]
fig, (ax_g, ax_s) = plt.subplots(2, 1, sharex=True)
for values, style, label in ((res.output_psd, "-", "Gyy (measured)"),
                             (res.coherent_psd, "--", "Gvv (coherent)"),
                             (res.noise_psd, ":", "Gnn (noise)")):
    ax_g.semilogx(freqs, 10 * np.log10(values[band]), style, label=label)
ax_g.set_ylabel("Spectral density [dB re 1/Hz]")
ax_s.semilogx(freqs, res.snr_db[band], label="Spectral SNR [dB]")
ax_s.axhline(10 * np.log10(0.64 / 0.25), color="r", linestyle="--",
             label="closed form 4.1 dB")
ax_s.set(xlabel="Frequency [Hz]", ylabel="SNR [dB]")
for ax in (ax_g, ax_s):
    ax.legend()
plt.show()

The coherence_bias field reports the small positive bias of the coherence estimate, b[γ̂²] ≈ (1-γ²)²/nd (Eq. 9.75), negligible once nd reaches a few hundred, and another reason to average generously before trusting a low coherence.

4. Fractional-octave smoothing

fractional_octave_smoothing averages a spectrum over a rectangular window of constant relative width: 1/n octave, [f·2^(-1/2n), f·2^(+1/2n)] around each frequency. This is the constant-percentage resolution bandwidth that Bendat & Piersol recommend for the spectra of resonant systems (Section 8.5.3), and the de facto standard for presenting loudspeaker and room responses. The average is always computed on power (amplitudes are squared first, dB levels converted and back), so band power is conserved rather than amplitude, and a flat spectrum passes through exactly unchanged.

from phonometry import fractional_octave_smoothing

smooth_psd = fractional_octave_smoothing(res.frequencies, res.psd, 3.0)
smooth_mag = fractional_octave_smoothing(freqs, np.abs(response), 6.0,
                                         domain="amplitude")  # an FRF |H|
smooth_db = fractional_octave_smoothing(freqs, levels, 3.0, domain="db")  # dB curve

A single spectral line with PSD ordinate P (units²/Hz) in a bin of width Δf smooths to the closed-form level P·Δf / (f₀·(2^{1/2n} - 2^{-1/2n})) over one kernel width, the oracle pinned in the tests.

5. Colored-noise generators

noise_signal produces Gaussian noise whose PSD follows Gxx(f) ∝ f^α exactly in expectation: seeded white noise is shaped in the frequency domain by the exact magnitude response (f/f_ref)^{α/2} bin by bin (a zero-phase filter applied circularly), so a measured slope deviates from the power law only by the random error of the spectral estimate, not by the piecewise or few-pole pink approximations whose slope ripples by fractions of a dB. The record is zero-mean and rescaled to the requested RMS exactly, and the same seed reproduces the same record bit for bit.

color	α	PSD slope
white	0	0 dB/octave
pink	-1	-3.01 dB/octave
red (Brownian)	-2	-6.02 dB/octave
blue	+1	+3.01 dB/octave
violet	+2	+6.02 dB/octave

from phonometry import noise_signal

pink = noise_signal(48000, 10.0, color="pink", seed=7)     # deterministic
white = noise_signal(48000, 10.0, color="white", rms=0.5, seed=7)

Measured over three decades (20 Hz – 20 kHz) with the estimator of section 1, the regression slope of each color lands within a few thousandths of a dB/octave of the exact value; the conformance suite pins the pink slope at -3.0116 against the exact -3.0103.

6. Choosing the window

Every estimator on this page accepts any window scipy.signal.get_window knows, but the choice is a quantified trade-off, not a preference. window_metrics computes the figures of merit Harris (1978) tabulated, for any taper and length, sampled DFT-even exactly as the Welch estimators apply it:

• ENBW (equivalent noise bandwidth, in bins): how much wider than one bin the effective analysis bandwidth is. This is the same number the PSD result reports as resolution_bandwidth (ENBW·fs/nperseg in Hz), and it enters directly in the tone/noise trade: a broadband noise floor read from a windowed spectrum sits 10·lg(ENBW) dB above the true density.
• Coherent gain: the DC gain Σw/N a bin-centered tone is scaled by.
• Scalloping loss: the worst-case attenuation of a tone that falls midway between two bins (3.92 dB for rectangular, 1.42 dB for Hann).
• Worst-case processing loss: scalloping plus 10·lg(ENBW), the worst-case reduction in output SNR for tone detection in white noise.
• Highest sidelobe and main-lobe -3 dB width: leakage floor versus resolution.

from phonometry import window_metrics

m = window_metrics("hann", 2048)
print(m.enbw_bins)              # 1.5, exactly
print(m.scalloping_loss_db)     # 1.42 dB
print(m.highest_sidelobe_db)    # -31.5 dB
m.plot()                        # window + spectrum with metrics marked

The closed forms anchor the tests: ENBW is exactly 1 for rectangular, 3/2 for Hann, 1987/1458 for Hamming and 1523/882 for Blackman (DFT-even sampling), and the rectangular scalloping loss is 20·lg(N·sin(π/2N)), the Dirichlet kernel evaluated half a bin off center.

Show the code for this figure

import matplotlib.pyplot as plt
import numpy as np
from phonometry import window_metrics

n, oversample = 1024, 256
fig, ax = plt.subplots(figsize=(10, 6.2))
for name in ("boxcar", "hann", "hamming", "blackman"):
    res = window_metrics(name, n)
    spectrum = np.abs(np.fft.rfft(res.taps, n=n * oversample))
    level = 20.0 * np.log10(spectrum / spectrum[0])
    bins = np.arange(level.size) / oversample
    shown = bins <= 16.0
    ax.plot(bins[shown], level[shown],
            label=(f"{name}: ENBW {res.enbw_bins:.2f} bins, "
                   f"sidelobe {res.highest_sidelobe_db:.1f} dB"))
ax.set_xlim(0.0, 16.0)
ax.set_ylim(-100.0, 5.0)
ax.set_xlabel("Frequency offset [DFT bins]")
ax.set_ylabel("Level re main lobe [dB]")
ax.legend(loc="upper right")
plt.tight_layout()
plt.show()

The Hann default of this module is the balanced choice: fast sidelobe falloff (-18 dB/octave) protects noise spectra from leakage, ENBW 1.5 costs only 1.76 dB against the rectangular window, and its overlap correlation at 50 % keeps nearly all segment information in the effective average count. Reach for a lower-sidelobe taper (Blackman, Kaiser with a high beta) when a weak tone must be found next to a strong one, and accept the wider main lobe; reach for the rectangular window only for self-windowing records (transients that decay inside the segment) or bin-centered synthesis.

7. Multitaper estimation for short records

Welch’s method buys stability with record length: each independent segment adds two degrees of freedom, so a record that only fits a couple of segments leaves an estimate barely better than a periodogram. multitaper_psd implements the alternative of Thomson (1982), as developed by Percival & Walden (1993, Chapter 7): the whole record is multiplied by K orthogonal discrete prolate spheroidal (Slepian) tapers - the sequences that concentrate the most spectral-window energy inside a chosen design band [-W, W] - and the K resulting eigenspectra are averaged:

$${\hat{S}}^{(mt)}(f)=\frac{1}{K}\sum _{k=0}^{K-1}{\hat{S}}_k(f),{\hat{S}}_k(f)=\Delta t|\sum _{t=1}^N{h}_{t,k}{x}_t{e}^{-i2\pi ft\Delta t}{|}^2.$$

Because the tapers are orthogonal the eigenspectra are nearly uncorrelated, so the average carries about 2K chi-square degrees of freedom from a single record - the same statistical machinery as the Welch result (random error, chi-square confidence interval), without segmenting. The half-bandwidth W = NW·fs/N is set through the duration x half-bandwidth product NW (default 4). 2W is the resolution of the estimate (reported as resolution_bandwidth), and only the tapers below the Shannon number 2·NW keep their spectral-window energy inside the design band - their concentrations λk are reported as eigenvalues, and the default taper count is K = 2·NW - 1, all tapers with near-unity concentration. Larger NW admits more tapers (lower variance) at the cost of resolution.

from phonometry import multitaper_psd

res = multitaper_psd(signal, fs)                 # NW = 4, K = 7, adaptive
print(res.degrees_of_freedom.mean())             # ~2K from one record
print(res.eigenvalues)                           # taper concentrations
res.plot()                                       # density with the CI band

By default the eigenspectra are combined with Thomson’s adaptive weights (P&W Eqs. 368a/370a, iterated to convergence). Each taper’s weight at each frequency balances the local spectrum against the broad-band leakage the taper could carry:

$$b_k(f)=\frac{S(f)}{{\lambda }_{k}S(f)+(1-{\lambda }_k){\sigma }^2\Delta t},{\hat{S}}^{(amt)}(f)=\frac{{\sum }_k{b}_k^2(f){\lambda }_k{\hat{S}}_k(f)}{{\sum }_k{b}_k^2(f){\lambda }_k},$$

so the leakier high-order tapers are downweighted exactly where the spectrum is locally weak, and nothing is lost where it is locally white (for white noise the weights are uniform). The price is bookkept honestly: the equivalent degrees of freedom become frequency dependent, ν(f) = 2·(Σk dk)²/Σk dk² with dk = b²k·λk (P&W Eq. 370b), and the confidence interval widens wherever leakage protection spent them. adaptive=False selects the plain eigenvalue-weighted average instead.

Calibration matches the Welch estimators exactly: no detrending, 'density' integrates to the signal power, and 'spectrum' reads A²/2 at the peak of a sinusoid of amplitude A (a tone’s power in 'density' scaling spreads over the 2W band). The Slepian tapers themselves come from scipy.signal.windows.dpss; their concentrations reproduce the quadruple-precision table of Percival & Walden (Table 382) to machine precision, which is the anchor oracle of the test suite.

Show the code for this figure

import matplotlib.pyplot as plt
import numpy as np
from phonometry import multitaper_psd, noise_signal

fs = 48000.0
x = noise_signal(fs, 8192 / fs, color="pink", seed=11)   # 171 ms record
single = multitaper_psd(x, fs, n_tapers=1, adaptive=False)
res = multitaper_psd(x, fs)                              # NW = 4, K = 7
band = (res.frequencies >= 20.0) & (res.frequencies <= 20000.0)
freqs = res.frequencies[band]

fig, ax = plt.subplots(figsize=(10, 6))
ax.semilogx(freqs, 10 * np.log10(single.psd[band]), color="gray",
            alpha=0.45, lw=0.7, label="Single Slepian taper (K = 1)")
ax.fill_between(freqs, 10 * np.log10(res.ci_lower[band]),
                10 * np.log10(res.ci_upper[band]), alpha=0.3,
                label="95 % chi-square confidence interval")
ax.semilogx(freqs, 10 * np.log10(res.psd[band]), lw=1.2,
            label="Multitaper estimate (K = 7, adaptive)")
ax.set_xlabel("Frequency [Hz]")
ax.set_ylabel("PSD [dB re 1/Hz]")
ax.legend()
plt.show()

Reach for multitaper_psd when the record is too short to segment (room impulse response tails, transient captures, single machine cycles) or when a high-dynamic-range spectrum needs leakage protection that a Hann-windowed Welch average cannot give; stay with power_spectral_density for long records, where segment averaging is cheaper than K full-length FFTs and the two estimators agree.

Relation to the H1/H2 estimators

The frequency-response estimators transfer_function and coherence, the two-microphone sound intensity probe and these estimators all share one Welch core (same taper, overlap policy and detrend-off calibration), so a PSD, a coherence and an H1 computed with the same segment length are mutually consistent bin by bin. The same cross-spectral matrix underlies multiple and partial coherence, which extends the ordinary coherence to several correlated inputs and one output.

What this guide covers

Covered. Bendat & Piersol’s Welch estimators (Random Data, 4th ed., 2010, Sections 5.2, 8.5, 9.1-9.2 and 11.5): power_spectral_density and cross_spectral_density with the effective number of averages, chi-square confidence intervals and the resolution-bias error; the coherent output spectrum and spectral SNR; the Harris (1978) window figures of merit of window_metrics; and the Thomson (1982) multitaper estimator of multitaper_psd, following Percival & Walden’s Chapters 7-8. fractional_octave_smoothing and noise_signal complete the toolbox.

Not covered. The frequency-response estimators transfer_function and coherence, and the sound-intensity probe, share this page’s Welch core but are documented on the electroacoustics and sound intensity pages, not here. So is multiple and partial coherence, on the MISO coherence page. This page implements Bendat & Piersol’s textbook estimators, not a certification standard, so it carries no clause numbers or compliance limits to check against.

References

• Bendat, J. S., & Piersol, A. G. (2010). Random data: Analysis and measurement procedures (4th ed.). Wiley. https://doi.org/10.1002/9781118032428Sections 5.2 and 8.5 (autospectra and their random/bias errors, chi-square intervals), 9.1-9.2 (cross-spectra, coherent output spectrum and their errors) and 11.5 (Welch processing, tapering and overlap). ISBN 978-0-470-24877-5.
• Harris, F. J. (1978). On the use of windows for harmonic analysis with the discrete Fourier transform. Proceedings of the IEEE, 66(1), 51-83. https://doi.org/10.1109/PROC.1978.10837The window figures of merit (Table 1): equivalent noise bandwidth, coherent gain, scalloping loss, worst-case processing loss, highest sidelobe level and main-lobe width, computed by window_metrics for any scipy taper.
• Percival, D. B., & Walden, A. T. (1993). Spectral analysis for physical applications: Multitaper and conventional univariate techniques. Cambridge University Press. https://doi.org/10.1017/CBO9780511622762Chapter 7 (multitaper estimation: eigenspectra, adaptive weighting, equivalent degrees of freedom) and Chapter 8 (computing the Slepian sequences); the Table 382 eigenvalues anchor the taper oracle in the test suite. ISBN 978-0-521-43541-3.
• Thomson, D. J. (1982). Spectrum estimation and harmonic analysis. Proceedings of the IEEE, 70(9), 1055-1096. https://doi.org/10.1109/PROC.1982.12433The multitaper method: Slepian tapers, eigenspectra and the adaptive weights implemented by multitaper_psd.
• Welch, P. D. (1967). The use of fast Fourier transform for the estimation of power spectra: A method based on time averaging over short, modified periodograms. IEEE Transactions on Audio and Electroacoustics, 15(2), 70-73. https://doi.org/10.1109/TAU.1967.1161901The overlapped-segment variance formula behind the effective number of averages (Bendat & Piersol Section 11.5.2.2, Ref. 11).