compute coherence in python

Question

I'm learning cross-spectrum and coherence. From what I understand, coherence is like the analogue of correlation in that you normalize the cross-spectrum by the product of individual power spectrum:

Here is my current python implementation

import numpy

def crossSpectrum(x,y):

    #-------------------Remove mean-------------------
    xp=x-numpy.mean(x)
    yp=y-numpy.mean(y)
    n=len(x)

    # Do FFT
    cfx=numpy.fft.fft(xp)/n
    cfy=numpy.fft.fft(yp)/n
    freq=numpy.fft.fftfreq(n)

    # Get cross spectrum
    cross=cfx.conj()*cfy

    return cross,freq


#-------------Main---------------------------------
if __name__=='__main__':

    x=numpy.linspace(-250,250,500)
    noise=numpy.random.random(len(x))
    y=10*numpy.sin(2*numpy.pi*x/10.)+5*numpy.sin(2*numpy.pi*x/5.)+\
            2*numpy.sin(2*numpy.pi*x/20.)+10
    y+=noise*10
    y2=5*numpy.sin(2*numpy.pi*x/10.)+5+noise*50

    p11,freq=crossSpectrum(y,y)
    p22,freq=crossSpectrum(y2,y2)
    p12,freq=crossSpectrum(y,y2)

    # coherence
    coh=numpy.abs(p12)**2/p11.real/p22.real
    print coh

And my computed coherence is an array of 1s. What am I doing wrong?

Also, sometimes the coherence plot has downward pointing spikes (like the output of scipy.signal.coherence, in other places that are pointing upward (e.g. here). I'm bit confused by the interpretation of coherence, shouldn't a larger coherence implies covariability between the 2 timeseries at that frequency?

Thanks in advance.

Answer 1

You should have been using the welch method. As an example attached code similar to yours (with some simplifications) with expected results.

import numpy 
from matplotlib.pyplot import plot, show, figure, ylim, xlabel, ylabel
def crossSpectrum(x, y, nperseg=1000):

#-------------------Remove mean-------------------
cross = numpy.zeros(nperseg, dtype='complex128')
for ind in range(x.size / nperseg):

    xp = x[ind * nperseg: (ind + 1)*nperseg] 
    yp = y[ind * nperseg: (ind + 1)*nperseg] 
    xp = xp - numpy.mean(xp)
    yp = yp - numpy.mean(xp)

    # Do FFT
    cfx = numpy.fft.fft(xp)
    cfy = numpy.fft.fft(yp)

    # Get cross spectrum
    cross += cfx.conj()*cfy
freq=numpy.fft.fftfreq(nperseg)
return cross,freq

#-------------Main---------------------------------
if __name__=='__main__':

x=numpy.linspace(-2500,2500,50000)
noise=numpy.random.random(len(x))
y=10*numpy.sin(2*numpy.pi*x)
y2=5*numpy.sin(2*numpy.pi*x)+5+noise*50

p11,freq=crossSpectrum(y,y)
p22,freq=crossSpectrum(y2,y2)
p12,freq=crossSpectrum(y,y2)

# coherence
coh=numpy.abs(p12)**2/p11.real/p22.real
plot(freq[freq > 0], coh[freq > 0])
xlabel('Normalized frequency')
ylabel('Coherence')

and visualization