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 = 60) → 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) – Convergence tolerance for Wilson algorithm iterations.

  • max_iterations (int, default=60) – Maximum number of iterations before stopping algorithm.

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
>>> from scipy.fft import fft
>>> # Create a simple 2x2 cross-spectral matrix
>>> n_times, n_freqs, n_signals = 10, 32, 2
>>> # Generate random coefficients for AR process
>>> coeffs = np.random.randn(n_times, n_freqs, n_signals, n_signals)
>>> cross_spec = np.matmul(coeffs, coeffs.conj().swapaxes(-1, -2))
>>>
>>> # Compute minimum phase decomposition
>>> min_phase = minimum_phase_decomposition(cross_spec)
>>>
>>> # Verify decomposition: should reconstruct original matrix
>>> reconstructed = np.matmul(min_phase, min_phase.conj().swapaxes(-1, -2))
>>> error = np.abs(reconstructed - cross_spec).max()
>>> print(f"Reconstruction error: {error:.2e}")

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.

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.