Bartlett周期图的Python实现

Question

我试图根据来自巴特利特法的描述在Python中实现周期图，并通过设置overlap=0，使用window='boxcar‘(矩形窗口)，将结果与Scipy的结果进行比较。然而，我的结果是由于某种规模因素的影响。有人能指出我的代码出了什么问题吗？谢谢

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

TL;DR：

密码没什么问题。但是welch返回功率谱密度，即功率谱乘以fs，它通过与2相乘来补偿一半的频谱。

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

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

1. 原始的N点数据段被分割成K(不重叠)数据段，每个数据段的长度为M。
2. 对于每一段，通过计算离散傅里叶变换(不除以M的DFT版本)来计算周期图，然后计算结果的平方大小并除以M。
3. 对K数据段的上述周期图的平均结果。

那么，我们应该更信任谁呢，维基百科还是枕木？我倾向于后者，但我们可以自己找出答案。根据Parseval定理，平方信号的积分应与平方FFT幅值上的积分相同。由于周期图是由平方FFT得到的，所以定理应该大致成立。

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

对于psd来说，这已经足够接近了，所以让我们假设它是正确的。但我不相信这种行为是如此明目张胆！让我们仔细看看文档，看看他们对scaling论点有什么看法(强调我的论点)：

在计算Pxx具有‘V2/ Hz’**单位的功率谱密度(“密度”)和计算Pxx具有V**2单位的功率谱(‘频谱’)之间选择，如果x以V测量，fs以赫兹测量。默认为“密度”

嗯哼！welch的结果是功率谱密度，这意味着它每赫兹有一个功率单位。然而，我们将它与信号功率进行比较。如果我们用采样率乘以psd2来去除1/Hz单位，那么它和psd是一样的。嗯，除了一个因子2。这个因子是用来补偿把一半的频谱减少的。如果我们设置return_onesided=False来得到全谱，那么这个因子就没有了。