Power Spectrum Fundamentals

• Power spectrum is a representation of how variance or power is distributed across frequencies, typically derived as the Fourier transform of the autocorrelation function or its analogs in various domains.
• It is estimated from finite data using methods like periodograms, Bartlett and Welch averaging, as well as interlacing techniques to minimize aliasing and shot noise effects.
• Its applications span cosmology, graph signal processing, and quantum systems, enabling precision parameter inference, robust noise reduction, and enhanced image and signal analysis.

The power spectrum is a representation of how variance or power is distributed across frequency, wavenumber, or an analogous spectral variable. In the classical wide-sense stationary setting it is the Fourier transform of the autocorrelation function, while in random-field theory it is the Fourier-space two-point statistic, in graph signal processing it is the squared graph Fourier transform, and in quantum spectral statistics it is the Fourier transform of the covariance of level displacements [(Ariananda et al., 2014); (Wang et al., 2024); (Riser et al., 2017); (Djima et al., 12 Mar 2025)]. Across these settings, the power spectrum serves as a compact summary of second-order structure, but its sufficiency, numerical estimation, and interpretation depend strongly on the underlying field, the sampling scheme, and the observational geometry (Neyrinck, 2011).

1. Formal definitions and domain-specific realizations

For a discrete-time wide-sense stationary process $x[n]$ with autocorrelation $R_x[k]=E\{x[n]x^*[n-k]\}$, the power spectral density is

$$P_x(\omega)=\sum_{k=-\infty}^{\infty} R_x[k]\,e^{-j\omega k},$$

and for a continuous-time wide-sense stationary process $x(t)$ with autocorrelation $R_x(\tau)$, it is

$$S_x(f)=\int_{-\infty}^{\infty} R_x(\tau)\,e^{-j2\pi f\tau}\,d\tau.$$

This is the Wiener–Khinchin relation in the forms used for compressive and FFT-based spectral estimation [(Ariananda et al., 2014); (Yang et al., 2019)].

In cosmology, the power spectrum is the Fourier-space two-point statistic of the overdensity field. With $$\delta(\mathbf{x})=\rho(\mathbf{x})/\langle\rho\rangle-1$$, the standard definition is

$$\langle \delta(\mathbf{k})\delta^*(\mathbf{k}')\rangle=(2\pi)^3\delta_D(\mathbf{k}-\mathbf{k}')P(k),$$

and in a periodic finite box this corresponds to $\langle |\delta(\mathbf{k})|^2\rangle=P(k)$ for the discrete Fourier modes, up to window and aliasing corrections [(Wang et al., 2024); (Neyrinck, 2011)].

For 2D imaging, the discrete power spectrum is the squared Fourier amplitude $P[u,v]=|\hat I[u,v]|^2$. Its isotropic or radial version aggregates over the radius $R_x[k]=E\{x[n]x^*[n-k]\}$0, commonly by mean or median within annular bins (Souza, 7 Apr 2026).

For graph signals, with normalized Laplacian $R_x[k]=E\{x[n]x^*[n-k]\}$1 and graph Fourier coefficients $R_x[k]=E\{x[n]x^*[n-k]\}$2, the power spectrum signature is the measure

$R_x[k]=E\{x[n]x^*[n-k]\}$3

or equivalently $R_x[k]=E\{x[n]x^*[n-k]\}$4, which is basis-invariant within degenerate eigenspaces (Djima et al., 12 Mar 2025).

For unfolded quantum spectra, the relevant “frequency” is the discrete index of ordered levels. If $R_x[k]=E\{x[n]x^*[n-k]\}$5 denotes the displacement from the mean position, then the power spectrum is

$R_x[k]=E\{x[n]x^*[n-k]\}$6

with Nyquist frequency $R_x[k]=E\{x[n]x^*[n-k]\}$7 for the discrete sequence (Riser et al., 2017).

Setting	Object	Power-spectrum form
WSS signal	Autocorrelation $R_x[k]=E\{x[n]x^*[n-k]\}$8	$R_x[k]=E\{x[n]x^*[n-k]\}$9
Cosmological field	Overdensity $$P_x(\omega)=\sum_{k=-\infty}^{\infty} R_x[k]\,e^{-j\omega k},$$0	$$P_x(\omega)=\sum_{k=-\infty}^{\infty} R_x[k]\,e^{-j\omega k},$$1
Quantum level sequence	Displacements $$P_x(\omega)=\sum_{k=-\infty}^{\infty} R_x[k]\,e^{-j\omega k},$$2	$$P_x(\omega)=\sum_{k=-\infty}^{\infty} R_x[k]\,e^{-j\omega k},$$3
Graph signal	GFT coefficients	$$P_x(\omega)=\sum_{k=-\infty}^{\infty} R_x[k]\,e^{-j\omega k},$$4
2D image	Fourier amplitude	$$P_x(\omega)=\sum_{k=-\infty}^{\infty} R_x[k]\,e^{-j\omega k},$$5

A central distinction is that the power spectrum completely characterizes Gaussian random fields at the two-point level, whereas in non-Gaussian regimes it remains informative but is no longer sufficient (Neyrinck, 2011).

2. Estimation from finite data

Finite-sample estimation usually proceeds through quadratic forms of Fourier coefficients. Standard nonparametric estimators for time series include the periodogram,

$$P_x(\omega)=\sum_{k=-\infty}^{\infty} R_x[k]\,e^{-j\omega k},$$6

Bartlett averaging over non-overlapping segments, and Welch averaging with overlapping windowed segments. The periodogram is asymptotically unbiased but inconsistent, while Bartlett and Welch trade frequency resolution for variance reduction (Yang et al., 2019).

In cosmological large-scale structure, a common bandpower estimator is

$$P_x(\omega)=\sum_{k=-\infty}^{\infty} R_x[k]\,e^{-j\omega k},$$7

and for Gaussian fields its covariance is diagonal,

$$P_x(\omega)=\sum_{k=-\infty}^{\infty} R_x[k]\,e^{-j\omega k},$$8

Late-time nonlinearity introduces trispectrum terms and off-diagonal covariance, degrading parameter constraints (Neyrinck, 2011).

For event-mode data, the relevant estimator is often the Rayleigh power

$$P_x(\omega)=\sum_{k=-\infty}^{\infty} R_x[k]\,e^{-j\omega k},$$9

or the equivalent Leahy-normalized form $x(t)$0, whose null expectation is $x(t)$1 for uniform phases (Caliandro et al., 2013).

For the CMB, with spherical-harmonic coefficients $x(t)$2, the standard full-sky estimator is

$x(t)$3

with cosmic variance

$x(t)$4

Because only one sky is observed, estimating a smooth underlying $x(t)$5 from a single realization is intrinsically a bias–variance problem (Paykari et al., 2012).

Several methods address low-sample or high-variance regimes by imposing structure beyond naive averaging. Factor-analysis spectral estimation assumes that multiple observed PSDs lie near a low-dimensional subspace spanned by a small number of source spectra; projecting individual multitaper estimates onto that learned subspace improves accuracy for short signals once enough examples are available (Andén et al., 2017). For the CMB, sparsity-based estimation in DCT and wavelet dictionaries yields a smooth estimate of the true underlying spectrum without requiring cosmological parameters; the paper reports that DCT recovers peaks and troughs accurately with only $x(t)$6 of coefficients and that the minimum normalized mean-square error occurs below $x(t)$7 retained coefficients for DCT, whereas WT needs $x(t)$8 (Paykari et al., 2012).

3. Discretization, aliasing, and observed spectra

FFT-based estimation replaces continuous fields by gridded ones, introducing mass-assignment windows and aliasing. In a periodic box with grid spacing $x(t)$9 and Nyquist frequency $R_x(\tau)$0, the measured spectrum is

$R_x(\tau)$1

Near $R_x(\tau)$2, high-frequency power is folded into the fundamental zone, biasing the estimate upward (Wang et al., 2024).

Interlacing mitigates this by averaging two grids shifted by $R_x(\tau)$3. The associated filter

$R_x(\tau)$4

cancels all odd alias images. The recommended estimator subtracts the exact interlaced Poisson shot noise and then deconvolves the mass-assignment window,

$R_x(\tau)$5

With this choice, percent-level accuracy is maintained essentially up to $R_x(\tau)$6; for example, for a down-sampled dark-matter sample with $R_x(\tau)$7, CIC with interlacing and exact shot-noise subtraction reaches $R_x(\tau)$8 error at $R_x(\tau)$9, whereas constant-shot-noise subtraction gives $$S_x(f)=\int_{-\infty}^{\infty} R_x(\tau)\,e^{-j2\pi f\tau}\,d\tau.$$0 errors at $$S_x(f)=\int_{-\infty}^{\infty} R_x(\tau)\,e^{-j2\pi f\tau}\,d\tau.$$1 (Wang et al., 2024).

Observed cosmological power spectra on the light cone are further modified by unequal-time structure. The observed 3D spectrum is constructed by folding the unequal-time angular spectrum along the line of sight. In the wide-depth limit it reduces to the equal-time theoretical spectrum, but unequal-time cross-correlations generate scale-dependent corrections that can amount to a few percent on large scales and high redshifts. These radial-mode effects also introduce anisotropies beyond standard redshift-space distortions, generating non-vanishing odd multiples and imaginary contributions (Raccanelli et al., 2023).

The same paper introduces the frequency-angular power spectrum,

$$S_x(f)=\int_{-\infty}^{\infty} R_x(\tau)\,e^{-j2\pi f\tau}\,d\tau.$$2

a dimensionless observable built from angular and redshift separations בלבד. Its defining advantage is that it does not depend on physical distances and is hence free of Alcock–Paczynski effects (Raccanelli et al., 2023).

4. Transformations, robust summaries, and generalized spectra

In cosmological density fields, non-Gaussian mode coupling motivates transforming the field before measuring its power spectrum. For a roughly lognormal field, the log transform

$$S_x(f)=\int_{-\infty}^{\infty} R_x(\tau)\,e^{-j2\pi f\tau}\,d\tau.$$3

reduces higher-order connected moments and strongly suppresses off-diagonal covariance in the transformed power spectrum. Empirically, this tightens cosmological parameter constraints by a factor of several, with up to a factor of $$S_x(f)=\int_{-\infty}^{\infty} R_x(\tau)\,e^{-j2\pi f\tau}\,d\tau.$$4 reduction in the best case for $$S_x(f)=\int_{-\infty}^{\infty} R_x(\tau)\,e^{-j2\pi f\tau}\,d\tau.$$5 when nonlinear scales up to $$S_x(f)=\int_{-\infty}^{\infty} R_x(\tau)\,e^{-j2\pi f\tau}\,d\tau.$$6 are included (Neyrinck, 2011). In the related Coyote Universe study, the power spectrum of the log-density $$S_x(f)=\int_{-\infty}^{\infty} R_x(\tau)\,e^{-j2\pi f\tau}\,d\tau.$$7 gives the tightest cosmological parameter error bars, yielding a factor of $$S_x(f)=\int_{-\infty}^{\infty} R_x(\tau)\,e^{-j2\pi f\tau}\,d\tau.$$8–$$S_x(f)=\int_{-\infty}^{\infty} R_x(\tau)\,e^{-j2\pi f\tau}\,d\tau.$$9 improvement over the conventional power spectrum in all five parameters tested, again with a factor of $$\delta(\mathbf{x})=\rho(\mathbf{x})/\langle\rho\rangle-1$$0 improvement for $$\delta(\mathbf{x})=\rho(\mathbf{x})/\langle\rho\rangle-1$$1 (Neyrinck, 2011).

For lognormal modeling, the exact real-space mapping is

$$\delta(\mathbf{x})=\rho(\mathbf{x})/\langle\rho\rangle-1$$2

with corresponding Fourier-space conversions $$\delta(\mathbf{x})=\rho(\mathbf{x})/\langle\rho\rangle-1$$3 and $$\delta(\mathbf{x})=\rho(\mathbf{x})/\langle\rho\rangle-1$$4. This formalism shows that the nonlinear increase in small-scale power is significantly smaller for the log-transformed spectrum, which fits the linear spectrum with less than $$\delta(\mathbf{x})=\rho(\mathbf{x})/\langle\rho\rangle-1$$5 error for redshifts down to $$\delta(\mathbf{x})=\rho(\mathbf{x})/\langle\rho\rangle-1$$6 and $$\delta(\mathbf{x})=\rho(\mathbf{x})/\langle\rho\rangle-1$$7; baryon acoustic oscillation smearing is at least three times weaker, though still present (Greiner et al., 2013).

The notion of power spectrum also generalizes beyond Euclidean Fourier analysis. On graphs, the power spectrum signature $$\delta(\mathbf{x})=\rho(\mathbf{x})/\langle\rho\rangle-1$$8 is invariant under graph automorphisms and stable under perturbations. The main stability theorem is

$$\delta(\mathbf{x})=\rho(\mathbf{x})/\langle\rho\rangle-1$$9

where $$\langle \delta(\mathbf{k})\delta^*(\mathbf{k}')\rangle=(2\pi)^3\delta_D(\mathbf{k}-\mathbf{k}')P(k),$$0 is the $$\langle \delta(\mathbf{k})\delta^*(\mathbf{k}')\rangle=(2\pi)^3\delta_D(\mathbf{k}-\mathbf{k}')P(k),$$1-Wasserstein distance. This makes the squared graph Fourier transform a permutation-invariant and basis-invariant descriptor of local geometry and symmetry (Djima et al., 12 Mar 2025).

In image analysis, radial power spectra aggregate $$\langle \delta(\mathbf{k})\delta^*(\mathbf{k}')\rangle=(2\pi)^3\delta_D(\mathbf{k}-\mathbf{k}')P(k),$$2 over angles. Median-based radial aggregation combined with Hann windowing reduces sensitivity to bright compact sources, masking artifacts, and localized features that can bias standard estimators. The slope in $$\langle \delta(\mathbf{k})\delta^*(\mathbf{k}')\rangle=(2\pi)^3\delta_D(\mathbf{k}-\mathbf{k}')P(k),$$3 versus $$\langle \delta(\mathbf{k})\delta^*(\mathbf{k}')\rangle=(2\pi)^3\delta_D(\mathbf{k}-\mathbf{k}')P(k),$$4 is then used as a low-dimensional summary of morphological complexity (Souza, 7 Apr 2026).

5. Applications across cosmology, sensing, and quantum systems

In observational cosmology, power spectra remain central because they interface directly with parameter inference pipelines. A simulation-based emulator for the redshift-space halo power spectrum uses 101 flat-geometry $$\langle \delta(\mathbf{k})\delta^*(\mathbf{k}')\rangle=(2\pi)^3\delta_D(\mathbf{k}-\mathbf{k}')P(k),$$5CDM cosmologies and a feed-forward neural network with two hidden layers of 200 units each to predict $$\langle \delta(\mathbf{k})\delta^*(\mathbf{k}')\rangle=(2\pi)^3\delta_D(\mathbf{k}-\mathbf{k}')P(k),$$6 on a $$\langle \delta(\mathbf{k})\delta^*(\mathbf{k}')\rangle=(2\pi)^3\delta_D(\mathbf{k}-\mathbf{k}')P(k),$$7 grid over $$\langle \delta(\mathbf{k})\delta^*(\mathbf{k}')\rangle=(2\pi)^3\delta_D(\mathbf{k}-\mathbf{k}')P(k),$$8. The reported accuracy is about $$\langle \delta(\mathbf{k})\delta^*(\mathbf{k}')\rangle=(2\pi)^3\delta_D(\mathbf{k}-\mathbf{k}')P(k),$$9 for the monopole and $\langle |\delta(\mathbf{k})|^2\rangle=P(k)$0 for the quadrupole for halos of $\langle |\delta(\mathbf{k})|^2\rangle=P(k)$1, sufficient for SDSS LOWZ- and CMASS-like applications; the emulator is designed to combine with halo occupation models, Finger-of-God damping, and Alcock–Paczynski distortions (Kobayashi et al., 2020).

For cosmic shear, a Bayesian hierarchical model jointly samples the shear field and tomographic $\langle |\delta(\mathbf{k})|^2\rangle=P(k)$2-, $\langle |\delta(\mathbf{k})|^2\rangle=P(k)$3-, and $\langle |\delta(\mathbf{k})|^2\rangle=P(k)$4-mode power spectra. The posterior factorizes into Gaussian map-conditionals and inverse-Wishart power-spectrum conditionals, and a messenger-field Gibbs sampler avoids large dense inversions in the presence of masks and anisotropic noise. In the simulated catalogue analyzed in the paper, the recovered posterior distributions reproduce the underlying simulated $\langle |\delta(\mathbf{k})|^2\rangle=P(k)$5- and $\langle |\delta(\mathbf{k})|^2\rangle=P(k)$6-mode spectra (Alsing et al., 2015).

In wideband sensing, compressive PSD estimation exploits the Wiener–Khinchin relation and multi-coset sampling. A cooperative architecture organizes sensors into groups, each measuring specific temporal lags; a fusion center reconstructs the autocorrelation vector by least squares once the union of circular difference sets covers all lags. The identifiability condition is

$\langle |\delta(\mathbf{k})|^2\rangle=P(k)$7

and non-overlapping circular Golomb rulers minimize the required number of groups (Ariananda et al., 2014). A related FFT-based compressed estimator builds an indicator sequence $\langle |\delta(\mathbf{k})|^2\rangle=P(k)$8 and a zero-padded sampled sequence $\langle |\delta(\mathbf{k})|^2\rangle=P(k)$9, computes lag sums through FFT/IFFT, and reconstructs the PSD with overall complexity $P[u,v]=|\hat I[u,v]|^2$0. For a $P[u,v]=|\hat I[u,v]|^2$1–$P[u,v]=|\hat I[u,v]|^2$2 GHz monitored band, $P[u,v]=|\hat I[u,v]|^2$3 branches, and $P[u,v]=|\hat I[u,v]|^2$4, the paper reports $P[u,v]=|\hat I[u,v]|^2$5 s runtime in MATLAB for $P[u,v]=|\hat I[u,v]|^2$6 kHz resolution and an FPGA estimate of $P[u,v]=|\hat I[u,v]|^2$7 ms for three $P[u,v]=|\hat I[u,v]|^2$8k-point FFTs (Yang et al., 2019).

Power spectra also underpin RF cartography. Spatial–frequency PSD maps $P[u,v]=|\hat I[u,v]|^2$9 are modeled as superpositions of known spectral templates $R_x[k]=E\{x[n]x^*[n-k]\}$00 with unknown spatial coefficients $R_x[k]=E\{x[n]x^*[n-k]\}$01, and learned from low-bit quantized power measurements via nonparametric and semiparametric kernel estimators. The resulting optimization problems take support-vector-machine-type dual forms with $R_x[k]=E\{x[n]x^*[n-k]\}$02-insensitive losses adapted to quantization intervals, and the paper develops both batch and online solvers (Romero et al., 2016). At the hardware level, a separate FPGA design sustains four simultaneous data streams, computes $R_x[k]=E\{x[n]x^*[n-k]\}$03-point FFTs, forms $R_x[k]=E\{x[n]x^*[n-k]\}$04, averages over $R_x[k]=E\{x[n]x^*[n-k]\}$05 blocks, and stores per-bin peaks for quick-look interference detection (Abhyankar et al., 2011).

In quantum-chaotic and complex-wave systems, the power spectrum of level fluctuations has acquired a distinct role as a diagnostic of long-range spectral correlations. For the unitary symmetry class, a nonperturbative random-matrix treatment yields a universal, parameter-free large-$R_x[k]=E\{x[n]x^*[n-k]\}$06 law expressed through a fifth Painlevé transcendent; the exact solution shows that the power spectrum is not determined by the spectral form factor alone, but depends on correlations of all orders (Riser et al., 2017, Riser et al., 2019). Dyson’s circular ensembles extend this program to $R_x[k]=E\{x[n]x^*[n-k]\}$07, leading to the universal small-frequency law

$R_x[k]=E\{x[n]x^*[n-k]\}$08

with finite-$R_x[k]=E\{x[n]x^*[n-k]\}$09 formulas tied to conditioned gap-probability generating functions, Fredholm determinants, and Painlevé VI/III$R_x[k]=E\{x[n]x^*[n-k]\}$10 systems (Forrester et al., 2024).

In quantum linear systems, the output field in the stationary regime is characterized by a power spectrum

$R_x[k]=E\{x[n]x^*[n-k]\}$11

Under pure stationary Gaussian input and global minimality, the transfer function is uniquely recoverable from the power spectrum, and equivalent systems are related by a symplectic transformation (Levitt et al., 2016).

6. Scope, limitations, and methodological boundaries

The power spectrum is a complete descriptor only for Gaussian random fields. In late-time large-scale structure, the density field becomes highly non-Gaussian on small scales, and the power spectrum loses power to describe the field because of substantial off-diagonal covariance generated by nonlinear mode coupling (Neyrinck, 2011). This is the main reason why transformations such as logarithmic Gaussianization can improve inference without changing the observable class.

Exact numerical recovery also depends on restrictive assumptions. The exact interlaced shot-noise separation near the Nyquist limit assumes Poisson sampling, periodic boundary conditions, and the specific mass-assignment windows used in the derivation; halos can be sub-Poisson on large scales, and real surveys require an additional survey-window treatment (Wang et al., 2024). The observed light-cone power spectrum and the frequency-angular power spectrum are derived in linear theory, so extensions to nonlinear bias, relativistic corrections, and full survey systematics remain necessary for precision applications (Raccanelli et al., 2023).

Across engineered sensing problems, information loss can be dominated by quantization, sparse spatial coverage, or kernel misspecification. In PSD-map learning from few-bit measurements, performance depends on interval widths, sensor density, and kernel choice; increasing the number of compressed measurements per sensor does not compensate for inadequate spatial coverage (Romero et al., 2016). In graph settings, the exact power spectrum signature requires eigendecomposition, while scalable approximations rely on Lanczos, polynomial filters, or quantile encodings of the spectral measure (Djima et al., 12 Mar 2025).

In quantum applications, a recurring misconception is that the power spectrum can be reduced to a two-point object such as the spectral form factor. The exact random-matrix analyses explicitly invalidate that simplification: the power spectrum involves generating functions of conditioned gap probabilities and hence the full hierarchy of spectral correlations (Riser et al., 2017, Riser et al., 2019). In quantum linear systems, the converse identifiability result is equally conditional: uniqueness of the transfer function from the power spectrum requires pure input and global minimality, whereas mixed inputs are not treated in the same theorem (Levitt et al., 2016).

A common implication across these literatures is that the power spectrum is best understood not as a single universal estimator, but as a family of second-order summaries whose adequacy depends on the interplay between stochastic assumptions, spectral domain, sampling operator, and the systematic structure of the measurement process.