scipy.signal.

lombscargle#

scipy.signal.lombscargle(x, y, freqs, precenter=False, normalize=False, *, weights=None, floating_mean=False)[原始碼]#

計算廣義 Lomb-Scargle 週期圖。

Lomb-Scargle 週期圖由 Lomb [1] 開發，並由 Scargle [2] 進一步擴展，用於尋找和測試時間採樣不均勻的微弱週期訊號的顯著性。此處使用的演算法基於 y(ω) = a*cos(ω*x) + b*sin(ω*x) + c 形式的加權最小平方擬合，其中針對每個頻率獨立計算擬合。此演算法由 Zechmeister 和 Kürster 開發，通過啟用個別樣本的加權並計算未知的 y 偏移（也稱為「浮動平均」模型）[3]，改進了 Lomb-Scargle 週期圖。有關更多詳細資訊和實際考量，請參閱關於 Lomb-Scargle 週期圖的優秀參考文獻 [4]。

當 normalize 為 False（或 “power”）（預設值）時，計算出的週期圖為未標準化的，對於振幅為 A 的諧波訊號，對於足夠大的 N，其值為 (A**2) * N/4。其中 N 是 x 或 y 的長度。

當 normalize 為 True（或 “normalize”）時，計算出的週期圖會根據資料在常數參考模型（在零處）附近的殘差進行標準化。

當 normalize 為 “amplitude” 時，計算出的週期圖是振幅和相位的複數表示。

輸入陣列應為實數浮點資料類型的 1 維陣列，在處理前會轉換為 float64 陣列。

參數:

xarray_like: 樣本時間。
yarray_like: 量測值。數值假定基準線為 y = 0。如果有可能存在 y 偏移，建議將 floating_mean 設定為 True。
freqsarray_like: 輸出週期圖的角頻率（例如，對於單位為 s 的 x，單位為 rad/s=2π/s）。頻率通常 >= 0，因為 -freq 處的任何峰值也將存在於 +freq 處。
precenterbool, optional: 若為 True，則通過減去平均值來預先居中量測值。這是一個舊版參數，如果 floating_mean 為 True，則不需要此參數。
normalizebool | str, optional: 計算標準化或複數（振幅 + 相位）週期圖。有效選項為：False/"power"、True/"normalize" 或 "amplitude"。
weightsarray_like, optional: 每個樣本的權重。權重必須為非負數。
floating_meanbool, optional: 若為 True，則獨立確定每個頻率的 y 偏移。否則，y 偏移假定為 0。

返回:

pgramarray_like: Lomb-Scargle 週期圖。

引發:

ValueError: 如果任何輸入陣列 x、y、freqs 或 weights 不是 1 維陣列，或者任何陣列的長度為零。或者，如果輸入陣列 x、y 和 weights 的形狀彼此不相同。
ValueError: 如果任何權重 < 0，或者權重總和 <= 0。
ValueError: 如果 normalize 參數不是允許的選項之一。

另請參閱

periodogram: 使用週期圖的功率譜密度
welch: Welch 方法的功率譜密度
csd: Welch 方法的交叉譜密度

註解

除非 floating_mean 為 True，否則使用的演算法不會自動考慮任何未知的 y 偏移。因此，對於大多數用例，如果有可能存在 y 偏移，建議將 floating_mean 設定為 True。如果 precenter 為 True，它會執行操作 y -= y.mean()。但是，precenter 是一個舊版參數，當 floating_mean 為 True 時不需要此參數。此外，precenter 移除的平均值不考慮樣本權重，也不會校正由於在峰值和/或谷值處始終遺漏觀測而導致的任何偏差。當 normalize 參數為 “amplitude” 時，對於 freqs 中低於 (2*pi)/(x.max() - x.min()) 的任何頻率，預測的振幅將趨於無窮大。「奈奎斯特頻率」限制（請參閱奈奎斯特-香農採樣定理）的概念通常不適用於不均勻採樣的資料。因此，對於不均勻採樣的資料，freqs 中的有效頻率通常可能遠高於預期。

參考文獻

[1]

N.R. Lomb “不均勻間隔資料的最小平方頻率分析”，《天文學和天體物理學》，第 39 卷，第 447-462 頁，1976 年 DOI:10.1007/bf00648343

[2]

J.D. Scargle “天文時間序列分析研究。II - 不均勻間隔資料的頻譜分析的統計方面”，《天體物理學雜誌》，第 263 卷，第 835-853 頁，1982 年 DOI:10.1086/160554

[3]

M. Zechmeister 和 M. Kürster，“廣義 Lomb-Scargle 週期圖。浮動平均和克卜勒週期圖的新形式體系”，《天文學和天體物理學》，第 496 卷，第 577-584 頁，2009 年 DOI:10.1051/0004-6361:200811296

[4]

J.T. VanderPlas，“理解 Lomb-Scargle 週期圖”，《天體物理學雜誌增刊系列》，第 236 卷，第 1 期，第 16 頁，2018 年 5 月 DOI:10.3847/1538-4365/aab766

範例

>>> import numpy as np
>>> rng = np.random.default_rng()

首先為訊號定義一些輸入參數

>>> A = 2.  # amplitude
>>> c = 2.  # offset
>>> w0 = 1.  # rad/sec
>>> nin = 150
>>> nout = 1002

隨機產生樣本時間

>>> x = rng.uniform(0, 10*np.pi, nin)

繪製選定時間的正弦波

>>> y = A * np.cos(w0*x) + c

定義要計算週期圖的頻率陣列

>>> w = np.linspace(0.25, 10, nout)

計算每個標準化選項的 Lomb-Scargle 週期圖

>>> from scipy.signal import lombscargle
>>> pgram_power = lombscargle(x, y, w, normalize=False)
>>> pgram_norm = lombscargle(x, y, w, normalize=True)
>>> pgram_amp = lombscargle(x, y, w, normalize='amplitude')
...
>>> pgram_power_f = lombscargle(x, y, w, normalize=False, floating_mean=True)
>>> pgram_norm_f = lombscargle(x, y, w, normalize=True, floating_mean=True)
>>> pgram_amp_f = lombscargle(x, y, w, normalize='amplitude', floating_mean=True)

現在繪製輸入資料的圖

>>> import matplotlib.pyplot as plt
>>> fig, (ax_t, ax_p, ax_n, ax_a) = plt.subplots(4, 1, figsize=(5, 6))
>>> ax_t.plot(x, y, 'b+')
>>> ax_t.set_xlabel('Time [s]')
>>> ax_t.set_ylabel('Amplitude')

然後繪製每個標準化選項以及 floating_mean=True 和 floating_mean=False 的週期圖

>>> ax_p.plot(w, pgram_power, label='default')
>>> ax_p.plot(w, pgram_power_f, label='floating_mean=True')
>>> ax_p.set_xlabel('Angular frequency [rad/s]')
>>> ax_p.set_ylabel('Power')
>>> ax_p.legend(prop={'size': 7})
...
>>> ax_n.plot(w, pgram_norm, label='default')
>>> ax_n.plot(w, pgram_norm_f, label='floating_mean=True')
>>> ax_n.set_xlabel('Angular frequency [rad/s]')
>>> ax_n.set_ylabel('Normalized')
>>> ax_n.legend(prop={'size': 7})
...
>>> ax_a.plot(w, np.abs(pgram_amp), label='default')
>>> ax_a.plot(w, np.abs(pgram_amp_f), label='floating_mean=True')
>>> ax_a.set_xlabel('Angular frequency [rad/s]')
>>> ax_a.set_ylabel('Amplitude')
>>> ax_a.legend(prop={'size': 7})
...
>>> plt.tight_layout()
>>> plt.show()

../../_images/scipy-signal-lombscargle-1.png