Bartlett周期图的Python实现

Question

我试图根据Bartlett方法的描述在Python中实现周期图，并将结果与​​Scipy的结果进行比较，通过设置overlap=0，使用window='boxcar'（矩形窗口）。然而，我的结果存在一定的偏差。有人能指出我的代码有什么问题吗？谢谢

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal


def my_bartlett_periodogram(x, fs, nperseg, nfft):    
    nsegments = len(x) // nperseg
    psd = np.zeros(nfft)
    for segment in x.reshape(nsegments, nperseg):
        psd += np.abs(np.fft.fft(segment))**2 / nfft
    psd[0] = 0   # important!!
    psd /= nsegments
    psd = psd[0 : nfft//2]
    freq = np.linspace(0, fs/2, nfft//2)
    return freq, psd

def plot_output(t, x, f1, psd1, f2, psd2):
    fig, axs = plt.subplots(3,1, figsize=(12,15))
    axs[0].plot(t[:300], x[:300])
    axs[1].plot(freq1, psd1)
    axs[2].plot(freq2, psd2)
    axs[0].set_title('Input (len=8192, fs=512)')
    axs[1].set_title('Bartlett Periodogram (nfft=512, zero-overlap, no-window)')
    axs[2].set_title('Scipy Periodogram (nfft=512, zero-overlap, no-window)')
    axs[0].set_xticks([])
    axs[2].set_xlabel('Freq (Hz)')
    plt.show()

# Run
fs = nfft = nperseg = 512
t = np.arange(8192) / fs
x = np.sin(2*np.pi*50*t) + np.sin(2*np.pi*100*t) + np.sin(2*np.pi*150*t)
freq1, psd1 = my_bartlett_periodogram(x, fs, nperseg, nfft)
freq2, psd2 = signal.welch(x, fs, nperseg=nperseg, nfft=nfft, window='boxcar', noverlap=0)
plot_output(t, x, freq1, psd1, freq2, psd2)

Answer 1

长话短说：

\n\n

代码没有任何问题。但welch返回功率谱密度，它是功率谱fs乘以 2 来补偿削减一半的频谱。

\n\n

为了补偿，psd2 * fs / 2应该与 非常相似psd。

\n\n\n\n

根据维基百科，计算psd似乎是正确的：

\n\n

\n

\n
1. 将原来的N点数据段分割成K个（不重叠）数据段，每个数据段长度为M
\n
2. 对于每个段，通过计算离散傅立叶变换（不除以 M 的 DFT 版本）来计算周期图，然后计算结果的平方幅值并将其除以 M。
\n
3. 对 K 个数据段的上述周期图结果进行平均。
\n

\n

\n\n

那么我们应该更信任谁，维基百科还是 scipy？我倾向于后者，但我们可以自己找出答案。根据Parseval 定理，平方信号的积分应与平方 FFT 幅度的积分相同。由于周期图是从平方 FFT 获得的，因此定理应该近似成立。

\n\n

print(np.mean(y**2))  # 1.499727698431174\nprint(np.mean(psd))  # (1.4999999999999991+0j)\nprint(np.mean(psd2))  # 0.0058365758754863788\n

\n\n

对于 来说这已经足够接近了psd，所以我们假设它是正确的。但我拒绝相信 scipy 应该如此公然错误！让我们仔细看看文档，看看他们对这个论点有什么看法scaling（强调我的）：

\n\n

\n

选择计算功率谱密度 (\xe2\x80\x98密度\xe2\x80\x99)（其中 Pxx 的单位为 ）V**2/Hz和计算功率谱（\xe2\x80\x98spectrum\xe2\x80\x99）（其中 Pxx 的单位为 ） V**2，如果 x 以 V 为单位测量，fs 以 Hz 为单位测量。默认为 \xe2\x80\x98密度\xe2\x80\x99

\n

\n\n

嗯！welch\ 的结果是功率谱密度，这意味着它的单位为每赫兹功率。然而，我们将其与信号功率进行了比较。如果我们乘以psd2采样率以去掉 1/Hz 单位，则与 相同psd。好吧，除了系数 2。这个系数是为了补偿削减一半的频谱。如果我们设置return_onesided=False为获得完整的频谱，那么该因素就消失了。

\n