Lomb-Scargle与FFT功率谱：在等间隔数据上崩溃的比较

Question

我正在尝试创建一些计算均匀和非均匀采样数据的功率谱的例程，使用Lomb-Scargle周期图(LSP)和FFT-Power spectrum。我遇到的问题是，在使用scipy中的LSP实现时，对于均匀采样的数据会崩溃。

下面的代码有效，并且产生几乎相同（并且正确）的输出，就我所知。然而，我被迫在Lomb-Scargle函数中插入一个hack来添加频率抖动，以便它们不会完全匹配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()

- Mike

---

不知道你的方法，你的代码在5处产生了一个漂亮的脉冲，这是你指定的频率。也许可以详细说明一下你的问题？ - deinonychusaur

问题在于Lomb-Scargle实现中，我必须随机扰动输入时间戳，以防止函数因除零错误而崩溃。我会在那里加强注释，方便未来读者理解。 - Mike

那么我的猜测是，由于它是为不均匀采样率设计的，因此具有完全均匀的采样率会导致崩溃——它会除以采样率的变化。如果是这种情况，您的回溯应该会提示。 - deinonychusaur

Answer 1

这是一个关于 lombscargle 的 bug。代码中包含一个 arctan 计算，实现为 atan(2 * cs / (cc - ss))，其中 cc 和 ss 取决于 x 和 freqs 的元素。对于某些输入，cc - ss 可能为 0。修复后的代码使用 atan2(2 * cs, cc - ss)，已经包含在 scipy 0.15.0 中。

- Warren Weckesser