简介面向信号处理与图像分析领域学习者的连续小波变换C源码包以单个cwt.cpp文件实现CWT核心算法适合正在研究时频局部化分析、希望快速掌握算法工程化写法的开发人员与高年级学生。代码覆盖小波基函数选择、尺度与位置参数遍历、小波系数计算等关键步骤并涉及标准库容器和复数运算通过计算结果可进一步绘制小波谱图用于观察语音、振动、生物医学等非平稳信号频率随时间的变化为去噪和特征提取提供基础。包内仅有1个cpp文件压缩包约2KB结构紧凑、便于阅读和移植到实际项目中。目前已有594人学习/下载。借助这份源码可把CWT理论公式落实为可运行的程序逻辑并结合Morlet小波、墨西哥帽等常见基函数的局部性质深入理解连续小波变换在实际信号分析中的实现要点与应用价值。1. 连续小波变换不是换张热力图那么简单做信号处理的人迟早会撞上连续小波变换CWT。它常以“时频热力图”的形式出现在论文里看起来只是把短时傅里叶变换STFT换了个配色于是不少人拿到pywt.cwt跑完就画图画完才发现纵轴上的“尺度”不知道怎么换算成赫兹改两个参数图就碎成雪花。这个工具的核心不是那张图而是用一组可变宽的窗去扫描信号低频段用长窗换频率分辨高频段用短窗换时间定位恰好匹配振动冲击、心电节律、地震波这类“缓变趋势加瞬态脉冲”共存的信号。这篇博文把尺度标定、母函数选择、参数陷阱和结果验证一次讲透结尾给一个轴承故障特征频率提取的完整示例新手能顺着跑通老手可以对照检查自己的参数习惯。2. 连续小波变换的数学骨架尺度、平移与频率映射2.1 CWT 与 STFT 的本质区别短时傅里叶变换固定窗长整个时频平面上每个点的分辨率标记完全一样窗短则频率分辨率差窗长则时间分辨率差这个矛盾在 STFT 里无法靠换参数消除。连续小波变换换了个思路用一族同源函数对信号做内积W(a, b) 1/√a ∫ x(t) · ψ*((t - b)/a) dt其中 a 是尺度b 是平移ψ 是小波母函数。a 变小ψ((t-b)/a) 被压缩等效窗变短适合抓高频瞬态a 变大窗被拉长频率分辨变细。和 STFT 不同CWT 在低频段看的是“很长一段时间里的频率结构”在高频段看的是“很短时间内发生的波形细节”这才是它处理非平稳信号时胜出的根本原因。母函数还必须满足小波允许条件即积分为零所以小波天然是振荡的不是随便拿个窗函数就能当母函数。这一点直接决定了后面选型和调参的方向允许条件保证 CWT 是可逆的也意味着系数的大小能横向比较。2.2 尺度怎么换算成频率大多数小波母函数有一个中心频率 fc可以理解为母函数在原始尺度下的主振荡频率。在采样率 fs 下伪频率的计算关系是f fc · fs / a其中伪频率的含义是小波带通滤波的中心频率而不是信号真实瞬时频率。PyWavelets 里可以直接查import pywt wavelet cmor1.5-1.0 for s in [1, 10, 100]: f_per_sample pywt.scale2frequency(wavelet, s) # 返回 fc / s单位是“周期/采样点”乘采样率才是赫兹 print(s, f_per_sample)这段代码的输出分别是 1.0、0.1、0.01乘上 fs20000 就是 20000 Hz、2000 Hz、200 Hz。逻辑上scale2frequency只负责算fc / a物理频率轴要自己乘采样率或者更省事直接让pywt.cwt的sampling_period参数替你算。常见小波的中心频率参考值如下小波中心频率 fc近似说明cmor1.5-1.01.00参数即带宽 1.5、中心频率 1.0morl约 0.813实数 Morlet系数包络带振荡mexh约 0.25墨西哥帽时间定位好频率选择性弱用的时候别手抄这些数直接在代码里pywt.scale2frequency(wavelet, 1)[0]取一下不同版本和不同实现存在细微差异。2.3 为什么工程上默认选复数小波实数小波morl、mexh对瞬态的响应是振荡包络时频图上会出现明暗相间的“梳齿”对脉冲起始时刻的定位很不友好。复数小波系数本身带实部和虚部取模|W(a,b)|得到的是平滑包络相位还能单独拿出来做瞬时频率估计。因此做故障诊断、心电分析时我一般默认用cmor这类复数小波只有纯瞬态检测任务才退回mexh。代价是复数系数计算量翻倍但对 10 万点量级的数据完全可忽略。3. 用 PyWavelets 跑通 CWT最小可复现代码3.1 pywt.cwt 的最小调用先把最简版本跑起来确认版本和返回结构。PyWavelets 1.x 的pywt.cwt返回(cfs, freqs)元组如果你的环境只返回一个数组说明版本太老先升级再继续。import numpy as np import pywt fs 1000 # 采样率 1000 Hz t np.arange(0, 2, 1 / fs) x np.cos(2 * np.pi * 10 * t) # 10 Hz 持续分量 x[500:520] np.exp(-np.arange(20) / 4) \ * np.cos(2 * np.pi * 300 * np.arange(20) / fs) # 短时冲击 scales np.arange(1, 200) # 尺度 1..199线性取 cfs, freqs pywt.cwt(x, scales, cmor1.5-1.0, sampling_period1 / fs) print(cfs.shape) # (199, 2000)复数系数矩阵 print(freqs[:3], freqs[-3:]) # 高频端在前低频端在后单位 Hzcfs的第一维对应scales的顺序不是频率升序freqs和scales一一对应第二维是时间。这里有个关键约定sampling_period1/fs传进去后返回的freqs直接就是物理赫兹这一步最容易漏。核心参数含义如下参数含义常见取值scales尺度序列决定频带覆盖范围np.arange或对数均匀序列wavelet小波名或小波对象cmor1.5-1.0sampling_period采样间隔单位秒决定 freqs 单位1 / fsmethod卷积实现方式默认auto3.2 用线性扫频校验频率轴光看形状不能证明标定正确。用一个瞬时频率已知的线性扫频信号做冒烟测试脊线和理论值对得上才算通过from scipy.signal import chirp T 4 t_chirp np.linspace(0, T, fs * T, endpointFalse) y chirp(t_chirp, f020, f180, t1T, methodlinear) scales np.arange(2, 400) cfs, freqs pywt.cwt(y, scales, morl, sampling_period1 / fs) ridge np.argmax(np.abs(cfs), axis0) # 每时刻最大能量对应的尺度行 f_inst freqs[ridge] # 沿脊线取伪频率 f_true 20 (80 - 20) * t_chirp / T # 理论瞬时频率 err np.abs(f_inst - f_true).max() print(最大频率偏差 Hz:, err)argmax(axis0)沿尺度轴找每个时刻能量最大的位置得到的就是瞬时频率脊线。干净扫频信号下误差应该只有几赫兹如果偏出两位数先检查尺度范围是否盖住 20~80 Hz再检查sampling_period是否写成了fs。注意脊线用argmax在噪声下会跳行实际项目里要对脊线做中值滤波或是加频率变化率约束。3.3 一次实验看懂幅值、相位与边界热力图的纵轴要换成freqs注意它是降序的绘图时用extent指定坐标范围import matplotlib.pyplot as plt plt.figure(figsize(10, 4)) plt.pcolormesh(t_chirp, freqs, np.abs(cfs), shadinggouraud) plt.yscale(log) plt.ylabel(Frequency (Hz)) plt.colorbar()对复数小波np.abs(cfs)是时频幅值谱np.angle(cfs)是相位谱。幅值谱看能量在哪相位谱在匀加速目标场景里可以用于提取瞬时频率变化率。这张图最明显的问题在两端卷积在信号边界截断会产生竖直的伪影带这就是锥形影响区COI第 6 章会给量化剔除方法。4. 连续小波变换的调参实战母函数、尺度序列与采样率4.1 母函数选型这张表贴在工位上选母函数没有绝对最优只有任务匹配母函数类型拿手场景要注意的坑cmorB-C复故障诊断、心电、语音基频B 太小频带窄瞬态响应拖尾morl实脊线轮廓、频带粗估包络振荡别直接取极大值定位mexh实突变点、奇异点检测高频衰减快尺度范围要放宽cgauP复边缘检测、相位分析阶数 P 越大振荡越剧烈gausP实求导特征近似中心频率随阶数变必须现场查询cmor的 B 和 C 是被问得最多的参数。cmorB-C中 B 是带宽C 是中心频率。B 越大母函数在频率域越宽时间定位越好但频率分辨越差C 决定小波的自然响应频率C 越高同样频带需要越小的尺度计算量随之上涨。我常用的起点是cmor1.5-1.0需要更细频率分辨就用cmor0.8-1.0更强调瞬态定位用cmor2.5-1.0。4.2 尺度序列别再用 np.arange 一套打天下线性尺度只适合频带很窄的场景。更通用的做法是对数均匀尺度让频率轴上的每个十倍频程都有相同的采样密度fmin, fmax 20.0, 500.0 # 目标频带单位 Hz wavelet cmor1.5-1.0 fc pywt.scale2frequency(wavelet, 1)[0] num 128 scales fc * fs / np.geomspace(fmax, fmin, num) cfs, freqs pywt.cwt(x, scales, wavelet, sampling_period1 / fs)逻辑尺度与频率成反比频率对数均匀分布自然对应尺度对数均匀分布。好处是低频端尺度间隔大但频带本身就窄点数不浪费高频端也不会因为线性取尺度而过密。num取 128 够可视化做脊线提取建议 192 到 256。fmax我一般压到0.4 * fs而不是 0.5给小波自身带宽留出余量避免 Nyquist 附近的畸变污染整个高频区。4.3 三个必踩的坑第一个坑是坐标轴方向。freqs降序排列直接 plot 会得到上下颠倒的热力图用extent或pcolormesh传参时注意纵轴范围写成[freqs[-1], freqs[0]]。第二个坑是尺度下限。scale1时小波只覆盖几个采样点系数基本是边缘伪影。建议scale_min不小于 2或者干脆用fmax 0.4 * fs反推最小尺度。第三个坑是内存。cfs是复数矩阵128 个尺度乘 100 万样本点一个complex128矩阵就要 2 GB。对策是两条路先按目标频带降采样scipy.signal.decimate找不到现成工具时再分段做 CWT段间留 10%~20% 重叠最后拼接时把 COI 段丢弃。提示对数尺度序列几乎都不是整数pywt 默认的methodauto会走 FFT 路径速度远快于逐尺度卷积如果尺度是等间隔整数直接np.arange反而更稳。4.4 怎么判断参数失败了热力图上出现竖直等间距条纹说明尺度上限太低没盖住低频缓变成分。高频区一片空白先查fmax是否被0.4 * fs卡掉了。整张图只有噪点没有结构通常是母函数带宽 B 太小频率域窄到只响应极窄的频带换大 B 再试。这些判断比看任何误差指标都快。5. CWT 故障特征频率提取从振动信号到故障频率5.1 合成一段带故障特征的振动信号用合成信号验证流程比直接上现场数据靠谱因为真值在手算法对不对一目了然import numpy as np from scipy.signal import butter, sosfiltfilt, hilbert import pywt fs 20000 T 1.0 t np.arange(0, T, 1 / fs) rng np.random.default_rng(42) fault_freq 50.0 # 故障特征频率就是要找回的目标 res_freq 3200.0 # 结构共振频带中心 imp np.zeros_like(t) for k in range(int(T * fault_freq)): idx int(round(k * fs / fault_freq)) rng.integers(-2, 3) imp[idx % len(t)] 1.0 tau np.arange(0, 0.002, 1 / fs) h np.exp(-tau * 1500.0) * np.cos(2 * np.pi * res_freq * tau) x np.convolve(imp, h)[:len(t)] x 0.25 * rng.standard_normal(len(t))这段模拟的是滚动轴承外圈故障的经典三要素周期冲击、共振衰减、白噪声。每个故障周期触发一次冲击冲击激发 3200 Hz 的共振衰减振荡。rng.integers(-2, 3)给冲击位置加 ±2 个采样点的抖动模拟轴承滚动体滑移避免频谱里出现完美校准的数学峰值。5.2 用 CWT 时频图提取瞬时能量脊wavelet cmor1.5-1.0 fmin, fmax 100.0, 6000.0 fc pywt.scale2frequency(wavelet, 1)[0] scales fc * fs / np.geomspace(fmax, fmin, 192) cfs, freqs pywt.cwt(x, scales, wavelet, sampling_period1 / fs) amp np.abs(cfs) ridge_idx np.argmax(amp, axis0) env_cwt amp[ridge_idx, np.arange(len(t))] env_cwt[:fs // 10] 0.0 # 丢弃左端 COI env_cwt[-fs // 10:] 0.0 # 丢弃右端 COI spec np.abs(np.fft.rfft(env_cwt - env_cwt.mean())) f_axis np.fft.rfftfreq(len(env_cwt), 1 / fs) top_idx np.argsort(spec)[-5:] print(峰值频率:, np.sort(f_axis[top_idx]))ridge_idx每时刻取能量最大的尺度行等价于自动选带。相比固定通带包络CWT 脊线不做任何频带假设共振频率随温度或负载漂移时仍然能抓到能量。对env_cwt做 FFT 后峰值应该出现在 50、100、150 Hz 等故障特征频率的整数倍处。如果只出基频不出谐波通常是冲击间隔抖动太大属正常现象。5.3 与固定带通包络谱对照经典做法是固定带通加希尔伯特包络代码很短sos butter(8, [res_freq - 400, res_freq 400], btypeband, fsfs, outputsos) y sosfiltfilt(sos, x) env_fix np.abs(hilbert(y))两种方法对照起来看方法需要先验频带对共振频率漂移额外信息固定带通 希尔伯特包络需要敏感漂移 10% 即失效无CWT 脊线包络不需要鲁棒脊线自动跟随热力图本身就是诊断证据固定带通的问题是共振频率一漂滤波器就落在带外包络幅度大幅衰减。CWT 脊线的好处不只是免调参时频图上任何一处异常都能倒回去定位发生时刻这在故障诊断报告里是硬需求。6. 收尾技巧给连续小波变换结果做三件自检6.1 用 icwt 重构检查信息是否丢失xr pywt.icwt(cfs, scales, wavelet, sampling_period1 / fs) keep slice(fs // 10, -fs // 10) # 掐掉边界段再算误差 nmse np.linalg.norm(x[keep] - xr[keep]) / np.linalg.norm(x[keep]) print(归一化重构误差:, nmse)pywt.icwt必须和pywt.cwt使用同一组scales和同一个wavelet。归一化重构误差超过 5%先怀疑尺度范围没盖住主要频带再检查是不是把 COI 段算进了误差。这一步能在三分钟内区分“算法问题”和“参数问题”。6.2 用单位脉冲量化 COIimp_test np.zeros(len(t)) imp_test[len(t) // 2] 1.0 cfs_imp, _ pywt.cwt(imp_test, scales, wavelet, sampling_period1 / fs) foot np.abs(cfs_imp[:, len(t) // 2:]) coi_half np.array([np.argmax(foot[i] foot[i].max() / np.e) for i in range(len(scales))]) coi_time coi_half / fs # 每个尺度对应的 e 折时间半径把单位脉冲放在时间中心脉冲响应的泄漏范围就是这个尺度下的最短可分辨时间。系数衰减到峰值 e 分之一的时间点之外都是受边界影响的不可信区域。如果某个尺度返回 0说明该尺度已经大于数据长度的一半整行系数都不能用。画图时把这个区域用半透明遮罩盖住审稿人和同事都不会再挑边界伪影的毛病。6.3 用叠加信号做消融校验最后一条经验把已知成分叠加起来验证整条链路。构造“20 Hz 正弦 80 Hz 线性扫频 300 Hz 衰减脉冲”跑完整套 CWT 和脊线提取确认三个成分都能在时频图上被辨认、脊线频率和预设一致再上真实数据。这个习惯能提前暴露尺度范围、母函数带宽和 COI 截断三者之间的匹配问题。实际操作中把fmax压到 0.4 倍采样率、尺度数开到 192、提取脊线前对幅值矩阵做一次 3×3 中值滤波绝大多数时频图上的雪花噪声都在这一步消失。本文还有配套的精品资源点击获取
