pyhctsa.operations.correlation.cyclostationarity_scan

pyhctsa.operations.correlation.cyclostationarity_scan(y, max_lag=20, n_smooth=16, min_cycles=4, max_ar_order=10, peak_level=0.05)

Scan for second-order cyclostationarity: periodic variation of the autocovariance.

Detects a series whose variability, or the correlation between nearby samples, rises and falls on a regular cycle, even though nothing periodic shows in the series itself: its mean stays constant and its power spectrum has no peak at that period. Examples are noise whose loudness swells and fades regularly, a process whose memory strengthens and weakens with a period, and seasonal series whose volatility follows the season.

Formally, a series is second-order cyclostationary (periodically correlated) when its autocovariance at some lag tau varies periodically in time, with cyclic frequency alpha: amplitude-modulated noise, dynamics whose coefficients vary with a period, or seasonal series [1], [2]. It need not have a periodic mean or any peak in its power spectrum, so spectral and autocorrelation features miss it. Its signature is Gardner’s spectral-line regeneration [1]: the lag-product series e_t e_{t+tau} of a cyclostationary series has a spectral line (a finite-strength sinusoidal component) at frequency alpha, while that of a stationary series has only a continuous spectrum. The cyclic autocorrelation R^alpha(tau) is that line’s amplitude.

The test for a line at one cyclic frequency and lag is a single-lag form of the Dandawate-Giannakis test [3], applied at every Fourier frequency:

  1. Prewhiten the series by an AR(p) model, p <= max_ar_order chosen by AIC (as in residual_autocorr_scan()), and z-score the residuals e. Prewhitening removes the spectral shape of a stationary series, so that a sharp spectral peak is not read as cyclostationarity.

  2. For each lag tau = 0..``max_lag``, form the mean-removed lag product e_t e_{t+tau} (truncated to even length n) and its periodogram I(alpha_k) at alpha_k = k/n.

  3. Compare each bin with the mean S(alpha_k) of the 2 * n_smooth bins on either side of it (skipping the adjacent bin). Under a locally flat continuum, I/S ~ F(2, 4 n_smooth) (F(1, 4 n_smooth) at alpha = 1/2), which gives a p-value for a line at that bin.

  4. Take the smallest p-value p_min over all M cells (alpha_k with at least min_cycles cycles in the record, and tau). The score -ln(M p_min) is approximately standard Gumbel for a stationary series (mean 0.58, 95th percentile 2.97) at any length, and grows with length for a cyclostationary one.

Deterministic phase-coherent oscillators (e.g. the Rossler or van der Pol flows) score highly. This is correct in the fraction-of-time sense of cyclostationarity [1], since their residual variance is locked to the phase of the oscillation. Smooth low-pass series that an AR(max_ar_order) model under-whitens give scores below the null mean. The cyclic frequency and lag are in samples, so the outputs depend on the sampling rate.

References

Parameters:
y : array-like

The input time series.

max_lag : int, optional

The largest lag tau of the lag products (default: 20). Capped at an eighth of the residual length.

n_smooth : int, optional

The number of periodogram bins on each side used to estimate the continuum at each cyclic frequency (default: 16).

min_cycles : int, optional

The lowest cyclic frequency scanned, in cycles per record (default: 4).

max_ar_order : int, optional

The largest order of the prewhitening AR model (default: 10).

peak_level : float, optional

The significance level at which the strongest line counts as a detected rhythm (default: 0.05). Below it, peak_alpha and peak_lag are NaN.

Returns:

  • max_score: -ln(M p_min) over all cyclic frequencies and lags.

  • lag0_score: the same over lag 0 only, i.e. spectral lines in the squared residuals (periodic variance).

  • peak_alpha: the cyclic frequency of the strongest line, in cycles per sample (its period is 1 / peak_alpha samples).

  • peak_lag: the lag of the strongest line.

peak_alpha and peak_lag locate a rhythm only when there is one: they are NaN unless max_score exceeds the Gumbel null’s 1 - peak_level quantile, -ln(-ln(1 - peak_level)) (2.97 at the default 0.05). Without a significant line, the strongest bin is just the largest noise fluctuation, at a position uniform over the cells scanned, so its location means nothing. At the default level about 5% of stationary series still get values.

Returns NaN for a constant series or fewer than 64 residuals.

Return type:

dict