API Reference

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

Spectral coefficients and physical wavenumbers of field sampled on grid. The coordinate system is the grid type — there is no coordinate guessing. The two backend axes compose freely:

• transform::AbstractSpectralBackend — the spectral math (DirectSumBackend (default), FFTBackend, NUFFTBackend, SHTBackend, NUFSHTBackend).
• execution::AbstractExecutionBackend — where/how it runs (SerialBackend, ThreadedBackend, GPUBackend{B}, DistributedBackend{Inner}, MPIBackend{Inner}, AutoBackend (default)).

Data model

field is an AbstractArray shaped (spatial…, batch…): the first ndims_spatial(grid) dims are spatial (and must equal spatial_size(grid)); every trailing dim is a batch dim (components, levels, time, ensemble — any number), carried through and preserved. Tensor-product grids take an (N_1,…,N_D, batch…) tensor; scattered grids take (N, batch…); structured spherical grids take (Nθ, Nφ, batch…). A Tuple of equal-shaped arrays is a convenience that 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.

Returns

(coeffs, ks_phys) — complex coefficients (ms…, batch…) (Cartesian) / (Nθ, Nφ, batch…) (spherical), and the physical wavenumber coordinates per spectral axis.

Example

using FlowFieldSpectra, FFTW
g = UniformCartesianGrid(; domain = (2π, 2π), n = (64, 64))
u = rand(64, 64, 3)                       # (Nx, Ny) with a 3-slice batch
coeffs, ks = calculate_spectrum(g, u, (64, 64); transform = FFTBackend())  # coeffs (64,64,3)

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

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

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

Inverse of calculate_spectrum: reconstruct field values at the grid points from coefficients coeffs ((ms…, batch…) Cartesian, (Nθ, Nφ, batch…) spherical). Returns an array (spatial…, batch…); real_output=true returns its real part. Uses the direct-sum inverse (works for any grid) with SerialBackend/ThreadedBackend.

plan_spectrum(grid, ::Type{T}, ms; transform=DirectSumBackend(), execution=AutoBackend(), n_transf=1, kwargs...)
plan_spectrum(transform, execution, grid, ::Type{T}, ms; n_transf=1, kwargs...)

Construct a reusable AbstractSpectralPlan for the transform×execution backend pair on grid at spectral resolution ms, transforming n_transf co-located fields/slices 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 number of batched transforms n_transf, 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.

AbstractGrid{FT, D}

Abstract supertype for all coordinate grids. FT is the coordinate element type and D the number of physical (spatial, transformed) dimensions. The grid type is the coordinate system, so backends dispatch on it rather than guessing from coordinate magnitudes.

AbstractCartesianGrid{FT, D} <: AbstractGrid{FT, D}

Cartesian grids in D dimensions (uniform / nonuniform tensor-product, or scattered).

AbstractSphericalGrid{FT} <: AbstractGrid{FT, 2}

Spherical $(\theta, \phi)$ grids (structured-quadrature or scattered).

UniformCartesianGrid(axes::NTuple{D,AbstractRange}; domain_size=nothing)
UniformCartesianGrid(axis::AbstractRange; domain_size=nothing)
UniformCartesianGrid(; domain, n)

Tensor-product Cartesian grid with uniform spacing along every axis — the FFT-eligible grid. Each axes[d] is a 1-D AbstractRange of length N_d (one value per grid line, not prod(N) per-point coordinates). A field on this grid is a D-dimensional tensor (N_1, …, N_D, batch…).

domain_size[d] is the physical period along axis d; when omitted it is step(axes[d]) * N_d (the periodic-domain convention, so range(0, L, N+1)[1:N] recovers L). The domain=…, n=… keyword form builds range(0, L_d, n_d + 1)[1:n_d] for each axis.

NonuniformCartesianGrid(axes::NTuple{D,AbstractVector}; domain_size=nothing)
NonuniformCartesianGrid(axis::AbstractVector; domain_size=nothing)

Tensor-product Cartesian grid with nonuniform spacing along one or more axes (e.g. Gaussian latitudes, a stretched vertical grid). Still fully gridded: each axes[d] is a 1-D coordinate vector of length N_d, and a field is a tensor (N_1, …, N_D, batch…). FFTW is not valid; use NUFFTBackend or DirectSumBackend. domain_size defaults to the per-axis coordinate span.

ScatteredCartesianGrid(coords::NTuple{D,AbstractVector}; domain_size=nothing)

Arbitrary scattered points in D-dimensional Cartesian space (a genuine point cloud — no product structure). coords[d] is the length-N vector of the d-th coordinate of every point; a field is (N, batch…). This per-axis-vector layout is what NUFFT/cuFINUFFT consume directly. Suitable for NUFFTBackend ($D \le 3$) and DirectSumBackend (any D). domain_size defaults to the coordinate bounding box.

StructuredSphericalGrid(θ, φ; weights=nothing, quad=ClenshawCurtis())

Structured spherical quadrature grid given by a colatitude axis θ (length Nθ) and a longitude axis φ (length Nφ) — a tensor-product (Nθ, Nφ) grid. A field is (Nθ, Nφ, batch…). weights are the per-colatitude quadrature weights (length Nθ) or nothing. Suitable for SHTBackend.

ScatteredSphericalGrid(θ, φ; weights=nothing)

Arbitrary scattered points $(\theta, \phi)$ on the sphere (length-N per-point vectors). A field is (N, batch…). Suitable for NUFSHTBackend and DirectSumBackend.

gauss_legendre_sphere(lmax; FT=Float64) -> StructuredSphericalGrid

Structured spherical grid with Gauss–Legendre colatitude nodes (Nθ = lmax+1) and uniform longitude (Nφ = 2·lmax+1), carrying the exact Gauss quadrature weights. On this grid the direct and GPU spherical-harmonic transforms are exact (the quadrature is exact to degree 2·lmax+1), matching FastSphericalHarmonics. Use this grid for SHTBackend × GPUBackend and for an exact DirectSumBackend spherical transform.

AbstractQuadrature

Quadrature scheme for a structured spherical grid. Dispatched on type, not a symbol.

ClenshawCurtis() — Clenshaw–Curtis latitude nodes (default).

GaussLegendre() — Gauss–Legendre latitude nodes.

Equiangular() — equiangular latitude nodes.

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

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|².

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

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.

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.

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…).

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).

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.

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).

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_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…).

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(; 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.

AbstractSpectralBackend

Abstract supertype for transform backends: which spectral math maps physical↔spectral coefficients. Orthogonal to AbstractExecutionBackend. Concrete subtypes dispatch the calculate_spectrum interface to different mathematical methods and third-party libraries.

DirectSumBackend <: AbstractSpectralBackend

Slow, dependency-free reference transform that computes the Discrete Fourier Transform (DFT) or Spherical Harmonic Transform (SHT) directly using $O(N \cdot M)$ direct summation. This backend is fully self-contained and requires no external packages to be loaded. It is the default transform.

Details

• Cartesian coordinates: Computes the exact Discrete Fourier Transform (DFT) at the target frequencies.
• Spherical coordinates: Computes SHT coefficients using a direct projection onto the Spherical Harmonic basis, with associated Legendre polynomials computed via a type-stable recurrence relation.
• Complexity: $O(N \cdot M)$ where $N$ is the number of spatial grid points and $M$ is the number of spectral modes.

FFTBackend <: AbstractSpectralBackend

Fast Fourier Transform backend for uniform Cartesian grids. Leverages FFTW.jl (via a package extension) to achieve optimal performance.

Requirements

To use this backend, you must import FFTW in your script:

using FFTW

Details

• Grid requirements: Grids must be uniform and rectilinear (Cartesian). Coordinates should represent grid axes rather than scattered point lists.
• Complexity: $O(N \log N)$ where $N$ is the number of grid points.

NUFFTBackend <: AbstractSpectralBackend

Non-Uniform Fast Fourier Transform (NUFFT) backend for non-uniform / scattered Cartesian grids. Leverages FINUFFT.jl (via a package extension). On a GPU execution backend (GPUBackend(CUDABackend())) this routes to cuFINUFFT for a fast device NUFFT.

Requirements

To use this backend, you must import FINUFFT in your script:

using FINUFFT

Details

• Grid requirements: Scattered / non-uniform Cartesian coordinates.
• Complexity: $O(N \log N + M \log(1/\epsilon))$ where $N$ is the number of points, $M$ is the number of modes, and $\epsilon$ is the target accuracy.
• Parameters: Supports passing an accuracy parameter eps (defaults to 1e-8).

SHTBackend <: AbstractSpectralBackend

Spherical Harmonic Transform backend for uniform / structured spherical grids. Leverages FastSphericalHarmonics.jl (via a package extension) for high-performance SHT on equiangular and Clenshaw-Curtis grids.

Requirements

To use this backend, you must import FastSphericalHarmonics in your script:

using FastSphericalHarmonics

Details

• Grid requirements: Latitude/longitude grids structured specifically for Clenshaw-Curtis quadrature nodes.
• Complexity: $O(L^3)$ or $O(L^2 \log L)$ where $L$ is the maximum spherical degree (lmax).

NUFSHTBackend <: AbstractSpectralBackend

Non-Uniform Fast Spherical Harmonic Transform (NUFSHT) backend for unstructured/scattered spherical grids. Leverages NUFSHT.jl (via a package extension).

Requirements

To use this backend, you must import NUFSHT in your script:

using NUFSHT

Details

• Grid requirements: Arbitrary scattered coordinates $(\theta, \phi)$ on the sphere.
• Complexity: $O(M \log M + N \log(1/\epsilon))$ using Double Fourier Sphere (DFS) folding and NUFFT techniques.
• Parameters: Supports solve::Bool to trigger an iterative CG solver (conjugate gradient) for recovering the spectral coefficients from scattered grid measurements.

AbstractExecutionBackend

Abstract supertype for execution backends — the local compute backends (SerialBackend, ThreadedBackend, GPUBackend) that say what one process computes on, and the distribution wrappers (DistributedBackend, MPIBackend), parametric over an inner local backend, that say how work is split across processes. Orthogonal to AbstractSpectralBackend.

Serial single-threaded CPU execution (always available, no extension needed).

Multithreaded CPU execution over the direct-sum outer loop / internal library threads (requires using OhMyThreads).

GPUBackend{B}

GPU execution on the KernelAbstractions backend object B (e.g. GPUBackend(CUDABackend()), GPUBackend(KernelAbstractions.CPU())). Requires using KernelAbstractions plus a vendor GPU package. On a CUDA device, FFTBackend/NUFFTBackend route to CUFFT/cuFINUFFT (fast); on any KA device DirectSumBackend (and every spherical transform) runs the portable KA direct-sum kernels.

DistributedBackend(inner = SerialBackend())

Multi-process execution via Distributed, each worker running inner locally. Requires using Distributed and workers started with addprocs() + @everywhere using FlowFieldSpectra. Parametric over the inner local backend, e.g. DistributedBackend(ThreadedBackend()) for multithreaded workers or DistributedBackend(GPUBackend(dev)) for one GPU per worker.

MPIBackend(inner = SerialBackend(); comm = nothing)

Multi-rank execution via MPI.jl, each rank running inner locally and partial coefficient buffers combined in place with MPI.Allreduce! (every rank ends with the full result). Requires using MPI and launching under mpiexec with MPI.Init() called. Not CPU-only: MPIBackend(GPUBackend(dev)) targets a multi-GPU cluster (one GPU per rank) and MPIBackend(ThreadedBackend()) is hybrid MPI+threads. comm === nothing ⇒ the MPI extension uses MPI.COMM_WORLD (the core never references MPI).

AutoBackend <: AbstractExecutionBackend

Select the best available local execution backend at call time: ThreadedBackend when the OhMyThreads extension is loaded and Threads.nthreads() > 1, else SerialBackend. Never resolves to GPU/Distributed/MPI — those require an explicit device/process context.

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.