braintools.analysis.spectral.band_power

Contents

braintools.analysis.spectral.band_power#

braintools.analysis.spectral.band_power(x, dt, band, nperseg=None, relative=False)#

Integrate the Welch power spectral density over a frequency band.

The PSD is estimated with Welch’s method [1] exactly as scipy.signal.welch with its defaults (periodic Hann window, 50 % overlap, mean detrending, density scaling, mean average), and integrated with the trapezoid rule over the frequency bins low <= f <= high.

Parameters:

  • x (Array | ndarray | bool | number | bool | int | float | complex | Quantity) – Real signal of shape (T, ...); time is axis 0. A brainunit.Quantity is accepted.

  • dt (float | Quantity) – Sampling interval (seconds if a plain number). Static.

  • band (Tuple[float, float]) – (low, high) in Hz with 0 < low < high < fs / 2. Static.

  • nperseg (int | None) – Welch segment length. Static. Defaults to min(256, T) (as SciPy).

  • relative (bool) – If True, divide by the power integrated over all frequencies.

Returns:

power – Band power of shape x.shape[1:]. Absolute power has unit U**2 when x has unit U; relative power is dimensionless and nan for an all-constant signal (zero total power).

Return type:

Array | Quantity

Raises:

ValueError – If a band edge is <= 0 or >= fs / 2, T < nperseg, or fewer than two frequency bins fall inside band (increase nperseg).

See also

spectrogram

Short-time Fourier power spectrogram.

multitaper_psd

Multitaper power spectral density.

braintools.analysis.spectral.power_spectral_density

Welch PSD.

Notes

Equivalent to f, p = scipy.signal.welch(x, fs, nperseg=nperseg, axis=0) followed by scipy.integrate.trapezoid(p[m], f[m], axis=0) with m = (f >= low) & (f <= high). Jittable, vmappable and differentiable with respect to x.

References

Examples

>>> import jax.numpy as jnp
>>> from braintools.analysis.spectral import band_power
>>> t = jnp.arange(4000) * 1e-3
>>> x = jnp.sin(2 * jnp.pi * 10 * t)
>>> round(float(band_power(x, dt=1e-3, band=(5.0, 15.0), nperseg=500)), 2)
0.5
>>> round(float(band_power(x, dt=1e-3, band=(5.0, 15.0), nperseg=500, relative=True)), 2)
1.0
Contents