braintools.analysis.spectral.spectrogram

Contents

braintools.analysis.spectral.spectrogram#

braintools.analysis.spectral.spectrogram(x, dt, nperseg, noverlap=None, window='hann', detrend='constant')#

Compute the short-time Fourier power spectrogram along time axis 0.

The signal is cut into segments of nperseg samples that overlap by noverlap samples. Each segment is mean-detrended, multiplied by the window and Fourier transformed. The one-sided power spectral density of every segment is returned, scaled as in scipy.signal.spectrogram with scaling='density' and mode='psd'.

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.

  • nperseg (int) – Segment length in samples. Static. Must not exceed T.

  • noverlap (int | None) – Overlap between segments in samples, 0 <= noverlap < nperseg. Static. Defaults to nperseg // 8 (the SciPy default).

  • window (str | tuple | Array | ndarray | bool | number | bool | int | float | complex | Quantity) – Window specification passed to scipy.signal.get_window (periodic windows, as in SciPy), or an explicit window of length nperseg.

  • detrend (str | bool | None) – 'constant' removes the mean of every segment (the SciPy default); False/None disables detrending.

Return type:

Tuple[Array, Array | Quantity, Array | Quantity]

Returns:

  • freqs (jax.Array) – Frequencies in Hz, shape (F,) with F = nperseg // 2 + 1.

  • times (jax.Array or brainunit.Quantity) – Segment centres, shape (S,) with S = (T - noverlap) // (nperseg - noverlap). Seconds as a plain array, or a Quantity in seconds if dt is a Quantity.

  • Sxx (jax.Array or brainunit.Quantity) – Power spectral density, shape (F, S, ...). If x is a Quantity with unit U, Sxx has unit U**2 * second (i.e. U**2/Hz).

Raises:

ValueError – If T < nperseg, noverlap is out of range, detrend is unknown or dt is not positive.

See also

band_power

Welch PSD integrated over a frequency band.

multitaper_psd

Multitaper power spectral density.

morlet_transform

Complex Morlet wavelet time-frequency transform.

Notes

The result equals scipy.signal.spectrogram(x, fs, window, nperseg, noverlap, axis=0) (shape (F, ..., S)) with the segment axis moved to position 1. Unlike SciPy, a signal shorter than nperseg raises instead of silently shrinking nperseg.

The function is jittable and vmappable (all shape-determining arguments are static) and differentiable with respect to x.

References

Examples

The spectrogram of a 50 Hz sinusoid peaks at 50 Hz in every segment (periodogram scaling as in Welch [1]):

>>> import brainunit as u
>>> import jax.numpy as jnp
>>> from braintools.analysis.spectral import spectrogram
>>> t = jnp.arange(1000) * 1e-3
>>> x = jnp.sin(2 * jnp.pi * 50 * t)
>>> freqs, times, Sxx = spectrogram(x, dt=1 * u.ms, nperseg=200, noverlap=100)
>>> Sxx.shape
(101, 9)
>>> float(freqs[jnp.argmax(Sxx[:, 0])])
50.0
>>> print(times[:2])
[0.1 0.2] s
Contents