# plasma_plots.spectral

*module*

Fourier and spectral diagnostics of labeled arrays, computed on demand.

Ported from Struphy's ``postprocessing-fft`` branch (``fft``, ``time_fft``, ``inverse_time_fft``,
``fwhm_window``, ``filter_time``), with the same conventions, plus further tools: explicit band
filters, spectral peaks, spectrograms, poloidal/toroidal mode decomposition, mode structure at a
frequency, cross-spectra and matrix-pencil fits of complex frequencies.

Conventions: every transform divides by the sample count ``N`` (numpy's ``norm="forward"``).
Frequencies are angular, ω = 2π f, in units inverse to the coordinate. The forward kernel is
numpy's ``exp(-i ω t)``, so ``exp(+i ω t)`` appears at positive ω and a mode
``exp(2π i m eta2)`` at ``m``. A right-moving wave ``exp(i (k x - ω t))`` therefore sits at
``(k, -ω)``. The dominant-band filter follows the TAE_example_Shrut workflow, with xarray
coordinates replacing dictionaries of post-processed snapshots.

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/spectral.py#L1-L1)

## plasma_plots.spectral.TimeFilterResult

*class* · *dataclass*

```python
class TimeFilterResult
```

The result of [`filter_time()`][plasma_plots.spectral.filter_time]: filtered field and reduced spectrum with the band.

A zero/constant signal has no oscillatory peak: ``has_peak=False``, indices -1, frequencies
NaN, and a zero filtered signal.

**Attributes**

- `filtered` (`xarray.DataArray`) — The real signal reconstructed from the selected band, on the input's grid, with attrs ``time_filter``, ``omega_min`` and ``pad_bins`` added.
- `spectrum` (`xarray.Dataset`) — The selected band for each retained dimension: ``power`` (summed over the reduced dimensions, over ``omega``), ``dominant_frequency``, ``idx_dominant``, ``idx_lo``, ``idx_hi``, ``omega_lo``, ``omega_hi`` and ``has_peak``. Its attrs record ``omega_min``, ``pad_bins`` and ``power_reduction_dims``.

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/spectral.py#L399-L419)

### plasma_plots.spectral.TimeFilterResult.filtered

*attribute* · *instance attribute*

```python
filtered: xr.DataArray
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/spectral.py#L418-L418)

### plasma_plots.spectral.TimeFilterResult.spectrum

*attribute* · *instance attribute*

```python
spectrum: xr.Dataset
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/spectral.py#L419-L419)

## plasma_plots.spectral.band_filter

*function*

```python
def band_filter(data: xr.DataArray, omega_lo: float, omega_hi: float, *, detrend: bool = False) -> xr.DataArray
```

Reconstruct only the frequencies in ``[omega_lo, omega_hi]`` (inclusive), at every point.

The explicit counterpart of [`filter_time()`][plasma_plots.spectral.filter_time], e.g. to separate two known modes, or to
keep a gap frequency read off a continuum plot. A rectangular band: finite records can
show leakage and edge ringing.

**Parameters**

- `data` (`xarray.DataArray`) — A real signal with a uniformly spaced ``t`` dimension.
- `omega_lo` (`float`) — The lowest angular frequency kept.
- `omega_hi` (`float`) — The highest angular frequency kept, at least ``omega_lo``.
- `detrend` (`bool`) (default: `False`) — Subtract the temporal mean first; the result then has zero mean even if the band includes DC. Default: False.

**Returns**

- (`xarray.DataArray`) — The filtered signal on the grid of ``data``, with attrs ``time_filter="band"``, ``omega_lo`` and ``omega_hi`` added.

**Raises**

- `ValueError` — If ``omega_lo`` exceeds ``omega_hi``.

> **See Also**
>
> [`filter_time()`][plasma_plots.spectral.filter_time] : Keep the dominant peak's band automatically.

**Examples**

```pycon
>>> slow = band_filter(phi, 0.5, 1.5)
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/spectral.py#L594-L642)

## plasma_plots.spectral.cross_spectrum

*function*

```python
def cross_spectrum(first: xr.DataArray, second: xr.DataArray, *, dims=None, detrend: bool = True, window: str | None = None) -> xr.Dataset
```

Compute the cross-spectrum of two real signals on the same time grid, with phase and coherence.

``cross = conj(F1) F2`` of their one-sided [`time_fft()`][plasma_plots.spectral.time_fft] coefficients, so ``phase``
(radians) is how far ``second`` leads ``first`` at each frequency: +π/2 for
``second = -sin(ω t)`` against ``first = cos(ω t)``. With ``dims`` (e.g. every spatial point
as an ensemble) the cross-spectrum is summed over them first, and the phase coherence
``|Σ cross| / Σ |cross|`` measures, between 0 and 1, how consistently the two signals keep
one phase across the ensemble, whatever their amplitude profiles (1: the same phase
everywhere, near 0: random phases, as for noise). Without averaging it would be 1 by
construction, and is left out.

**Parameters**

- `first` (`xarray.DataArray`) — The reference signal, real, with a uniformly spaced ``t`` dimension.
- `second` (`xarray.DataArray`) — The other real signal, on exactly the same coordinates.
- `dims` (`str or sequence of str`) (default: `None`) — Dimensions to sum the cross-spectrum over, as an ensemble. Default: none (no ``coherence``).
- `detrend` (`bool`) (default: `True`) — Subtract each signal's temporal mean first. Default: True.
- `window` (`(None, 'hann')`) (default: `None`) — Window applied to both signals before transforming. Default: None.

**Returns**

- (`xarray.Dataset`) — Over ``omega`` and the dimensions not summed: ``cross`` (complex), ``magnitude``, ``phase`` (in rad) and, with ``dims``, ``coherence`` (NaN where both signals vanish).

**Raises**

- `ValueError` — If the coordinates of the two signals differ, or either is complex or non-uniform in
time.

> **See Also**
>
> [`plasma_plots.spectral_plots.plot_cross_spectrum()`][plasma_plots.spectral_plots.plot_cross_spectrum] : Draw magnitude, coherence and phase.

**Examples**

```pycon
>>> cross = cross_spectrum(phi, density, dims=["eta1", "eta2", "eta3"])
>>> cross.phase.sel(omega=1.0, method="nearest")
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/spectral.py#L1189-L1263)

## plasma_plots.spectral.drop_periodic_endpoint

*function*

```python
def drop_periodic_endpoint(data: xr.DataArray, dim: str, *, period: float = 1.0) -> xr.DataArray
```

Drop the last sample along ``dim`` if it repeats the first one period later.

Struphy's logical grids often include both ends of a periodic direction (``eta2 = 0`` and
``eta2 = 1``); a Fourier transform must see each point once. Arrays without a duplicate come
back unchanged.

**Parameters**

- `data` (`xarray.DataArray`) — The array.
- `dim` (`str`) — The periodic dimension.
- `period` (`float`) (default: `1.0`) — The period in the coordinate of ``dim``. Default: 1.0.

**Returns**

- (`xarray.DataArray`) — ``data`` without its last sample along ``dim`` if the last coordinate equals the first plus ``period``, otherwise ``data`` itself.

**Raises**

- `ValueError` — If ``dim`` is not a dimension of ``data``.

**Examples**

```pycon
>>> fft(
...     drop_periodic_endpoint(phi.isel(t=-1, eta1=0, eta3=0), "eta2"),
...     dim="eta2",
... )
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/spectral.py#L546-L591)

## plasma_plots.spectral.fft

*function*

```python
def fft(data: xr.DataArray, *, dim: str, detrend: bool = False, window: str | None = None) -> xr.DataArray
```

Compute the two-sided, shifted FFT along a named uniform coordinate, normalized by N.

Frequencies are angular (2π times cycles per unit of the supplied coordinate) and run from
negative to positive (``fftshift``). With numpy's kernel ``exp(-i ω t)``, a right-moving wave
``exp(i (k x - ω t))`` sits at ``(k, -ω)``. Remove duplicate endpoints of periodic spatial
grids before calling this function (see [`drop_periodic_endpoint()`][plasma_plots.spectral.drop_periodic_endpoint]); no endpoint is
dropped automatically. Windowed powers refer to the windowed signal, without
amplitude/energy compensation. Coefficient phases are relative to the first sample,
recorded as ``sample_origin``.

**Parameters**

- `data` (`xarray.DataArray`) — The signal, real or complex, with finite values. ``dim`` must be a dimension with a strictly increasing, uniformly spaced one-dimensional coordinate.
- `dim` (`str`) — The dimension to transform.
- `detrend` (`bool`) (default: `False`) — Subtract the mean along ``dim`` first. Default: False.
- `window` (`(None, 'hann')`) (default: `None`) — Multiply by a periodic Hann window (see [`hann()`][plasma_plots.spectral.hann]) first. Default: None (boxcar).

**Returns**

- (`xarray.DataArray`) — Complex coefficients named ``coefficients``, ``fft(data) / N``, with ``dim`` replaced by ``omega`` (for ``dim="t"``) or ``k_<dim>`` (otherwise). Coordinates that depend on ``dim`` are dropped. The input's attrs are kept and extended by ``transform_dim``, ``n_samples``, ``sample_spacing``, ``sample_origin``, ``frequency_resolution`` (2π / (N dt)), ``nyquist_frequency`` (π / dt), ``normalization`` (``"forward"``), ``window``, ``detrend`` and ``label``.

**Raises**

- `TypeError` — If ``data`` is not a DataArray.
- `ValueError` — If ``dim`` has no uniform, increasing numeric coordinate with at least two samples, the
values are not finite, the frequency coordinate already exists, or ``detrend``/``window``
are invalid.

> **See Also**
>
> [`time_fft()`][plasma_plots.spectral.time_fft] : The one-sided transform of a real signal in time, with power.
> [`drop_periodic_endpoint()`][plasma_plots.spectral.drop_periodic_endpoint] : Remove a duplicate periodic endpoint first.

**Examples**

```pycon
>>> coefficients = fft(phi.isel(t=-1, eta2=0, eta3=0), dim="eta1")
>>> power = abs(coefficients) ** 2
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/spectral.py#L132-L201)

## plasma_plots.spectral.filter_time

*function*

```python
def filter_time(data: xr.DataArray, *, dims=None, omega_min: float = 1e-08, pad_bins: int = 0) -> TimeFilterResult
```

Keep the dominant peak's FWHM frequency band and reconstruct the real signal.

The [`time_fft()`][plasma_plots.spectral.time_fft] power is summed over ``dims`` to select a shared band. Each remaining
coordinate (e.g. each component) gets its own band, which is applied at every spatial
point. The sum is unweighted, not a physical energy integral. The band is the half-power
band of the strongest bin at ``omega >= omega_min`` (see [`fwhm_window()`][plasma_plots.spectral.fwhm_window]). Frequencies
below ``omega_min`` (including DC) are always removed, even with padding. This
rectangular-bin filter uses no taper; finite records can exhibit spectral leakage and edge
ringing. Source data is not mutated.

**Parameters**

- `data` (`xarray.DataArray`) — A real signal with a uniformly spaced ``t`` dimension.
- `dims` (`str or sequence of str`) (default: `None`) — Non-time dimensions to sum the power over before choosing the band. Default: all dimensions except ``t`` and ``component``. Pass ``dims=()`` for independent filtering at each point.
- `omega_min` (`float`) (default: `1e-08`) — Finite, positive lowest frequency considered, which excludes DC. Default: 1e-8.
- `pad_bins` (`int`) (default: `0`) — Nonnegative number of extra bins on each side of the band. Default: 0.

**Returns**

- (`TimeFilterResult`) — The ``filtered`` signal and the ``spectrum`` with the selected band.

**Raises**

- `ValueError` — If ``omega_min`` is not finite and positive or exceeds every bin, ``pad_bins`` is not a
nonnegative integer, or ``dims`` does not name distinct non-time dimensions.

> **See Also**
>
> [`band_filter()`][plasma_plots.spectral.band_filter] : Keep an explicit frequency band instead.
> [`plasma_plots.spectral_plots.plot_filtered()`][plasma_plots.spectral_plots.plot_filtered] : Compare the signal with its reconstruction.

**Examples**

```pycon
>>> result = filter_time(phi, pad_bins=1)
>>> result.spectrum.dominant_frequency.item()
>>> result.filtered.isel(t=-1)
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/spectral.py#L422-L538)

## plasma_plots.spectral.fwhm_window

*function*

```python
def fwhm_window(power, idx_peak: int, idx_min: int = 0, pad_bins: int = 0)
```

Find the inclusive, contiguous half-power band about a peak, padded and clamped.

Starting at ``idx_peak``, the band grows to each side while the power stays at or above half
the peak power, never below ``idx_min``; it is then widened by ``pad_bins`` on each side and
clamped to ``[idx_min, len(power) - 1]``.

**Parameters**

- `power` (`array_like`) — A finite, nonnegative one-dimensional power spectrum.
- `idx_peak` (`int`) — The bin of the peak; its power must be positive.
- `idx_min` (`int`) (default: `0`) — The lowest bin the band may include, at most ``idx_peak``. Default: 0.
- `pad_bins` (`int`) (default: `0`) — Nonnegative number of extra bins on each side. Default: 0.

**Returns**

- (`(int, int)`) — The first and last bin of the band, both included.

**Raises**

- `ValueError` — If ``power`` is not a finite, nonnegative 1-D array, the indices are not integers, or
the peak, minimum bin or padding are invalid.

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/spectral.py#L347-L396)

## plasma_plots.spectral.hann

*function*

```python
def hann(n: int) -> np.ndarray
```

Return the periodic Hann window of length ``n``.

The same as ``scipy.signal.windows.hann(n, sym=False)``: ``0.5 - 0.5 cos(2π j / n)`` for
``j = 0, ..., n - 1``.

**Parameters**

- `n` (`int`) — The number of samples.

**Returns**

- (`numpy.ndarray`) — The window, of length ``n``, starting at 0.

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/spectral.py#L63-L79)

## plasma_plots.spectral.inverse_time_fft

*function*

```python
def inverse_time_fft(coefficients: xr.DataArray, template: xr.DataArray) -> xr.DataArray
```

Invert forward-normalized rFFT coefficients, using a template's length and coordinates.

The original length is required to distinguish odd and even sample counts. For
windowed/detrended coefficients this reconstructs the processed signal; it does not undo the
window or add the mean back.

**Parameters**

- `coefficients` (`xarray.DataArray`) — One-sided coefficients as from [`time_fft()`][plasma_plots.spectral.time_fft] (``rfft / N``), over ``omega`` and the template's other dimensions, in any order. Possibly modified, e.g. with bins set to zero.
- `template` (`xarray.DataArray`) — A real signal on the original time grid, which gives the sample count, the times, the other coordinates and the attrs of the result.

**Returns**

- (`xarray.DataArray`) — A copy of ``template`` with the reconstructed real values.

**Raises**

- `ValueError` — If the dimensions, the frequency grid or another coordinate do not match the template,
or the coefficients are not finite.

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/spectral.py#L285-L344)

## plasma_plots.spectral.matrix_pencil

*function*

```python
def matrix_pencil(data: xr.DataArray, *, n_modes: int = 1, pencil: int | None = None, detrend: bool = False) -> xr.Dataset
```

Fit frequencies and growth rates of a sum of exponentially growing or damped oscillations.

Fits ``f(t) = Σ_j a_j exp((γ_j + i ω_j)(t - t0))`` with the matrix-pencil method (Hua &
Sarkar, 1990): an SVD of the Hankel matrix of the samples, whose signal subspace shifts by
one sample as multiplication by ``exp((γ + i ω) dt)``. Unlike an FFT peak, this is not
limited to the bin spacing 2π/T: a clean record shorter than one period can still give the
frequency, together with the growth (γ > 0) or damping rate.

For a real signal, ``n_modes`` counts real oscillations (each a conjugate pair), and one
extra real exponential is fitted to absorb an offset or slow trend; oscillations are
returned first, then a non-oscillating component if there are fewer than ``n_modes``
oscillations.

**Parameters**

- `data` (`xarray.DataArray`) — A ``(t,)`` series, real or complex, uniformly spaced in time.
- `n_modes` (`int`) (default: `1`) — The number of modes to fit and return: real oscillations for a real signal, complex exponentials for a complex one. The fit needs well over ``2 * n_modes`` samples. Default: 1.
- `pencil` (`int`) (default: `None`) — The pencil parameter (Hankel matrix width minus one), which trades noise robustness against resolution. Default: ``N // 2``.
- `detrend` (`bool`) (default: `False`) — Subtract the mean first. Default: False.

**Returns**

- (`xarray.Dataset`) — Along ``mode``, strongest first (oscillations first for real input): ``omega`` (≥ 0 for real input), ``gamma``, the real ``amplitude`` (``2|a|`` for a conjugate pair) and ``phase`` (rad) at ``t0``. The attrs hold ``label``, ``sample_origin`` (``t0``, the first time), the relative rms ``residual`` of the reconstruction, ``pencil`` and the leading ``singular_values``.

**Raises**

- `ValueError` — If ``data`` is not a ``(t,)`` series, or has too few samples for ``n_modes`` with this
``pencil``.

> **See Also**
>
> [`pencil_reconstruction()`][plasma_plots.spectral.pencil_reconstruction] : The fitted signal at any times.
> [`plasma_plots.spectral_plots.plot_pencil_fit()`][plasma_plots.spectral_plots.plot_pencil_fit] : Draw the fit and the complex frequencies.

**Examples**

```pycon
>>> fit = matrix_pencil(energy.sel(t=slice(0.0, 5.0)), n_modes=2)
>>> fit.omega.values, fit.gamma.values
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/spectral.py#L1266-L1393)

## plasma_plots.spectral.mode_amplitudes

*function*

```python
def mode_amplitudes(modes: xr.DataArray, *, top: int | None = None, real: bool = True, relative: bool = False) -> xr.DataArray
```

Stack the mode amplitudes from [`mode_spectrum()`][plasma_plots.spectral.mode_spectrum] along one labeled ``mode`` dimension.

**Parameters**

- `modes` (`xarray.DataArray`) — The output of [`mode_spectrum()`][plasma_plots.spectral.mode_spectrum] (complex, with ``mode_names`` in its attrs, or with ``m``/``n`` dimensions).
- `top` (`int`) (default: `None`) — Keep only the ``top`` modes with the largest peak amplitude over every other dimension, strongest first. Default: all modes.
- `real` (`bool`) (default: `True`) — The field is real: each ``(m, n)`` is combined with its conjugate ``(-m, -n)``. Only the half with the first nonzero mode number positive is kept, and its amplitude doubled, so a field ``A cos(...)`` gives ``A``. A mode without a twin on the grid (the mean, or the Nyquist mode of an even grid) is kept as it is. Default: True.
- `relative` (`bool`) (default: `False`) — Divide by the amplitude of the mean (the mode with all numbers zero), which is then left out: e.g. density perturbations relative to the background density, as growth plots of an instability often show. NaN where the mean vanishes. Default: False.

**Returns**

- (`xarray.DataArray`) — ``amplitude``: the real amplitudes over ``mode`` and the remaining dimensions (e.g. ``t``). Coordinates on ``mode`` give each mode's numbers (``m``, ``n``) and a label such as ``"(10, -1)"``. The attrs hold ``label`` and ``mode_names``.

**Raises**

- `ValueError` — If ``modes`` has no mode numbers, or ``relative=True`` and there is no mean mode.

> **See Also**
>
> [`plasma_plots.spectral_plots.plot_mode_amplitudes()`][plasma_plots.spectral_plots.plot_mode_amplitudes] : Amplitudes over time, with growth rates.

**Examples**

```pycon
>>> amplitudes = mode_amplitudes(mode_spectrum(phi.isel(eta1=8)), top=4)
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/spectral.py#L1006-L1099)

## plasma_plots.spectral.mode_spectrum

*function*

```python
def mode_spectrum(data: xr.DataArray, *, dims=None, names=('m', 'n'), periods=None, scale=None) -> xr.DataArray
```

Compute complex Fourier amplitudes over integer mode numbers along periodic directions.

For a torus with θ = 2π eta2 and φ = 2π eta3, the default gives coefficients over poloidal
``m`` and toroidal ``n``, as functions of every remaining dimension, e.g.
``(t, eta1, m, n)``. On GVEC's angles (see [`plasma_plots.gvec.from_gvec()`][plasma_plots.gvec.from_gvec]) the default
uses their periods, 2π and 2π/nfp, and gives the full-torus ``n``, a multiple of nfp. A duplicate periodic endpoint is dropped first (see
[`drop_periodic_endpoint()`][plasma_plots.spectral.drop_periodic_endpoint]), and each transform is normalized by N (see [`fft()`][plasma_plots.spectral.fft]).
The mode ``exp(2π i (m eta2 + n eta3))`` appears at ``(m, n)``, so a real field ``cos(...)``
of amplitude ``A`` has ``A/2`` at ``(m, n)`` and at ``(-m, -n)``; see
[`mode_amplitudes()`][plasma_plots.spectral.mode_amplitudes]. A sector of the torus (``tor_period`` in Struphy) counts ``n`` per
sector; multiply by the number of sectors (or use ``scale``) for the full-torus mode number.

**Parameters**

- `data` (`xarray.DataArray`) — The field, sampled uniformly over one full period along each of ``dims``.
- `dims` (`str or sequence of str`) (default: `None`) — The periodic dimensions to transform. Default: the two angles of the logical dimensions (see [`plasma_plots.arrays.logical_dims()`][plasma_plots.arrays.logical_dims]), ``("eta2", "eta3")`` for Struphy.
- `names` (`str or sequence of str`) (default: `('m', 'n')`) — The name of the mode number of each dimension, one per dimension. Default: ``("m", "n")``.
- `periods` (`float or sequence of float`) (default: `None`) — Each direction's period in its coordinate: one number for all, or one per dimension. Default: the coordinate's ``period`` attribute (see [`plasma_plots.arrays.angle_period()`][plasma_plots.arrays.angle_period]), else 1.0.
- `scale` (`int or sequence of int`) (default: `None`) — Multiplies the mode numbers (one number, or one per dimension; cast to integers), e.g. ``scale=(1, 6)`` labels a sixth of a torus (Struphy's ``tor_period=6``) with full-torus toroidal mode numbers. Default: 2π over the ``period`` attribute of an angle (nfp for GVEC's toroidal angle), else 1.

**Returns**

- (`xarray.DataArray`) — ``modes``: complex amplitudes with each of ``dims`` replaced by its integer mode number (``names``). The attrs hold ``label``, ``mode_dims``, ``mode_names`` and the run's provenance.

**Raises**

- `ValueError` — If ``dims``, ``names``, ``periods`` and ``scale`` differ in length, or a dimension does
not sample one full period uniformly.

> **See Also**
>
> [`mode_amplitudes()`][plasma_plots.spectral.mode_amplitudes] : Real amplitudes of the modes along one ``mode`` dimension.
> [`mode_structure()`][plasma_plots.spectral.mode_structure] : Transform the complex amplitude at one frequency for its harmonics.

**Examples**

```pycon
>>> modes = mode_spectrum(phi)
>>> abs(modes.sel(m=2, n=-1)).isel(t=-1).plot()
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/spectral.py#L894-L1003)

## plasma_plots.spectral.mode_structure

*function*

```python
def mode_structure(data: xr.DataArray, omega: float, *, window: str | None = 'hann', detrend: bool = True) -> xr.DataArray
```

Compute the complex amplitude of the oscillation at an exact frequency ``omega``, everywhere.

``A(x) = 2 Σ_t w(t) f(t, x) exp(-i ω (t - t0)) / Σ_t w(t)``, so a field
``a(x) cos(ω t + φ(x))`` gives ``a exp(i φ)``: ``abs()`` is the eigenfunction's amplitude
and ``np.angle()`` its phase, e.g. at a frequency from [`spectral_peaks()`][plasma_plots.spectral.spectral_peaks]
(``omega_refined``), not limited to FFT bins. The Hann window suppresses leakage from other
frequencies; the record should still span a few periods of their difference. For the radial
profile of each poloidal harmonic, transform the result over the angles:
``mode_spectrum(mode_structure(field, omega))``. With the numpy sign convention, a wave
``cos(m θ + n φ - ω t)`` then appears at ``(-m, -n)``.

**Parameters**

- `data` (`xarray.DataArray`) — A real field with a uniformly spaced ``t`` dimension.
- `omega` (`float`) — The angular frequency ω.
- `window` (`('hann', None)`) (default: `"hann"`) — The weights ``w(t)``: a periodic Hann window, or ``None`` for uniform weights. Default: ``"hann"``.
- `detrend` (`bool`) (default: `True`) — Subtract the temporal mean at every point first. Default: True.

**Returns**

- (`xarray.DataArray`) — ``mode_structure``: the complex amplitude over every dimension of ``data`` but ``t``. The attrs hold ``label``, ``omega`` and ``sample_origin`` (``t0``, the first time, which the phases refer to).

**Raises**

- `ValueError` — If ``data`` is complex (apply this before [`mode_spectrum()`][plasma_plots.spectral.mode_spectrum]), ``t`` is not uniform,
or ``window`` is invalid.

> **See Also**
>
> [`plasma_plots.spectral_plots.plot_mode_profiles()`][plasma_plots.spectral_plots.plot_mode_profiles] : Draw the harmonics' radial profiles.

**Examples**

```pycon
>>> omega = float(
...     spectral_peaks(phi.isel(eta1=8, eta2=0, eta3=0)).omega_refined[0]
... )
>>> harmonics = mode_spectrum(mode_structure(phi, omega))
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/spectral.py#L1102-L1186)

## plasma_plots.spectral.pencil_reconstruction

*function*

```python
def pencil_reconstruction(fit: xr.Dataset, t) -> xr.DataArray
```

Evaluate the real signal described by a [`matrix_pencil()`][plasma_plots.spectral.matrix_pencil] fit of a real series.

The sum over modes of ``amplitude exp(gamma (t - t0)) cos(omega (t - t0) + phase)``.

**Parameters**

- `fit` (`xarray.Dataset`) — The result of [`matrix_pencil()`][plasma_plots.spectral.matrix_pencil] for a real series.
- `t` (`array_like`) — The times to evaluate at.

**Returns**

- (`xarray.DataArray`) — ``fit`` over ``t``.

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/spectral.py#L1396-L1429)

## plasma_plots.spectral.spectral_peaks

*function*

```python
def spectral_peaks(data, *, n_peaks: int = 3, dims=None, omega_min: float = 1e-08, rel_height: float = 0.001, detrend: bool = True, window: str | None = None) -> xr.Dataset
```

Find the ``n_peaks`` strongest local maxima of a power spectrum, with sub-bin frequencies.

Only maxima at ``omega >= omega_min`` and above ``rel_height`` times the largest of them
count; a rising last bin counts as a peak at the Nyquist edge. Each peak's frequency is
refined below the bin spacing by a parabolic fit to the log power of the peak and its two
neighbors, which locates an off-bin frequency to a fraction of a bin (the shift is clipped
to half a bin).

**Parameters**

- `data` (`xarray.DataArray or xarray.Dataset`) — A time series with a ``t`` dimension (transformed by [`time_fft()`][plasma_plots.spectral.time_fft] with ``detrend`` and ``window``), a [`time_fft()`][plasma_plots.spectral.time_fft] Dataset (its ``power`` is used), or an array over ``omega``: power, or complex coefficients whose squared magnitude is used.
- `n_peaks` (`int`) (default: `3`) — The number of peaks to return at most. Default: 3.
- `dims` (`str or sequence of str`) (default: `None`) — Dimensions to sum the power over. Default: every dimension but ``omega``. Only ``omega`` may remain.
- `omega_min` (`float`) (default: `1e-08`) — The lowest frequency a peak may have, which excludes DC. Default: 1e-8.
- `rel_height` (`float`) (default: `0.001`) — The weakest peak kept, relative to the strongest. Default: 1e-3.
- `detrend` (`bool or int`) (default: `True`) — For a time series: ``True`` subtracts the mean, ``False`` nothing, and an integer is the degree of a least-squares polynomial in ``t`` removed at every point first. An energy such as LinearMHD's ``en_U`` oscillates at twice the wave frequency around a slow trend, so its peaks with ``detrend=2`` sit at ``2 * omega``. Ignored for spectra. Default: True.
- `window` (`(None, 'hann')`) (default: `None`) — Window applied before transforming a time series. Default: None.

**Returns**

- (`xarray.Dataset`) — Along ``peak``, strongest first: ``omega`` (the peak bin), ``omega_refined`` (the sub-bin frequency), ``power``, and the half-power band ``omega_lo``/``omega_hi`` (see [`fwhm_window()`][plasma_plots.spectral.fwhm_window]). The attrs hold the ``frequency_resolution`` (bin spacing).

**Raises**

- `ValueError` — If dimensions other than ``omega`` remain after the sum over ``dims``.

> **See Also**
>
> [`matrix_pencil()`][plasma_plots.spectral.matrix_pencil] : Frequencies and growth rates beyond the bin resolution.
> [`mode_structure()`][plasma_plots.spectral.mode_structure] : The eigenfunction at a peak's frequency.

**Examples**

```pycon
>>> peaks = spectral_peaks(phi.isel(eta2=0, eta3=0), n_peaks=2)
>>> peaks.omega_refined.values
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/spectral.py#L687-L799)

## plasma_plots.spectral.spectrogram

*function*

```python
def spectrogram(data: xr.DataArray, *, length: int | float, step: int | float | None = None, detrend: bool = True, window: str | None = 'hann') -> xr.DataArray
```

Compute short-time power spectra: [`time_fft()`][plasma_plots.spectral.time_fft] power in sliding windows along ``t``.

For following a frequency that drifts (a chirping mode), or telling a persistent
oscillation from an initial transient. The frequency resolution is 2π / ``length``.

**Parameters**

- `data` (`xarray.DataArray`) — A real signal with a uniformly spaced ``t`` dimension.
- `length` (`int or float`) — The window length: a sample count (integer, 4 to the number of samples) or a time span (float).
- `step` (`int or float`) (default: `None`) — The shift between windows: a sample count (integer) or a time span (float). Default: a quarter of ``length``.
- `detrend` (`bool`) (default: `True`) — Subtract each window's mean. Default: True.
- `window` (`('hann', None)`) (default: `"hann"`) — The taper of each window, not compensated for. Default: ``"hann"``.

**Returns**

- (`xarray.DataArray`) — ``power`` over ``(t, omega, ...)``, where ``t`` is each window's center and ``...`` the other dimensions of ``data``. The attrs hold ``label``, ``window_length``, ``window_step`` (both in time units) and ``frequency_resolution``.

**Raises**

- `ValueError` — If ``length`` spans fewer than 4 or more than all samples, or ``step`` less than one.

> **See Also**
>
> [`plasma_plots.spectral_plots.plot_spectrogram()`][plasma_plots.spectral_plots.plot_spectrogram] : Draw the result.

**Examples**

```pycon
>>> power = spectrogram(phi.isel(eta1=8, eta2=0, eta3=0), length=10.0)
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/spectral.py#L802-L891)

## plasma_plots.spectral.time_fft

*function*

```python
def time_fft(data: xr.DataArray, *, detrend: bool = False, window: str | None = None) -> xr.Dataset
```

Compute the one-sided time FFT, with complex coefficients and mean-square power per bin.

Replaces ``t`` by ``omega`` (ω ≥ 0) and preserves other dimensions and coordinates.
``coefficients = rfft(data) / N``. ``power`` doubles the positive-frequency bins, except the
even-N Nyquist bin: the DC and Nyquist amplitudes are never doubled. Its sum over ω equals
the temporal mean square of the input (after any mean subtraction/windowing), not a PSD per
unit frequency. A Hann window is not compensated for.

Sampling uses the saved times, including their units, not the simulation dt. Bin spacing is
2π / (N dt); padding is not used to claim extra resolution. Coordinates depending on ``t``
are dropped; time-independent mapped coordinates and provenance are retained. Computation
eagerly loads the selected array. Phases are relative to the first saved time, recorded as
``sample_origin``.

**Parameters**

- `data` (`xarray.DataArray`) — A real signal with a ``t`` dimension, uniformly spaced in time (select a uniform interval first if the saved times are not).
- `detrend` (`bool`) (default: `False`) — Subtract the temporal mean first. Default: False.
- `window` (`(None, 'hann')`) (default: `None`) — Multiply by a periodic Hann window (see [`hann()`][plasma_plots.spectral.hann]) first. Default: None (boxcar).

**Returns**

- (`xarray.Dataset`) — ``coefficients`` (complex) and ``power`` (real, with ``units`` the square of the input's), both over ``omega`` and the other dimensions of ``data``. The attrs are those of the coefficients (see [`fft()`][plasma_plots.spectral.fft]): ``n_samples``, ``sample_spacing``, ``sample_origin``, ``frequency_resolution``, ``nyquist_frequency``, ``window``, ...

**Raises**

- `ValueError` — If the signal is complex (use [`fft()`][plasma_plots.spectral.fft]), not finite, or ``t`` is not uniformly
spaced.

> **See Also**
>
> [`fft()`][plasma_plots.spectral.fft] : The two-sided transform along any dimension.
> [`inverse_time_fft()`][plasma_plots.spectral.inverse_time_fft] : Back from the coefficients to a signal.
> [`spectral_peaks()`][plasma_plots.spectral.spectral_peaks] : The strongest peaks of the power.

**Examples**

```pycon
>>> spectrum = time_fft(phi.isel(eta2=0, eta3=0), detrend=True)
>>> spectrum.power.sum("eta1").plot()
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/spectral.py#L204-L282)

## plasma_plots.spectral.trace_branch

*function*

```python
def trace_branch(spectrum: xr.DataArray, theory, *, window: float = 0.2, k_range: tuple[float, float] | None = None, threshold: float = 0.001) -> xr.Dataset
```

Measure the frequency of a dispersion branch near a theory curve, at every ``k``.

At each non-negative ``k`` (in ``k_range``) the power of waves travelling either way is
searched for its maximum within ``omega_theory (1 ± window)`` at positive ω, and the peak
refined below the bin spacing with a parabola through its log power. With numpy's sign
convention a right-moving wave sits at ``(k, -ω)``, i.e. mirrored at ``(-k, +ω)``, so the
power at ``-k`` is added to the power at ``+k``. Unlike
[`fit_dispersion_branches()`][plasma_plots.analysis.fit_dispersion_branches], the branch may be curved (e.g. a
whistler or Bohm-Gross branch).

**Parameters**

- `spectrum` (`xarray.DataArray`) — An ``(omega, k)`` power spectrum, e.g. from ``array.plasma.analysis.dispersion()``.
- `theory` (`callable`) — The expected branch ``omega(k)``, applied to an array of ``k``; of a complex frequency (as [`plasma_plots.theory`][plasma_plots.theory] returns), the real part is used.
- `window` (`float`) (default: `0.2`) — The relative half-width of the search window about the theory. Default: 0.2.
- `k_range` (`(float, float)`) (default: `None`) — The ``k`` interval to trace (negative ``k`` are never traced). Default: every ``k >= 0``.
- `threshold` (`float`) (default: `0.001`) — Maxima weaker than this times the strongest one found count as no wave. Default: 1e-3.

**Returns**

- (`xarray.Dataset`) — Over ``k``: ``omega`` (measured), ``omega_theory`` and ``relative_error``. A ``k`` is NaN where the window holds no local maximum (only the flank of a peak outside it), or where that maximum is weaker than ``threshold`` times the strongest one found (no wave at that ``k``). The attrs hold ``window`` and ``frequency_resolution``.

**Raises**

- `ValueError` — If ``spectrum`` lacks an ``omega`` or ``k`` dimension.

**Examples**

```pycon
>>> branch = trace_branch(
...     spectrum, lambda k: np.sqrt(1 + 3 * k**2), k_range=(0.0, 2.0)
... )
>>> branch.relative_error.plot()
```

[View source](https://github.com/max-models/plasma-plots/blob/devel/src/plasma_plots/spectral.py#L1432-L1539)