功率谱、Welch 方法和 window 函数

Power spectrum, Welch method and window functions

我需要使用 Welch 方法查找信号的功率谱。然后使用切比雪夫、凯撒和高斯windows,找到频率分离方面最好的window。

我编写了一个函数来产生带噪声的正弦信号并对这个信号使用 Welch 方法,还生成了所有这三个 window 函数。

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

def sig_noise(f,snr,n):
    fs=1000
    dt=1/fs
    tk=n*dt
    A=1
    B=10**(-snr/20)
    t=np.arange(0,tk,dt)
    y=(np.random.rand(n)-0.5)
    z=A*np.sin(2*np.pi*f*t)+B*y
    return z,t

[s,t]=sig_noise(140,-10,64)

f, Pxx_den = sg.welch(s,nperseg=64,noverlap=None) # Welch method, parameters are taken from book

plt.plot(f,Pxx_den)
plt.xlabel("Frequency")
plt.grid()

plt.figure(figsize=(10,10))

chebyshev=sg.chebwin=(64,40) # Chebyshev window

kaiser=sg.kaiser(64,0) # Kaiser window

gauss=sg.gaussian(64,np.std((f,Pxx_den),ddof=0)) # Gauss window

# plots for check

plt.subplot(221)
plt.grid()
plt.title("Gauss")
plt.plot(gauss)
plt.subplot(222)
plt.grid()
plt.title("Chebyshev")
plt.plot(chebyshev)
plt.subplot(223)
plt.grid()
plt.title("Kaiser")
plt.plot(kaiser)

f, Pxx_den = sg.welch(s,nperseg=64,noverlap=None,window=chebyshev)

最后一行代码不起作用,因为我想使用例如 chebyshev window。未知 window 类型是我得到的答案。谁能告诉我哪里错了?

问题在于以下行中的额外 = 字符:

chebyshev=sg.chebwin=(64,40) # Chebyshev window
#                  ^^^
#                  extra "=" character

这导致 chebyshev 等于元组 (64,40) 而不是预期的数组(遵循与其他 windows 类似的方法)。稍后当给定元组作为 windows 参数时,scipy.signal.welch tries to construct the window by calling scipy.signal.get_window 用于由元组的第一个元素指定的 window 类型。显然 64 (您的情况下元组的第一个元素)不是有效的 window 类型。

因此,解决方法是简单地更正拼写错误:

chebyshev=sg.chebwin(64,40) # Chebyshev window

或者使用适当构造的元组调用 welch

f, Pxx_den = sg.welch(s,nperseg=64,noverlap=None,window=('chebwin',40)