功率谱、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)
我需要使用 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)