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 偏移（也称为“浮动均值”模型）改进了 Lomb-Scargle 周期图 [3]。有关更多详细信息和实际考虑因素，请参阅关于 Lomb-Scargle 周期图的优秀参考资料 [4]。当 normalize 为 False（或“power”）（默认值）时，计算的周期图是未标准化的，对于幅度为 A 的谐波信号，在 N 足够大的情况下，其值为
(A**2) * N/4。其中 N 是 x 或 y 的长度。当 normalize 为 True（或“normalize”）时，计算的周期图通过数据围绕恒定参考模型（在零处）的残差进行标准化。
当 normalize 为“amplitude”时，计算的周期图是幅度和相位的复数表示。
输入数组应为实数浮点数据类型的 1-D 数组，在处理之前将其转换为 float64 数组。
- 参数:
- xarray_like
采样时间。
- yarray_like
测量值。假设值的基线为
y = 0。如果可能存在 y 偏移，建议将 floating_mean 设置为 True。- freqsarray_like
输出周期图的角频率（例如，对于单位为 s 的 x，单位为 rad/s=2π/s）。频率通常 >= 0，因为
-freq处的任何峰值也会存在于+freq处。- precenterbool, 可选
如果为 True，则通过减去平均值来预先居中测量值。这是一个遗留参数，如果 floating_mean 为 True，则不必要。
- normalizebool | str, 可选
计算标准化或复数（幅度 + 相位）周期图。有效选项为：
False/"power"、True/"normalize"或"amplitude"。- weightsarray_like, 可选
每个样本的权重。权重必须为非负数。
- floating_meanbool, 可选
如果为 True，则独立确定每个频率的 y 偏移。否则，y 偏移假定为 0。
- 返回:
- pgramarray_like
Lomb-Scargle 周期图。
- 引发:
- ValueError
如果任何输入数组 x、y、freqs 或 weights 不是 1D，或者如果任何一个数组的长度为零。或者，如果输入数组 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 的周期图
>>> 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()
