[FieldTrip] Spectral factorization: Wilson-Burg algorithm does not converge

[vpapenm at uni-koeln.de]

Dear all,

I use sfactorization_wilson.m (in my own pipeline outside field trip  
actually) to estimate the spectral Granger causality of my time series  
(local field potentials from deep brain stimulation electrodes). As  
input I use a 3D spectral matrix S with dimensions 2x2x32 calculated  
from wavelet transformation at frequencies f=[0:1:31] Hz. It consists  
of time-averaged auto- (Sxx, Syy) and cross-spectral densities (Sxy,  
conj(Sxy)), obtained from smoothing the wavelet matrix.

My problem is that the implemented Wilson-Burg algorithm quite often  
does not converge at all and I wonder if anybody else has that  
problem. The non-convergence leads to a new spectral matrix  
Snew=psi*psi', which has non-zero imaginary parts or which is simply  
too different from the original matrix S.

I also wonder why the only criteria to stop the iteration is when the  
increment step psi_old-psi becomes small. Why is there no criterion to  
stop the iteration if the new spectral matrix is close enough to the  
original one (e.g. S-Snew<tol)?

At the moment I have tried two different approaches to estimate the  
power at f=0 Hz: 1) I use the power at f=0Hz from FFT [which is mean(  
sum(x).^2/N*dt ); for the auto-spectra, e.g.] and 2) I use a linear  
extrapolation of the PSD to get P(f=0). However, both methods do not  
lead to sufficient convergence of the algorithm.

I would be very grateful if someone found the time to answer.

Best regards,
Mitch