我试图创建一些例程来计算均匀和不均匀采样数据的功率谱,使用Lomb Scargle周期图(LSP)和FFT功率谱。我遇到的问题是,当在scipy中使用LSP实现时,我会遇到均匀采样数据崩溃的情况。在
就我所知,下面的代码可以工作,并生成几乎相同(且正确)的输出。然而,我被迫在Lomb Scargle函数中插入一个kludge来增加频率抖动,因此它们并不完全匹配FFT。当我注释掉那一行时,我得到一个除以零的错误。在
这是scipy中Lomb Scargle实现的问题,还是我不应该对均匀采样的数据使用它??提前谢谢。在
import numpy as np
import scipy.signal as sp
import matplotlib.pyplot as plt
def one_sided_fft(t,x):
full_amplitude_spectrum = np.abs(np.fft.fft(x))/x.size
full_freqs = np.fft.fftfreq(x.size, np.mean(np.ediff1d(t)))
oneinds = np.where(full_freqs >=0.0)
one_sided_freqs = full_freqs[oneinds]
one_sided_amplitude_spectrum=2*full_amplitude_spectrum[oneinds]
return one_sided_freqs, one_sided_amplitude_spectrum
def power_spectrum(t,x):
onef, oneamps = one_sided_fft(t,x)
return onef, oneamps**2
def lomb_scargle_pspec(t, x):
tstep = np.mean(np.ediff1d(t))
freqs = np.fft.fftfreq(x.size, tstep)
idxx = np.argsort(freqs)
one_sided_freqs = freqs[idxx]
one_sided_freqs = one_sided_freqs[one_sided_freqs>0]
#KLUDGE TO KEEP PERIODOGRAM FROM CRASHING
one_sided_freqs = one_sided_freqs+0.00001*np.random.random(one_sided_freqs.size)
#THE FOLLOWING LINE CRASHES WITHOUT THE KLUDGE
pgram = sp.lombscargle(t, x, one_sided_freqs*2*np.pi)
return one_sided_freqs, (pgram/(t.size/4))
if __name__ == "__main__":
#Sample data
fs = 100.0
fund_freq=5
ampl = 0.4
t = np.arange(0,10,1/fs)
x = ampl*np.cos(2*np.pi*fund_freq*t)
#power spectrum calculations
powerf, powerspec = power_spectrum(t,x)
lsf, lspspec = lomb_scargle_pspec(t,x)
#plotting
fig, (ax0, ax1, ax2)= plt.subplots(nrows=3)
fig.tight_layout()
ax0.plot(t, x)
ax0.set_title('Input Data, '+str(fund_freq)+' Hz, '+
'Amplitude: '+str(ampl)+
' Fs = '+str(fs)+' Hz')
ax0.set_ylabel('Volts')
ax0.set_xlabel('Time[s]')
ax1.plot(powerf, powerspec)
ax1.set_title('FFT-based Power Spectrum')
ax1.set_ylabel('Volts**2')
ax1.set_xlabel('Freq[Hz]')
ax2.plot(lsf, lspspec)
ax2.set_title('Lomb-Scargle Power Spectrum')
ax2.set_ylabel('Volts**2')
ax2.set_xlabel('Freq[Hz]')
plt.show()
它是一个bug in ^{} 。代码包含一个实现为}依赖于}的元素。对于某些输入,
atan(2 * cs / (cc - ss))
的arctan计算,其中cc
和{x
和{cc - ss
可以是0。使用atan2(2 * cs, cc - ss)
的修复代码包含在scipy0.15.0中。在相关问题 更多 >
编程相关推荐