如何匹配同一数据集的 lomb-scargle 和 FFT 图？

Question

我正在做一些工作，将一段时间内某些气体浓度的插值 fft 与相同数据的 lomb-scargle 周期图进行比较，其中采样不均匀。我正在使用 scipy 的 fft 函数来计算傅立叶变换，然后对其模数求平方，得到我认为的功率谱密度，单位为十亿分之一（ppb）平方。

我可以让 lomb-scargle 图与 FFT 几乎完全匹配，但幅度不同，FFT 功率谱密度总是更高，即使我认为 lomb-scargle 功率是功率谱密度.现在我使用的 lomb 代码：http://www.astropython.org/snippet/2010/9/Fast-Lomb-Scargle-algorithm，对数据集进行归一化，去掉平均值并除以数据方差的 2 倍，因此我以相同的方式对 FFT 数据进行了归一化，但幅度仍然不匹配。

因此我做了更多的研究，发现归一化的 lomb-scargle 幂可以没有单位，因此我无法匹配。这引出了两个问题：

1. 归一化 lim-scargle perioogram 的功率谱密度是什么单位（如果有）？

2. 如何在幅度和模式方面将我的 fft 图与我的 lomb-scargle 图进行匹配？

谢谢。

Answer 1

自从提出和回答这个问题以来，AstroPy 项目已经获得了 Lomb-Scargle 方法，这个问题在文档中得到了解决：http://docs.astropy.org/en/stable/stats/lombscargle.html#psd-normalization-unnormalized

简而言之，您可以计算傅立叶周期图并将其与 astropy Lomb-Scargle 周期图进行比较，如下所示

import numpy as np
from astropy.stats import LombScargle

def fourier_periodogram(t, y):
    N = len(t)
    frequency = np.fft.fftfreq(N, t[1] - t[0])
    y_fft = np.fft.fft(y)
    positive = (frequency > 0)
    return frequency[positive], (1. / N) * abs(y_fft[positive]) ** 2

t = np.arange(100)
y = np.random.randn(100)
frequency, PSD_fourier = fourier_periodogram(t, y)
PSD_LS = LombScargle(t, y).power(frequency, normalization='psd')

np.allclose(PSD_fourier, PSD_LS)
# True

由于 AstroPy 是天文学中常用的工具，我认为这可能比基于上述代码 sn-p 的答案更有用。

Answer 2

一个级数的傅里叶变换的平方模被定义为能谱密度 (ESD)。您需要将 ESD 除以系列的长度，以转换为 功率谱密度 (PSD) 的估计值。

单位

PSD 的单位是 [units]**2/[frequency]，其中 [units] 代表原始系列的单位。

标准化

为了检查正确的归一化，可以对白噪声的 PSD（具有已知方差）进行数值积分。如果积分谱等于序列的方差，则归一化是正确的。不过，因子 2（太低）并不是不正确的，它可能表明 PSD 被归一化为 双面；在这种情况下，只需乘以 2，您就有了一个正确标准化的单面 PSD。

使用 numpy，randn 函数生成高斯分布的伪随机数。例如

10 * np.random.randn(1, 100)

生成一个 1×100 数组，均值 = 0，方差 = 100。如果采样频率是 1-Hz，则 单面 PSD 理论上将在 200 个单位**2/Hz 处平坦，从 [0,0.5] Hz 开始；因此，积分谱将为 10，等于系列的方差。

更新

我修改了您链接的 python 代码中包含的示例，以演示长度为 20、方差为 1、采样频率为 10 的正态分布系列的标准化：

import numpy
import lomb
numpy.random.seed(999)
nd = 20
fs = 10
x = numpy.arange(nd)
y = numpy.random.randn(nd)
fx, fy, nout, jmax, prob = lomb.fasper(x, y, 1., fs)
fNy = fx[-1]
fy = fy/fs
Si = numpy.mean(fy)*fNy
print fNy, Si, Si*2

这给了我：

5.26315789474 0.482185882163 0.964371764327

它向您展示了一些东西：

1. 要求的“奈奎斯特”频率实际上是采样频率。
2. 结果需要除以采样频率。
3. 输出是针对双面 PSD 进行归一化的，因此乘以 2 会使积分频谱接近 1。