Spectral connectivity from wavelet transform

[mgm248]

Hello,

I am trying to implement the spectral connectivity measures (specifically granger causality) using Fourier coefficients obtained from wavelet transforms of my data. My issue is I only get NaNs back from the function. Things go awry on the line

minnimum_phase_factor = (minimum_phase_decomposition(cross_spectral_matrix[
                            ..., pair_indices, pair_indices.T]))

Which hits the error numpy.linalg.LinAlgError: Singular matrix, probably on this line

linear_predictor = _get_linear_predictor( minimum_phase_factor, cross_spectral_matrix, I) (Although the command seems to run fine for a good number of iterations before it doesn't).

Have you come across this error before? Any idea what could be going wrong?

[edeno]

Well, this is a failure of the optimization algorithm to find minimum phase matrix square root of the cross spectral density.

A couple questions:

1. Are you passing the Connectivity class both positive and negative frequencies from the wavelet transform?
2. Are you computing this for multiple signals? It might be possible to remove the offending signal.
3. Have you computed the power spectrum for your signals to make sure the power spectrums make sense (this might help diagnose whatever problem the offending signal has)?

[mgm248]

I think 1 is my problem -- I am only passing positive frequencies from the wavelet transform. I obtain these using a package (MNE), where I don't have the option to return negative frequencies from the transform. So it sounds like I'll need to do things differently than that. Sorry, I hadn't realized this requirement; thanks for your help!

[edeno]

If your signal is real (as in real number valued), the negative fourier coefficients should be conjugate symmetric to the positive fourier coefficients. This means you can just take the complex conjugate of the positive fourier coefficients and use them as your negative fourier coefficients.

[mgm248]

Ah, that's really helpful, thanks!

[mgm248]

I'm still not having much luck. I'm wondering if the problem might be with the size of the cross-spectral matrix that gets input into the minimum_phase_decomposition. With wavelet decomposition, I input a matrix that has 1400 points in time (one for each sample), and 200 frequencies (including the negative coefficients). Maybe the default tolerance is designed for much smaller matrices. If I set the tolerance to be much higher, I get better results; a tolerance of 5000 gets 1231 of 1400 timepoints to converge, with the max_iterations still at 60. Maybe upping the tolerance that much is a bad idea, so I could increase the max iterations and try to do better, but does it seem like I am going down a dangerous path, or is this ok to do?

[edeno]

4000 is definitely too high. I would be more confident in something like 1e-2 or 1e-3, although of course you are subject to more error. In my experience not everything converges using this algorithm, although I haven't had time to think about what controls that convergence.