在 python 中计算连贯性答案

【问题标题】：compute coherence in python在 python 中计算连贯性
【发布时间】：2018-07-10 06:20:54
【问题描述】：

我正在学习跨谱和连贯性。据我了解，相干性就像相关性的类似物，因为您可以通过单个功率谱的乘积对交叉谱进行归一化：

这是我当前的 python 实现

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

我计算出的连贯性是一个 1 的数组。我做错了什么？

此外，有时连贯图有向下指向的尖峰（如scipy.signal.coherence 的输出，在其他指向上方的地方（例如here）。我对连贯性的解释有点困惑，不应该更大的相干性意味着该频率的 2 个时间序列之间存在协变？

提前致谢。

【问题讨论】：

标签： python signal-processing fft

【解决方案1】：

您应该一直在使用 welch 方法。作为示例，附加的代码与您的类似（经过一些简化），具有预期的结果。

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')

和可视化

【讨论】：