巴特利特周期图的Python实现

Question

巴特利特周期图的Python实现

7

我正在尝试基于Bartlett's method的描述，在Python中实现周期图，并通过将重叠设置为0，使用窗口='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)

- Scoodood

我刚刚发现，我的计算中的直流分量必须设置为0，否则当输入中有直流偏移时，输出将不匹配。 - Scoodood

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- MB-F · Accepted Answer

简而言之：

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

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

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

原始的N点数据段被分成K个（不重叠的）数据段，每个长度为M。

对于每个数据段，通过计算离散傅里叶变换（DFT版本，不除以M），然后计算结果的平方幅值并将其除以M来计算出周期图。

对K个数据段的周期图结果进行平均。

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

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

对于psd来说，这已经足够接近了，因此让我们假设它是正确的。但我拒绝相信scipy应该如此明显地错误！让我们更仔细地查看文档，看看他们对scaling参数有何说明（强调我的）：

在计算功率谱密度（'density'）和功率谱（'spectrum'）之间进行选择，其中Pxx的单位分别为V**2/Hz和V**2，如果x以V为单位，fs以Hz为单位。默认为“density”。

嗯！welch的结果是功率谱密度，意味着其单位为每赫兹的功率。然而，我们将其与信号功率进行了比较。如果我们将psd2乘以采样率以消除1/Hz单位，则与psd相同。好吧，除了一个因子2以外。这个因子旨在补偿截去一半频谱的影响。如果我们设置return_onesided=False以获取完整频谱，则该因子被消除。