[Python] python信号处理绘制信号频谱

python信号处理绘制信号频谱：scipy.signal.welch

scipy.signal.welch 是 Python 中用于计算信号功率谱密度（PSD）的常用函数，采用 Welch 方法实现。这种方法通过将信号分段、加窗和平均来减少频谱估计的方差。以下是详细说明：

一、函数介绍

f, Pxx = scipy.signal.welch(
    x,                   # 输入信号
    fs=1.0,              # 采样频率
    window='hann',       # 窗函数
    nperseg=None,        # 每段长度
    noverlap=None,       # 重叠点数
    nfft=None,           # FFT长度
    detrend='constant',  # 去趋势方法
    return_onesided=True,# 是否返回单边谱
    scaling='density',   # 缩放类型
    axis=-1,             # 计算轴
    average='mean'       # 平均方法
)

二、核心参数详解

参数	默认值	说明
x	-	输入信号（1D数组或2D数组）
fs	1.0	采样频率（Hz），影响频率坐标
window	‘hann’	窗函数类型：‘hann’（汉宁窗）、‘hamming’（汉明窗）、‘blackman’（布莱克曼窗）、‘boxcar’（矩形窗）等
nperseg	None	每段长度（点数）。None时取256或len(x)较小值
noverlap	None	段间重叠点数。None时为nperseg // 2（50%重叠）
nfft	None	FFT长度。None时等于nperseg；大于nperseg时进行零填充
detrend	‘constant’	去趋势方法：‘constant’（去均值）、‘linear’（去线性趋势）、False（不去趋势）
scaling	‘density’	输出缩放：‘density’（PSD, V²/Hz）、‘spectrum’（功率谱, V²）
average	‘mean’	平均方法：‘mean’（算术平均）、‘median’（中位数平均）

三、返回值

返回值	说明
f	频率数组（Hz）
Pxx	功率谱密度数组（单位取决于scaling参数）

四、算法原理

1. 分段：将长度为N的信号分成K段

   • 每段长度：nperseg
   • 段数：$K=\frac{N-\text{noverlap}}{\text{nperseg}-\text{noverlap}}$

2. 处理每段：

   for segment in segments:
       # 1. 去趋势 (detrend)
       # 2. 加窗 (apply window)
       # 3. 计算FFT
       # 4. 计算周期图: |FFT|²
   

3. 平均：对所有段的周期图进行平均
   $$P_{xx}=\frac{\text{avg}(|FF{T}_k{|}^2)}{f_s\cdot S_2}$$
   其中 $S_2={\sum }_{n=0}^{N-1}w[n{]}^2$ 是窗函数的能量补偿因子

五、关键特性

1. 方差减小：通过分段平均降低随机噪声影响
2. 频率分辨率：$\Delta f=\frac{f_s}{\text{nperseg}}$
   • 增加nperseg ➔ 提高分辨率
   • 增加noverlap ➔ 增加段数（降低方差）
3. 偏置-方差权衡：
   • 长分段：低偏置、高方差
   • 短分段：高偏置、低方差

六、完整示例

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

# 生成测试信号
fs = 1000  # 采样率
t = np.arange(0, 1, 1/fs)
x = np.sin(2*np.pi*50*t) + 0.5*np.sin(2*np.pi*120*t)  # 50Hz + 120Hz
x += 0.8 * np.random.randn(len(t))  # 添加高斯噪声

# 计算功率谱密度
f, Pxx = signal.welch(
    x,
    fs=fs,
    window='hann',      # 汉宁窗（默认）
    nperseg=1024,       # 每段1024点
    noverlap=512,       # 50%重叠
    nfft=2048,          # 2048点FFT（零填充）
    detrend='constant', # 去除均值
    scaling='density'   # 功率谱密度
)

# 绘制结果
plt.figure(figsize=(10, 5))
plt.semilogy(f, Pxx)  # 对数坐标
plt.title('Power Spectral Density (Welch Method)')
plt.xlabel('Frequency [Hz]')
plt.ylabel('PSD [V²/Hz]')
plt.xlim(0, 200)  # 限制显示范围
plt.grid(True)
plt.show()

绘制得到频谱图如下：

七、应用场景推荐配置

场景	推荐参数
平稳信号	nperseg=1024, noverlap=512
低信噪比信号	增加noverlap（如75%）
高分辨率需求	增大nperseg（牺牲计算速度）
瞬态信号分析	减小nperseg，使用detrend='linear'
精确峰值测量	使用window='flattop'（平坦窗）

八、常见问题解决

1. 频率泄露：

   • 改用旁瓣衰减更好的窗（如window='blackman'）
   • 增加nperseg提高分辨率

2. 弱信号被淹没：

   # 使用中位数平均增强弱信号可见性
   f, Pxx = welch(..., average='median')
   

3. 直流偏移影响：

   # 确保去趋势
   f, Pxx = welch(..., detrend='constant')
   

4. 频率坐标错误：

   • 检查fs是否设置正确
   • 确认return_onesided=True（实信号）

九、与FFT方法的对比

特性	welch	直接FFT
方差	低	高
计算量	中高	低
频率分辨率	可调节	固定
抗噪能力	强	弱
适合信号类型	长时记录、噪声信号	短信号、瞬态信号

Welch 方法特别适合分析实际工程信号，在噪声环境下能提供更稳定的频谱估计结果。

---

研究学习不易，点赞易。
工作生活不易，收藏易，点收藏不迷茫 ：）

---