python - 似乎无法在 scipy 中正确使用 welch 函数

Question

我正在尝试通过使用Python的scipy包来学习DSP。我已经测量了来自机器的一些 200Hz 信号。现在我想检查信号的频谱密度。这是绘制出来的信号:

正如您所看到的，该信号显然不是低频信号。但是，我仍然得到如下所示的 PSD 图。这是我一直在使用的代码:

def plot_and_save(sample,i=""):
filepath = sample + "/" + "merged_data.csv"
df = pd.read_csv(filepath)

sns.set_style("darkgrid")
sns.set_colodef plot_and_save(sample,i=""):r_codes("dark")
fig1, axes = plt.subplots(3,2)
fig1.figsize = [1920,1080]

axes[0,0].plot(df["Time"], df["Voltage"], marker="",color="blue")
axes[0,0].set_title("Plot of Voltage against time")
axes[0,0].set_xlabel("Time (us)", ha= "right", x=1.0)
axes[0,0].set_ylabel("Voltage (V)")

axes[0,1].plot(df["Time"], df["Current"], marker="",color="red")
axes[0,1].set_title("Plot of Current against time")
axes[0,1].set_xlabel("Time (us)", ha= "right",x=1.0)
axes[0,1].set_ylabel("Current (A)")

sampling_freq = 1e6
vol_f, vol_psd = welch(df["Voltage"].to_numpy(), fs=sampling_freq)
axes[1,0].plot(vol_f, vol_psd, marker="",color="green")
axes[1,0].set_title("Plot of Voltage's Power Spectral Density")
axes[1,0].set_xlabel("Frequency", ha= "right",x=1.0)
axfilepath = sample + "/" + "merged_data.csv"
df = pd.read_csv(filepath)

sns.set_style("darkgrid")
sns.set_color_codes("dark")
fig1, axes = plt.subplots(3,2)
fig1.figsize = [1920,1080]

axes[0,0].plot(df["Time"], df["Voltage"], marker="",color="blue")
axes[0,0].set_title("Plot of Voltage against time")
axes[0,0].set_xlabel("Time (us)", ha= "right", x=1.0)
axes[0,0].set_ylabel("Voltage (V)")

axes[0,1].plot(df["Time"], df["Current"], marker="",color="red")
axes[0,1].set_title("Plot of Current against time")
axes[0,1].set_xlabel("Time (us)", ha= "right",x=1.0)
axes[0,1].set_ylabel("Current (A)")

sampling_freq = 1e6
vol_f, vol_psd = welch(df["Voltage"].to_numpy(), fs=sampling_freq)
axes[1,0].plot(vol_f, vol_psd, marker="",color="green")
axes[1,0].set_title("Plot of Voltage's Power Spectral Density")
axes[1,0].set_xlabel("Frequency", ha= "right",x=1.0)
axes[1,0].set_ylabel("PSD")

cur_f, cur_psd = welch(df["Current"], fs=sampling_freq)
axes[1,1].plot(cur_f, cur_psd, marker="",color="yellow")
axes[1,1].set_title("Plot of Current's Power Spectral Density")
axes[1,1].set_xlabel("Frequency", ha= "right",x=1.0)
axes[1,1].set_ylabel("PSD")

axes[2,0].plot(df["Time"], df["Current"]*df["Voltage"], marker="",color="red")def plot_and_save(sample,i=""):
axes[2,0].set_title("Plot of power against time")
axes[2,0].set_xlabel("Time (us)", ha= "right",x=1.0)
axes[2,0].set_ylabel("Watt")

sns.scatterplot(x="Current", y="Voltage", data=df, ax=axes[2,1], color = 'r', edgecolor="None", linewidth=0)
axes[2,1].set_title("Plot of Voltage against Current")
axes[2,1].set_xlabel("Current (A)", ha="right", x=1.0)
axes[2,1].set_ylabel("Voltage (V)")

fig1.suptitle("Sample number: " + str(i+5),fontsize=16)
# plt.savefig("formatted/"+date+"/figure/"+str(i+5)def plot_and_save(sample,i=""):".svg",format="svg")
plt.show()es[1,0].set_ylabel("PSD")

cur_f, cur_psd = welch(df["Current"], fs=sampling_freq)
axes[1,1].plot(cur_f, cur_psd, marker="",color="yellow")
axes[1,1].set_title("Plot of Current's Power Spectral Density")
axes[1,1].set_xlabel("Frequency", ha= "right",x=1.0)
axes[1,1].set_ylabel("PSD")

axes[2,0].plot(df["Time"], df["Current"]*df["Voltage"], marker="",color="red")
axes[2,0].set_title("Plot of power against time")
axes[2,0].set_xlabel("Time (us)", ha= "right",x=1.0)
axes[2,0].set_ylabel("Watt")

sns.scatterplot(x="Current", y="Voltage", data=df, ax=axes[2,1], color = 'r', edgecolor="None", linewidth=0)
axes[2,1].set_title("Plot of Voltage against Current")
axes[2,1].set_xlabel("Current (A)", ha="right", x=1.0)
axes[2,1].set_ylabel("Voltage (V)")

fig1.suptitle("Sample number: " + str(i+5),fontsize=16)
# plt.savefig("formatted/"+date+"/figure/"+str(i+5)".svg",format="svg")
plt.show()

非常感谢您的帮助!

Answer 1

这是计算 PSD 的最小工作示例。我认为你的问题可能是你将 nperseg 选项保留为默认值，根据 documentation 为 256 点。您想要做的是在每个片段中包含足够的样本，以确保包含信号的主要特征。我的粗略指导至少是周期的 1.5-2.0 倍。

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import welch

# Generate signal
freq = 1.0 # Hz
t = np.linspace(0, 10, 1000)
y = np.sin(2.*np.pi*freq*t) + 0.2*np.random.random(len(t))

# Compute no. samples per segment
nSegments = 3
overlap = 0.5

nPerSeg = np.round(len(y)//nSegments / overlap)
if nSegments == 1:
    nPerSeg = len(y)
nOverlap = np.round(overlap * nPerSeg)

# Compute sampling rate.
dt = t[1] - t[0]

# Compute PSD
f, psd = welch(y, fs=1./dt,
    window="hamm", nperseg=nPerSeg, noverlap=nOverlap, nfft=None, detrend="constant",
    return_onesided=True, scaling="density")

# Plot
plt.figure()
plt.plot(t, y)

plt.figure()
plt.plot(f/freq, psd)
plt.xscale("log")
plt.yscale("log")
plt.grid()
plt.xlabel("f/freq")
plt.ylabel("PSD of y")

plt.show()

预期结果显示主信号频率处的峰值: