2

我正在做一些工作，将一段时间内某些气体浓度的插值 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 图相匹配？

谢谢你。

4

2 回答 2

6

一系列傅里叶变换的平方模被定义为能谱密度(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。
于 2013-02-18T19:52:19.227 回答
2

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

简而言之，您可以计算傅立叶周期图并将其与 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 是天文学中常用的工具，我认为这可能比基于上述代码片段的答案更有用。

于 2017-03-09T19:36:29.677 回答