Сложная кросс-спектральная плотность

Question

mlab.csd от matplotlib: http://matplotlib.org/api/mlab_api.html#matplotlib.mlab.csd можно использовать для получения real value cross spectral density. Если я хочу получить информацию о фазе из спектральной плотности, мне нужен расчет csd, который возвращает комплексные значения. Есть ли это?

Answer 1

Это обсуждается, например, в этом ответе: https://stackoverflow.com/a/29306730/3920342

Если вы используете csd библиотеки mlab, вы получите комплексные значения, чтобы вы могли рассчитать фазовые углы (и вещественную когерентность). В следующем коде 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") 

Вот реальная и мнимая часть поперечного спектральной плотности (!):

Этот код взят из вопроса How to use the cross-spectral density to calculate the phase shift of two related signals для создания двух сигналов s1 и 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