spectral_connectivity.minimum_phase_decomposition.minimum_phase_decomposition#

minimum_phase_decomposition(cross_spectral_matrix: ndarray[tuple[int, ...], dtype[complexfloating]], tolerance: float = 1e-08, max_iterations: int = 500, *, _warn_on_failure: bool = True) → ndarray[tuple[int, ...], dtype[complexfloating]][source]#

Compute minimum phase decomposition using Wilson algorithm.

Finds a minimum phase matrix square root G of the cross-spectral density matrix S such that S = G G^H, where all poles of G are inside the unit circle. This decomposition is essential for computing directed connectivity measures like spectral Granger causality.

Parameters:

cross_spectral_matrix (NDArray[complexfloating],) – shape (n_time_samples, …, n_fft_samples, n_signals, n_signals) Cross-spectral density matrix to be decomposed. Must be Hermitian positive semidefinite for each frequency.
tolerance (float, default=1e-8) – Relative convergence tolerance for Wilson algorithm iterations.
max_iterations (int, default=500) – Maximum number of iterations before stopping algorithm. Near-singular cross-spectral matrices (highly correlated channels) can need several hundred iterations to reach the relative tolerance; the loop returns early as soon as every sub-spectrum has converged.

Returns:

minimum_phase_factor – shape (n_time_samples, …, n_fft_samples, n_signals, n_signals) Minimum phase square root of cross_spectral_matrix. All eigenvalues have negative real parts (minimum phase property).

Return type:

NDArray[complexfloating],

Examples

>>> import numpy as np
>>> rng = np.random.default_rng(0)
>>> n_times, n_freqs, n_signals = 1, 32, 2
>>> # A valid cross-spectral matrix (Hermitian positive definite). Here it is
>>> # constant across frequency (a white process), which factors exactly.
>>> a = rng.standard_normal((n_signals, n_signals))
>>> spd = a @ a.T + n_signals * np.eye(n_signals)
>>> cross_spec = np.tile(spd, (n_times, n_freqs, 1, 1))
>>> min_phase = minimum_phase_decomposition(cross_spec)
>>> # The factor reconstructs the input: G @ G^H == S.
>>> reconstructed = np.matmul(min_phase, min_phase.conj().swapaxes(-1, -2))
>>> error = np.abs(reconstructed - cross_spec).max()
>>> bool(error < 1e-6)
True

Notes

The Wilson algorithm iteratively refines an initial guess using the “plus” operator (causal projection) until convergence. The algorithm may not converge for all time points; warnings are issued when the maximum iteration count is reached.

When the spectrum is exactly conjugate-symmetric, S(-f) == conj(S(f)) bit for bit (as SciPy’s FFT gives for real-valued signals), every iterate keeps that symmetry, so the iteration runs on the non-negative frequencies with real FFTs and mirrors the rest. Spectra symmetric only up to rounding (possible with other FFT backends such as CuPy’s), or not at all (complex-valued signals), take the two-sided iteration, which does about twice the arithmetic per iteration.

Each update is built from a square root L of S (S = L Lᴴ), computed once, which lets sub-spectra with near-collinear channels converge. An indefinite or non-finite S has no such square root and is returned as NaN with the non-convergence warning.

Convergence of the iterate does not by itself guarantee G Gᴴ ≈ S. The factorization assumes the cross-spectrum is resolved finely enough in frequency that the corresponding autocovariance decays within the analysis window; an under-resolved (aliased) spectrum can satisfy the successive- iterate convergence test yet reconstruct S poorly, silently biasing the directed-connectivity measures built on it. Use minimum_phase_reconstruction_error() to check factorization quality explicitly, and a longer FFT (larger n_fft_samples / n_time_samples_per_window) if the error is large.

References

[1]

Wilson, G. T. (1972). The factorization of matricial spectral densities. SIAM Journal on Applied Mathematics, 23(4), 420-426.

[2]

Dhamala, M., Rangarajan, G., & Ding, M. (2008). Analyzing information flow in brain networks with nonparametric Granger causality. NeuroImage, 41(2), 354-362.