如何匹配同一数据集的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，在处理数据时会先减去平均值，然后再除以数据的方差的两倍，因此我也以同样的方式对FFT数据进行了归一化，但幅度仍然不匹配。

因此，我做了一些额外的研究，发现归一化的Lomb-Scargle功率可能是无单位的，所以我无法让这两个图匹配。这让我有了两个问题：

1. 归一化的Lomb-Scargle周期图的功率谱密度有什么单位（如果有的话）？

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是天文学中常用的工具，所以我觉得这可能比上面提到的代码片段的答案更有用。

Answer 2

傅里叶变换的平方模被称为能量谱密度（ESD）。要把能量谱密度转换为功率谱密度（PSD），你需要把ESD除以序列的长度。

单位

功率谱密度的单位是[单位]²/[频率]，其中[单位]代表你原始序列的单位。

归一化

为了检查归一化是否正确，可以对已知方差的白噪声的PSD进行数值积分。如果积分后的谱等于序列的方差，那么归一化就是正确的。不过，如果结果是2（太低了），也不是错误，这可能表示PSD是归一化为双侧的；在这种情况下，只需乘以2，就得到了正确归一化的单侧PSD。

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

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

这会生成一个1行100列的数组，均值为0，方差为100。如果采样频率是1赫兹，那么理论上，单侧PSD在[0,0.5]赫兹范围内会是200单位²/赫兹；因此，积分后的谱将是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。