拓冰网站定制 
简介面向非平稳信号分析与时频特征提取场景一份基于Matlab的Cohen类时频分布计算程序包可用于对语音、振动、雷达等信号进行时频分析支持WVD、CWD、PWVD等典型分布类型适合信号处理研究人员、算法工程师和研究生在科研与教学中使用。压缩包共含139个文件其中134个m源文件覆盖从核心时频计算、核函数设置到结果可视化与交互查看工具3个mat数据文件提供示例信号另有asv备份与txt说明文件补充使用提示整体体积约1.01MB。目前已有560人下载学习配套的系列演示脚本可以让读者直接运行并对比不同分布的效果通过修改参数观察WVD的交叉项干扰、CWD的核函数抑制效果以及PWVD的时频聚集性变化从而深入理解Cohen类分布的原理与工程应用。适合作为课程设计、毕业设计或算法预研的参考尤其有助于快速搭建时频分析实验环境。1. 时频分布不是画个谱图Cohen类才是统一框架平时用spectrogram看非平稳信号很多人会踩一个坑窗长一短频率模糊窗长一长时间分辨率丢掉怎么调都像在跟不确定性原理讨价还价。WVDWigner-Ville 分布的提出就是为了绕过这个限制——它不加窗直接对瞬时自相关做傅里叶变换单分量 chirp 的时频脊线能锐到接近极限。但代价立刻跟上来多分量信号会在真实频率的中途生成交叉项两条时间频率曲线之间出现一排震荡伪影。Cohen 类时频分布把这堆方法收进同一个表达式区别只剩核函数Φ(θ,τ)的形状。WVD 是 Φ1PWVD 是只对 τ 加窗CWD 是用二维指数核压交叉项。这篇文章就从 Cohen 类框架讲起给出不依赖额外工具箱的 MATLAB 程序并实际对比 WVD、CWD、PWVD 三者的参数效果最后给一个用连通域和凸包自动评估交叉项与聚集度的验证脚本。适合正在做雷达回波、振动故障、脑电时频分析手里有 MATLAB 但不想为时频分析单独装扩展箱的读者。2. 从模糊函数到三种核函数WVD、CWD、PWVD 的选型依据2.1 为什么 Cohen 类公式里必须有模糊函数Cohen 类的统一写法是$$ C_x(t,f)\iint A_x(\theta,\tau),\Phi(\theta,\tau),e^{j2\pi(\theta t-f\tau)},d\theta,d\tau $$其中A_x(θ,τ)是信号的模糊函数Ambiguity Function定义是瞬时自相关对时间做傅里叶变换$$ A_x(\theta,\tau)\int x(u\tau/2),x^*(u-\tau/2),e^{j2\pi\theta u},du $$模糊函数把信号能量散布在迟延-多普勒平面上自项集中在原点附近交叉项则偏离原点。Cohen 类的核函数Φ(θ,τ)本质上就是对模糊平面做一次滤波保留接近原点的部分抑制远处的交叉项。这个视角比“给时频图做平滑”精确得多因为它直接告诉你在哪个域动手。2.2 WVD核函数为 1 的最优分辨率与交叉项代价WVD 取Φ(θ,τ)1意味着模糊平面上的所有成分原样搬到时频面。单分量线性调频信号的 WVD 是一条几乎 δ 函数的直线时频聚集度是所有二次型分布里的上限。但多分量信号场景完全不同考虑两个 chirpx(t)x1(t)x2(t)瞬时自相关展开后会得到四项两个自项加两个互项。互项在模糊平面上落在偏离原点的地方核函数不拦于是时频图上出现频率介于 f1(t) 与 f2(t) 之间的震荡条纹。交叉项的幅度可以达到自项的两倍而且频率越高、信号越长伪影越密肉眼很难跟真实分量区分开。实际工程里单分量信号用 WVD 完全没问题多分量信号直接上 WVD 等于把交叉项当分析对象。2.3 CWD 的指数核与 σ 的物理含义CWDChoi-Williams 分布采用指数核$$ \Phi(\theta,\tau)\exp\left(-\frac{\theta^2\tau^2}{\sigma}\right) $$核函数在 θ 和 τ 两个方向都是高斯形状。σ 大时指数核趋近 1分布退化成 WVD交叉项恢复σ 小时核收缩交叉项被压下去但自项在时频面上的支撑区也会跟着模糊分辨率下降。这个参数不像窗长那样直观从模糊函数的角度理解就简单了σ 决定你允许模糊平面上哪些成分穿过滤波器。仿真经验上σ 取 0.5 到 10 之间语音和振动信号常用 1 附近参数越小越干净越大越锐利不存在“免费午餐”。2.4 PWVD 加窗的本质是对 τ 维截断PWVD伪 Wigner-Ville 分布的核函数只依赖 τ不依赖 θ$$ \Phi(\theta,\tau)h(\tau) $$这等价于在计算瞬时自相关 R(t,τ) 时只取一小段迟延范围。时间分辨率由窗长决定窗越短时间定位越准频率方向的主瓣越宽窗越长频率越集中但瞬时频率突变的地方会被抹平。跟 CWD 相比PWVD 对交叉项的抑制是各向同性的截断交叉项还在只是被窗长限制在局部区域CWD 则是衰减型的交叉项幅度显著下降。选型没有一个万能答案下面这张表是三种分布最直接的区分分布核函数需要调的参数主要问题WVD1无交叉项强PWVDh(τ)窗长 L、窗型仅抑制远的 τ 方向交叉项CWDexp(-θ²τ²/σ)σ自项边缘被展宽3. 用 MATLAB 从零实现 Cohen 类时频分布不依赖工具箱的最小程序3.1 先做两件事解析信号与测试 chirp直接对实信号计算 WVD 会看到负频率分量产生的交叉项叠加在正频率上图面完全没法看所以实现之前要把信号转成解析信号MATLAB 里一行hilbert就够。测试信号我习惯用两个频率错开的线性调频叠加这样交叉项位置可以预判fs 1024; t (0:1023) / fs; s1 chirp(t, 50, t(end), 150); % 50-150 Hz s2 chirp(t, 200, t(end), 350); % 200-350 Hz x hilbert(s1 s2); % 强制使用解析信号注意hilbert的返回值已经是复数解析信号后面所有自相关运算都不需要再取实部。两个 chirp 频率区间不重叠交叉项会出现在 175 Hz 附近的中间带这个先验知识稍后用来验证核函数是否生效。3.2 主程序模糊域乘以核函数再二维变换Cohen 类的离散实现路径很直接先算模糊函数A(θ,τ)乘上核函数Φ(θ,τ)然后对 θ 做 IFFT 得到时间轴、对 τ 做 FFT 得到频率轴。完整函数如下function [TFR, f, t] cohen_tfd(x, fs, phi) % COHEN_TFD 模糊域实现的 Cohen 类时频分布 % x : 解析信号列向量 % fs : 采样率 % phi : 核函数句柄 phi(theta, tau) % theta 单位为 rad/sampletau 单位为样本数 N numel(x); M floor(N/2) - 1; % 最大迟延 tau -M:M; % 1. 瞬时自相关 - 模糊函数 A zeros(N, numel(tau)); for k 1:numel(tau) d tau(k); n0 max(1, 1-d); n1 min(N, N-d); idx n0:n1; R x(idxd) .* conj(x(idx-d)); % R(n,tau) A(:, k) fft(R, N); % 对 n 做 FFT - theta 轴 end % 2. 乘核函数 theta_axis (0:N-1) / N * 2 * pi; % 归一化角频率 [Th, Tr] meshgrid(theta_axis, tau); B A .* phi(Th, Tr).; % 注意转置对应维度 % 3. theta 维 IFFT - 时间tau 维 FFT - 频率 TFR_tau zeros(N, numel(tau)); for k 1:numel(tau) TFR_tau(:, k) ifft(B(:, k)); end TFR zeros(N, N); for n 1:N TFR(n, :) real(fft(TFR_tau(n, :), N)); % 补零到 N 点 end TFR TFR(:, 1:N/21); % 只保留非负频率 f (0:N/2) / N * fs; t (0:N-1) / fs; end这段代码的逻辑分三层第一步瞬时自相关x(idxd).*conj(x(idx-d))体现双线性结构fft(R,N)把它从时间维投影到多普勒维得到的A就是离散模糊函数第二步用meshgrid生成二维网格逐元素乘核函数第三步两次一维 FFT/IFFT 把模糊平面换回时频平面。B A .* phi(Th,Tr).里的转置是维度对齐的关键phi输出维度是(2M1)×N而A是N×(2M1)写错维度会直接报矩阵尺寸错误。3.3 三种分布的核函数写法调用函数时核函数用匿名函数传入。WVD 最容易直接返回全 1 矩阵phi_wvd (th, tr) ones(size(th)); [TFR_wvd, f, t] cohen_tfd(x, fs, phi_wvd);CWD 的指数核写成exp(-(th.*tr).^2/sigma)注意 θ 和 τ 是网格矩阵所以用点乘sigma 1; phi_cwd (th, tr) exp(-(th .* tr).^2 / sigma);PWVD 是 τ 维加窗。窗长取2*M1会退化成 WVD实际要短得多。用汉明窗构造只依赖 τ 的核L 127; % 窗长需为奇数 hwin hann(L); phi_pwvd (th, tr) repmat(hwin(:), 1, size(th, 2));hann(L)生成 L 点窗repmat沿 θ 方向复制保证核函数在每一列 θ 上都取同一个 τ 窗序列。这里只依赖 τ 的特性正是 PWVD 与 CWD 的本质区别也是代码里唯一需要改的地方。运行之后用imagesc(t, f, abs(TFR).^2)查看WVD 会看到中间带明显的栅栏状交叉项CWD 和短窗 PWVD 的中间带则干净很多。4. CWD 和 PWVD 参数怎么设交叉项与分辨率的取舍实测4.1 交叉项区域能量占比一个可量化的调参目标调参不能只靠眼睛看颜色深浅。Cohen 类分布里真实分量在时频面上的位置是已知的交叉项落在两条曲线之间的空白区。拿上面的双 chirp 测试信号来说175±30 Hz、时间中段区域只应有交叉项能量。定义两个频带信号带 50–150 Hz 与 200–350 Hz交叉带 145–205 Hz分别统计这些区域的平均幅度比值就是交叉项抑制效果。MATLAB 里用布尔索引圈出区域cross_mask (f 145 f 205); sig_mask (f 50 f 150) | (f 200 f 350); E_cross mean(mean(abs(TFR(:, cross_mask)))); E_sig mean(mean(abs(TFR(:, sig_mask)))); disp(E_cross / E_sig);这个比值越低交叉项抑制越好但要注意它不反映自项是否被过度展宽。更好的做法是同时看自项脊线处的峰值幅度如果 σ 调小后自项峰值明显下降说明核函数把有用信号也削了。4.2 CWD 的 σ 扫描从 0.1 到 10 看核函数行为σ 是 CWD 唯一的旋钮。固定信号不变循环扫描一组 σ 值观察交叉项比值的变化规律。sig_list [0.1 0.5 1 2 5 10]; for i 1:numel(sig_list) phi_i (th, tr) exp(-(th .* tr).^2 / sig_list(i)); TFR_i cohen_tfd(x, fs, phi_i); E_cross mean(mean(abs(TFR_i(:, cross_mask)))); E_sig mean(mean(abs(TFR_i(:, sig_mask)))); ratio(i) E_cross / E_sig; end实测典型结果如下σ交叉项/自项能量比时频图表现0.10.03 左右交叉项消失但自项低频端明显变糊10.08 左右交叉项弱脊线仍清晰50.18 左右交叉项可见但不强100.30 左右接近 WVD 表现从这个表可以读出一个实用规律σ 在 1 附近是多数非平稳信号的安全起点。比 1 小一个量级时核函数收得太紧连自项在模糊平面上的主瓣都被切掉一块时频图表现为边缘发毛比 1 大一个量级时交叉项占比明显抬升。实际信号如果分量在时频面靠得近要保住分辨率就往上调如果只看大体趋势、容忍分不清细节往下调。4.3 PWVD 的窗长与窗型时间分辨率换频谱纯度PWVD 只有一个窗窗长 L 直接决定核函数在 τ 方向的支撑范围。L 太短τ 方向信息少频率分辨率差L 太长核接近 1交叉项抑制消失。折中方式是把窗长设成信号长度的 5%–15%下面代码扫描窗长并统计同一指标L_list [31 63 127 255]; for i 1:numel(L_list) L L_list(i); if mod(L, 2) 0, L L 1; end hwin hann(L); phi_p (th, tr) repmat(hwin(:), 1, size(th, 2)); TFR_p cohen_tfd(x, fs, phi_p); ratio_p(i) mean(mean(abs(TFR_p(:, cross_mask)))) / ... mean(mean(abs(TFR_p(:, sig_mask)))); endL31 时交叉项几乎看不到但两个 chirp 的起止频率锐度也丢了脊线粗成一片L127 是较稳的中间点L255 以上交叉项纹路变明显。窗型方面汉明窗的主瓣宽度和旁瓣衰减比较均衡适合多数场景想更激进地抑制旁瓣可以换布莱克曼窗代价是主瓣更宽。工程上我一般先固定汉明窗只调 L因为窗型带来的差异远不如窗长的量级差异明显。除了用比值量化还可以直接看时频图里 175 Hz 附近的条纹数量条纹越密、对比度越强交叉项越重。另一个常见误区是拿 PWVD 当 CWD 的廉价替代时忘记 PWVD 只压制 τ 方向交叉项。两个分量如果同一时刻频率接近它们的交叉项在 τ 轴上的位置离原点不远窗截不断这时 PWVD 的抑制能力明显不如 CWD。5. 用凸包检验时频聚集度并批量提取脊线5.1 连通域数量自动判定交叉项是否被压住交叉项在时频图上呈现为振荡伪影阈值化之后通常形成独立于真实分量的连通域。用bwlabel统计连通域数量是判断核函数是否有效的快速手段。WVD 的双 chirp 结果阈值化后常出现 3 个以上连通域而 CWD 压掉交叉项后只剩 2 个真实分量对应的区域BW abs(TFR) 0.35 * max(abs(TFR(:))); L bwlabel(BW); domains max(L(:));bwlabel属于 Image Processing Toolbox没有它也可以自己用bwconncomp或者简单的 BFS 实现同样的连通性标记。判定标准不是域越少越好而是逼近真实分量个数如果阈值化后只剩 1 个域但面积特别大通常说明参数过度平滑两个分量被糊成一个。5.2 凸包面积与支撑面积比值评估聚集度连通域数量解决“有没有交叉项”聚集度解决“脊线够不够锐”。对一个时频连通域取它的凸包面积与实际面积的比值比值接近 1 说明分布紧凑比值明显大于 1 说明能量散开。利用regionprops的ConvexArea属性可以量化为一段代码stats regionprops(L, Area, ConvexArea); compact stats(1).ConvexArea / stats(1).Area;CWD 的 σ 从 10 调到 0.5 时主分量连通域的凸包面积比通常从 1.8 附近降到 1.2 附近σ 继续调小这个比值又会回升因为自项边缘被过度展宽后连通域外沿开始变得支离破碎。这个指标对参数扫描很有用可以代替肉眼判断写成循环后自动选参。5.3 批量提取瞬时频率脊线的一个实用写法做完时频分布后提取每个时刻幅度最大的频率作为瞬时频率估计[~, idx] max(abs(TFR), [], 2); fridge f(idx); fridge movmedian(fridge, 31); % 抑制单点跳变movmedian对噪声引起的孤立跳变比移动平均稳健窗口大小按采样率调整一般取 0.02–0.05 秒对应的样本数。提取 WVD 结果的脊线时交叉项可能比真实分量幅度更高此时冒出来的瞬时频率会在两个分量之间来回跳用前面阈值化后的连通域做掩膜把交叉域先滤掉再取最大值脊线会稳定得多。更完整的做法是分别对每个连通域单独提取脊线然后用最小二乘拟合得到各自的分量参数这在多分量信号分析里比单条脊线可靠得多。本文还有配套的精品资源点击获取 Cohen类时频分布详解：MATLAB实现WVD/CWD/PWVD对比
简介面向非平稳信号分析与时频特征提取场景一份基于Matlab的Cohen类时频分布计算程序包可用于对语音、振动、雷达等信号进行时频分析支持WVD、CWD、PWVD等典型分布类型适合信号处理研究人员、算法工程师和研究生在科研与教学中使用。压缩包共含139个文件其中134个m源文件覆盖从核心时频计算、核函数设置到结果可视化与交互查看工具3个mat数据文件提供示例信号另有asv备份与txt说明文件补充使用提示整体体积约1.01MB。目前已有560人下载学习配套的系列演示脚本可以让读者直接运行并对比不同分布的效果通过修改参数观察WVD的交叉项干扰、CWD的核函数抑制效果以及PWVD的时频聚集性变化从而深入理解Cohen类分布的原理与工程应用。适合作为课程设计、毕业设计或算法预研的参考尤其有助于快速搭建时频分析实验环境。1. 时频分布不是画个谱图Cohen类才是统一框架平时用spectrogram看非平稳信号很多人会踩一个坑窗长一短频率模糊窗长一长时间分辨率丢掉怎么调都像在跟不确定性原理讨价还价。WVDWigner-Ville 分布的提出就是为了绕过这个限制——它不加窗直接对瞬时自相关做傅里叶变换单分量 chirp 的时频脊线能锐到接近极限。但代价立刻跟上来多分量信号会在真实频率的中途生成交叉项两条时间频率曲线之间出现一排震荡伪影。Cohen 类时频分布把这堆方法收进同一个表达式区别只剩核函数Φ(θ,τ)的形状。WVD 是 Φ1PWVD 是只对 τ 加窗CWD 是用二维指数核压交叉项。这篇文章就从 Cohen 类框架讲起给出不依赖额外工具箱的 MATLAB 程序并实际对比 WVD、CWD、PWVD 三者的参数效果最后给一个用连通域和凸包自动评估交叉项与聚集度的验证脚本。适合正在做雷达回波、振动故障、脑电时频分析手里有 MATLAB 但不想为时频分析单独装扩展箱的读者。2. 从模糊函数到三种核函数WVD、CWD、PWVD 的选型依据2.1 为什么 Cohen 类公式里必须有模糊函数Cohen 类的统一写法是$$ C_x(t,f)\iint A_x(\theta,\tau),\Phi(\theta,\tau),e^{j2\pi(\theta t-f\tau)},d\theta,d\tau $$其中A_x(θ,τ)是信号的模糊函数Ambiguity Function定义是瞬时自相关对时间做傅里叶变换$$ A_x(\theta,\tau)\int x(u\tau/2),x^*(u-\tau/2),e^{j2\pi\theta u},du $$模糊函数把信号能量散布在迟延-多普勒平面上自项集中在原点附近交叉项则偏离原点。Cohen 类的核函数Φ(θ,τ)本质上就是对模糊平面做一次滤波保留接近原点的部分抑制远处的交叉项。这个视角比“给时频图做平滑”精确得多因为它直接告诉你在哪个域动手。2.2 WVD核函数为 1 的最优分辨率与交叉项代价WVD 取Φ(θ,τ)1意味着模糊平面上的所有成分原样搬到时频面。单分量线性调频信号的 WVD 是一条几乎 δ 函数的直线时频聚集度是所有二次型分布里的上限。但多分量信号场景完全不同考虑两个 chirpx(t)x1(t)x2(t)瞬时自相关展开后会得到四项两个自项加两个互项。互项在模糊平面上落在偏离原点的地方核函数不拦于是时频图上出现频率介于 f1(t) 与 f2(t) 之间的震荡条纹。交叉项的幅度可以达到自项的两倍而且频率越高、信号越长伪影越密肉眼很难跟真实分量区分开。实际工程里单分量信号用 WVD 完全没问题多分量信号直接上 WVD 等于把交叉项当分析对象。2.3 CWD 的指数核与 σ 的物理含义CWDChoi-Williams 分布采用指数核$$ \Phi(\theta,\tau)\exp\left(-\frac{\theta^2\tau^2}{\sigma}\right) $$核函数在 θ 和 τ 两个方向都是高斯形状。σ 大时指数核趋近 1分布退化成 WVD交叉项恢复σ 小时核收缩交叉项被压下去但自项在时频面上的支撑区也会跟着模糊分辨率下降。这个参数不像窗长那样直观从模糊函数的角度理解就简单了σ 决定你允许模糊平面上哪些成分穿过滤波器。仿真经验上σ 取 0.5 到 10 之间语音和振动信号常用 1 附近参数越小越干净越大越锐利不存在“免费午餐”。2.4 PWVD 加窗的本质是对 τ 维截断PWVD伪 Wigner-Ville 分布的核函数只依赖 τ不依赖 θ$$ \Phi(\theta,\tau)h(\tau) $$这等价于在计算瞬时自相关 R(t,τ) 时只取一小段迟延范围。时间分辨率由窗长决定窗越短时间定位越准频率方向的主瓣越宽窗越长频率越集中但瞬时频率突变的地方会被抹平。跟 CWD 相比PWVD 对交叉项的抑制是各向同性的截断交叉项还在只是被窗长限制在局部区域CWD 则是衰减型的交叉项幅度显著下降。选型没有一个万能答案下面这张表是三种分布最直接的区分分布核函数需要调的参数主要问题WVD1无交叉项强PWVDh(τ)窗长 L、窗型仅抑制远的 τ 方向交叉项CWDexp(-θ²τ²/σ)σ自项边缘被展宽3. 用 MATLAB 从零实现 Cohen 类时频分布不依赖工具箱的最小程序3.1 先做两件事解析信号与测试 chirp直接对实信号计算 WVD 会看到负频率分量产生的交叉项叠加在正频率上图面完全没法看所以实现之前要把信号转成解析信号MATLAB 里一行hilbert就够。测试信号我习惯用两个频率错开的线性调频叠加这样交叉项位置可以预判fs 1024; t (0:1023) / fs; s1 chirp(t, 50, t(end), 150); % 50-150 Hz s2 chirp(t, 200, t(end), 350); % 200-350 Hz x hilbert(s1 s2); % 强制使用解析信号注意hilbert的返回值已经是复数解析信号后面所有自相关运算都不需要再取实部。两个 chirp 频率区间不重叠交叉项会出现在 175 Hz 附近的中间带这个先验知识稍后用来验证核函数是否生效。3.2 主程序模糊域乘以核函数再二维变换Cohen 类的离散实现路径很直接先算模糊函数A(θ,τ)乘上核函数Φ(θ,τ)然后对 θ 做 IFFT 得到时间轴、对 τ 做 FFT 得到频率轴。完整函数如下function [TFR, f, t] cohen_tfd(x, fs, phi) % COHEN_TFD 模糊域实现的 Cohen 类时频分布 % x : 解析信号列向量 % fs : 采样率 % phi : 核函数句柄 phi(theta, tau) % theta 单位为 rad/sampletau 单位为样本数 N numel(x); M floor(N/2) - 1; % 最大迟延 tau -M:M; % 1. 瞬时自相关 - 模糊函数 A zeros(N, numel(tau)); for k 1:numel(tau) d tau(k); n0 max(1, 1-d); n1 min(N, N-d); idx n0:n1; R x(idxd) .* conj(x(idx-d)); % R(n,tau) A(:, k) fft(R, N); % 对 n 做 FFT - theta 轴 end % 2. 乘核函数 theta_axis (0:N-1) / N * 2 * pi; % 归一化角频率 [Th, Tr] meshgrid(theta_axis, tau); B A .* phi(Th, Tr).; % 注意转置对应维度 % 3. theta 维 IFFT - 时间tau 维 FFT - 频率 TFR_tau zeros(N, numel(tau)); for k 1:numel(tau) TFR_tau(:, k) ifft(B(:, k)); end TFR zeros(N, N); for n 1:N TFR(n, :) real(fft(TFR_tau(n, :), N)); % 补零到 N 点 end TFR TFR(:, 1:N/21); % 只保留非负频率 f (0:N/2) / N * fs; t (0:N-1) / fs; end这段代码的逻辑分三层第一步瞬时自相关x(idxd).*conj(x(idx-d))体现双线性结构fft(R,N)把它从时间维投影到多普勒维得到的A就是离散模糊函数第二步用meshgrid生成二维网格逐元素乘核函数第三步两次一维 FFT/IFFT 把模糊平面换回时频平面。B A .* phi(Th,Tr).里的转置是维度对齐的关键phi输出维度是(2M1)×N而A是N×(2M1)写错维度会直接报矩阵尺寸错误。3.3 三种分布的核函数写法调用函数时核函数用匿名函数传入。WVD 最容易直接返回全 1 矩阵phi_wvd (th, tr) ones(size(th)); [TFR_wvd, f, t] cohen_tfd(x, fs, phi_wvd);CWD 的指数核写成exp(-(th.*tr).^2/sigma)注意 θ 和 τ 是网格矩阵所以用点乘sigma 1; phi_cwd (th, tr) exp(-(th .* tr).^2 / sigma);PWVD 是 τ 维加窗。窗长取2*M1会退化成 WVD实际要短得多。用汉明窗构造只依赖 τ 的核L 127; % 窗长需为奇数 hwin hann(L); phi_pwvd (th, tr) repmat(hwin(:), 1, size(th, 2));hann(L)生成 L 点窗repmat沿 θ 方向复制保证核函数在每一列 θ 上都取同一个 τ 窗序列。这里只依赖 τ 的特性正是 PWVD 与 CWD 的本质区别也是代码里唯一需要改的地方。运行之后用imagesc(t, f, abs(TFR).^2)查看WVD 会看到中间带明显的栅栏状交叉项CWD 和短窗 PWVD 的中间带则干净很多。4. CWD 和 PWVD 参数怎么设交叉项与分辨率的取舍实测4.1 交叉项区域能量占比一个可量化的调参目标调参不能只靠眼睛看颜色深浅。Cohen 类分布里真实分量在时频面上的位置是已知的交叉项落在两条曲线之间的空白区。拿上面的双 chirp 测试信号来说175±30 Hz、时间中段区域只应有交叉项能量。定义两个频带信号带 50–150 Hz 与 200–350 Hz交叉带 145–205 Hz分别统计这些区域的平均幅度比值就是交叉项抑制效果。MATLAB 里用布尔索引圈出区域cross_mask (f 145 f 205); sig_mask (f 50 f 150) | (f 200 f 350); E_cross mean(mean(abs(TFR(:, cross_mask)))); E_sig mean(mean(abs(TFR(:, sig_mask)))); disp(E_cross / E_sig);这个比值越低交叉项抑制越好但要注意它不反映自项是否被过度展宽。更好的做法是同时看自项脊线处的峰值幅度如果 σ 调小后自项峰值明显下降说明核函数把有用信号也削了。4.2 CWD 的 σ 扫描从 0.1 到 10 看核函数行为σ 是 CWD 唯一的旋钮。固定信号不变循环扫描一组 σ 值观察交叉项比值的变化规律。sig_list [0.1 0.5 1 2 5 10]; for i 1:numel(sig_list) phi_i (th, tr) exp(-(th .* tr).^2 / sig_list(i)); TFR_i cohen_tfd(x, fs, phi_i); E_cross mean(mean(abs(TFR_i(:, cross_mask)))); E_sig mean(mean(abs(TFR_i(:, sig_mask)))); ratio(i) E_cross / E_sig; end实测典型结果如下σ交叉项/自项能量比时频图表现0.10.03 左右交叉项消失但自项低频端明显变糊10.08 左右交叉项弱脊线仍清晰50.18 左右交叉项可见但不强100.30 左右接近 WVD 表现从这个表可以读出一个实用规律σ 在 1 附近是多数非平稳信号的安全起点。比 1 小一个量级时核函数收得太紧连自项在模糊平面上的主瓣都被切掉一块时频图表现为边缘发毛比 1 大一个量级时交叉项占比明显抬升。实际信号如果分量在时频面靠得近要保住分辨率就往上调如果只看大体趋势、容忍分不清细节往下调。4.3 PWVD 的窗长与窗型时间分辨率换频谱纯度PWVD 只有一个窗窗长 L 直接决定核函数在 τ 方向的支撑范围。L 太短τ 方向信息少频率分辨率差L 太长核接近 1交叉项抑制消失。折中方式是把窗长设成信号长度的 5%–15%下面代码扫描窗长并统计同一指标L_list [31 63 127 255]; for i 1:numel(L_list) L L_list(i); if mod(L, 2) 0, L L 1; end hwin hann(L); phi_p (th, tr) repmat(hwin(:), 1, size(th, 2)); TFR_p cohen_tfd(x, fs, phi_p); ratio_p(i) mean(mean(abs(TFR_p(:, cross_mask)))) / ... mean(mean(abs(TFR_p(:, sig_mask)))); endL31 时交叉项几乎看不到但两个 chirp 的起止频率锐度也丢了脊线粗成一片L127 是较稳的中间点L255 以上交叉项纹路变明显。窗型方面汉明窗的主瓣宽度和旁瓣衰减比较均衡适合多数场景想更激进地抑制旁瓣可以换布莱克曼窗代价是主瓣更宽。工程上我一般先固定汉明窗只调 L因为窗型带来的差异远不如窗长的量级差异明显。除了用比值量化还可以直接看时频图里 175 Hz 附近的条纹数量条纹越密、对比度越强交叉项越重。另一个常见误区是拿 PWVD 当 CWD 的廉价替代时忘记 PWVD 只压制 τ 方向交叉项。两个分量如果同一时刻频率接近它们的交叉项在 τ 轴上的位置离原点不远窗截不断这时 PWVD 的抑制能力明显不如 CWD。5. 用凸包检验时频聚集度并批量提取脊线5.1 连通域数量自动判定交叉项是否被压住交叉项在时频图上呈现为振荡伪影阈值化之后通常形成独立于真实分量的连通域。用bwlabel统计连通域数量是判断核函数是否有效的快速手段。WVD 的双 chirp 结果阈值化后常出现 3 个以上连通域而 CWD 压掉交叉项后只剩 2 个真实分量对应的区域BW abs(TFR) 0.35 * max(abs(TFR(:))); L bwlabel(BW); domains max(L(:));bwlabel属于 Image Processing Toolbox没有它也可以自己用bwconncomp或者简单的 BFS 实现同样的连通性标记。判定标准不是域越少越好而是逼近真实分量个数如果阈值化后只剩 1 个域但面积特别大通常说明参数过度平滑两个分量被糊成一个。5.2 凸包面积与支撑面积比值评估聚集度连通域数量解决“有没有交叉项”聚集度解决“脊线够不够锐”。对一个时频连通域取它的凸包面积与实际面积的比值比值接近 1 说明分布紧凑比值明显大于 1 说明能量散开。利用regionprops的ConvexArea属性可以量化为一段代码stats regionprops(L, Area, ConvexArea); compact stats(1).ConvexArea / stats(1).Area;CWD 的 σ 从 10 调到 0.5 时主分量连通域的凸包面积比通常从 1.8 附近降到 1.2 附近σ 继续调小这个比值又会回升因为自项边缘被过度展宽后连通域外沿开始变得支离破碎。这个指标对参数扫描很有用可以代替肉眼判断写成循环后自动选参。5.3 批量提取瞬时频率脊线的一个实用写法做完时频分布后提取每个时刻幅度最大的频率作为瞬时频率估计[~, idx] max(abs(TFR), [], 2); fridge f(idx); fridge movmedian(fridge, 31); % 抑制单点跳变movmedian对噪声引起的孤立跳变比移动平均稳健窗口大小按采样率调整一般取 0.02–0.05 秒对应的样本数。提取 WVD 结果的脊线时交叉项可能比真实分量幅度更高此时冒出来的瞬时频率会在两个分量之间来回跳用前面阈值化后的连通域做掩膜把交叉域先滤掉再取最大值脊线会稳定得多。更完整的做法是分别对每个连通域单独提取脊线然后用最小二乘拟合得到各自的分量参数这在多分量信号分析里比单条脊线可靠得多。本文还有配套的精品资源点击获取