信息处理 - 理解周期图 - 吾爱随笔录

我正在尝试通过遵循此处的 astropy Lomb-scargle 周期图教程来使用周期图来判断信号何时是周期性的。

http://docs.astropy.org/en/stable/stats/lombscargle.html

我模拟了一些数据，一个是正弦曲线（周期 = 200），另一个是斜高斯曲线（即单个瞬态事件）。希望周期图能够为周期性对象挑选周期，并为瞬态给出一个周期，这意味着窗口中发生了单个事件。

不幸的是，结果根本没有意义。我在最后附上了代码和下面生成的数字。由于我的模拟中存在随机噪声，每个结果都不同，我在下面提供了两个结果示例。我使用上面链接中概述的相同方法，我们LombScargle.model()在最佳拟合频率上使用函数frequency[argmax[power]]

红线是我模拟数据的真实函数。绿色是周期图中的最佳拟合。右侧图是周期图中的 PSD。

示例 1

在这里，我们可以看到挑选出的正弦曲线（右上图）的最佳拟合频率是 0.105（即 9-10 天的周期），这并不接近于 0.05 的频率，我期望的周期为 200 天，但是当我将这个 0.105 的最佳拟合频率提供给 lomb-scargle 模型拟合器时，会得出一条周期为 200 天的完美匹配的周期曲线？

这根本不符合逻辑。

示例 2

在这里我再次运行代码，这一次结果被切换了？它适合瞬态的一个非常大的周期，因此我可以自信地说这是一个单一的瞬态事件，但正弦拟合很糟糕。最佳频率仍然是 0.105（周期 = 10），但 lomb-scargle 模型拟合器覆盖了似乎有 60 天周期的东西，这是错误的吗？

如果我做错了什么，我可以得到澄清吗？有人告诉我，周期图实际上是用于像这样的不均匀采样数据的事实上的工具……结果似乎有一半时间很糟糕。

澄清我的问题

能否解释在第一个图中，我输入 astropy lomb-scargle 模型拟合器的最佳拟合频率 0.105 以某种方式创建了频率为 0.05 的正确匹配正弦曲线？解释是什么？
当我只期望 1 时，为什么示例 1 中的周期图的右上角有 5 个强峰？中间两个接近0.05的真实值（分别为0.045和0.055）

这是用于模拟和绘制数据的短代码

import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as sc
from astropy.stats import LombScargle
import math


#simulate parameters
data_range = [i for i in range(1,1001)]
number_of_samples = 50
gauss_skew = sc.skewnorm.pdf
skew = -10
period = 200
location = data_range[int(len(data_range)/2)]


y1= [(2 * (1. + np.sin(2. * np. pi * x/period)) + np.random.normal(loc =0.0, scale = 0.5)) for x in data_range]
errors1 = [np.random.normal(loc = 0.0, scale = 1) for x in data_range]
y2 = [(1000* (gauss_skew(x,skew,loc=location ,scale = 50)) + np.random.normal(loc =0.0, scale = 1)) for x in data_range]
errors2 = [np.random.normal(loc = 0.0, scale = 1) for x in data_range]

sample_rate = int(len(data_range)/number_of_samples)# To thin the data a bit
y1 = y1[::sample_rate]
y2 = y2[::sample_rate]
errors1 = errors1[::sample_rate]
errors2 = errors2[::sample_rate]
x1= x2 =data_range[::sample_rate]

truth1 = [2* (1. + np.sin(2. * np. pi * x/period)) for x in data_range]
truth2 = [1000 * (gauss_skew(x,-10,loc=location ,scale = 50)) for x in data_range]

truths = [truth1,truth2]
x= [x1,x2]
y=[y1,y2]
errors = [errors1,errors2]

fig,ax  = plt.subplots(nrows=2,ncols=2)

ax[0][0].errorbar(x1,y1,yerr=errors1,fmt='o')
ax[0][0].set_xlabel('time')
ax[1][0].errorbar(x2,y2,yerr=errors2,fmt='o')
ax[1][0].set_xlabel('time')
ax[0][1].set_xlabel('frequency')
ax[1][1].set_xlabel('frequency')


for i in range(0,2):
    frequency, power = LombScargle(x[i],y[i],errors[i]).autopower()
    #Get the best fit frequency as in the tutorial
    best_frequency = frequency[np.argmax(power)]
    print('best frequency:',best_frequency)
    t_fit = np.linspace(x[i][0], math.floor(x[i][-1]),num =50)

    #Fit the best fit frequency
    #plot the best best model based on the best fit
    y_fit = LombScargle(x[i], y[i], errors[i]).model(t_fit, best_frequency)
    ax[i][0].plot(t_fit,y_fit,'g')
    ax[i][0].plot(data_range,truths[i],'r')

    #Plot the PSD
    ax[i][1].plot(frequency,power)
    ax[i][1].axvline(x=best_frequency,color='black', ls='--')
plt.show()