2010-11-16 42 views

Répondre

1

Peut-être que cette Matlab toolbox aidera; il est assez facile de traduire Matlab en Python, en général.

3

Le meilleur package pour cela dans la terre python est http://pypi.python.org/pypi/nitime

Il a plusieurs estimateurs de cohérence, mais je ne regarde pas très attentivement ceux-ci. C'est un paquetage pour la neuro-imagerie, mais les algorithmes n'utilisent que numpy et scipy, intentionnellement, pour pouvoir être utilisés par d'autres applications.

1

Voici une fonction qui repose sur la fonction scipy.spectrogram (version scipy> 0.17) et calcule la bicohérence entre deux signaux.

Définition Hagihira 2001 et 2007. Voir Hayashi Wikipedia-bicoherence

Hope this helps.

Cordialement,

def compute_bicoherence(s1, s2, rate, nperseg=1024, noverlap=512): 
    """ Compute the bicoherence between two signals of the same lengths s1 and s2 
    using the function scipy.signal.spectrogram 
    """ 
    from scipy import signal 
    import numpy 
    # compute the stft 
    f1, t1, spec_s1 = signal.spectrogram(s1, fs = rate, nperseg = nperseg, noverlap = noverlap, mode = 'complex',) 
    f2, t2, spec_s2 = signal.spectrogram(s2, fs = rate, nperseg = nperseg, noverlap = noverlap, mode = 'complex') 

    # transpose (f, t) -> (t, f) 
    spec_s1 = numpy.transpose(spec_s1, [1, 0]) 
    spec_s2 = numpy.transpose(spec_s2, [1, 0]) 

    # compute the bicoherence 
    arg = numpy.arange(f1.size/2) 
    sumarg = arg[:, None] + arg[None, :] 
    num = numpy.abs(
     numpy.mean(spec_s1[:, arg, None] * spec_s1[:, None, arg] * numpy.conjugate(spec_s2[:, sumarg]), 
     axis = 0) 
     ) ** 2 
    denum = numpy.mean(
     numpy.abs(spec_s1[:, arg, None] * spec_s1[:, None, arg]) ** 2, axis = 0) * numpy.mean(
      numpy.abs(numpy.conjugate(spec_s2[:, sumarg])) ** 2, 
      axis = 0) 
    bicoh = num/denum 
    return f1[arg], bicoh 

# exemple of use and display 
freqs, bicoh = compute_bicoherence(s1, s2, rate) 
f = plt.figure(figsize = (9, 9)) 
plt.pcolormesh(freqs, freqs, bicoh, 
    # cmap = 'inferno' 
    ) 
plt.colorbar() 
plt.clim(0, 0.5) 
plt.show()