2016-03-08 1 views
3

У меня есть два сигналаПочему кросс-спектры отличаются в mlab и scipy.signal?

import numpy as np 
import matplotlib.pyplot as plt 
from matplotlib import mlab 
import mpld3 
from scipy import signal 

mpld3.enable_notebook() 

nfft = 256 
dt = 0.01 
t = np.arange(0, 30, dt) 
nse1 = np.random.randn(len(t)) * 0.1    # white noise 1 
nse2 = np.random.randn(len(t)) * 0.1    # white noise 2 



# two signals with a coherent part and a random part 
s1 = np.sin(2*np.pi*1*t) + nse1 
s2 = np.sin(2*np.pi*1*t+np.pi) + nse2 

plt.plot(s1, 'r', s2, 'g') 
plt.show() 

Я хотел бы получить когерентность

cxy, fcoh = plt.cohere(s1, s2, nfft, 1./dt) 
fcoh,cxy = signal.coherence(s1,s2, nfft=nfft, fs=1./dt) 
plt.hold(True) 
plt.plot(fcoh, cxy) 
#plt.xlim(0, 5) 
plt.show() 

Coherence

и фазовый сдвиг

(csd, f) = mlab.csd(s1, s2, NFFT=nfft, Fs=1./dt) 
fig = plt.figure() 
angle = np.angle(csd,deg =False) 
angle[angle<-np.pi/2] += 2*np.pi 

plt.plot(f, angle, 'g') 
plt.hold(True) 

(f, csd) = signal.csd(s1, s2, fs=1./dt, nfft=nfft) 

angle = np.angle(csd,deg =False) 
angle[angle<-np.pi/2] += 2*np.pi 

plt.plot(f, angle,'r') 


#plt.xlim(0,5) 
plt.show() 

enter image description here

Я попытался использовать scipy и mlab. Может ли кто-нибудь объяснить, почему я получаю разные результаты?

ответ

4

Поскольку две функции имеют разные значения по умолчанию для некоторых параметров.

Например, если вы передаете к plt.cohere() вариант noverlap=128 вы получите почти идеальный матч с numpy.signal() решения: enter image description here

Помимо для небольшого несовпадения при 0 частоте Гц, и мы действительно не заботятся о когерентности компонентов DC мы? Готов поспорить, что если вы углубитесь в документацию, вы обнаружите еще одну небольшую причуду в стандартных значениях этих двух.