pyhctsa.operations.spectral.rhythmicity

pyhctsa.operations.spectral.rhythmicity(y, m=5, tau=1.5, n_freq=40, f_low=0.02, f_high=0.4, smooth_win=5, min_cycles=40)

Phase-persistence (“rhythmicity”) spectrum of a time series.

A power spectrum cannot tell a sustained oscillation apart from a train of brief, high-power bursts: both deposit energy at the same frequency. This operation instead asks, at each frequency, how predictable the phase is a fixed number of cycles into the future: high for an oscillation that keeps ticking (sustained, “rhythmic”), low for one that only flickers on intermittently (transient, “bursty”).

The method is the lagged-angle vector index (LAVI) of Karvat et al. (2026), a phase-only refinement of the lagged coherence of Fransen et al. (2015):

  1. A complex Morlet wavelet transform at each analysis frequency f (m cycles per wavelet) gives the instantaneous phase phi(t, f). The transform is an FFT-domain convolution.

  2. LAVI(f) = |mean_t exp(i [phi(t, f) - phi(t + L, f)])|, with the lag L = round(tau / f) samples, a fixed tau cycles. This is the mean resultant length of the lagged phase differences, in [0, 1]; it is 1 when the phase advances by a perfectly repeatable amount every tau cycles (a stable rhythm). A coefficient with |x(t, f)| < 1e-6 max_t |x(t, f)| (the maximum over the edge-free samples at that frequency; under 1e-12 of the peak power) has no meaningful phase, typically where a spiky or skewed series is (nearly) constant across the wavelet’s support, down to coefficients that are pure round-off. Any pair (t, t + L) involving one is left out of the mean, so that the output does not depend on round-off.

  3. The spectral null. Part of LAVI is set by the power spectrum alone: a Gaussian process with the same spectrum (equivalently, the mean over phase-randomized surrogates) has a lag-L correlation rho(f) of its wavelet output, computed here exactly from the series’ own periodogram, and its expected LAVI is then the mean resultant length of the phase difference of a circular complex Gaussian pair, LAVI0(f) = (pi/4) rho 2F1(1/2, 1/2; 2; rho**2) = (E(rho**2) - (1 - rho**2) K(rho**2)) / rho, with K and E the complete elliptic integrals. This includes a floor that the method has even for white noise: the wavelet overlaps itself across the lag, so rho = exp(-pi**2 tau**2 / m**2) (0.41 for the defaults) and LAVI0 = 0.33. The excess, LAVI(f) - LAVI0(f), is the part of the rhythmicity spectrum that the power spectrum does not determine: about zero for any linear Gaussian process, positive for a sustained oscillation, negative for bursts. This replaces the surrogate “noise ribbon” of Karvat et al., deterministically and at the cost of one extra FFT per frequency.

The amplitude-weighted lagged coherence is deliberately not reported: up to edge effects it is the normalized lag-L autocorrelation of the wavelet output, a fixed function of the periodogram (exactly rho(f) above).

The series is treated as evenly sampled with unit spacing: frequencies are in cycles per sample (0 to 0.5) and reported angular frequencies are 2*pi*f, as in spectral_summaries(). LAVI is biased upward when few cycles are available, so a frequency is used only when the pairs used (free of edge effects, with a defined phase) span at least min_cycles cycles, i.e. number of pairs times f is at least min_cycles; the output is NaN if fewer than 10 frequencies qualify (for the defaults, N below about 240).

Parameters:
y : array-like

The input time series.

m : float, optional

Width of the Morlet wavelet, in cycles. Larger values give sharper frequency resolution and blunter time resolution. Default is 5, as in Karvat et al.

tau : float, optional

The lag at which phase persistence is measured, in cycles. Default is 1.5, as in Karvat et al., who report that the relative shape of the rhythmicity spectrum is nearly independent of this choice.

n_freq : int, optional

Number of log-spaced analysis frequencies. Default is 40.

f_low : float, optional

Lowest analysis frequency, in cycles per sample. Default is 0.02.

f_high : float, optional

Highest analysis frequency, in cycles per sample; must be below the Nyquist frequency of 0.5. Default is 0.4.

smooth_win : int, optional

Width, in frequency bins (forced odd), of the moving-average (Daniell) spectral window applied before counting rhythmic bands; 1 disables it. Only lavi_numPeaks is taken from the smoothed spectrum; every other statistic uses the raw spectrum. Default is 5.

min_cycles : float, optional

The minimum number of cycles of f that the pairs used must span for f to be used. Default is 40, which keeps the upward bias of LAVI to about 0.02.

Returns:

Statistics over the frequencies that qualify. NaN if the series is constant or fewer than 10 frequencies qualify. For the rhythmicity spectrum LAVI(f):

  • lavi_mean, lavi_median, lavi_std, lavi_iqr, lavi_max, lavi_min: level and spread of rhythmicity across frequency (lavi_median is the baseline).

  • lavi_maxOnMedian: peak rhythmicity relative to the baseline.

  • lavi_freqAtMax, lavi_freqAtMin: angular frequency (2*pi*f) of the most sustained and most transient band.

  • lavi_meanAboveMedian: mean amount by which above-baseline frequencies exceed the baseline.

  • lavi_fracSustained, lavi_fracTransient: fraction of frequencies more than 2 scaled MAD above / below the baseline.

  • lavi_ac1: lag-1 autocorrelation of the spectrum across frequency.

  • lavi_numPeaks: number of local maxima of the smoothed spectrum that exceed the baseline by at least one scaled MAD (distinct rhythmic bands).

  • lavi_atPeakPower: LAVI at the frequency with the most wavelet power.

For the excess over the spectral null, LAVI(f) - LAVI0(f):

  • ex_mean, ex_std, ex_max, ex_min: level and spread of the rhythmicity that the power spectrum does not account for.

  • ex_freqAtMax, ex_freqAtMin: angular frequency of the largest and smallest excess.

  • ex_atPeakPower: the excess at the frequency with the most wavelet power; is the dominant oscillation more sustained (> 0) or more bursty (< 0) than its spectrum implies?

Return type:

dict or float

Notes

Scaled MAD is 1.4826 times the median absolute deviation from the median. The wavelet at frequency f is a Gaussian-windowed complex exponential with envelope standard deviation m / (2 pi f) samples, truncated at 3.5 standard deviations, made zero-mean and normalized to unit energy.

For a broadband series with no rhythm, the frequency locations (*_freqAtMax) tend to fall at the lowest qualifying frequency, where a small residual upward bias of LAVI (about 0.01 to 0.02 at min_cycles = 40) remains; that frequency depends on N.

References

Karvat, G., et al. (2026). Universal rhythmic architecture uncovers two modes of neural dynamics. Nature Communications, 17.

Fransen, J. W., van Ede, F., & Maris, E. (2015). Identifying neuronal oscillations using rhythmicity. NeuroImage, 118, 256-267. doi:10.1016/j.neuroimage.2015.06.003