Python实战：用PyWavelets库实现连续小波变换（CWT）信号分析
信号处理领域里，时频分析一直是个让人又爱又恨的话题。传统傅里叶变换就像个固执的老学究，非要你把整个乐章听完才肯告诉你用了哪些音符。而连续小波变换（CWT）则像个敏锐的音乐侦探，能准确指出第三小节第二拍那个走调的音符。今天我们就用Python的PyWavelets库，来解开这个时频分析的魔法。
1. 环境配置与数据准备
工欲善其事，必先利其器。在开始CWT冒险之前，我们需要准备好Python环境和示例数据。推荐使用Anaconda创建专属环境：
conda create -n wavelet python=3.9 conda activate wavelet pip install pywavelets numpy matplotlib scipy对于实际工程信号，采样率的选择至关重要。假设我们分析一个包含多种频率成分的合成信号：
import numpy as np import matplotlib.pyplot as plt # 生成测试信号 fs = 1000 # 采样率1kHz t = np.linspace(0, 1, fs, endpoint=False) signal = (np.sin(2*np.pi*10*t) * (t>0.2) * (t<0.4) + np.sin(2*np.pi*50*t) * (t>0.6) * (t<0.8)) noise = 0.5 * np.random.randn(len(t)) noisy_signal = signal + noise # 可视化信号 plt.figure(figsize=(10,4)) plt.plot(t, noisy_signal) plt.title('原始含噪信号') plt.xlabel('时间(s)') plt.ylabel('幅值') plt.tight_layout()这个信号包含两个时段性的正弦波（10Hz和50Hz），并添加了高斯白噪声。典型的工业场景中，这种瞬时出现的频率成分很常见，比如机械故障诊断中的冲击信号。
2. PyWavelets核心参数解析
PyWavelets提供了多种小波基函数，我们需要理解几个关键参数：
小波类型选择矩阵
| 小波类型 | 复数/实数 | 适用场景 | 参数范围 | 时间分辨率 | 频率分辨率 |
|---|---|---|---|---|---|
| Morlet | 复数 | 通用分析 | m=6-20 | 中等 | 优秀 |
| Paul | 复数 | 瞬态检测 | order=4-40 | 优秀 | 一般 |
| DOG | 实数 | 边缘检测 | order=2-8 | 良好 | 良好 |
尺度参数设计技巧
- 最小尺度：至少2倍采样间隔
- 最大尺度：不超过信号长度的1/4
- 尺度数量：通常30-100个，对数分布更合理
import pywt # 尺度计算函数 def get_scales(min_scale, max_scale, num_scales): return np.logspace(np.log2(min_scale), np.log2(max_scale), num=num_scales, base=2) scales = get_scales(2, 100, 50)提示：实际项目中，建议先用少量尺度测试内存占用，PyWavelets的cwt函数会预先分配 scales×n 的矩阵
3. 完整CWT分析流程
现在我们把所有组件组装起来，实现端到端的CWT分析：
# 执行CWT计算 coef, freqs = pywt.cwt(noisy_signal, scales, 'morl', sampling_period=1/fs) # 时频图绘制 plt.figure(figsize=(10,6)) plt.imshow(np.abs(coef), extent=[0, 1, 1, 100], cmap='jet', aspect='auto', vmax=abs(coef).max(), vmin=-abs(coef).max()) plt.colorbar(label='幅值') plt.title('Morlet小波时频分析') plt.ylabel('频率(Hz)') plt.xlabel('时间(s)') # 标记理论频率位置 plt.axhline(y=10, color='white', linestyle='--', alpha=0.5) plt.axhline(y=50, color='white', linestyle='--', alpha=0.5)这段代码会生成彩色的时频热力图，其中：
- 横轴表示时间
- 纵轴表示频率（线性坐标）
- 颜色深浅表示该时频点的能量强度
常见问题排查指南
- 内存不足：尝试分块处理或减少尺度数量
- 频率显示不准：检查sampling_period参数
- 边缘效应：考虑使用padding模式
- 结果不稳定：尝试不同小波基
4. 工业级优化技巧
在实际工程应用中，我们还需要考虑以下高级技巧：
内存优化方案
# 分块处理大数据 chunk_size = 10000 for i in range(0, len(signal), chunk_size): chunk = signal[i:i+chunk_size] coef, _ = pywt.cwt(chunk, scales, 'morl') # 处理或保存当前块结果并行计算加速
from joblib import Parallel, delayed def process_scale(s): return pywt.cwt(signal, [s], 'morl')[0] results = Parallel(n_jobs=4)(delayed(process_scale)(s) for s in scales) coef = np.concatenate(results, axis=0)自动化特征提取
# 寻找时频图中的局部极值 from skimage.feature import peak_local_max peaks = peak_local_max(np.abs(coef), min_distance=5, threshold_abs=0.2) for peak in peaks: print(f"在时间 {t[peak[1]]:.2f}s 检测到 {freqs[peak[0]]:.1f}Hz 成分")5. 典型应用场景解析
让我们看几个实际案例，理解CWT如何解决具体问题：
案例1：轴承故障诊断
# 模拟轴承故障信号（冲击响应） bearing_signal = np.zeros(2000) for i in range(100, 1900, 150): bearing_signal[i:i+20] = np.exp(-0.5*(np.arange(20)-10)**2) * np.sin(2*np.pi*30*np.arange(20)/1000) # CWT分析 scales = get_scales(5, 50, 40) coef, _ = pywt.cwt(bearing_signal, scales, 'gaus8')案例2：ECG信号分析
from scipy.misc import electrocardiogram ecg = electrocardiogram()[1000:3000] # 使用Mexican hat小波检测R波 coef, _ = pywt.cwt(ecg, np.arange(10,30), 'mexh') r_peaks = np.argmax(coef[15], axis=0) # 选择合适尺度案例3：语音信号处理
import soundfile as sf audio, rate = sf.read('speech.wav') # 使用复数小波分析共振峰 scales = 1/(np.linspace(80, 300, 50)*2*np.pi/rate) coef, _ = pywt.cwt(audio[:8000], scales, 'cmor1.5-1.0')6. 性能对比与进阶方向
当处理超长信号时，我们需要权衡计算精度和效率：
算法性能对比表
| 方法 | 时间复杂度 | 空间复杂度 | 适合场景 | 精度 |
|---|---|---|---|---|
| 直接CWT | O(N²) | O(N²) | 短信号精确分析 | 高 |
| 分块CWT | O(kN) | O(k√N) | 长信号处理 | 中 |
| 多级DWT | O(N) | O(N) | 实时系统 | 低 |
| STFT | O(NlogN) | O(N) | 平稳信号分析 | 中 |
前沿扩展方向
- 结合深度学习：用CNN处理时频图
- 自适应小波选择：根据信号特性自动调整参数
- GPU加速：使用CuPy替代NumPy
- 在线分析：滑动窗口实时处理
在最近的一个电机监测项目中，我们使用Morlet小波（m=12）成功捕捉到了轴承早期磨损产生的7.5Hz特征频率，比传统FFT方法提前两周发现了潜在故障。关键是要根据具体场景反复调整尺度范围和小波参数，有时候微调一个参数就能让特征从噪声中浮现出来。
