API Reference

calculate_spectrum(grid, field, ms::Tuple;
                   transform=AutoSpectralBackend(), execution=AutoBackend(), kwargs...)

Spectral coefficients and physical wavenumbers of field sampled on grid (a FlowGeometries grid). The grid's architecture × geometry is the coordinate system — no coordinate guessing. The two backend axes compose freely:

• transform::SpectralBackends.AbstractSpectralBackend — the spectral math (AutoSpectralBackend (default), FFTSpectralBackend, NUFFTSpectralBackend, FSHTSpectralBackend, NUFSHTSpectralBackend, DirectSumSpectralBackend).

  AutoSpectralBackend picks the fastest transform this grid admits among the extensions loaded, and only ever substitutes one computing the same coefficients: an FFT on a uniform Cartesian grid (asked for its own length per axis), the hybrid FFT/NUFFT composite where a grid is uniform in some directions and stretched in others, a NUFFT on a nonuniform or scattered Cartesian grid, the FastSphericalHarmonics analysis on its own Clenshaw–Curtis grid, and NUFSHT on any other spherical grid. With no fast transform loaded it warns once and runs the direct sum, an O(∏ms · N) / O(L³) correctness reference; name DirectSumSpectralBackend() explicitly to ask for it.

• execution::ComputationalBackends.AbstractExecutionBackend — where/how it runs (SerialBackend, ThreadedBackend, GPUBackend, DistributedBackend, MPIBackend, AutoBackend (default)).

Data model

field is an AbstractArray shaped (spatial…, batch…): the first ndims(grid) dims are spatial (and must equal size(grid)); every trailing dim is a batch dim (components, levels, time, ensemble — any number), carried through and preserved. Structured grids take an (N_1,…,N_D, batch…) tensor (spherical: (nlon, nlat, batch…)); unstructured grids take (N, batch…). A Tuple of equal-shaped arrays stacks them along a new trailing batch axis.

ms is the spectral resolution: Cartesian (m_1,…,m_D); spherical (Nθ, Nφ) with lmax = Nθ-1.

Preprocessing

preprocess::Preprocess detrends, tapers and zero-pads the field first (see preprocess_field); omitting it transforms the field as given, copying nothing. A taper is scaled to unit mean square, so Parseval holds on the returned coefficients with no correction applied afterwards. pad > 1 lengthens the grid's axes and scales ms by the same factor, so the wavenumber spacing narrows.

Returns

(coeffs, ks_phys) — the coefficients and the physical wavenumber coordinates per spectral axis.

The field's element type picks the coefficient layout on both geometries. Cartesian: a real field gives the rfft-packed half (m_1÷2+1, m_2…, batch…), a complex field the full native cube (ms…, batch…); both complex-valued. Spherical: (Nθ, Nφ, batch…), REAL for a real field (the basis is the real spherical harmonics, so the coefficients are real) and complex for a complex one — see sph_coeff_type.

calculate_spectrum!(coeffs, plan::HybridPlan, field) -> ks_phys

Execute a prebuilt hybrid plan in place. field is (N_1…N_D, batch…); coeffs is the packed half for a real field and the full native spectrum for a complex one. The FFTW plan, the per-axis NUFFT plans and the working arrays are reused across calls.

calculate_spectrum!(coeffs, grid, field, ms; transform=DirectSumSpectralBackend(), execution=AutoBackend(), kwargs...)

In-place calculate_spectrum: write coefficients (ms…, batch…) into the preallocated coeffs and return ks_phys. Supported for transform=DirectSumSpectralBackend() with SerialBackend/ThreadedBackend; for the library transforms build a reusable plan with plan_spectrum and call calculate_spectrum!(coeffs, plan, field).

calculate_spectrum!(coeffs, plan::DirectSumSphericalPlan, field) -> ks

Fill preallocated coeffs (Nθ, Nφ, batch…) with the spherical spectrum of field (spatial…, batch…), reusing plan's nodes / ring table / quadrature / Legendre tables.

calculate_spectrum!(coeffs, plan::DirectSumCartesianPlan, field) -> ks

Fill preallocated coeffs (packed (m₁÷2+1, m₂…, batch…) for a real field, full native for a complex one) with the Cartesian spectrum of field (spatial…, batch…), reusing plan's DFT matrices, working arrays and twin storage.

synthesize(grid, coeffs, ms::Tuple; transform=AutoSpectralBackend(), execution=AutoBackend(),
           real_output=true, iflag=1)

Inverse of calculate_spectrum: reconstruct field values at the grid points from its coefficients. real_output=true consumes the packed half a real Cartesian field transforms to ((m_1÷2+1, m_2…, batch…)) and writes a real array directly, folding each stored mode's conjugate partner into the sum; real_output=false consumes the full native spectrum (ms…, batch…) and returns it complex. Spherical coefficients are (Nθ, Nφ, batch…) either way. Returns an array (spatial…, batch…) shaped as the forward's input.

transform selects the inverse the same way calculate_spectrum selects the forward, and AutoSpectralBackend resolves it by the same rules; the direct-sum inverse serves any grid.

plan_spectrum(grid, ::Type{T}, ms; transform=AutoSpectralBackend(), execution=AutoBackend(), batch=(), kwargs...)
plan_spectrum(transform, execution, grid, ::Type{T}, ms; batch=(), kwargs...)

Construct a reusable AbstractSpectralPlan for the transform×execution backend pair on grid at spectral resolution ms, transforming a field with trailing batch shape batch of element type T in one batched execution. The keyword form resolves execution and forwards to the canonical positional form implemented by the backend extensions. Requires the transform's extension to be loaded.

Execute a plan with calculate_spectrum!(coeffs, plan, fields).

AbstractSpectralPlan

Supertype for reusable transform plans. A plan is tied to the fixed geometry of a problem — the grid coordinates, the spectral resolution ms, the trailing batch shape, and the element type — but not to the field values. Build a plan once with plan_spectrum and reuse it across many fields / batch slices / time steps via calculate_spectrum!, avoiding repeated FFTW/FINUFFT plan construction and point sorting.

Concrete plan types are defined in the backend extensions (e.g. the FFTW and FINUFFT extensions); this module only declares the shared interface.

sph_mode_index(l, m) -> CartesianIndex

CartesianIndex of degree l, order m in the (lmax+1, 2lmax+1) coefficient array.

sph_coeff_type(::Type{T}, ::Type{FT}) -> Type

Element type of a spherical coefficient array for a field of element type T on a grid of element type FT: FT for a real field and Complex{FT} for a complex one.

The basis is the real spherical harmonics, so a real field's coefficients are real and a real array holds them exactly. This mirrors the Cartesian rule, where a real field's coefficients are the packed Hermitian half and a complex field's the full native cube.

coefficient_size(plan) -> Dims

Shape of the coefficient array calculate_spectrum!(coeffs, plan, field) fills, batch axes included.

A real Cartesian field's is the packed half (ms[1]÷2+1, ms[2:D]…, batch…) and a complex one's the full native (ms…, batch…); a spherical plan's is (lmax+1, 2lmax+1, batch…).

coefficient_type(plan) -> Type

Element type of that array. Cartesian coefficients are complex; a spherical plan's follow the field, real for a real one, because the basis is the real spherical harmonics.

An allocation needs this as well as coefficient_size, and allocate_coefficients takes both.

wavenumbers(plan) -> Tuple

The physical wavenumber axes the matching calculate_spectrum returns as its second value, in the same packed native order as the coefficients: rfftfreq-like on a halved axis, fftfreq on a full one, and (0:lmax, -lmax:lmax) for a spherical plan.

A halved axis carries the Nyquist twin the transform attached, so the axes a reduction needs come from here as well.

allocate_coefficients(plan) -> Array

A zeroed coefficient array of this plan's own size and element type, ready for calculate_spectrum!(coeffs, plan, field).

AbstractSynthesisPlan

Supertype for reusable inverse-transform plans, the counterpart to AbstractSpectralPlan. A synthesis plan is tied to the grid, the spectral resolution ms, the trailing batch shape, the element type, and whether the output is real — everything except the coefficient values.

Build one with plan_synthesis and reuse it across many coefficient sets via synthesize!, so a snapshot series on one grid builds its backward transform once. The forward and inverse are separate objects, so a caller that only inverts never builds a forward plan.

plan_synthesis(grid, ::Type{T}, ms; transform=AutoSpectralBackend(), execution=AutoBackend(),
               batch=(), real_output=true, iflag=1, kwargs...)
plan_synthesis(transform, execution, grid, ::Type{T}, ms; batch=(), real_output=true, kwargs...)

Construct a reusable AbstractSynthesisPlan inverting an ms-resolution spectrum on grid back to a field with trailing batch shape batch. T is the FIELD's element type, which fixes both the coefficient layout the plan consumes (the packed half for a real field, the full native cube for a complex one) and the output it writes.

Execute with synthesize!(out, plan, coeffs).

synthesize!(out, plan::AbstractSynthesisPlan, coeffs; ks=nothing) -> out

Fill preallocated out with the inverse transform of coeffs, reusing everything plan holds. out is field_size-shaped of field_type; allocate_field provides one.

ks is the wavenumber tuple the forward returned. A plan holds what the GRID fixes, while the Nyquist twin a packed inverse needs on a nonuniformly-sampled grid is a functional of the coefficients, so it arrives with them on ks[1]. A uniform grid needs none.

field_size(plan::AbstractSynthesisPlan) -> Dims

Shape of the field synthesize! writes: the grid's own spatial shape followed by the batch axes, so a round trip returns the shape the forward transform consumed.

field_type(plan::AbstractSynthesisPlan) -> Type

Element type of that field — real when the plan was built for a real field, complex otherwise.

allocate_field(plan) -> Array

A zeroed field array of this plan's own size and element type, ready for synthesize!(out, plan, coeffs).

unpacked(coeffs, ns::NTuple{D,Int}, ks = nothing) -> full

Full native-order complex spectrum (ns…, batch…) reconstructed from an rfft-packed half coeffs (ns[1]÷2+1, ns[2:D]…, batch…) via C[k] = conj(C[-k]). Passing the ks that calculate_spectrum returned supplies the halved axis's NyquistTwin, which the k₁ < 0 rows need wherever an even axis d ≥ 2 sits at −n_d/2: negating that index lands on +n_d/2, off the native axis. With ks omitted, those entries take the negated-index value, which matches the twin under uniform sampling.

unpacked!(full, coeffs, ns::NTuple{D,Int}, ks = nothing) -> full

Fill full (ns…, batch…) with the complete native-order spectrum of the rfft-packed half coeffs (ns[1]÷2+1, ns[2:D]…, batch…), by Hermitian symmetry over the spectral dims 1:D.

Each mode is either stored outright or the conjugate of its negated-frequency partner, whose per-axis index is j ↦ 1 for j == 1 and n-j+2 otherwise. Where an even axis d ≥ 2 sits at −ns[d]/2 that negation lands off the native axis, so the partner comes from the halved axis's NyquistTwin, carried on ks; on a nonuniformly sampled grid ks is required for an exact result. full and coeffs share rank, and the batch dims ride along.

isotropic_spectrum(ks_phys::Tuple, coeffs; num_bins=0, dims=(), convention=SpectralConvention())

1D radially-integrated (isotropic) energy spectrum of coeffs (ms…, batch…), binning over the D = length(ks_phys) spectral dims and preserving all batch dims → (k_bins, E) with E of shape (num_bins, batch…). dims (absolute batch-dim index/indices > D) folds those axes into the energy (e.g. dims=D+1 sums a vector-component axis into a single kinetic-energy spectrum). E(k) = (1/2dk) Σ_{|k|∈bin} Σ_{folded} |C|².

convention.scaling selects DensityScaling() (divide by the bin width, so Σ E·dk recovers the variance — the default) or PowerScaling() (per-bin power). A radial bin covers every direction, so its Σ already includes each mode's −k partner and convention.sided is fixed at OneSided() here; see Normalization.SpectralConvention.

isotropic_spectrum!(E, k_bins, ks_phys, coeffs; num_bins=0, convention=SpectralConvention())

In-place, allocation-free isotropic spectrum that preserves all batch dims (no folding): fills preallocated E (shape (num_bins, batch…)) and k_bins (length num_bins). Reusable across a time loop with zero steady-state heap traffic. convention is read as in isotropic_spectrum.

transect_spectrum(ks_phys::Tuple, coeffs, dims::Tuple)

Integrate the spectral energy density ½|C|² along the spectral dimensions dims (1-indexed, ⊆ 1:D), scaling by their wavenumber spacing. Returns (ks_reduced, E_reduced) where E_reduced has the kept spectral dims followed by the batch dims.

Σ E_reduced / ∏dk over the kept axes is the field's folded Parseval total, whichever axes are kept: a kept full axis carries both signs of its wavenumber outright, and a kept halved axis reports |k₁| with the −k₁ energy folded in. So a kept full axis's entries at ±k are equal for a real field's spectrum and differ for a complex field's, where the sign carries the propagation direction.

This is the one reduction with no radial cutoff, so it reaches the −N_d/2 modes whose k₁ < 0 partners are the +N_d/2 twins. It reads those from the halved axis's Packing.NyquistTwin, which the transform attaches whenever index negation cannot reach them.

transect_spectrum!(E_reduced, ks_phys, coeffs, dims) -> nothing

In-place, allocation-free transect_spectrum: fills preallocated E_reduced (kept spectral dims + batch dims) with the dims-integrated ½|C|² density.

spherical_energy_spectrum(coeffs; lmax=size(coeffs,1)-1)

Degree energy spectrum E(ℓ, batch…) = ½ Σ_{m=-ℓ}^{ℓ} |C_ℓ^m|² of spherical-harmonic coefficients (Nθ, Nφ, batch…). Returns (0:lmax, E_l), E_l of shape (lmax+1, batch…).

Reads real coefficients (what a real field transforms to, the real spherical harmonics being real) and complex ones alike; E_l carries the real type either way.

spherical_energy_spectrum!(E_l, coeffs; lmax=size(coeffs,1)-1) -> nothing

In-place spherical_energy_spectrum: fills preallocated E_l (shape (lmax+1, batch…)).

anisotropic_spectrum(ks_phys::Tuple, coeffs; num_k_bins=0, num_θ_bins=16, dims=())

Anisotropy-resolved 2D energy spectrum E(k, θ, batch…) for a 2D field: bin ½|C|² by wavenumber magnitude and polar angle, preserving batch. Integrating over θ recovers the isotropic spectrum.

cross_spectrum(ks_phys::Tuple, coeffs_f, coeffs_g; num_bins=0, dims=())

Radially-binned cross-spectrum S_fg(k, batch…) = ½ Σ_{|k|∈bin} Σ_{folded} f̂ conj(ĝ); coeffs_f, coeffs_g share shape (ms…, batch…). Real part → co-spectrum, negative imag part → quad spectrum.

cospectrum(ks, cf, cg; …) — Re S_fg(k, batch…) (in-phase, flux-carrying part).

quadspectrum(ks, cf, cg; …) — -Im S_fg(k, batch…) (90°-out-of-phase part).

welch_power_spectrum(ks_phys::Tuple, coeffs; num_bins=0)

Variance-reduced (Welch / ensemble-averaged) isotropic power spectrum. The trailing batch dims of coeffs (ms…, realization…) index independent segments/realizations whose periodograms are averaged before radial binning. Returns (k_bins, E_k).

welch_power_spectrum!(E_k, k_bins, ks_phys, coeffs; num_bins=0) -> nothing

In-place, allocation-free welch_power_spectrum: fills preallocated E_k and k_bins (both length num_bins). Reusable across a loop with zero steady-state heap traffic.

coherence_spectrum(ks_phys::Tuple, cf, cg; num_bins=0) -> (k_bins, coherence², phase)

Magnitude-squared coherence $\gamma^2(k) = |S_{fg}|^2 / (S_{ff} S_{gg})$ and phase between two fields whose coefficients cf, cg share (ms…, realization…). Cross/auto spectra are averaged over the realization batch and over the modes in each radial bin before the ratio is formed.

coherence_spectrum!(coherence², phase, k_bins, ks_phys, cf, cg; num_bins=0) -> nothing

In-place coherence_spectrum: fills preallocated coherence², phase, k_bins (each length num_bins). Uses O(num_bins) internal scratch for the complex cross-spectrum accumulator.

lomb_scargle(t, y, freqs; center=true) -> Vector

Lomb–Scargle periodogram of an irregularly-sampled 1D series y at sample times/locations t, evaluated at the (strictly positive) frequencies freqs. This is the standard estimator for gappy / non-uniformly sampled records (moorings, drifters, satellite tracks, astronomical time series) where an FFT cannot be applied directly.

With center=true (default) the sample mean is removed first. The classic time-shift τ is chosen per frequency to make the cosine and sine bases orthogonal, giving a periodogram that is invariant to time translation:

\[P(f) = \tfrac12\left[ \frac{(\sum_j y_j\cos\omega(t_j-\tau))^2}{\sum_j\cos^2\omega(t_j-\tau)} + \frac{(\sum_j y_j\sin\omega(t_j-\tau))^2}{\sum_j\sin^2\omega(t_j-\tau)} \right], \quad \omega = 2\pi f .\]

This is the direct $O(N\,M)$ reference implementation (N samples, M frequencies); it needs no FFT. freqs must be positive (the f=0 term is undefined).

lomb_scargle!(P, t, y, freqs; center=true) -> P

In-place lomb_scargle: writes the periodogram into preallocated P (length length(freqs)). Allocation-free — reuse P across many series that share freqs.

spectral_divergence(ks_phys::Tuple, coeffs) -> AbstractArray

Spectral divergence $\widehat{\nabla\cdot u} = i\sum_d k_d \hat u_d$ of a D-component vector field with coefficients (ms…, D, extra…). Returns (ms…, 1, extra…). Defined for D = 1, 2, 3.

spectral_divergence!(out, ks_phys::Tuple, coeffs) -> out

In-place spectral_divergence: writes the divergence into preallocated out of shape (ms…, 1, extra…). Allocation-free — reuse out across a batch/time loop.

spectral_vorticity(ks_phys::Tuple, coeffs) -> AbstractArray

Spectral vorticity $\hat\omega = i\,k \times \hat u$ of a vector field with coefficients (ms…, D, extra…): D = 2 → scalar out-of-plane vorticity (ms…, 1, extra…); D = 3 → 3-component vorticity (ms…, 3, extra…).

spectral_vorticity!(out, ks_phys::Tuple, coeffs) -> out

In-place spectral_vorticity: writes the vorticity into preallocated out — shape (ms…, 1, extra…) for D = 2, (ms…, 3, extra…) for D = 3. Allocation-free.

compensate(k_bins, E_k, p) -> AbstractArray

Compensated spectrum $k^p E(k)$ (e.g. p = 5/3 Kolmogorov plateau, p = 2 for Z(k)=k²E(k)). E_k may be (num_bins,) or (num_bins, batch…); k_bins broadcasts along the wavenumber axis.

band_energy(k_bins, E_k, k1, k2) -> Real

Energy integrated over the wavenumber band $[k_1, k_2]$ via the trapezoidal rule over the bins whose centers fall in the band. E_k is a 1D spectrum (num_bins,).

preprocess_field(grid, field, spec::Preprocess) -> (grid′, field′, ms_scale)

field detrended, tapered and zero-padded per spec, with the grid and mode-count scaling the padding implies. Returns (grid, field, 1) unchanged for a spec that is the identity, so an unpreprocessed call copies nothing.

The caller's array is never mutated: a spec with anything to do works on a copy.

spec.window tapers each spatial axis (the tensor product of Preprocessing.axis_taper) and spec.pad > 1 extends each axis with zeros. Both read the grid's axes, so both apply to a structured grid; a node cloud has no axis to taper or extend along, and asking for either raises. spec.detrend needs no axes for Demean, so that applies to any grid.

Preprocessing lives outside plan_spectrum: a plan is fixed by the grid and the resolution while a taper is an estimator choice, and multitaper reuses ONE plan across K differently tapered copies of the same field transformed as a batch.

Preprocess(; detrend=Demean(), window=NoWindow(), pad=1.0)

Preprocessing applied to a field (per spectral axis) before transforming. Fields are typed (not symbols) so downstream code dispatches at compile time.

• detrend::AbstractDetrend: Demean() (default), NoDetrend(), or LinearDetrend().
• window::AbstractWindow: NoWindow() (default), Hann(), Hamming(), Blackman(), Tukey(α).
• pad::Float64: zero-padding factor (≥ 1); 1.0 means none.

Window power/amplitude corrections for variance preservation are applied by the normalization layer via window_correction.

AbstractWindow

Supertype for apodization tapers applied per spectral axis before transforming. Concrete windows dispatch window_function!. Reduces spectral leakage for non-periodic data; only meaningful on uniform axes.

NoWindow() — rectangular window (all ones); no apodization.

Hann() — Hann (raised-cosine) taper.

Hamming() — Hamming taper.

Blackman() — Blackman taper.

Tukey(alpha=0.5)

Tukey (tapered-cosine) window with taper fraction alpha ∈ [0,1]; alpha=0 is rectangular and alpha=1 is Hann.

AbstractDetrend

Supertype for detrending operations applied before transforming. Concrete subtypes dispatch detrend!.

NoDetrend() — leave the data unchanged.

Demean() — subtract the mean (remove the DC component). The default.

LinearDetrend() — subtract the least-squares linear trend.

dpss(N::Integer, NW::Real, K::Integer = floor(Int, 2NW) - 1; T = Float64) -> Matrix{T}

Discrete prolate spheroidal sequences (Slepian tapers): the K length-N sequences with maximal spectral concentration in the half-bandwidth W = NW/N, returned as an N×K orthonormal matrix whose column k is the order-(k-1) taper. NW is the time–bandwidth product (typical 2.5–4); K ≈ 2·NW − 1 tapers are usefully concentrated.

Multitaper power spectral estimation reuses the ensemble machinery: apply each taper to the (demeaned) signal, transform the K tapered copies as a batch, and average their periodograms with [welch_power_spectrum]. The tapers are the eigenvectors of a symmetric tridiagonal matrix (Percival & Walden), computed here with LinearAlgebra.eigen.

SpectralConvention(; sided=OneSided(), scaling=DensityScaling(), parseval_check=false)

Convention governing how spectral coefficients become reported spectra, so the package never silently guesses normalization. Fields are typed for compile-time dispatch.

• sided::AbstractSidedness: OneSided() (default) or TwoSided().
• scaling::AbstractScaling: DensityScaling() (default) or PowerScaling().
• parseval_check::Bool: when true, callers assert Σ E·Δk ≈ Var(field) (after mean removal) as a correctness self-test.

The variance-preservation property — ∫ E(k) dk = Var(f) after demeaning — is the invariant the test suite enforces across every backend and grid.

AbstractSidedness

Whether a spectrum keeps both signs of wavenumber (TwoSided) or folds negatives onto positives (OneSided). Dispatches sided_factor.

OneSided() — fold negative wavenumbers onto positives (doubles interior bins). The usual convention for real fields.

TwoSided() — keep ± wavenumbers (no folding).

AbstractScaling

Whether a reduced spectrum is reported as a spectral density (DensityScaling, divided by the bin width so ∫E dk recovers variance) or as per-bin/per-mode PowerScaling.

DensityScaling() — spectral density (per unit wavenumber); ∫E dk = Var(f).

PowerScaling() — per-bin/per-mode power (no dk division).

TransformProblem{NS, B}

Shape of one transform: the leading NS spatial array dims (matching the grid) and the trailing B batch dims. Built from (grid, field) via TransformProblem; drives buffer ranks and the coefficient output size.

FINUFFTBackend()
NonuniformFFTsBackend()

The concrete NUFFT providers (<: SpectralBackends.AbstractNUFFTSpectralBackend). A provider is a library choice, not spectral math, so these live in FlowFieldSpectra, not SpectralBackends. They are symmetric — neither is a default: FINUFFTBackend() uses FINUFFT (using FINUFFT; GPU via cuFINUFFT), NonuniformFFTsBackend() uses NonuniformFFTs (using NonuniformFFTs; real-data fast path, half memory). Both may be loaded at once. The abstract SpectralBackends.NUFFTSpectralBackend() selects no provider — pass one of these.

See FINUFFTBackend — the NonuniformFFTs.jl NUFFT provider.

plot_spectrum(ks_phys::Tuple, coeffs; title="Energy Spectrum", kwargs...)

Plot a 1D isotropic / 2D Cartesian / spherical-degree energy spectrum. Requires using CairoMakie.

compare_spectra(spectra_list; labels, kwargs...)

Overlay multiple 1D energy spectra. Requires using CairoMakie.

compare_spectral_analysis(true_coeffs, approx_coeffs; kwargs...)

Coefficient comparison + error maps. Requires using CairoMakie.