企业官网建设流程全解析

简介：面向非平稳信号分析与时频特征提取场景，一份基于Matlab的Cohen类时频分布计算程序包，可用于对语音、振动、雷达等信号进行时频分析，支持WVD、CWD、PWVD等典型分布类型，适合信号处理研究人员、算法工程师和研究生在科研与教学中使用。压缩包共含139个文件，其中134个m源文件覆盖从核心时频计算、核函数设置到结果可视化与交互查看工具，3个mat数据文件提供示例信号，另有asv备份与txt说明文件补充使用提示，整体体积约1.01MB。目前已有560人下载学习，配套的系列演示脚本可以让读者直接运行并对比不同分布的效果，通过修改参数观察WVD的交叉项干扰、CWD的核函数抑制效果以及PWVD的时频聚集性变化，从而深入理解Cohen类分布的原理与工程应用。适合作为课程设计、毕业设计或算法预研的参考，尤其有助于快速搭建时频分析实验环境。

1. 时频分布不是画个谱图，Cohen类才是统一框架

平时用spectrogram看非平稳信号，很多人会踩一个坑：窗长一短频率模糊，窗长一长时间分辨率丢掉，怎么调都像在跟不确定性原理讨价还价。WVD（Wigner-Ville 分布）的提出就是为了绕过这个限制——它不加窗，直接对瞬时自相关做傅里叶变换，单分量 chirp 的时频脊线能锐到接近极限。但代价立刻跟上来：多分量信号会在真实频率的中途生成交叉项，两条时间频率曲线之间出现一排震荡伪影。Cohen 类时频分布把这堆方法收进同一个表达式，区别只剩核函数Φ(θ,τ)的形状。WVD 是 Φ=1，PWVD 是只对 τ 加窗，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 是一条几乎 δ 函数的直线，时频聚集度是所有二次型分布里的上限。但多分量信号场景完全不同，考虑两个 chirp：x(t)=x1(t)+x2(t)，瞬时自相关展开后会得到四项，两个自项加两个互项。互项在模糊平面上落在偏离原点的地方，核函数不拦，于是时频图上出现频率介于 f1(t) 与 f2(t) 之间的震荡条纹。交叉项的幅度可以达到自项的两倍，而且频率越高、信号越长，伪影越密，肉眼很难跟真实分量区分开。实际工程里，单分量信号用 WVD 完全没问题，多分量信号直接上 WVD 等于把交叉项当分析对象。

2.3 CWD 的指数核与 σ 的物理含义

CWD（Choi-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 则是衰减型的，交叉项幅度显著下降。选型没有一个万能答案，下面这张表是三种分布最直接的区分：

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/sample，tau 单位为样本数 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(idx+d) .* 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/2+1); % 只保留非负频率 f = (0:N/2) / N * fs; t = (0:N-1) / fs; end

这段代码的逻辑分三层：第一步瞬时自相关x(idx+d).*conj(x(idx-d))体现双线性结构，fft(R,N)把它从时间维投影到多普勒维，得到的A就是离散模糊函数；第二步用meshgrid生成二维网格，逐元素乘核函数；第三步两次一维 FFT/IFFT 把模糊平面换回时频平面。B = A .* phi(Th,Tr).'里的转置是维度对齐的关键，phi输出维度是(2M+1)×N，而A是N×(2M+1)，写错维度会直接报矩阵尺寸错误。

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*M+1会退化成 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

实测典型结果如下：

从这个表可以读出一个实用规律：σ 在 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)))); end

L=31 时交叉项几乎看不到，但两个 chirp 的起止频率锐度也丢了，脊线粗成一片；L=127 是较稳的中间点；L=255 以上交叉项纹路变明显。窗型方面，汉明窗的主瓣宽度和旁瓣衰减比较均衡，适合多数场景；想更激进地抑制旁瓣可以换布莱克曼窗，代价是主瓣更宽。工程上我一般先固定汉明窗，只调 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 结果的脊线时，交叉项可能比真实分量幅度更高，此时冒出来的瞬时频率会在两个分量之间来回跳；用前面阈值化后的连通域做掩膜，把交叉域先滤掉再取最大值，脊线会稳定得多。更完整的做法是分别对每个连通域单独提取脊线，然后用最小二乘拟合得到各自的分量参数，这在多分量信号分析里比单条脊线可靠得多。

本文还有配套的精品资源，点击获取