第1章：Python信号处理工具链

各位同学好，我是老张。在信号处理这个行当摸爬滚打了十几年，我越来越觉得Python就是咱们的瑞士军刀。今天咱们就来聊聊这套工具链的核心——NumPy、SciPy和Matplotlib。说白了，这三样东西就是信号处理的"三驾马车"。

1.1 NumPy基础与数组操作

NumPy，全称Numerical Python。我个人习惯叫它"数组引擎"。为什么这么说？因为信号处理本质上就是跟数组打交道——采样点就是数组元素，时域信号就是一维数组，频谱图就是二维数组。

创建数组的几种方式

import numpy as np

# 从列表创建
signal = np.array([0.1, 0.5, 0.8, 1.0, 0.6])

# 生成等间隔采样点
t = np.linspace(0, 1, 1000)  # 0到1秒，1000个点

# 生成正弦波
f = 50  # 频率50Hz
x = np.sin(2 * np.pi * f * t)

# 全零数组（常用于初始化缓冲区）
buffer = np.zeros(1024)

我在项目中遇到过一个问题：用Python列表做实时信号处理，结果延迟高得离谱。换成NumPy数组后，速度提升了近20倍。你想想看，这就是向量化操作的力量。

数组索引与切片

嗯，这里要注意。NumPy的索引方式跟Python列表有区别，但更强大。

# 一维信号切片
signal = np.sin(np.linspace(0, 2*np.pi, 100))
first_half = signal[:50]  # 前50个采样点

# 二维数组（比如频谱图）
spectrogram = np.random.randn(256, 128)
row_10 = spectrogram[9, :]  # 第10行所有列
col_20 = spectrogram[:, 19]  # 第20列所有行

# 花式索引
indices = [0, 10, 20, 30]
selected = signal[indices]

数组运算

说白了，NumPy的运算就是"逐元素操作"。这跟MATLAB很像，但更灵活。

# 向量化运算
x = np.array([1, 2, 3, 4])
y = np.array([5, 6, 7, 8])

z = x + y      # 逐元素相加
z = x * y      # 逐元素相乘
z = np.sqrt(x) # 开平方

# 广播机制
a = np.array([[1, 2, 3], [4, 5, 6]])
b = np.array([10, 20, 30])
c = a + b  # b被广播到每一行

1.2 SciPy信号模块(signal)入门

SciPy的signal模块，是我做信号处理时最常用的工具箱。它封装了从滤波器设计到频谱分析的几乎所有功能。

滤波器设计

我记得刚入行时，设计滤波器还要手算系数。现在用SciPy，几行代码搞定。

from scipy import signal

# 设计一个低通滤波器
fs = 1000  # 采样率
cutoff = 100  # 截止频率
order = 4

b, a = signal.butter(order, cutoff/(fs/2), 'low')

# 应用滤波器
filtered_signal = signal.filtfilt(b, a, raw_signal)

频谱分析

频谱分析是信号处理的"眼睛"。SciPy提供了多种工具。

# 快速傅里叶变换
freqs = np.fft.fftfreq(len(signal), 1/fs)
spectrum = np.fft.fft(signal)

# 功率谱密度
f, Pxx = signal.welch(signal, fs, nperseg=256)

# 短时傅里叶变换
f, t, Sxx = signal.spectrogram(signal, fs)

信号生成

做算法验证时，我经常需要生成测试信号。SciPy提供了丰富的波形函数。

# 生成测试信号
t = np.linspace(0, 1, 1000)

# 正弦波 + 噪声
clean = signal.chirp(t, f0=10, t1=1, f1=100)
noisy = clean + 0.5 * np.random.randn(len(t))

# 方波、锯齿波
square_wave = signal.square(2 * np.pi * 5 * t)
sawtooth = signal.sawtooth(2 * np.pi * 5 * t)

1.3 Matplotlib可视化基础

做信号处理，光有数据不行，得看得见。Matplotlib就是咱们的"眼睛"。

基本绘图流程

import matplotlib.pyplot as plt

# 创建图形
plt.figure(figsize=(10, 6))

# 绘制时域信号
plt.subplot(2, 1, 1)
plt.plot(t, signal)
plt.title('时域波形')
plt.xlabel('时间 (s)')
plt.ylabel('幅度')

# 绘制频谱
plt.subplot(2, 1, 2)
plt.plot(freqs[:len(freqs)//2], np.abs(spectrum[:len(spectrum)//2]))
plt.title('频谱')
plt.xlabel('频率 (Hz)')
plt.ylabel('幅度')

plt.tight_layout()
plt.show()

高级可视化技巧

有时候我们需要同时观察时域和频域信息，这时候就得用点"花活"。

# 频谱图（时频分析）
plt.figure(figsize=(12, 4))
plt.specgram(signal, Fs=fs, NFFT=256, noverlap=128)
plt.colorbar(label='强度 (dB)')
plt.title('频谱图')
plt.xlabel('时间 (s)')
plt.ylabel('频率 (Hz)')

# 三维可视化
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(111, projection='3d')
X, Y = np.meshgrid(t, f)
ax.plot_surface(X, Y, Sxx, cmap='viridis')

知识体系结构图

下面这张图，是我花时间整理的。它展示了本章三个工具之间的关系。

实战小练习

光说不练假把式。我建议你动手试试这个练习：

1. 用NumPy生成一个50Hz正弦波，叠加20Hz和80Hz的谐波
2. 用SciPy设计一个带通滤波器，保留40-60Hz成分
3. 用Matplotlib绘制滤波前后的时域波形和频谱对比

好了，这一章的内容就到这里。记住：NumPy是基础，SciPy是工具，Matplotlib是眼睛。三者缺一不可。下一章咱们会深入探讨时域分析的具体技术，到时候见。