2014-12-11 1 views
1

mlab.csd from matplotlib : http://matplotlib.org/api/mlab_api.html#matplotlib.mlab.csd실제 값은 교차 스펙트럼 밀도입니다. 스펙트럼 밀도에서 위상 정보를 얻으려면 복잡한 값을 반환하는 csd 계산이 필요합니다. 하나 있습니까?복잡한 교차 스펙트럼 밀도

+1

(!) 나는 아마 그것을 직접 계산할 것이다. Matplotlib은 편의를 위해 이러한 것들을 포함하고 있지만, 데이터를 가지고 정확히 무슨 일이 일어나는지 직접 알지 못한다. – Ajean

답변

0

이 대답은 : https://stackoverflow.com/a/29306730/3920342

mlab 라이브러리의 csd를 사용하면 복잡한 값을 얻을 수 있으므로 위상 각 (및 실수 값 일관성)을 계산할 수 있습니다. 다음 코드에서 s1과 s2는 (시간 영역에서) 상관 될 두 개의 신호를 포함합니다. 여기

from matplotlib import mlab 

# First create power sectral densities for normalization 
(ps1, f) = mlab.psd(s1, Fs=1./dt, scale_by_freq=False) 
(ps2, f) = mlab.psd(s2, Fs=1./dt, scale_by_freq=False) 
plt.plot(f, ps1) 
plt.plot(f, ps2) 

# Then calculate cross spectral density 
(csd, f) = mlab.csd(s1, s2, NFFT=256, Fs=1./dt,sides='default', scale_by_freq=False) 
fig = plt.figure() 
ax1 = fig.add_subplot(1, 2, 1) 
# Normalize cross spectral absolute values by auto power spectral density 
ax1.plot(f, np.absolute(csd)**2/(ps1 * ps2)) 
ax2 = fig.add_subplot(1, 2, 2) 
angle = np.angle(csd, deg=True) 
angle[angle<-90] += 360 
ax2.plot(f, angle) 

# zoom in on frequency with maximum coherence 
ax1.set_xlim(9, 11) 
ax1.set_ylim(0, 1e-0) 
ax1.set_title("Cross spectral density: Coherence") 
ax2.set_xlim(9, 11) 
ax2.set_ylim(0, 90) 
ax2.set_title("Cross spectral density: Phase angle") 

크로스 스펙트럼 밀도의 실수 부 및 허수 부 : real and imaginary part of the cross spectral density

이 코드는 두 개의 신호 (S1)를 만드는 문제 How to use the cross-spectral density to calculate the phase shift of two related signals에서 촬영하고, S2 :

""" 
Compute the coherence of two signals 
""" 
import numpy as np 
import matplotlib.pyplot as plt 

# make a little extra space between the subplots 
plt.subplots_adjust(wspace=0.5) 

nfft = 256 
dt = 0.01 
t = np.arange(0, 30, dt) 
nse1 = np.random.randn(len(t))     # white noise 1 
nse2 = np.random.randn(len(t))     # white noise 2 
r = np.exp(-t/0.05) 

cnse1 = np.convolve(nse1, r, mode='same')*dt # colored noise 1 
cnse2 = np.convolve(nse2, r, mode='same')*dt # colored noise 2 

# two signals with a coherent part and a random part 
s1 = 0.01*np.sin(2*np.pi*10*t) + cnse1 
s2 = 0.01*np.sin(2*np.pi*10*t) + cnse2