资讯动态

同步压缩小波变换原理与Python实现:从时频分析到脊线提取

发布时间:2026/9/15 1:59:37 来源:尧图企业网站定制
简介同步压缩小波变换SST程序包面向信号处理与时频分析研究人群适合需要处理非平稳信号、语音、心电或金融数据的工程师和研究生。它将小波多分辨率特性与频率重分配机制结合有效解决传统小波变换在非线性信号中出现的频率混叠与时频模糊问题。压缩包共9个文件以6个Matlab源程序.m为主辅以2个.mat数据文件与1个txt说明文档整体约31KB源码涵盖SST示例、线性小波变换、多带宽检测、小波变换封装及多分量SST主程序覆盖从系数分解、频率重分配到同步压缩谱输出的完整流程。已有1021人浏览学习使用者可据此快速搭建实验环境理解SST核心步骤并将算法与数据结构直接迁移到自身项目中提高非平稳信号处理的效率与可解释性为相关研究与工程落地提供可复用基础。1. 同步压缩小波变换程序为什么值得自己调一遍同步压缩小波变换Synchrosqueezed Wavelet TransformSSWT是近几年时频分析里被反复翻牌的处理工具。它把小波变换的系数沿频率方向“挤”回真实瞬时频率位置解决了一根时变频率线上能量发散的问题。标题里把“变换”写成“变化”搜索引擎会同时命中这两种写法实际指的都是同一个东西。对做振动、生物电、雷达、通信信号处理的人来说这通常不是看一篇科普而是要在本地跑出一张干净的时频图并提取瞬时频率。这篇内容从原理、最小可运行程序到参数调优和脊线提取按我实际做信号的顺序一步一步写。新手可以直接抄代码熟手可以重点看参数边界和逆变换验证手段。2. 同步压缩小波变换的原理频带重排如何拿回时频分辨率2.1 从连续小波变换看时频模糊连续小波变换把信号 (x(t)) 与小波母函数 (\psi) 的伸缩、平移做内积[ W_x(a,b)\frac{1}{a}\int x(t)\psi^*\left(\frac{t-b}{a}\right)dt ]其中 (a) 是尺度对应频率的倒数(b) 是平移对应时间。小波基的等效带宽随中心频率变化低频处频率分辨率高、时间分辨率低高频反过来。这种自适应分辨率是它比短时傅里叶变换更适合非平稳信号的原因。问题是小波系数的能量并不集中在真实的瞬时频率 (\omega(b)) 附近而是发散在一个频带内。例如一个线性调频信号在 (b0.5s) 处的瞬时频率是 100Hz但小波系数会分布在 90Hz 到 110Hz 之间看起来像一条胖带。胖带意味着时频图的频率定位精度不够直接找峰值去估计瞬时频率会产生系统的偏斜。同步压缩的基本思路就是既然小波系数在数学上满足一个相位关系——对纯正弦分量(W_x(a,b)) 对 (b) 求偏导后除以 (W_x) 再取虚部就能得到一个局部的瞬时频率估计那么就可以把这个频率值作为新的横坐标把小波系数累积进去。2.2 同步压缩沿频率方向压缩小波系数先定义候选同步压缩频率。对小波系数 (W_x(a,b))瞬时频率估计为[ \omega_x(a,b)-\frac{\partial_b W_x(a,b)}{2\pi i W_x(a,b)} ]当 (W_x(a,b)\neq 0) 时这个公式对单一分量的解析信号给出准确频率。对多分量信号只要分量够稀疏而且小波母函数有足够好的频域局部性(\omega_x(a,b)) 仍然近似等于每个分量的瞬时频率。接下来做重排。把尺度轴划分为离散区间 (a_j)频率轴划分为离散区间 (\omega_l)。对每个 ((a_j,b))如果计算出的 (\omega_x(a_j,b)) 落在频率区间 ((\omega_l-\frac{1}{2}\Delta\omega, \omega_l\frac{1}{2}\Delta\omega))就把 (W_x(a_j,b),a_j^{-1/2}) 累加到 (T_x(\omega_l,b))[ T_x(\omega_l,b)\sum_{a_j: |\omega_x(a_j,b)-\omega_l| \Delta\omega/2} W_x(a_j,b)a_j^{-1/2} ]这一步把原本散布在尺度方向的系数压缩到频率方向的窄带里。理想情况下一个纯正弦分量在时频图上被压缩成一条沿时间方向的亮线线宽只受频率离散化限制。注意这里压缩的是小波系数的模和相位而不是像重排reassignment那样同时改变时间方向的位置。同步压缩只改变频率坐标保留时间坐标因此支持逆变换原始信号可以从 (T_x) 重建出来这是它相对经典重排方法最大的优势。2.3 与重排方法的区别经典重排方法把时频图的能量同时滑向质心时间和频率两个方向都“重定位”图像非常锐利但不可逆。同步压缩只重排频率方向以损失一定锐利度为代价换取可逆性。这个特性在需要从时频表示中重构模态、滤波降噪的场景里特别有用。另一个区别是瞬时频率的估计。重排法不一定依赖相位导数而同步压缩法的核心就是 (\omega_x(a,b))所以它对信号的解析性有要求。实际使用前通常要把信号做 Hilbert 变换得到解析信号或者直接用复数小波母函数比如 Morlet 小波。ssqueezepy 里的cwt默认采用复数 Morlet 小波就是为了同时获取幅值和相位。从实现角度看同步压缩还有一个隐藏的超参数压缩的阶数。一阶同步压缩只估计瞬时频率 (\omega_x)二阶同步压缩会额外修正频率的变化率对变频速率高的信号如急速扫频效果更好。代价是计算量增加且对噪声更敏感。这个在后面的参数表中会再谈。3. 本地跑通同步压缩小波变换的最小程序3.1 安装 ssqueezepy 与依赖Python 生态里做同步压缩最省事的库是ssqueezepy它把连续小波、同步压缩、逆变换、时频图绘制都封装好了。安装命令pip install ssqueezepy matplotlib numpy依赖里需要numpy、scipy和matplotlib。装完以后验证一下 import 是否正常python -c import ssqueezepy; print(ssqueezepy.__version__)ssqueezepy的底层 CWT 是自己实现的不需要额外安装 PyWavelets。它会根据信号长度自动选择一个合适的尺度数量但如果你传入很大的nv计算时间会明显增加后面会用具体案例说明。3.2 生成多分量仿真信号我先构造一个能体现同步压缩效果的信号前半段是 50Hz 正弦后半段跳到 120Hz并在 0.4s 到 0.6s 加入一个线性调频成分。这样既有瞬时频率突变又有连续扫频适合对比算法效果。import numpy as np fs 1000 T 1.0 N 1024 t np.linspace(0, T, N, endpointFalse) x np.zeros_like(t) x[:N//2] np.cos(2*np.pi*50 * t[:N//2]) x[N//2:] np.cos(2*np.pi*120 * t[N//2:]) x np.cos(2*np.pi*(50 100*t) * t) # 0~100Hz 线性调频作为干扰这里混合了三个分量。注意第三个分量的瞬时频率是 (50200t)在中间时刻达到 150Hz会和 120Hz 分量混叠。代码注释里写清楚每个分量的频率范围方便后面观察同步压缩的选择性。3.3 一行计算同步压缩小波变换核心调用只需要一行from ssqueezepy import ssq_cwt Tx, CWT, ssq_freqs, cwt_freqs ssq_cwt(x, waveletmorlet, nv16)Tx是同步压缩后的时频系数形状为(len(ssq_freqs), len(x))CWT是压缩前的小波系数ssq_freqs是同步压缩频率轴单位 Hzcwt_freqs是小波变换的离散频率轴。绘制时频图时我一般把系数取模并做对数压缩import matplotlib.pyplot as plt def plot_tf(coefficients, freqs, t, title): plt.figure(figsize(8, 4)) extent [t[0], t[-1], freqs[0], freqs[-1]] plt.imshow(np.abs(coefficients), aspectauto, extentextent, originlower, cmapturbo) plt.colorbar(labelmagnitude) plt.xlabel(Time (s)) plt.ylabel(Frequency (Hz)) plt.title(title) plt.tight_layout() plot_tf(CWT, cwt_freqs, t, CWT) plot_tf(Tx, ssq_freqs, t, SSWT) plt.show()这段代码的关键在于extent必须用对应的频率轴否则图像会上下颠倒或错位。originlower让频率从小到大沿 y 轴向上排列。3.4 关键参数表ssq_cwt的参数并不算多但每个都对结果有实质影响。我经常用的几个如下参数默认值作用常见取值waveletmorlet小波母函数复数小波才能提相位morlet需要时用bumpnv16每倍频程内的语音数voices per octave8/16/32越大频带越细n_octavesauto尺度轴覆盖的倍频程数自动也能手动给整数l1_normTrue尺度归一化方式L1 范数更适合频率定位通常保持 Trueorder1同步压缩阶数1 或 2变频率信号用 2threshold0.1小波系数的模相对阈值小于阈值不参与压缩0.05~0.2噪声大调高nv是最关键的参数。举例来说如果nv8每倍频程 8 个尺度一个 8 倍频程的信号就有 (8\times864) 个尺度nv32时尺度数到 256压缩后的频率分辨率更高但计算量和内存占用也翻倍。对 1 秒钟、1024 点的信号nv16是速度和精度的平衡点。4. 参数调整与信号分离实战4.1 用小波母函数控制时频聚焦Morlet 小波的形状由它的中心频率和带宽参数决定。ssqueezepy里的morlet用小波基类管理可以通过eps和res等参数间接控制。在实际业务里大部分场景直接用默认 Morlet 就够唯一的例外是希望得到更窄的频带或者更强的频率分离能力。bump小波在频域内有紧支撑同步压缩效果更干净因为它没有 Morlet 长尾巴带来的频谱泄漏。但 bump 小波在低频处的时间定位更差对瞬态冲击信号反而模糊。我在做轴承故障数据时通常会对比两种母函数from ssqueezepy import ssq_cwt Tx_bump, _, fs_bump, _ ssq_cwt(x, waveletbump, nv16)比较fs_bump和morlet的频率轴可以发现bump 在高频段的频率分布更稀疏。因此如果你关注的是 1kHz 以上的高频分量用 bump 未必划算反而 Morlet 更均匀。4.2 压缩阈值和噪声抑制同步压缩的噪声表现有个特点小波系数中的噪声幅度通常小于信号成分因此设置一个幅度阈值把低于阈值的系数在压缩前置零能显著减少时频图上的散点。ssq_cwt的threshold参数是相对值它会统计小波系数的峰值然后把低于峰值threshold倍的系数当成噪声。默认值 0.1对干净信号表现很好对信噪比低于 10dB 的信号我一般调到 0.2 到 0.3把背景噪声压掉但调高了也会削弱弱信号的分量。一个更稳妥的做法是保留原始 CWT不做阈值而是在同步压缩后应用数据依赖的掩膜Tx_soft Tx.copy() max_mag np.abs(Tx_soft).max() Tx_soft[np.abs(Tx_soft) 0.05 * max_mag] 0这样做的理由同步压缩已经将能量集中到窄带噪声如果不在那一带就被压缩掉在频带内的噪声用统一阈值处理后时频图更利于脊线提取。4.3 对比 CWT 与 SSWT 的时频图用上面构造的信号跑一遍你会看到 CWT 时频图里 50Hz 频带明显比 SSWT 宽尤其是边界处有模糊拖尾。120Hz 分量与线性调频分量交叉的区域CWT 图上两条线黏在一起SSWT 图上则能看到清晰的分叉。如果看不到这个效果先检查两个问题第一nv是否太小小于 8 时频率离散化不够第二信号是否包含不可忽略的直流分量直流会让最低频率附近出现一条虚假高频带需要先x x - np.mean(x)去直流。还有一个常见的坑是横轴单位ssq_cwt默认把采样率当成 1 来算频率轴是归一化频率。如果你的信号采样率是 1000Hz必须传入fsfs参数Tx, CWT, ssq_freqs, cwt_freqs ssq_cwt(x, waveletmorlet, nv16, fsfs)原来的ssq_freqs的取值范围只有 0 到最大尺度对应的归一化频率乘以fs才是真实频率。我经常用这种信号做回归验证如果一个算法在合成信号上都不能把两个相近频率分开那它到了实测信号上也不会变好。同步压缩最典型的失败场景是两个分量的频率比小于 (2^{1/nv})这时尺度轴上的采样间隔不够压缩后依然是一条粗带。要分开它们只能提高nv比如 32 或 64或者先带通滤波把两个分量分离再分别分析。5. 进阶同步压缩逆变换与脊线提取5.1 用逆变换验证参数是否合适同步压缩保留相位信息所以理论上能够从 (T_x) 重建原始信号。ssqueezepy提供了issq_cwt函数from ssqueezepy import issq_cwt x_rec issq_cwt(Tx, ssq_freqs, cwt_freqs, l1_normTrue, order1)重建后的信号和原信号的误差能直观反映同步压缩过程中丢了多少信息。如果只是观察频带参数怎么调都可以但如果后续要做滤波或者反变换就要关注x_rec与x的相关系数。相关系数低于 0.9 时说明压缩丢掉了太多细节常见的应对是增大nv或使用order2的二阶同步压缩。二阶同步压缩对快速扫频信号的重建更准但不适合强噪声环境。有一个经验法则是先用一阶跑通再用信号重建误差和时频图锐利度交叉验证两种方式同时变好才确定参数不要只看时频图是否漂亮。5.2 脊线提取估计瞬时频率很多场景下最终产出不是时频图而是瞬时频率曲线。同步压缩后的时频矩阵里每个时刻能量最大处就对应主分量的瞬时频率。提取脊线可以按下面几步做def extract_ridge(Tx, ssq_freqs): Tx_mag np.abs(Tx) ridge_idx np.argmax(Tx_mag, axis0) ridge_freq ssq_freqs[ridge_idx] return ridge_freq ridge_freq extract_ridge(Tx, ssq_freqs)但直接取最大值的缺点是如果某个时刻主分量和干扰分量强度相近脊线会跳变。更稳健的做法是加一个频率连续性惩罚让相邻时刻的频率不会突变。实现时可以对每一列的能量峰用动态规划惩罚项系数一般在 10~50 之间惩罚越大脊线越平滑。把脊线叠加在时频图上同时画出理论瞬时频率就能直观评估整个处理流程。同步压缩小波变换程序的价值就在这里用出来了它能把时频图从“一团雾”变成“几条线”而几条线就足够支撑故障判据、调制识别、模态分离等后续判断。拿到任何新信号我建议都先跑一遍第 3 节的最小程序再决定要不要上二阶压缩和动态规划脊线而不是一上来就堆参数。本文还有配套的精品资源点击获取

读完文章,也想定制专属网站?

尧图设计师 24 小时内与您沟通定制方案

免费获取报价