Bartlett周期图的Python实现

我试图根据 Bartlett's method 的描述在 Python 中实现周期图，并将结果与 Scipy 的结果进行比较，通过设置重叠=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)

最佳答案

特尔；博士:

代码没有任何问题。但是 welch 返回功率谱密度，即功率谱乘以 fs 并且它通过乘以 2 来补偿切除一半的频谱。

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

根据 Wikipedia psd 的计算似乎是正确的:

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

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

这对于 psd 来说已经足够接近了，所以让我们假设它是正确的。但我拒绝相信 scipy 应该如此公然错误!让我们仔细看看 documentation，看看他们对 scaling 参数有什么看法(重点是我的):

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

关于Bartlett周期图的Python实现，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/46882312/

PSD

Bartlett周期图的Python实现