Scipy周期图术语混淆答案

【问题标题】：Scipy periodogram terminology confusionScipy周期图术语混淆
【发布时间】：2014-04-15 19:32:04
【问题描述】：

我对 scipy.signal.periodogram 中使用的术语感到困惑，即：

缩放：{“密度”，“光谱”}，可选在计算功率谱密度（“密度”）之间进行选择其中Pxx 的单位为 V*2/Hz，如果 x 以 V 为单位测量并计算如果x 是，则Pxx 的单位为V*2 的功率谱（'spectrum'）以 V 为单位测量。默认为“密度”

（见：http://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.periodogram.html）

1) 一些测试表明，选项“密度”的结果取决于信号和窗口长度以及采样频率（随着信号长度的增加而增加）。怎么会？我会说不应该依赖于这些东西的正是密度。如果我采用更长的信号，我应该得到更准确的估计，而不是不同的结果。更不用说对窗口长度的依赖也很令人惊讶。结果在无限信号的极限处发散，这可能是能量的特征，但不是功率的特征。当长度增加时，周期图不应该收敛到实际的理论 PSD 吗？如果，那么，我应该在 signal.periodogram 方法之外执行另一个归一化吗？

2) 相反，我看到替代选项“频谱”给出了我之前所说的功率谱密度，也就是说，它给出的结果与窗口段和窗口长度无关，并且与理论计算一致。例如，对于 Asin(2PIft)，两侧的解决方案在 -f 和 f 处产生两个峰，每个峰的高度为 0.25*A^2。

关于这个主题有很多文献，但我的印象是也有很多不兼容的术语，所以我会感谢任何澄清。直截了当的问题是如何解释这些选项及其单位。（我习惯于看到标有“功率谱密度”的 V^2/Hz）。

【问题讨论】：

标签： scipy signal-processing

【解决方案1】：

让我们使用一个名为 data 的实数数组，它的长度为 N，采样频率为 fs。我们称时间 bin dt=1/fs，并且 T = N * dt。在频域中，频域df = 1/T = fs/N。

功率谱 PS（scaling='spectrum' in scipy.periodogram）计算如下：

import numpy as np
import scipy.fft as fft
dft = fft.fft(data)
PS = np.abs(dft)**2 / N ** 2

它的单位是 V^2。可以这样理解。类推连续傅里叶变换，信号的能量E为：

E := np.sum(data**2) * dt = 1/N * np.sum(np.abs(dft)**2) * dt

（根据 Parseval 定理）。信号的功率P是总能量E除以信号持续时间T：

P := E/T = 1/N**2 * np.sum(np.abs(dft)**2)

功率 P 仅取决于离散傅里叶变换 (DFT) 和样本数 N。不直接取决于采样频率 fs 或信号持续时间 T。因此，每个频率通道的功率，即功率谱 SP，为由上式给出：

PS = np.abs(dft)**2 / N ** 2

对于功率谱密度PSD（scaling='density' in scipy.periodogram），需要将PS除以DFT的频率bin，df：

PSD := PS/df = PS * N * dt = PS * N / fs

因此：

PSD = np.abs(dft)**2 / N * dt

单位为 V^2/Hz = V^2 * s，现在取决于采样频率。这样，在频率范围内对 PSD 进行积分得到的结果与对 PS 的各个值求和相同。

这应该解释您在更改窗口、采样频率、持续时间时看到的关系。

【讨论】：

【解决方案2】：

scipy.signal.peridogram 使用 scipy.signal.welch 函数，重叠为 0。因此，缩放类似于韦尔奇函数、密度或谱提供的缩放。
在密度缩放的情况下，幅度将随窗口长度而变化，因为窗口越长，频率分辨率越高，例如\Delta_f 较小。由于估计的密度是平均密度，\Delta_f 越小，在平均时考虑的零能量就越少。
正如您所提到的，频谱缩放是对频谱上的能量密度进行积分以产生能量。因此，零值上的积分不会影响最终值。

【讨论】：

【解决方案3】：

傅立叶变换实际上需要无限时间序列中的有限能量（如衰减）。因此，如果您只是通过“复制”使您的时间序列样本更长，那么能量将是无限的，持续时间是无限的。
我的主要困惑在于 scipy.signal.periodogram 的“频谱”选项，即使时间序列变长，它似乎也会创建一个恒定的能谱。通常，0.5*A^2=S(f)*delta_f，其中 S(f) 是功率密度谱。 S(f)*delta_f，如果A是常数，则表示能量是常数。但是当使用更长的时间序列时，delta_f（即增量频率）会根据 FFT 过程相应减少。例如，100s 时间序列将导致 delta_f=0.01Hz，而 1000s 时间序列将导致 delta_f=0.001Hz。代表密度的 S(f) 会相应变化。

【讨论】：