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

Python实战用PyWavelets库实现连续小波变换CWT信号分析信号处理领域里时频分析一直是个让人又爱又恨的话题。传统傅里叶变换就像个固执的老学究非要你把整个乐章听完才肯告诉你用了哪些音符。而连续小波变换CWT则像个敏锐的音乐侦探能准确指出第三小节第二拍那个走调的音符。今天我们就用Python的PyWavelets库来解开这个时频分析的魔法。1. 环境配置与数据准备工欲善其事必先利其器。在开始CWT冒险之前我们需要准备好Python环境和示例数据。推荐使用Anaconda创建专属环境conda create -n wavelet python3.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, endpointFalse) signal (np.sin(2*np.pi*10*t) * (t0.2) * (t0.4) np.sin(2*np.pi*50*t) * (t0.6) * (t0.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复数通用分析m6-20中等优秀Paul复数瞬态检测order4-40优秀一般DOG实数边缘检测order2-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), numnum_scales, base2) scales get_scales(2, 100, 50)提示实际项目中建议先用少量尺度测试内存占用PyWavelets的cwt函数会预先分配 scales×n 的矩阵3. 完整CWT分析流程现在我们把所有组件组装起来实现端到端的CWT分析# 执行CWT计算 coef, freqs pywt.cwt(noisy_signal, scales, morl, sampling_period1/fs) # 时频图绘制 plt.figure(figsize(10,6)) plt.imshow(np.abs(coef), extent[0, 1, 1, 100], cmapjet, aspectauto, vmaxabs(coef).max(), vmin-abs(coef).max()) plt.colorbar(label幅值) plt.title(Morlet小波时频分析) plt.ylabel(频率(Hz)) plt.xlabel(时间(s)) # 标记理论频率位置 plt.axhline(y10, colorwhite, linestyle--, alpha0.5) plt.axhline(y50, colorwhite, linestyle--, alpha0.5)这段代码会生成彩色的时频热力图其中横轴表示时间纵轴表示频率线性坐标颜色深浅表示该时频点的能量强度常见问题排查指南内存不足尝试分块处理或减少尺度数量频率显示不准检查sampling_period参数边缘效应考虑使用padding模式结果不稳定尝试不同小波基4. 工业级优化技巧在实际工程应用中我们还需要考虑以下高级技巧内存优化方案# 分块处理大数据 chunk_size 10000 for i in range(0, len(signal), chunk_size): chunk signal[i:ichunk_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_jobs4)(delayed(process_scale)(s) for s in scales) coef np.concatenate(results, axis0)自动化特征提取# 寻找时频图中的局部极值 from skimage.feature import peak_local_max peaks peak_local_max(np.abs(coef), min_distance5, threshold_abs0.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:i20] 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)案例2ECG信号分析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], axis0) # 选择合适尺度案例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. 性能对比与进阶方向当处理超长信号时我们需要权衡计算精度和效率算法性能对比表方法时间复杂度空间复杂度适合场景精度直接CWTO(N²)O(N²)短信号精确分析高分块CWTO(kN)O(k√N)长信号处理中多级DWTO(N)O(N)实时系统低STFTO(NlogN)O(N)平稳信号分析中前沿扩展方向结合深度学习用CNN处理时频图自适应小波选择根据信号特性自动调整参数GPU加速使用CuPy替代NumPy在线分析滑动窗口实时处理在最近的一个电机监测项目中我们使用Morlet小波m12成功捕捉到了轴承早期磨损产生的7.5Hz特征频率比传统FFT方法提前两周发现了潜在故障。关键是要根据具体场景反复调整尺度范围和小波参数有时候微调一个参数就能让特征从噪声中浮现出来。