Power Spectral Density and Nonstationarity

A power spectral density compresses the second-order statistics of a stationary random process into frequency. Detector data are neither perfectly stationary nor perfectly Gaussian, so a PSD is always an estimate tied to an interval, window, averaging rule, and quality selection. This page defines the normalization, derives a practical estimator, and shows how leakage and drift can invalidate apparently precise spectra.

Required background. Probability and Inference for Gravity supplies expectation, covariance, random processes, and sampling distributions.

Helpful background. Gravitational Radiation and Waveforms fixes strain and Fourier conventions. Computational Relativity supplies discretization, convergence, and reproducible numerical workflows.

One-sided spectra and their units

Use the Fourier pair

$$\widetilde x(f)=\int_{-\infty}^{\infty}x(t)e^{2\pi i f t}\,dt, \qquad x(t)=\int_{-\infty}^{\infty}\widetilde x(f)e^{-2\pi i f t}\,df.$$

For a real, zero-mean, wide-sense stationary process $n(t)$, the covariance depends only on lag:

$$C_n(\tau)=\mathbb E[n(t)n(t+\tau)].$$

Its one-sided power spectral density $S_n(f)$ for $f>0$ is defined by

$$\mathbb E\!\left[ \widetilde n(f)\widetilde n^*(f') \right] =\frac12S_n(f)\delta(f-f').$$

Equivalently, the Wiener–Khinchin relation is

$$C_n(\tau)=\int_0^\infty S_n(f)\cos(2\pi f\tau)\,df, \qquad C_n(0)=\int_0^\infty S_n(f)\,df.$$

For dimensionless strain, $S_n$ has units $\mathrm{Hz}^{-1}$ and the amplitude spectral density $\sqrt{S_n}$ has units $\mathrm{Hz}^{-1/2}$. A plotted ASD and a PSD contain the same information but different numerical values and units. The factor of two disappears for a two-sided PSD; convention must therefore accompany every quoted normalization.

Stationarity diagonalizes the covariance in the ideal infinite-duration Fourier basis. It does not imply Gaussianity. Gaussianity specifies the full distribution from the covariance; stationarity only says that second-order statistics do not change under time translation.

A windowed periodogram

Take a segment $x_k(t)$ of duration $T$ and multiply it by a real window $w(t)$. Define

$$X_k(f)=\int_0^T w(t)x_k(t)e^{2\pi i f t}\,dt.$$

For positive frequencies away from DC and, for sampled data, Nyquist, a noise-power-normalized one-sided periodogram is

$$\widehat S_k(f) =\frac{2\left\lvert X_k(f)\right\rvert^2} {\displaystyle\int_0^T w^2(t)\,dt}.$$

The denominator is the window’s noise-power normalization. Using $(\int w\,dt)^2$ instead would calibrate a bin-centered deterministic tone, not broadband noise. For samples $x_j$ separated by $\Delta t$ and an unnormalized discrete transform $X_m=\sum_{j=0}^{N-1}w_jx_je^{2\pi ijm/N}$, the corresponding formula is

$$\widehat S_m=\frac{2\Delta t}{\sum_jw_j^2}\lvert X_m\rvert^2.$$

The factor of two is omitted for DC and the Nyquist bin when those bins are retained. This explicit $\Delta t$ supplies the correct PSD units.

Welch’s estimator averages $K$ overlapping segment periodograms,

$$\widehat S_{\rm W}(f)=\frac1K\sum_{k=1}^{K}\widehat S_k(f).$$

When each segment’s coefficient is circular complex Gaussian with variance set by the local PSD, independent segments give $2K\widehat S_{\rm W}/S_n\sim\chi^2_{2K}$ and fractional standard deviation $1/\sqrt K$. This is exact in the ideal independent-bin model and approximate for finite tapered records whose spectrum varies across the window. Overlapping segments correlate the averages, while tapering correlates neighboring frequency bins; the effective degrees of freedom must reflect those dependencies. A median or trimmed mean is more resistant to glitches, but it has a distribution-dependent bias correction and should not be called an ordinary Welch mean. The original averaged-periodogram construction is given in Welch 1967.

Resolution, windowing, and leakage

Finite duration replaces the true spectrum by a convolution with a spectral window. Define the two-sided density $P_n(f)=S_n(\lvert f\rvert)/2$ away from DC, $W(\nu)=\int_0^T w(t)e^{2\pi i\nu t}\,dt$, and

$$K_w(\nu)=\frac{\lvert W(\nu)\rvert^2}{\int_0^T w^2(t)\,dt}, \qquad \int_{-\infty}^{\infty}K_w(\nu)\,d\nu=1.$$

The expectation at positive frequency is then

$$\mathbb E[\widehat S_k(f)] =2\int_{-\infty}^{\infty}K_w(f-f')P_n(f')\,df'.$$

Writing the convolution in terms of $P_n$ accounts for both frequency signs of a real process. Near a sharp line or a spectral edge, the estimator therefore averages neighboring spectral power instead of measuring a point value of $S_n$.

A rectangular window has a narrow main lobe but high sidelobes. A tapered window suppresses distant leakage at the cost of a wider main lobe and fewer effectively weighted samples. The equivalent noise bandwidth

$$B_{\rm ENBW} =\frac{\displaystyle\int_0^T w^2(t)\,dt} {\left(\displaystyle\int_0^T w(t)\,dt\right)^2}$$

quantifies the noise bandwidth of a tone-calibrated bin. It has units of hertz and exceeds $1/T$ for common tapers. Window properties and normalization conventions are tabulated in Harris 1978.

A narrow spectral line generally does not land exactly on a discrete Fourier bin. Without tapering, its power leaks through sidelobes and can imitate broadband excess. Increasing $T$ narrows the frequency spacing $\Delta f=1/T$, but also assumes stationarity over a longer interval. There is no universally optimal segment length: it balances resolution, estimator variance, and temporal drift.

Local stationarity and time–frequency structure

Real strain data can be approximately stationary over a selected interval without being stationary over an observing run. Slow drift changes the broadband floor and line amplitudes; short glitches add localized excess; wandering lines violate a fixed-frequency model. A spectrogram exposes these effects by plotting $\widehat S(f;t_k)$ over successive segments.

Useful diagnostics compare rather than merely display:

• spectra from neighboring time blocks with the same estimator;
• mean and robust-median averages, including their bias correction;
• cumulative or band-limited power as a function of time;
• standardized segment coefficients $X_k(f)/\sqrt{[S_n(f)/2]\int w^2dt}$;
• results under several segment durations, tapers, overlaps, and quality masks.

If neighboring spectra differ beyond sampling uncertainty, a single stationary likelihood is misspecified. Possible responses include shortening the analyzed interval, using a time-dependent PSD, marginalizing over spectral uncertainty, or excluding a known transition. Quiet-looking data selected after inspecting the candidate can create selection bias, so the rule for choosing off-source data must be declared.

The statistical idealizations and their departures in interferometer data are discussed in LIGO Scientific Collaboration and Virgo Collaboration 2020, §§3–5, PDF.

A reproducible spectral benchmark

Start with a synthetic time series containing colored Gaussian noise, one off-bin sinusoidal line, a short chirp, and a transient impulse. Compute spectra with a rectangular and a tapered window; vary $T$, overlap, and averaging rule. The benchmark should recover the injected broadband PSD away from the line, show the expected line leakage, and reveal how the chirp and impulse contaminate averages.

Then apply the same recorded workflow to a declared public strain interval. Pin the GWOSC dataset and release, detector, GPS interval, file hash, sample rate, quality mask, calibration caveat, software environment, segment duration, taper and normalization, overlap, averaging rule, and output frequency grid. The public interface is frozen to the GWOSC API documentation available on 2026-08-18.

Refine sample rate, $T$, taper, gating choice, and—when supplied—calibration-response samples. Verify Parseval consistency between time-domain variance and integrated one-sided PSD after applying the same band limits. Report failures caused by edge transients, insufficient anti-alias filtering, nonstationarity, or line leakage. This is a benchmark specification, not a claim that a particular public interval was processed for this page.

Common pitfalls

Reading a PSD as an instantaneous property. It summarizes an estimator over selected data. Short transients and drift can be hidden by a long average.

Comparing curves with unmatched estimators. Different tapers, segment durations, robust averages, resolutions, and quality selections can change both lines and broadband levels.

Using the candidate itself to estimate noise without checking contamination. A strong signal or glitch can bias its own PSD and therefore its weight in a likelihood.

Exercises

1. Periodogram uncertainty. An average uses $K=64$ independent Gaussian-noise segments. Estimate the fractional standard deviation of each PSD bin. What changes if 50% overlap makes the effective number of independent averages $K_{\rm eff}=40$?

Solution

For $K$ independent averages, $\operatorname{sd}(\widehat S)/S=1/\sqrt K=0.125$, or $12.5\%$. Replacing $K$ by the effective value gives

$$\frac{1}{\sqrt{40}}\simeq0.158,$$

or $15.8\%$. Counting overlapped segments as independent would understate the uncertainty.

2. Parseval normalization. A dimensionless strain record has variance $4.0\times10^{-42}$ after ideal band limitation to $20$–$500\,\mathrm{Hz}$. What must the one-sided PSD integrate to over that band?

Solution

For the declared one-sided convention,

$$\int_{20\,\mathrm{Hz}}^{500\,\mathrm{Hz}}S_n(f)\,df =4.0\times10^{-42}.$$

The PSD has units $\mathrm{Hz}^{-1}$, so the integral is dimensionless squared. A factor-of-two discrepancy suggests mixing one- and two-sided conventions; another discrepancy can arise if the time-domain variance and spectral integral use different filters or windows.

References

• Gravitational Wave Open Science Center. “API Documentation.” Version available 2026-08-18. Official documentation.
• Harris, Fredric J. “On the Use of Windows for Harmonic Analysis with the Discrete Fourier Transform.” Proceedings of the IEEE 66 (1978): 51–83. DOI.
• LIGO Scientific Collaboration and Virgo Collaboration. “A Guide to LIGO–Virgo Detector Noise and Extraction of Transient Gravitational-Wave Signals.” Classical and Quantum Gravity 37 (2020): 055002. DOI. Open PDF.
• Welch, Peter D. “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 (1967): 70–73. DOI.