Python实战：用PyWavelets库实现连续小波变换（CWT）信号分析

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提供了多种小波基函数，我们需要理解几个关键参数：

小波类型选择矩阵

尺度参数设计技巧

• 最小尺度：至少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)

这段代码会生成彩色的时频热力图，其中：

• 横轴表示时间
• 纵轴表示频率（线性坐标）
• 颜色深浅表示该时频点的能量强度

常见问题排查指南

1. 内存不足：尝试分块处理或减少尺度数量
2. 频率显示不准：检查sampling_period参数
3. 边缘效应：考虑使用padding模式
4. 结果不稳定：尝试不同小波基

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. 性能对比与进阶方向

当处理超长信号时，我们需要权衡计算精度和效率：

算法性能对比表

前沿扩展方向

• 结合深度学习：用CNN处理时频图
• 自适应小波选择：根据信号特性自动调整参数
• GPU加速：使用CuPy替代NumPy
• 在线分析：滑动窗口实时处理

在最近的一个电机监测项目中，我们使用Morlet小波（m=12）成功捕捉到了轴承早期磨损产生的7.5Hz特征频率，比传统FFT方法提前两周发现了潜在故障。关键是要根据具体场景反复调整尺度范围和小波参数，有时候微调一个参数就能让特征从噪声中浮现出来。