简介压缩包内是一份连续小波变换CWT的C实现源码面向信号处理学习者和需要开展时频分析的开发者提供从理论到代码落地的参考范例。包内仅含1个cpp文件体积约2KB代码精简聚焦小波基函数选择、尺度与位置参数设置、小波系数计算等核心环节单文件设计便于快速阅读、修改和移植。已有594人学习下载适合作为理解CWT算法细节的入门素材。通过运行和调试这段代码可以直观观察不同尺度下信号的时频响应理解连续小波变换对非平稳信号的局部化分析能力代码中关于卷积计算和尺度遍历的思路还能为后续去噪、边缘检测等应用提供扩展基础。实现中还涉及时频谱图的可视化处理思路可帮助使用者将抽象的小波系数转化为直观的频率分布视图进一步提升对信号特征的把握能力。1. 连续小波变换不是“更聪明的傅里叶变换”连续小波变换不是把信号拆成正弦波的叠加——它把同一个被称为母小波的有限长波形做缩放和平移再去和信号做内积。傅里叶变换只能告诉你频率存在说不清频率出现在哪个时刻短时傅里叶变换虽然引入了时间窗口窗口一旦固定低频和高频的顾此失彼就无法避免。连续小波变换CWT在低频处用宽窗口换频率分辨率高频处用窄窗口换时间分辨率因此特别适合 chirp 扫描、瞬态冲击和机电振动这类时变频率信号。如果你要做时频图、提取瞬时频率或者从强噪声里抓短时脉冲这就是选 CWT 而不是其他工具的理由。下面从系数定义、参数设置到验证手段把这条链路讲清楚。2. 连续小波变换的核心母小波、尺度与系数2.1 母小波为什么必须能量归一化CWT 的定义式是W_f(a, b) (1/√a) ∫ f(t) ψ*((t − b)/a) dt其中 a 是尺度b 是平移量ψ* 是母小波的共轭。式中的 1/√a 就是能量归一化项。假设母小波 ψ 的能量为 1那么缩放后的 ψ((t−b)/a) 能量会变成 a除以 √a 后能量回到 1这样不同尺度下得到的系数才有可比性。如果不做这一步大尺度对应的波形被拉长内积结果天然偏大时频图会整体向低频倾斜高频细节被压制。母小波还需要满足容许条件即它的均值为零∫ψ(t)dt 0。这个条件保证 CWT 在零频率处没有响应直流分量和缓慢变化的趋势不会污染高频段结果。morl 这种高斯包络复指数波形本身均值不为零PyWavelets 在实现里做了去均值修正修正会微调中心频率。所以做精确的频率换算时不要用教科书上的理论中心频率而应该以 pywt.scale2frequency 的返回值为准。提示CWT 系数不是信号某一频率分量的幅值而是信号与伸展后母小波的投影值。系数大小既取决于信号分量幅度也取决于信号与母小波的波形相似程度。2.2 a尺度与b位移时频格子怎么移动a 变小母小波在时域上被压缩频谱变宽适合捕捉高频瞬态a 变大波形被拉长频带收窄适合观察低频持续振荡。b 决定分析窗口的中心位置b 连续移动就得到时间轴。CWT 与 STFT 最大的差别就在这儿STFT 的窗长固定一旦选定 nperseg全频段的时间分辨率都一致CWT 的窗长随尺度自动伸缩把时间和频率分辨率的选择权交给了 a。海森堡不确定原理在这套体系里表现为时间分辨率 Δt 与频率分辨率 Δf 的乘积有下界。CWT 没有绕过这个下界只是在每个尺度上做了不同的分配高频处 Δt 小、Δf 大低频处反过来。这也是为什么 CWT 时频图里同一段 chirp 信号的高频段看起来“更细”低频段“更宽”这是物理约束不是绘图问题。代码层面可以先确认一下中心频率import pywt # 计算 morl 小波在 scale1 时的中心频率单位周期/采样点 fc pywt.scale2frequency(morl, 1, precision10) print(fc) # 约 0.8125 周期/采样点precision 参数控制中心频率积分时的计算精度。默认精度在批量换算大量尺度时误差会累积用 precision10 更稳。这个 fc 乘上采样率 fs才是 scale1 对应的实际频率单位是 Hz。2.3 系数矩阵的行列与相位信息怎么读pywt.cwt 返回的 coefs 是二维数组行对应 scales 里的每个尺度列对应时间采样点。配套返回的 freqs 与 scales 逐位对应单位由 sampling_period 决定。读系数矩阵时先搞清楚行列含义再去取幅值或相位否则很容易把尺度轴当成频率轴直接用。coefs 元素索引含义物理含义读取方式coefs[i, j]i 行 第 i 个尺度j 列 第 j 个采样点尺度 i、时刻 j 处的复数投影np.abs(coefs[i, j]) 为强度coefs.real同行列与偶对称母小波匹配的分量复小波下可用coefs.realcoefs.imag同行列与奇对称形态匹配的分量np.imag(coefs) 可用freqs[i]与尺度 i 一一对应伪频率单位与 sampling_period 设定有关有了相位信息才能沿尺度方向做瞬时频率插值。实小波如 mexh的 coefs 是实数不存在合法的相位通道这时强行 np.angle(coefs) 得到的只有 0 或 π画出来是黑白噪声没有任何物理意义。这是一个新手高频踩坑点后面 4.6 会说清楚复小波和实小波在相位上的分工。3. 用Python实现连续小波变换的最小可跑代码3.1 构造一个带时变频率的测试信号先造一个能验证 CWT 效果的信号频率随时间线性变化的 chirp叠加一段固定频率正弦。这样时频图上既能看到连续的频率斜坡又能看到局部干扰对脊线的影响。import numpy as np fs 1000 # 采样率 1 kHz t np.arange(0, 2, 1/fs) # 线性 chirp2 秒内瞬时频率从 5 Hz 升到 40 Hz f_start, f_end 5.0, 40.0 inst_freq f_start (f_end - f_start) * t / t[-1] phase 2 * np.pi * np.cumsum(inst_freq) / fs x np.sin(phase) # 1.2 秒到 1.6 秒之间注入一段 20 Hz 正弦模拟叠加干扰 x[1200:1600] 0.5 * np.sin(2 * np.pi * 20 * t[1200:1600]) # 归一化到 [-1, 1]避免幅度量级影响后续系数阈值 x x / np.max(np.abs(x))这里直接写 np.sin(2 * np.pi * inst_freq * t) 是错误的因为瞬时频率是相位的导数必须对频率做积分得到相位cumsum(inst_freq) / fs 等价于做数值积分。用错公式后信号的实际瞬时频率会偏移后面验证脊线轨迹时会对不上。叠加的 20 Hz 正弦幅度取 0.5是为了让它在时频图上比主 chirp 暗一个量级方便观察小波系数幅值对能量的响应。3.2 用PyWavelets的cwt函数做变换import pywt # 尺度从 2^0 到 2^6对数取 100 个点 scales np.logspace(0, 6, num100, base2) # sampling_period1/fs让返回的 freqs 为实际 Hz coefs, freqs pywt.cwt(x, scales, waveletmorl, sampling_period1/fs) print(coefs.shape) # (100, 2000) print(freqs[0], freqs[-1])scales 不用 np.arange 是因为尺度与频率大致成反比线性间隔会让低频段在频率轴上越来越稀疏高频段又密集冗余。np.logspace(0, 6, base2) 表示每个倍频程采样点固定时频图在频率轴上视觉分布更均匀。num100 对大多数分析够用再大只是让图更平滑内存和计算时间会成倍增加。waveletmorl 是默认复小波coefs 是复数数组。sampling_period1/fs 必须传否则 freqs 的单位是周期/采样点画图时纵轴会差 fs 倍。如果漏了这个参数也可以事后手动换算freq_hz pywt.scale2frequency(morl, scale) * fs。3.3 把尺度轴换算成频率轴再画时频图import matplotlib.pyplot as plt plt.figure(figsize(10, 5)) # extent 的 y 顺序是从 freqs[-1] 到 freqs[0] # 因为 imshow 第 0 行画在图像顶部第 0 行对应最小尺度、最大频率 plt.imshow(np.abs(coefs), aspectauto, cmapturbo, extent[0, t[-1], freqs[-1], freqs[0]]) plt.ylim(0, 60) plt.colorbar(label|CWT coef|) plt.xlabel(Time [s]) plt.ylabel(Frequency [Hz]) plt.title(CWT scalogram) plt.show()extent 的 y 参数写成 [freqs[-1], freqs[0]] 而不是 [freqs[0], freqs[-1]]是因为 imshow 默认第 0 行放在最上面而 coefs 第 0 行对应最小尺度、最高频率。如果发现频率轴上下颠倒先检查这里。plt.ylim(0, 60) 把眼光集中在有用频段否则 scale1 对应的 812.5 Hz 会把整个图压成一条细带。4. 连续小波变换的6个必调参数与高频坑4.1 小波基选择morl、paul、mexh与复/实取舍小波基决定系数矩阵的形态也决定 scales 的取值范围。PyWavelets 里最常见的四个小波基差别很大小波名实/复中心频率采样周期1特点与典型场景morl复约 0.8125通用时频分析频率分辨率和相位信息平衡paul复m4 时约 0.75频域更集中脊线更窄适合瞬时频率提取mexh实约 0.25墨西哥帽二阶导数特征适合奇异点检测gaus8实约 0.089高斯八阶导数消失矩高适合冲击检测中心频率直接决定尺度 1 对应的伪频率因此切换小波基后原来的 scales 数组必须重新换算。Paul 小波对高频细节的分辨能力更强但时域波形尾部有振荡对孤立脉冲响应不干净做脉冲检测时不如 gaus8。4.2 scales数组对数间隔、线性间隔与Nyquist边界scales 数组是最容易出问题的参数。常见错误是用 np.arange(1, 100, 1) 这种线性序列。线性间隔看似覆盖范围宽实际在低频端稀疏、高频端密集画出来的时频图高频处挤成一片低频处又连不成脊线。推荐做法# 目标频率范围 5 Hz 到 100 Hz target_f_min, target_f_max 5.0, 100.0 # fc_cycles 是母小波中心频率fs 是采样率 scales fc * fs / np.linspace(target_f_max, target_f_min, 80)这里的 fc 是 2.2 用 scale2frequency 算出的中心频率周期/采样点。scales 是从小到大排列所以分母从 target_f_max 到 target_f_min 递减对应尺度从小到大递增。80 个尺度在视觉上足够再多不会带来额外信息。还有一个 Nyquist 边界坑。当 fs1000 Hz 时morl 小波 scale1 对应 812.5 Hz已经超过 500 Hz 的奈奎斯特频率。直接用 scale 从 1 开始的数组会在高频段画出混叠响应。所以 scales 下限要取 fc * fs / (fs/2) 以上对 morl 大约是 1.63安全起见从 2 开始。4.3 采样率与伪频率的换算关系伪频率公式是f_actual fc / (scale * dt) fc * fs / scale实际工作中往往是先有目标频率范围再反推采样率或尺度范围。给定目标最高频率 f_max 和采样率 fsscales 下限应不小于 fc * fs / f_max。给定目标最低频率 f_min 和最大尺度可用长度scales 上限受到信号长度限制不能无限大。def scale_for_freq(target_hz, fs, fc_cycles): 把目标频率换算成对应尺度 return fc_cycles * fs / target_hz # 例fs2000morl 小波想看到 10 Hz 成分 print(scale_for_freq(10.0, 2000.0, 0.8125)) # 162.5尺度 162.5 意味着母小波被拉长到 162 个采样周期以上对应时域支撑大约 162/1016 秒的分析窗。如果信号本身只有 2 秒长这个尺度下限根本没有足够数据支撑时频图低频部分全是边界伪迹。这就是 4.4 要处理的问题。4.4 边界效应与锥形区域COICWT 内核是有限支撑的信号两端会截断低频大尺度下截断影响更明显。时频图两边会出现向低频发散的亮色区域这就是锥形区域COI。处理办法是延拓再变换最后裁回原长度# 取最大尺度对应支撑长度的 2 倍作为延拓长度 max_scale scales[-1] pad_len int(2 * max_scale) xp np.pad(x, (pad_len, pad_len), modereflect) coefs_pad, _ pywt.cwt(xp, scales, morl, sampling_period1/fs) coefs_trim coefs_pad[:, pad_len:-pad_len]反射延拓比零填充好零填充等于在两边强行加阶跃会在所有尺度上引入高频响应。pad_len 取最大尺度的 1 到 3 倍具体看母小波的衰减速度。裁剪后的 coefs_trim 长度与原始信号一致但边界区系数仍不算数分析时至少舍弃时间轴前后 5% 的区域。4.5 模极大值定位奇异点与脉冲用模极大值可以定位脉冲和突变点。基本步骤是先取幅度矩阵沿尺度方向找局部极大值再用阈值过滤弱响应from scipy.signal import argrelmax mag np.abs(coefs_trim) thr 0.2 * mag.max() ridge_rows, ridge_cols argrelmax(mag, axis0, order2) ridge_matrix np.zeros_like(mag, dtypebool) ridge_matrix[ridge_rows, ridge_cols] mag[ridge_rows, ridge_cols] thrargrelmax 在 axis0 上找局部极大值order2 表示与相邻两个点比较抑制单点噪声。返回的是索引元组两个一维数组分别对应行和列直接解包成 ridge_rows、ridge_cols。阈值过滤后剩下的是能量足够的脊线片段。想要连成完整脊线可以用 scipy.ndimage.label 做连通域标记再按长度筛选而不是把散点直接连成线。4.6 实小波与复小波在相位信息上的差别复小波同时给出幅度和相位。对 morl 的 coefs 逐点取 np.angle再沿时间轴做 np.unwrap就得到连续相位轨迹相位求导除以 2π 就是瞬时频率# ridge_idx 是沿尺度方向取得极大值的行索引 inst_phase np.unwrap(np.angle(coefs_trim[ridge_idx, :])) inst_freq_est np.diff(inst_phase) / (2 * np.pi * dt)dt 是采样周期 1/fs。这个办法信噪比好但因为复小波的高斯包络会给相位引入固定偏移使用前要先把中心频率对应的相位参考扣除。实小波 coefs 是实数取 angle 得到 0 或 π求导全是窄脉冲没有任何物理意义。脉冲检测、边缘提取这类只关心位置的任务可以用实小波瞬时频率和相位类任务必须选复小波。拿实小波结果硬算相位是看到无意义锯齿最常见的来源。5. 验证连续小波变换结果的两个手段5.1 用逆变换icwt验证重构误差rec pywt.icwt(coefs, scales, morl, sampling_period1/fs) rel_err np.sqrt(np.mean((x - rec) ** 2)) / np.sqrt(np.mean(x ** 2)) print(rel_err)icwt 是 CWT 的逆变换重构误差能反映参数选得是否合理。实际操作会发现低频端误差偏大原因是 scales 数组没有覆盖足够低的频率范围信号在低频段的信息没被完整捕获。增大最大尺度、或选频域衰减更慢的小波基能降低误差。注意 icwt 对 scales 有单调递增要求且首尾尺度对应的伪频率要覆盖主要信号带宽随意给个不完整的 scales 数组重构结果会差到无法使用。5.2 用STFT交叉验证瞬时频率轨迹from scipy.signal import stft f_st, t_st, Zxx stft(x, fs, nperseg256, noverlap224) plt.pcolormesh(t_st, f_st, np.abs(Zxx))CWT 的结果不直观时用 STFT 做交叉验证最直接。STFT 的 nperseg256 在 fs1000 Hz 下的频率分辨率约 3.9 Hz时间分辨率约 0.128 秒chirp 主脊线在图上是一条略宽的亮带。CWT 在低频端能给出更细的频率分辨率高频端时间分辨率更好两条脊线的主轨迹应基本重合。对比时只看主脊线穿过每个频率的时刻是否一致背景噪声形态不同是正常的。5.3 先看脊线再看误差验证CWT结果的一个习惯ridge_col np.argmax(np.abs(coefs_trim), axis0) ridge_freq freqs[ridge_col] # 舍弃前后各 10% 的边界样本避免 COI 干扰 margin int(len(ridge_freq) * 0.1) mae np.mean(np.abs(ridge_freq[margin:-margin] - inst_freq[margin:-margin])) print(mae)验证 CWT 结果是否正确第一步不是看重构误差而是把脊线频率和理论瞬时频率画在同一张图上。如果脊线轨迹偏移超过一个频带宽度说明 scales 范围选窄了或者小波基选错了如果脊线在某段突然断裂大概率是边界效应在干扰。MAE 数值只作参考更可靠的判断是看最大偏差出现在哪个时间区段——出现在两端是延拓不够出现在中间说明信号本身有突变。先画脊线再做定量误差计算每次跑 CWT 都重复这个习惯能省下大量隐蔽的参数调试时间。本文还有配套的精品资源点击获取 连续小波变换Python实战：从时频分析原理到参数调优
简介压缩包内是一份连续小波变换CWT的C实现源码面向信号处理学习者和需要开展时频分析的开发者提供从理论到代码落地的参考范例。包内仅含1个cpp文件体积约2KB代码精简聚焦小波基函数选择、尺度与位置参数设置、小波系数计算等核心环节单文件设计便于快速阅读、修改和移植。已有594人学习下载适合作为理解CWT算法细节的入门素材。通过运行和调试这段代码可以直观观察不同尺度下信号的时频响应理解连续小波变换对非平稳信号的局部化分析能力代码中关于卷积计算和尺度遍历的思路还能为后续去噪、边缘检测等应用提供扩展基础。实现中还涉及时频谱图的可视化处理思路可帮助使用者将抽象的小波系数转化为直观的频率分布视图进一步提升对信号特征的把握能力。1. 连续小波变换不是“更聪明的傅里叶变换”连续小波变换不是把信号拆成正弦波的叠加——它把同一个被称为母小波的有限长波形做缩放和平移再去和信号做内积。傅里叶变换只能告诉你频率存在说不清频率出现在哪个时刻短时傅里叶变换虽然引入了时间窗口窗口一旦固定低频和高频的顾此失彼就无法避免。连续小波变换CWT在低频处用宽窗口换频率分辨率高频处用窄窗口换时间分辨率因此特别适合 chirp 扫描、瞬态冲击和机电振动这类时变频率信号。如果你要做时频图、提取瞬时频率或者从强噪声里抓短时脉冲这就是选 CWT 而不是其他工具的理由。下面从系数定义、参数设置到验证手段把这条链路讲清楚。2. 连续小波变换的核心母小波、尺度与系数2.1 母小波为什么必须能量归一化CWT 的定义式是W_f(a, b) (1/√a) ∫ f(t) ψ*((t − b)/a) dt其中 a 是尺度b 是平移量ψ* 是母小波的共轭。式中的 1/√a 就是能量归一化项。假设母小波 ψ 的能量为 1那么缩放后的 ψ((t−b)/a) 能量会变成 a除以 √a 后能量回到 1这样不同尺度下得到的系数才有可比性。如果不做这一步大尺度对应的波形被拉长内积结果天然偏大时频图会整体向低频倾斜高频细节被压制。母小波还需要满足容许条件即它的均值为零∫ψ(t)dt 0。这个条件保证 CWT 在零频率处没有响应直流分量和缓慢变化的趋势不会污染高频段结果。morl 这种高斯包络复指数波形本身均值不为零PyWavelets 在实现里做了去均值修正修正会微调中心频率。所以做精确的频率换算时不要用教科书上的理论中心频率而应该以 pywt.scale2frequency 的返回值为准。提示CWT 系数不是信号某一频率分量的幅值而是信号与伸展后母小波的投影值。系数大小既取决于信号分量幅度也取决于信号与母小波的波形相似程度。2.2 a尺度与b位移时频格子怎么移动a 变小母小波在时域上被压缩频谱变宽适合捕捉高频瞬态a 变大波形被拉长频带收窄适合观察低频持续振荡。b 决定分析窗口的中心位置b 连续移动就得到时间轴。CWT 与 STFT 最大的差别就在这儿STFT 的窗长固定一旦选定 nperseg全频段的时间分辨率都一致CWT 的窗长随尺度自动伸缩把时间和频率分辨率的选择权交给了 a。海森堡不确定原理在这套体系里表现为时间分辨率 Δt 与频率分辨率 Δf 的乘积有下界。CWT 没有绕过这个下界只是在每个尺度上做了不同的分配高频处 Δt 小、Δf 大低频处反过来。这也是为什么 CWT 时频图里同一段 chirp 信号的高频段看起来“更细”低频段“更宽”这是物理约束不是绘图问题。代码层面可以先确认一下中心频率import pywt # 计算 morl 小波在 scale1 时的中心频率单位周期/采样点 fc pywt.scale2frequency(morl, 1, precision10) print(fc) # 约 0.8125 周期/采样点precision 参数控制中心频率积分时的计算精度。默认精度在批量换算大量尺度时误差会累积用 precision10 更稳。这个 fc 乘上采样率 fs才是 scale1 对应的实际频率单位是 Hz。2.3 系数矩阵的行列与相位信息怎么读pywt.cwt 返回的 coefs 是二维数组行对应 scales 里的每个尺度列对应时间采样点。配套返回的 freqs 与 scales 逐位对应单位由 sampling_period 决定。读系数矩阵时先搞清楚行列含义再去取幅值或相位否则很容易把尺度轴当成频率轴直接用。coefs 元素索引含义物理含义读取方式coefs[i, j]i 行 第 i 个尺度j 列 第 j 个采样点尺度 i、时刻 j 处的复数投影np.abs(coefs[i, j]) 为强度coefs.real同行列与偶对称母小波匹配的分量复小波下可用coefs.realcoefs.imag同行列与奇对称形态匹配的分量np.imag(coefs) 可用freqs[i]与尺度 i 一一对应伪频率单位与 sampling_period 设定有关有了相位信息才能沿尺度方向做瞬时频率插值。实小波如 mexh的 coefs 是实数不存在合法的相位通道这时强行 np.angle(coefs) 得到的只有 0 或 π画出来是黑白噪声没有任何物理意义。这是一个新手高频踩坑点后面 4.6 会说清楚复小波和实小波在相位上的分工。3. 用Python实现连续小波变换的最小可跑代码3.1 构造一个带时变频率的测试信号先造一个能验证 CWT 效果的信号频率随时间线性变化的 chirp叠加一段固定频率正弦。这样时频图上既能看到连续的频率斜坡又能看到局部干扰对脊线的影响。import numpy as np fs 1000 # 采样率 1 kHz t np.arange(0, 2, 1/fs) # 线性 chirp2 秒内瞬时频率从 5 Hz 升到 40 Hz f_start, f_end 5.0, 40.0 inst_freq f_start (f_end - f_start) * t / t[-1] phase 2 * np.pi * np.cumsum(inst_freq) / fs x np.sin(phase) # 1.2 秒到 1.6 秒之间注入一段 20 Hz 正弦模拟叠加干扰 x[1200:1600] 0.5 * np.sin(2 * np.pi * 20 * t[1200:1600]) # 归一化到 [-1, 1]避免幅度量级影响后续系数阈值 x x / np.max(np.abs(x))这里直接写 np.sin(2 * np.pi * inst_freq * t) 是错误的因为瞬时频率是相位的导数必须对频率做积分得到相位cumsum(inst_freq) / fs 等价于做数值积分。用错公式后信号的实际瞬时频率会偏移后面验证脊线轨迹时会对不上。叠加的 20 Hz 正弦幅度取 0.5是为了让它在时频图上比主 chirp 暗一个量级方便观察小波系数幅值对能量的响应。3.2 用PyWavelets的cwt函数做变换import pywt # 尺度从 2^0 到 2^6对数取 100 个点 scales np.logspace(0, 6, num100, base2) # sampling_period1/fs让返回的 freqs 为实际 Hz coefs, freqs pywt.cwt(x, scales, waveletmorl, sampling_period1/fs) print(coefs.shape) # (100, 2000) print(freqs[0], freqs[-1])scales 不用 np.arange 是因为尺度与频率大致成反比线性间隔会让低频段在频率轴上越来越稀疏高频段又密集冗余。np.logspace(0, 6, base2) 表示每个倍频程采样点固定时频图在频率轴上视觉分布更均匀。num100 对大多数分析够用再大只是让图更平滑内存和计算时间会成倍增加。waveletmorl 是默认复小波coefs 是复数数组。sampling_period1/fs 必须传否则 freqs 的单位是周期/采样点画图时纵轴会差 fs 倍。如果漏了这个参数也可以事后手动换算freq_hz pywt.scale2frequency(morl, scale) * fs。3.3 把尺度轴换算成频率轴再画时频图import matplotlib.pyplot as plt plt.figure(figsize(10, 5)) # extent 的 y 顺序是从 freqs[-1] 到 freqs[0] # 因为 imshow 第 0 行画在图像顶部第 0 行对应最小尺度、最大频率 plt.imshow(np.abs(coefs), aspectauto, cmapturbo, extent[0, t[-1], freqs[-1], freqs[0]]) plt.ylim(0, 60) plt.colorbar(label|CWT coef|) plt.xlabel(Time [s]) plt.ylabel(Frequency [Hz]) plt.title(CWT scalogram) plt.show()extent 的 y 参数写成 [freqs[-1], freqs[0]] 而不是 [freqs[0], freqs[-1]]是因为 imshow 默认第 0 行放在最上面而 coefs 第 0 行对应最小尺度、最高频率。如果发现频率轴上下颠倒先检查这里。plt.ylim(0, 60) 把眼光集中在有用频段否则 scale1 对应的 812.5 Hz 会把整个图压成一条细带。4. 连续小波变换的6个必调参数与高频坑4.1 小波基选择morl、paul、mexh与复/实取舍小波基决定系数矩阵的形态也决定 scales 的取值范围。PyWavelets 里最常见的四个小波基差别很大小波名实/复中心频率采样周期1特点与典型场景morl复约 0.8125通用时频分析频率分辨率和相位信息平衡paul复m4 时约 0.75频域更集中脊线更窄适合瞬时频率提取mexh实约 0.25墨西哥帽二阶导数特征适合奇异点检测gaus8实约 0.089高斯八阶导数消失矩高适合冲击检测中心频率直接决定尺度 1 对应的伪频率因此切换小波基后原来的 scales 数组必须重新换算。Paul 小波对高频细节的分辨能力更强但时域波形尾部有振荡对孤立脉冲响应不干净做脉冲检测时不如 gaus8。4.2 scales数组对数间隔、线性间隔与Nyquist边界scales 数组是最容易出问题的参数。常见错误是用 np.arange(1, 100, 1) 这种线性序列。线性间隔看似覆盖范围宽实际在低频端稀疏、高频端密集画出来的时频图高频处挤成一片低频处又连不成脊线。推荐做法# 目标频率范围 5 Hz 到 100 Hz target_f_min, target_f_max 5.0, 100.0 # fc_cycles 是母小波中心频率fs 是采样率 scales fc * fs / np.linspace(target_f_max, target_f_min, 80)这里的 fc 是 2.2 用 scale2frequency 算出的中心频率周期/采样点。scales 是从小到大排列所以分母从 target_f_max 到 target_f_min 递减对应尺度从小到大递增。80 个尺度在视觉上足够再多不会带来额外信息。还有一个 Nyquist 边界坑。当 fs1000 Hz 时morl 小波 scale1 对应 812.5 Hz已经超过 500 Hz 的奈奎斯特频率。直接用 scale 从 1 开始的数组会在高频段画出混叠响应。所以 scales 下限要取 fc * fs / (fs/2) 以上对 morl 大约是 1.63安全起见从 2 开始。4.3 采样率与伪频率的换算关系伪频率公式是f_actual fc / (scale * dt) fc * fs / scale实际工作中往往是先有目标频率范围再反推采样率或尺度范围。给定目标最高频率 f_max 和采样率 fsscales 下限应不小于 fc * fs / f_max。给定目标最低频率 f_min 和最大尺度可用长度scales 上限受到信号长度限制不能无限大。def scale_for_freq(target_hz, fs, fc_cycles): 把目标频率换算成对应尺度 return fc_cycles * fs / target_hz # 例fs2000morl 小波想看到 10 Hz 成分 print(scale_for_freq(10.0, 2000.0, 0.8125)) # 162.5尺度 162.5 意味着母小波被拉长到 162 个采样周期以上对应时域支撑大约 162/1016 秒的分析窗。如果信号本身只有 2 秒长这个尺度下限根本没有足够数据支撑时频图低频部分全是边界伪迹。这就是 4.4 要处理的问题。4.4 边界效应与锥形区域COICWT 内核是有限支撑的信号两端会截断低频大尺度下截断影响更明显。时频图两边会出现向低频发散的亮色区域这就是锥形区域COI。处理办法是延拓再变换最后裁回原长度# 取最大尺度对应支撑长度的 2 倍作为延拓长度 max_scale scales[-1] pad_len int(2 * max_scale) xp np.pad(x, (pad_len, pad_len), modereflect) coefs_pad, _ pywt.cwt(xp, scales, morl, sampling_period1/fs) coefs_trim coefs_pad[:, pad_len:-pad_len]反射延拓比零填充好零填充等于在两边强行加阶跃会在所有尺度上引入高频响应。pad_len 取最大尺度的 1 到 3 倍具体看母小波的衰减速度。裁剪后的 coefs_trim 长度与原始信号一致但边界区系数仍不算数分析时至少舍弃时间轴前后 5% 的区域。4.5 模极大值定位奇异点与脉冲用模极大值可以定位脉冲和突变点。基本步骤是先取幅度矩阵沿尺度方向找局部极大值再用阈值过滤弱响应from scipy.signal import argrelmax mag np.abs(coefs_trim) thr 0.2 * mag.max() ridge_rows, ridge_cols argrelmax(mag, axis0, order2) ridge_matrix np.zeros_like(mag, dtypebool) ridge_matrix[ridge_rows, ridge_cols] mag[ridge_rows, ridge_cols] thrargrelmax 在 axis0 上找局部极大值order2 表示与相邻两个点比较抑制单点噪声。返回的是索引元组两个一维数组分别对应行和列直接解包成 ridge_rows、ridge_cols。阈值过滤后剩下的是能量足够的脊线片段。想要连成完整脊线可以用 scipy.ndimage.label 做连通域标记再按长度筛选而不是把散点直接连成线。4.6 实小波与复小波在相位信息上的差别复小波同时给出幅度和相位。对 morl 的 coefs 逐点取 np.angle再沿时间轴做 np.unwrap就得到连续相位轨迹相位求导除以 2π 就是瞬时频率# ridge_idx 是沿尺度方向取得极大值的行索引 inst_phase np.unwrap(np.angle(coefs_trim[ridge_idx, :])) inst_freq_est np.diff(inst_phase) / (2 * np.pi * dt)dt 是采样周期 1/fs。这个办法信噪比好但因为复小波的高斯包络会给相位引入固定偏移使用前要先把中心频率对应的相位参考扣除。实小波 coefs 是实数取 angle 得到 0 或 π求导全是窄脉冲没有任何物理意义。脉冲检测、边缘提取这类只关心位置的任务可以用实小波瞬时频率和相位类任务必须选复小波。拿实小波结果硬算相位是看到无意义锯齿最常见的来源。5. 验证连续小波变换结果的两个手段5.1 用逆变换icwt验证重构误差rec pywt.icwt(coefs, scales, morl, sampling_period1/fs) rel_err np.sqrt(np.mean((x - rec) ** 2)) / np.sqrt(np.mean(x ** 2)) print(rel_err)icwt 是 CWT 的逆变换重构误差能反映参数选得是否合理。实际操作会发现低频端误差偏大原因是 scales 数组没有覆盖足够低的频率范围信号在低频段的信息没被完整捕获。增大最大尺度、或选频域衰减更慢的小波基能降低误差。注意 icwt 对 scales 有单调递增要求且首尾尺度对应的伪频率要覆盖主要信号带宽随意给个不完整的 scales 数组重构结果会差到无法使用。5.2 用STFT交叉验证瞬时频率轨迹from scipy.signal import stft f_st, t_st, Zxx stft(x, fs, nperseg256, noverlap224) plt.pcolormesh(t_st, f_st, np.abs(Zxx))CWT 的结果不直观时用 STFT 做交叉验证最直接。STFT 的 nperseg256 在 fs1000 Hz 下的频率分辨率约 3.9 Hz时间分辨率约 0.128 秒chirp 主脊线在图上是一条略宽的亮带。CWT 在低频端能给出更细的频率分辨率高频端时间分辨率更好两条脊线的主轨迹应基本重合。对比时只看主脊线穿过每个频率的时刻是否一致背景噪声形态不同是正常的。5.3 先看脊线再看误差验证CWT结果的一个习惯ridge_col np.argmax(np.abs(coefs_trim), axis0) ridge_freq freqs[ridge_col] # 舍弃前后各 10% 的边界样本避免 COI 干扰 margin int(len(ridge_freq) * 0.1) mae np.mean(np.abs(ridge_freq[margin:-margin] - inst_freq[margin:-margin])) print(mae)验证 CWT 结果是否正确第一步不是看重构误差而是把脊线频率和理论瞬时频率画在同一张图上。如果脊线轨迹偏移超过一个频带宽度说明 scales 范围选窄了或者小波基选错了如果脊线在某段突然断裂大概率是边界效应在干扰。MAE 数值只作参考更可靠的判断是看最大偏差出现在哪个时间区段——出现在两端是延拓不够出现在中间说明信号本身有突变。先画脊线再做定量误差计算每次跑 CWT 都重复这个习惯能省下大量隐蔽的参数调试时间。本文还有配套的精品资源点击获取