pyhctsa.operations.spectral.doubled_phase_coherence

pyhctsa.operations.spectral.doubled_phase_coherence(y, harmonic=2, m=5, n_freq=40, f_low=0.02, f_high=0.3, min_cycles=40, pad=4)

Record-long coherence of the doubled phase of each frequency band.

Does an oscillation keep its timing across the whole record, up to a sign? Squaring a band signal X(t) = A(t) exp(i phi(t)) doubles its phase, which makes it blind to a sign flip (phi -> phi + pi). An oscillation whose phase is kept only modulo pi (polarity reversals, bursts of alternating or random sign on one underlying clock, a binary phase-shift-keyed carrier) is incoherent as X but coherent as X**2: X**2 is then a single phasor at twice the carrier frequency, running through the whole record. An oscillation whose phase is re-drawn at random (bursts that reset) stays incoherent in both, and so does any linear Gaussian process, whose band signal is circular (its phase is uniform and independent of its amplitude). The power spectrum cannot see this: flipping the sign of each burst, or re-drawing its phase, leaves the spectrum’s expectation unchanged.

The method, at each of n_freq log-spaced frequencies f from f_low to f_high:

  1. A complex Morlet wavelet transform (m cycles, as in rhythmicity()) gives the band signal X(t, f); only its n_e edge-free samples are used, and f is used only if they span at least min_cycles cycles.

  2. The coherence of X**h (h = harmonic, 2 by default) is its amplitude-weighted resultant length at the best demodulation frequency, C(f) = max_nu |sum_t X(t)**h exp(-2 pi i nu t)| / sum_t |X(t)|**h, with nu searched on a DFT zero-padded to pad * n_e points, over the band where the Gaussian null below gives X**h non-negligible power (at least 1% of its peak; around h f). C is in [0, 1]: 1 when X**h is a single steady phasor (a sinusoid, or a sinusoid of randomly flipping sign at h = 2), and about 1/sqrt(K) for K bursts of independent phase.

  3. The spectral null. Part of C is noise: for a finite record C is positive even for a circular process, by an amount that falls as N**-1/2 and depends on the spectrum. A Gaussian process with the series’ own (zero-padded) periodogram has E X**h = 0 and, by Isserlis’ theorem for circular complex Gaussians (E[X_t**h conj(X_s)**h] = h! c**h), E|S(nu)|**2 = h! sum_k (n_e - |k|) c(k)**h exp(-2 pi i nu k), with S the DFT of X**h and c(k) = E X(t + k) X*(t) the lag covariance of the band signal, computed exactly from the periodogram as in rhythmicity(). Treating the |S|**2 at the n_e-point (unpadded) DFT frequencies in the band as independent exponentials with these means, their expected maximum (one quadrature) gives E0[C(f)**2], including the upward bias of maximizing over nu. The excess C(f)**2 - E0[C(f)**2] is about zero for any linear Gaussian process, whatever its spectrum or length, and positive when the h-fold phase stays coherent for longer than the wavelet resolves.

For h = 2, C is the sample circularity quotient |E z**2| / E|z|**2 (Ollila 2008) of the band signal demodulated at nu / 2 (z = X exp(-i pi nu t)), maximized over nu: 0 for a circular (proper) variable, 1 for one confined to a line through the origin. In the language of cyclostationarity it is the normalized conjugate cyclic autocorrelation of the band signal at lag 0 and cycle frequency nu (Gardner 1991; Gardner et al. 2006), the feature by which squaring recovers the carrier of a binary phase-shift-keyed signal. Doubling an angle to make it insensitive to theta -> theta + pi is the standard treatment of axial data (Mardia & Jupp 2000). Expanding |sum_t X_t**2 e_t|**2 over pairs of times shows that C**2 is also the amplitude-weighted lagged coherence of Fransen et al. (2015), applied to X**2 and pooled over every lag rather than taken at one. At h = 1 the numerator is a smoothed periodogram ordinate, so h = 1 adds nothing to the power spectrum; h >= 2 is a statistic of order 2h, and only h = 2 is in the default feature set. The general-h null is the h-th power of the one at h = 1: for a circular Gaussian pair with lag correlation rho, the Cauchy-Schwarz-normalized coherence of (X_t**h, X_s**h) is rho**h. The null and its expected-maximum correction are this implementation’s own derivation, in the spirit of the analytic LAVI null of rhythmicity() (Karvat et al. 2026): they replace a surrogate distribution by its expectation.

The series is treated as evenly sampled with unit spacing: frequencies are in cycles per sample and reported angular frequencies are 2*pi*f, as in rhythmicity().

Parameters:
y : array-like

The input time series.

harmonic : int, optional

The power h to which the band signal is raised (the phase is multiplied by h). Default is 2 (coherence modulo pi); 1 gives a spectral quantity, kept as a check.

m : float, optional

Width of the Morlet wavelet, in cycles. Default is 5, as in rhythmicity().

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. It must satisfy f_high (1 + 3/m) <= 0.5, so that the wavelet’s passband (spectral standard deviation f/m) stays three standard deviations below the Nyquist frequency. Otherwise the wavelet also passes the negative-frequency image of the real series, X is no longer analytic, E X**2 is no longer zero, and X**2 carries a spurious line even for white noise (at f = 0.4 with m = 5 the excess of white noise grows with N, to about 0.04 at N = 16384, instead of staying at 0). Default is 0.3, the largest round value allowed for m = 5 (limit 0.3125), rather than the 0.4 of rhythmicity().

min_cycles : float, optional

A frequency is used only if the edge-free band signal spans at least this many of its cycles. Default is 40, as in rhythmicity().

pad : int, optional

Zero-padding factor of the DFT of X**h on which the demodulation frequency nu is searched. Default is 4: the maximum of an unpadded DFT can miss a line between two bins and lose up to 36% of C (a pure sinusoid can give 0.64 instead of 1); padding by 4 caps this scalloping loss at 2.6%.

Returns:

Statistics over the frequencies that qualify. NaN if the input is constant or not finite, or if fewer than 10 frequencies qualify (with the defaults, N below 286). For the coherence spectrum C(f):

  • coh_max, coh_freqAtMax: the largest coherence and its angular frequency 2*pi*f.

  • coh_mean: the mean coherence over frequency.

  • coh_atPeakPower: the coherence at the frequency with the most wavelet power.

These raw values include the finite-record noise floor, which falls as N**-1/2 (by a factor of about 2 from N = 500 to 10000 on noise), so they are excluded from the default feature set. For the excess over the spectral null, C(f)**2 - E0[C(f)**2]:

  • ex_max: the largest excess, the most sign-invariant-coherent band. It keeps a mild downward drift with N on noise (about 0.1 at N = 500 to 0.03 at N = 10000), being the maximum of about 40 zero-mean estimates whose spread shrinks with N.

  • ex_mean: the mean excess over frequency.

  • ex_atPeakPower: the excess at the frequency with the most wavelet power: does the dominant oscillation keep its timing (modulo pi at h = 2) across the record?

Typical ex_atPeakPower at N = 4096: about 0 for linear Gaussian processes and the Lorenz system; 0.7 for a sustained noisy sinusoid; 0.8 for bursts of random sign on one clock; 0.3 for bursts with independent phases; 0.15 for Student-t (2 degrees of freedom) noise.

Return type:

dict or float

Notes

  • The excess reads “the doubled phase is coherent for longer than the wavelet resolves”. A single long burst satisfies this within itself, so bursts with independent phases also give positive values; they fall below bursts that keep a common clock only when the record holds many bursts (the coherence of K independent bursts is about 1/sqrt(K)), and with about 3 bursts or fewer per record the two are indistinguishable.

  • Heavy tails give small positive values: each large isolated event makes a short, self-coherent transient in X**2.

  • The null treats the in-band DFT ordinates of X**h as independent exponentials (the DFT of a long stationary record is asymptotically complex Gaussian and independent across frequencies). This is approximate when the band holds few ordinates (short series, low f), and it does not hold for strongly non-Gaussian amplitudes, by design: that departure is part of what the excess reports.

  • About N = 1000 samples or more are needed for reliable values; with the defaults the output is defined from N = 286.

  • The cost is about three times that of rhythmicity() (a 4-times padded FFT and the expected-maximum quadrature at each frequency).

References

Ollila, E. (2008). On the circularity of a complex random variable. IEEE Signal Processing Letters, 15, 841-844. doi:10.1109/LSP.2008.2005050

Gardner, W. A. (1991). Exploitation of spectral redundancy in cyclostationary signals. IEEE Signal Processing Magazine, 8(2), 14-36. doi:10.1109/79.81007

Gardner, W. A., Napolitano, A., & Paura, L. (2006). Cyclostationarity: half a century of research. Signal Processing, 86(4), 639-697. doi:10.1016/j.sigpro.2005.06.016

Mardia, K. V., & Jupp, P. E. (2000). Directional Statistics. Wiley.

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

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