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

发布于:2025-06-04 ⋅ 阅读:(22) ⋅ 点赞:(0)

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 = N − noverlap nperseg − noverlap K = \frac{N - \text{noverlap}}{\text{nperseg} - \text{noverlap}} K=npersegnoverlapNnoverlap
  2. 处理每段

    for segment in segments:
        # 1. 去趋势 (detrend)
        # 2. 加窗 (apply window)
        # 3. 计算FFT
        # 4. 计算周期图: |FFT|²
    
  3. 平均:对所有段的周期图进行平均
    P x x = avg ( ∣ F F T k ∣ 2 ) f s ⋅ S 2 P_{xx} = \frac{\text{avg}(|FFT_k|^2)}{f_s \cdot S_2} Pxx=fsS2avg(FFTk2)
    其中 S 2 = ∑ n = 0 N − 1 w [ n ] 2 S_2 = \sum_{n=0}^{N-1} w[n]^2 S2=n=0N1w[n]2 是窗函数的能量补偿因子

五、关键特性

  1. 方差减小:通过分段平均降低随机噪声影响
  2. 频率分辨率 Δ f = f s nperseg \Delta f = \frac{f_s}{\text{nperseg}} Δf=npersegfs
    • 增加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 方法特别适合分析实际工程信号,在噪声环境下能提供更稳定的频谱估计结果。


研究学习不易,点赞易。
工作生活不易,收藏易,点收藏不迷茫 :)



网站公告

今日签到

点亮在社区的每一天
去签到