SciPy 信号处理解析：滤波、时频分析与噪声去除

在科学计算、工程信号、语音识别与通信等领域，信号处理 是理解和分析时间序列数据的核心方法。Python 的 SciPy 提供了功能强大的信号处理模块 scipy.signal，能够实现卷积、滤波、时频分析、系统建模等多种操作，为科研与工程应用提供高效工具。

---

1. 信号处理概览与 SciPy 模块结构

信号处理旨在从观测数据中提取有用信息，包括噪声抑制、特征提取以及系统响应分析。SciPy 的信号处理模块主要功能包括：

功能分类	核心函数	应用方向
卷积与相关性	convolve, correlate	信号平滑、模板匹配
滤波器设计	butter, firwin, cheby1	去噪、带通滤波
频率分析	freqz, periodogram, welch	频谱特征提取
时频变换	stft, cwt	瞬态信号与非平稳信号分析
系统分析	lfilter, lti	差分方程系统模拟

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

# 设置Seaborn绘图主题
sns.set_theme(style="whitegrid", font="SimHei", rc={"axes.unicode_minus": False})

---

2. 卷积与相关性：信号平滑与匹配

卷积用于描述系统对输入信号的响应，相关性用于信号匹配与相似性分析，其数学定义为：
$$y\left[n\right]=\left(x\ast h\right)\left[n\right]=\sum _{k=-\infty }^{+\infty }x\left[k\right]\ast h\left[n-k\right]$$

$$r_{xy}\left[n\right]=\sum _{k=-\infty }^{+\infty }x\left[k\right],y\left[k+n\right]$$

通过卷积实现信号平滑，可以有效抑制高频噪声：

# 参数设置与信号构造
np.random.seed(42)
t = np.linspace(0, 2, 400)

# 原始信号：低频 + 高频 + 噪声
x = np.sin(2 * np.pi * 3 * t) + 0.5 * np.sin(2 * np.pi * 20 * t) + 0.2 * np.random.randn(len(t))

# 卷积平滑核
window_size = 20
h = np.ones(window_size) / window_size

# 卷积平滑
y_conv = signal.convolve(x, h, mode='same')

# 可视化卷积效果
plt.figure(figsize=(10, 5))

# 原始信号
plt.plot(t, x, color='steelblue', label='原始信号（含高频噪声）')

# 卷积平滑信号
plt.plot(t, y_conv, color='darkorange', linewidth=2.2, label=f'卷积平滑信号 (窗口={window_size})')

# 标题与标签
plt.title('卷积平滑效果对比：高频成分被抑制', fontsize=14, weight='bold')
plt.xlabel('时间 t')
plt.ylabel('幅度')

# 可视化高频抑制效果
plt.xlim(0, 1.5)  # 局部高频区域
plt.ylim(np.min(x[:150]-0.2), np.max(x[:150])+0.2)

# 网格与图例
plt.grid(True, linestyle='--', alpha=0.8)
plt.legend(loc='lower right', fontsize=11)

plt.tight_layout()
plt.show()

说明： 卷积平滑显著降低信号的高频噪声，是基础的降噪方法。

---

3. 差分方程与线性系统分析

线性时不变系统（LTI）可用差分方程表示：
$$y[n]=\sum _{i=0}^M{b}_i,x[n-i]-\sum _{j=1}^N{a}_j,y[n-j]$$

$$H(z)=\frac{Y(z)}{X(z)}=\frac{b_0+b_1{z}^{-1}+\cdots +b_M{z}^{-M}}{1+a_1{z}^{-1}+\cdots +a_N{z}^{-N}}$$

SciPy 的 lfilter 可以模拟系统响应：

# 设置随机种子保证可重复性
np.random.seed(42)

# 定义系统系数（差分方程 y[n] = 0.1x[n] + 0.1x[n-1] + 0.1x[n-2] - 0.9y[n-1]）
b = [0.1, 0.1, 0.1]   # 分子系数（输入权重）
a = [1, -0.9]         # 分母系数（反馈项）

# 生成输入信号（白噪声）
x = np.random.randn(100)

# 系统输出信号
y = signal.lfilter(b, a, x)

# 可视化
plt.figure(figsize=(10, 5))
plt.plot(x, color='skyblue', linestyle='--', label='输入信号 x[n]')
plt.plot(y, color='orange', linewidth=2, label='输出信号 y[n]')

# 添加标题与说明
plt.title('差分方程系统响应可视化', fontsize=14, weight='bold')
plt.xlabel('样本点 n', fontsize=12)
plt.ylabel('幅值', fontsize=12)

# 美化图表
plt.grid(True, linestyle='--', alpha=0.8)
plt.legend(frameon=True, fontsize=11)
plt.tight_layout()
plt.show()

说明： 差分方程模型可以实现滤波器、反馈系统及线性系统的模拟。

---

4. 频率响应与滤波特性

数字滤波器的频率响应定义为：
$$H(e^{j\omega })=\sum _{n=0}^{N-1}h[n]e^{-j\omega n}$$

$$|H(e^{j\omega })|=\sqrt{(\text{Re}H{)}^2+(\text{Im}H{)}^2}$$

可用于分析滤波器在不同频率下的增益和相位特性：

b = [0.1, 0.1, 0.1]
a = [1, -0.9]

# 频率响应计算
w, h = signal.freqz(b, a)

# 可视化
fig, ax = plt.subplots(2, 1, figsize=(10, 6), sharex=True)

# 幅度响应（dB）
ax[0].plot(w/np.pi, 20*np.log10(np.abs(h)), color='darkorange', linewidth=2)
ax[0].set_title('数字滤波器频率响应', fontsize=14, weight='bold')
ax[0].set_ylabel('幅度 (dB)')
ax[0].grid(True, linestyle='--', alpha=0.4)

# 高亮通带（示例，归一化频率0~0.3）
ax[0].axvspan(0, 0.3, color='skyblue', alpha=0.2, label='通带')
ax[0].legend(loc='lower right')

# 相位响应
ax[1].plot(w/np.pi, np.angle(h), color='seagreen', linewidth=2)
ax[1].set_xlabel('归一化频率 (×π rad/sample)')
ax[1].set_ylabel('相位 (rad)')
ax[1].set_title('滤波器相位响应', fontsize=13)
ax[1].grid(True, linestyle='--', alpha=0.4)

plt.tight_layout()
plt.show()

说明： 幅频与相位响应是滤波器设计与性能评估的重要依据。

---

5. FIR 与 IIR 滤波器设计

5.1 FIR 滤波器

有限冲激响应滤波器（FIR）特点为线性相位和绝对稳定：
$$y[n]=\sum _{k=0}^{N-1}{b}_k,x[n-k]$$

# FIR 低通滤波器设计
fs, cutoff = 500, 50  # 采样率与截止频率
numtaps = 101         # 滤波器阶数
b = signal.firwin(numtaps, cutoff/(fs/2))  # 归一化截止频率

# 频率响应
w, h = signal.freqz(b)

# 可视化
plt.figure(figsize=(10, 5))

# 幅度响应
plt.plot(w*fs/(2*np.pi), 20*np.log10(np.abs(h)), color='darkorange', linewidth=2, label='幅度响应')

# 标注截止频率
plt.axvline(cutoff, color='skyblue', linestyle='--', linewidth=1.5, label=f'截止频率 {cutoff} Hz')

# 标题与标签
plt.title('FIR 低通滤波器频率响应', fontsize=14, weight='bold')
plt.xlabel('频率 (Hz)')
plt.ylabel('幅度 (dB)')
plt.grid(True, linestyle='--', alpha=0.8)
plt.legend()
plt.tight_layout()
plt.show()

5.2 IIR 滤波器

无限冲激响应滤波器（IIR）常用 Butterworth 设计，特点为阶数低、幅频平滑：
$$H(z)=\frac{b_0+b_1{z}^{-1}}{1+a_1{z}^{-1}+a_2{z}^{-2}}$$

# Butterworth IIR 滤波器设计
order = 4        # 滤波器阶数
cutoff = 0.2     # 数字频率 (归一化 0~1 对应 Nyquist)
b, a = signal.butter(order, cutoff)

# 频率响应
w, h = signal.freqz(b, a)

# 可视化
fig, ax = plt.subplots(2, 1, figsize=(10, 6), sharex=True)

# 幅度响应
ax[0].plot(w/np.pi, 20*np.log10(np.abs(h)), color='darkorange', linewidth=2)
ax[0].set_title(f'Butterworth IIR 滤波器频率响应 (阶数={order})', fontsize=14, weight='bold')
ax[0].set_ylabel('幅度 (dB)')
ax[0].grid(True, linestyle='--', alpha=0.4)

# 标注截止频率
ax[0].axvline(cutoff, color='skyblue', linestyle='--', label=f'截止频率 {cutoff*0.5:.2f} × fs/2')
ax[0].legend()

# 相位响应
ax[1].plot(w/np.pi, np.angle(h), color='seagreen', linewidth=2)
ax[1].set_xlabel('归一化频率 (×π rad/sample)')
ax[1].set_ylabel('相位 (rad)')
ax[1].set_title('相位响应', fontsize=13)
ax[1].grid(True, linestyle='--', alpha=0.4)

plt.tight_layout()
plt.show()

说明： FIR/IIR 滤波器可根据应用需求灵活选择。

---

6. 短时傅里叶变换（STFT）

STFT 用滑动窗口分析信号的时变频率：
$$X(t,\omega )=\int x(\tau ),w(\tau -t),e^{-j\omega \tau }d\tau$$

# 信号构造
fs = 500
t = np.linspace(0, 2, fs*2)

# 0-1s 高频 20Hz，1-2s 低频 5Hz
x = np.sin(2*np.pi*20*t*(t<1)) + np.sin(2*np.pi*5*t*(t>=1))

# STFT 计算
f, t_spec, Zxx = signal.stft(x, fs, nperseg=128)

# 可视化
plt.figure(figsize=(10,5))

# 对数幅度显示
plt.pcolormesh(t_spec, f, 20*np.log10(np.abs(Zxx)+1e-6), 
               shading='gouraud', cmap='viridis')  # 'viridis' 更科学，1e-6 避免 log(0)

plt.title('短时傅里叶变换 (STFT)', fontsize=14, weight='bold')
plt.xlabel('时间 [s]', fontsize=12)
plt.ylabel('频率 [Hz]', fontsize=12)
plt.ylim(0, 50)  # 聚焦低频部分
plt.colorbar(label='幅度 [dB]')

plt.tight_layout()
plt.show()

说明： STFT 能有效捕获非平稳信号的频率变化。

---

7. 连续小波变换（CWT）

小波变换提供比傅里叶更好的时频局部化能力：
$$W(a,b)=\int x(t),{\psi }^{\ast }\left(\frac{t-b}{a}\right),dt$$

其中 $a$ 为尺度参数，$b$ 为平移参数。

from scipy.signal import cwt, ricker

# 信号构造
t = np.linspace(0, 1, 500)
x = np.sin(2*np.pi*5*t) + 0.5*np.sin(2*np.pi*20*t)

# CWT 计算
widths = np.arange(1, 50)
cwt_matrix = cwt(x, ricker, widths)

# 可视化
plt.figure(figsize=(10,5))

# 对数幅度显示，增强低幅值可见性
plt.imshow(np.log1p(np.abs(cwt_matrix)), 
           extent=[t[0], t[-1], widths[0], widths[-1]],
           cmap='magma',  # 更适合科学可视化
           aspect='auto', 
           origin='lower')

# 美化标题和坐标轴
plt.title('连续小波变换 (CWT) 时频图', fontsize=16, weight='bold')
plt.xlabel('时间 (s)', fontsize=12)
plt.ylabel('尺度 (Scale)', fontsize=12)

# 添加色条并标注
cbar = plt.colorbar(label='对数幅度 [log(1+|CWT|)]')
cbar.ax.tick_params(labelsize=10)

# 取消网格
plt.grid(False)

plt.tight_layout()
plt.show()

说明： 小尺度对应高频，CWT 对瞬态信号分析更为敏感。

---

8. 实战应用：噪声去除与信号恢复

使用滤波器去噪，恢复信号主体：

# 信号构造
t = np.linspace(0, 1, 500)
x_clean = np.sin(2*np.pi*10*t)
x_noisy = x_clean + 0.5*np.random.randn(len(t))

# Butterworth 滤波器设计
b, a = signal.butter(4, 0.2)  # 4阶低通
x_filtered = signal.filtfilt(b, a, x_noisy)

# 可视化
plt.figure(figsize=(10,5))

# 含噪信号
plt.plot(t, x_noisy, color='steelblue', alpha=0.4, label='含噪信号')

# 滤波后信号
plt.plot(t, x_filtered, color='darkorange', linewidth=2, label='滤波后信号')

# 标题与标签
plt.title('Butterworth 滤波去噪效果', fontsize=14, weight='bold')
plt.xlabel('时间 [s]', fontsize=12)
plt.ylabel('幅值', fontsize=12)

# 网格和图例
plt.grid(True, linestyle='--', alpha=0.8)
plt.legend(loc='upper right', fontsize=11)

plt.tight_layout()
plt.show()

说明： 低通滤波器有效保留信号主体，同时去除高频噪声。

---

9. 总结

SciPy 的 signal 模块提供从卷积、相关性、滤波器设计到时频分析的完整工具链，可实现信号平滑、系统响应模拟、频率特性分析及噪声去除。通过 FIR/IIR 滤波器、STFT 与 CWT 等方法，可高效处理非平稳信号，既支撑科研实验，也满足工程应用，实现信号处理从理论到实践的无缝衔接。