我在学习跨谱和连贯性。据我所知,相干性就像是相关的类比,因为你通过单个功率谱的乘积来规范化互谱:
这是我当前的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)。我对连贯性的解释有点困惑,难道较大的连贯性不应该意味着在该频率下两个时间序列之间的协变量吗?提前谢谢。
最佳答案
你应该用韦尔奇方法。作为一个示例,附加了与您的代码类似的代码(有一些简化),并得到了预期的结果。
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')
以及形象化