基于MATLAB的Cohen类时频分布实现：WVD、CWD与PWVD核心解析

简介面向信号处理与时频分析方向的MATLAB程序资源包用于计算Cohen类时频分布支持WVDWigner-Ville分布、CWDChoi-Williams分布、PWVD伪Wigner-Ville分布等常用类型可满足非平稳信号分析中的不同需求适合科研人员、工程师及高年级本科生直接调用或二次开发。资源共139个文件以134个.m脚本为核心涵盖核心算法实现、多组演示程序及可视化辅助工具另配有少量.mat数据文件、.asv自动备份文件和.txt说明文本压缩包整体约1.01MB体量小巧便于快速获取与部署。目前已有560人学习下载具备一定参考热度。通过随包提供的实例演示可直观对比WVD、CWD、PWVD在时频聚集性与交叉项抑制方面的差异结合可视化界面工具还能交互式调节参数、观察不同窗函数和核函数对时频分布的影响有助于深入理解时频分析原理并迁移到实际信号处理任务中。1. 为什么自己写Cohen类时频分布而不是用现成函数当你面对一段非平稳信号比如语音、机械振动或脑电频谱分析只能告诉你信号里有哪几个频率却说不清这些频率是什么时候出现的。时频分布把时间与频率联合起来让瞬态、调频和跳变在同一张图上呈现。MATLAB里确实有spectrogram也能调用一些工具箱里的时频函数但那些函数大多把算法固定死了想看 Choi-Williams 分布 CWD或者想比较 WVD 和 PWVD 的交叉项抑制效果现成接口常常不给核函数的控制权。下面这个程序能指定 Cohen 类分布类型在 WVD、CWD、PWVD 之间切换。核心思路是把三种分布统一到同一个 Cohen 类框架下通过切换核函数决定输出哪一种。读完你不仅能跑出这三张时频图还能直接调 σ、窗长这些参数理解每一次改动在时频平面上发生了什么。适合刚接触时频分析的研究生也适合需要把核函数实验性改造的工程师。如果你还没装好 MATLAB先按常规的 matlab 下载安装教程准备好环境再继续。2. Cohen类时频分布的核函数WVD、CWD、PWVD的差别在哪2.1 从WVD到Cohen类核函数控制交叉项2.1.1 WVD的直觉与交叉项Wigner-Ville分布的定义是[ W(t,f)\int_{-\infty}^{\infty} z\left(t\frac{\tau}{2}\right)z^*\left(t-\frac{\tau}{2}\right)e^{-j2\pi f\tau}d\tau ]其中 (z(t)) 是解析信号。这个定义可以理解为在时间 (t) 附近把信号的左半段和右半段做相关再做傅里叶变换。因为积分里同时包含两段时间所以它对线性调频信号能给出教科书般的“倾斜直线”时频聚焦性远好于短时傅里叶变换。代价是双线性运算引入了交叉项如果信号包含两个分量比如两个正弦波时频平面上会出现一个频率在两个真分量之间位置的虚假振荡项幅度甚至在真实项之上。交叉项不是噪声而是双线性结构的数学产物。两个分量在时频平面上的几何中心位置必然出现交叉项它不会单独存在也无法用普通滤波去掉。Cohen类分布的出现就是为了在你可接受的时频分辨率损失范围内换掉这些交叉项。2.1.2 Cohen类统一公式Cohen类把WVD推广为[ C(t,f)\iint \phi(\xi,\tau) z\left(t-\xi\frac{\tau}{2}\right)z^*\left(t-\xi-\frac{\tau}{2}\right)e^{-j2\pi f\tau}d\xi d\tau ](\phi(\xi,\tau)) 就是核函数。当 (\phi(\xi,\tau)\delta(\xi)) 时退化为WVD。核函数本质上是对瞬时相关函数在时间 (\xi) 和时延 (\tau) 两个维度做加权平均或平滑。平滑会抹掉交叉项也会让自项的聚焦边缘变钝这就是整个Cohen类设计的权衡。WVD、CWD、PWVD三者可以用同一种 MATLAB 程序主干实现先计算瞬时相关矩阵再按核函数修改这个矩阵最后沿时延维做 FFT。差别只在核上这决定了程序结构可以高度统一。2.2 CWD与PWVD的核设计2.2.1 Choi-Williams核的指数衰减CWD的核在模糊域定义为[ \Phi_{\text{CWD}}(\theta,\tau)e^{-\theta^2\tau^2/\sigma} ](\theta) 是频率偏移(\tau) 是时延。(\sigma) 越大指数衰减越慢核越接近1分布越像WVD交叉项保留越多(\sigma) 越小平滑越强交叉项被压掉但自项会被展宽。实际经验是 (\sigma) 在0.5到10之间调整信号分量密集时用小 (\sigma)分量稀疏时用大 (\sigma)。注意这个核在模糊域是一个二维指数衰减函数所以要实现CWD通常需要对瞬时相关矩阵的时间维做 FFT得到模糊域乘上核再逆 FFT 回来。这也是后文代码里 CWD 分支单独处理的原因。2.2.2 PWVD的时间窗思路伪Wigner-Ville分布的做法更直接在时延方向加一个窗函数 (w(\tau))截断积分范围。[ \text{PWVD}(t,f)\int w(\tau)z\left(t\frac{\tau}{2}\right)z^*\left(t-\frac{\tau}{2}\right)e^{-j2\pi f\tau}d\tau ]等效于在瞬时相关矩阵的每一行乘以窗函数。窗长决定频率分辨率窗越长频率分辨率越高但时间方向上更多的非平稳信息被“平均”进来实际上增加了时间方向的模糊。所以PWVD没有完全消除交叉项只是把交叉项的干涉图样限制在窗内。代码里 PWVD 和 WVD 的差别就是一行 (R \cdot w) 的乘法正好体现了核函数可指定的意义。3. MATLAB程序实现统一Cohen类框架计算时频分布3.1 函数接口设计与参数表函数命名为cohen_tfd接口设计成function [TFR, f] cohen_tfd(x, fs, type, varargin)x是单通道信号fs是采样率type在wvd、cwd、pwvd之间选择。可选参数用inputParser统一解析避免调用时记参数顺序。参数名作用于含义默认值sigmacwd模糊域核的指数衰减系数1winlenpwvd时延方向窗长度会被强制为奇数127这个表格在代码里对应到addParameter。把sigma和winlen分开是因为两者物理含义完全不同sigma控制二维平滑力度winlen只控制时延方向截断。新手容易把它们混在一起调接口上分开能减少误操作。3.2 核函数离散化从连续公式到矩阵操作Cohen类计算的第一步是构建瞬时相关矩阵 (R)。设解析信号为 (z(n))(n1..N)定义时延索引 (m-M..M)则[ R(n,m)z(nm)z^*(n-m) ]这一步用双重循环写最直观边界用if判断。第二步按type对 (R) 做变换wvd(R) 保持原样。pwvd将窗函数 (w) 补零到 (R) 的列数然后每行乘以 (w)。cwd对 (R) 沿时间维做 FFT 得到模糊域 (A(\theta,m))乘以核 (\Phie^{-\theta^2m^2/\sigma})再逆 FFT。第三步对每个时刻 (n)将 (R(n,:)) 做 FFT 并取幅值平方得到时间-频率矩阵TFR。频率轴为f linspace(0, fs/2, nfft/2 1);这里有一个关键点CWD 的核是在模糊域相乘不是在时域卷积。很多简化实现会用二维高斯窗对 (R) 做卷积那是平滑伪WVD的近似严格说不是CWD。下面给的代码走的是模糊域路线和 Choi-Williams 原始定义一致所以当typecwd时交叉项抑制效果能真正体现。3.3 完整代码与调用示例完整函数如下保存为cohen_tfd.mfunction [TFR, f] cohen_tfd(x, fs, type, varargin) % cohen_tfd 计算Cohen类时频分布支持WVD、CWD、PWVD % 输入 % x 单通道信号列向量 % fs 采样率单位Hz % type wvd / cwd / pwvd % 可选参数 % sigma CWD核参数默认1 % winlen PWVD窗长默认127需为奇数 % 输出 % TFR 时频矩阵行为时间列为频率 % f 频率轴 x x(:); % 列向量 N length(x); z hilbert(x); % 解析信号抑制负频率伪影 % 解析可选参数 p inputParser; addParameter(p, sigma, 1); addParameter(p, winlen, min(127, N)); parse(p, varargin{:}); sigma p.Results.sigma; winlen p.Results.winlen; % 时延维设置 M floor((N-1)/2); % 最大时延采样数 m -M:M; nm length(m); winlen min(winlen, nm); if mod(winlen, 2) 0 winlen winlen - 1; % 强制奇数且不超过 nm end % --- 第一步瞬时相关矩阵 R(n, m) --- R zeros(N, nm); for n 1:N for id 1:nm mm m(id); % nmm 和 n-mm 必须同时落在 [1, N] 内 if (n mm 1) (n mm N) ... (n - mm 1) (n - mm N) R(n, id) z(n mm) * conj(z(n - mm)); else R(n, id) 0; end end end % --- 第二步按类型应用核函数 --- switch lower(type) case wvd % 核为1直接使用R case pwvd % 窗函数补零中心对准 m 0 w hamming(winlen); pad_left M - floor(winlen/2); w_pad zeros(1, nm); w_pad(pad_left (1:winlen)) w; R R .* w_pad; case cwd % 模糊域乘Choi-Williams核 A fft(R, N, 1); % 时间维FFT到模糊域 A fftshift(A, 1); % 将零频移到中间 theta ((0:N-1) - floor(N/2)). / N * 2 * pi; Phi exp(-(theta.^2) * (m.^2) / sigma); % 核矩阵 N x nm A A .* Phi; A ifftshift(A, 1); R ifft(A, N, 1); % 回到时间-时延域 otherwise error(type只支持 wvd / cwd / pwvd); end % --- 第三步沿时延维FFT得到时频分布 --- nfft 2^nextpow2(nm); TFR zeros(N, nfft/2 1); for n 1:N spec fft(R(n, :), nfft); TFR(n, :) abs(spec(1:nfft/2 1)).^2; % 功率谱显示更稳定 end f linspace(0, fs/2, nfft/2 1); end代码逻辑说明瞬时相关矩阵 (R) 的行索引是时间 (n)列索引是时延 (m)。双重循环里的边界判断保证所有索引都落在信号范围内越界部分补零对应信号两端用零延拓。PWVD 分支中hamming(winlen)是 MATLAB 自带窗补零后与 (R) 逐行相乘等效于只截取中心窗内的相关函数。CWD 分支的fftshift和ifftshift配对使用是为了把 FFT 输出的 (0) 到 (2\pi) 频率轴转换成 (-\pi) 到 (\pi)否则核函数在周期边界上不连续时频图上会出现明显的横向条纹。调用示例fs 1000; t 0:1/fs:1; x chirp(t, 50, t(end), 250, quadratic); % 50Hz到250Hz二次调频 x x sin(2*pi*80*t); % 叠加80Hz正弦 [TFR, f] cohen_tfd(x, fs, cwd, sigma, 1); imagesc(t, f, TFR); axis xy; xlabel(时间 (s)); ylabel(频率 (Hz));这里imagesc的横轴直接用时间向量纵轴用f。axis xy把 y 轴方向反转让频率低的显示在下方符合大多数人对频谱图的阅读习惯。运行后能看到一条弯曲的调频轨迹和一条 80Hz 的水平线CWD 的交叉项应该比直接用 WVD 淡得多。4. 参数调优与边界效应让核函数真正工作起来4.1 CWD核参数σ平滑强度与自项保留的平衡sigma是 CWD 唯一需要调的参数。它在模糊域核 (\Phie^{-\theta^2m^2/\sigma}) 里控制衰减速度。sigma越大(\Phi) 在高 (\theta) 区域也接近1模糊域被保存的信息越多输出越接近 WVDsigma越小核把高 (\theta) 成分全部压掉交叉项消失但自项的瞬时频率边缘会被展宽。在代码里改一个数就能看到变化[TFR1, f] cohen_tfd(x, fs, cwd, sigma, 0.2); [TFR2, ~] cohen_tfd(x, fs, cwd, sigma, 10);sigma取0.2时图面上交叉项区域几乎变成一片均匀底色真实分量的轨迹边缘也明显变粗。sigma取10时交叉项重新出现但轨迹线细得像用 WVD 画的。工程上一般先用sigma1跑一遍若交叉项干扰判读就往0.5方向调若能接受少量交叉项而需要更锐利的时频峰就往2到5方向调。注意sigma不是越大越好超过10后 CWD 基本失去平滑意义退化成带有轻微模糊的 WVD。4.2 PWVD窗长时频分辨率的此消彼长PWVD 的窗长winlen直接影响时延维长度 (nm)。窗越长时延方向的截断越小频率分辨率越高但每个时刻会混入更多邻域时间的信息窗越短时间定位越好但频率分辨率下降。这个权衡可以做成表格winlen频率分辨率时间定位交叉项表现31低轨迹粗好瞬变清楚局部振荡仍明显127中中交叉项集中在短窗内511高差强非平稳信号会糊交叉项幅度变大代码中winlen被强制为奇数是因为窗函数的中心要放在时延 (m0) 处。如果给偶数程序会-1改成奇数你看到的实际生效窗长可能与期望差一个点。建议在脚本里加一行disp(winlen)确认最终值。4.3 三个容易踩的坑第一个坑是忘了解析信号。直接用原始信号 (x) 计算 (R)会在频率轴零频附近出现镜像分量原因是实信号的频谱关于 0Hz 对称双线性结构会产生额外交叉项。代码里用hilbert获取解析信号目的就是去掉负频率只保留单边谱。第二个坑是边界效应。(R) 矩阵在信号两端是被零填充的所以 TFR 图的左右边缘会出现暗色条纹这是数据缺失造成的伪迹不是真实时频特性。处理方法是先给信号做两端延拓再调用本函数。第三个坑是输出取绝对值。WVD 和 CWD 在自项位置并不总是正值取实部才能看到负峰结构但负值在显示上不直观。上面的代码用abs(...).^2显示功率适合看能量分布如果你要分析 WVD 的负值特性把最后一行的输出改为real(spec(1:...))即可。5. 用LFM叠加正弦验证三种分布的实际效果5.1 构造验证信号并运行程序验证信号用线性调频加正弦叠加线性调频覆盖 50Hz 到 250Hz正弦固定在 140Hz采样率 1000Hz时长1秒。fs 1000; t 0:1/fs:1; x chirp(t, 50, 1, 250); % 线性调频 x x 0.8 * sin(2*pi*140*t); types {wvd, pwvd, cwd}; names {WVD, PWVD, CWD}; for k 1:3 [TFR, f] cohen_tfd(x, fs, types{k}, ... sigma, 1, winlen, 151); figure; imagesc(t, f, TFR); axis xy; xlabel(时间 (s)); ylabel(频率 (Hz)); title(names{k}); end5.2 读图交叉项、自项形状与核的作用WVD 图上你能找到一条斜线和一个水平线但在这两条线的中间位置会出现一个“第三个分量”频率在两者的动态中点附近这就是交叉项。交叉项以规则条纹出现不是随机的所以需要用核函数去平滑。PWVD 图上条纹变成短促的局部振荡原因是窗函数把时延方向截断了但交叉项并没有消失只是在每个时刻呈现为短时波动。CWD 图的交叉项区域最干净两条真分量的轨迹也略微模糊。这正好对应三种核的设计目标WVD 追求锐度PWVD 用窗换取局部时间精度CWD 用二维平滑换取全局清洁。5.3 与spectrogram的数值对照验证频率轴spectrogram是 MATLAB 自带的短时傅里叶变换它的频率轴和我们的f轴不一定完全一致但峰值位置应当对应。做一个数值对照[S, ~, ~] spectrogram(x, hamming(128), 64, 256, fs); [~, idx1] max(sum(abs(S), 2)); % 平均谱峰值索引 [~, idx2] max(sum(TFR, 1)); f1 idx1 / 256 * fs / 2; % spectrogram的峰值频率估算 f2 f(idx2); % 本程序的峰值频率估算这个技巧用来检查频率轴的单位和零点是否一致。如果f2与f1差了一大截问题通常出在nfft与时延列数的对应关系上。参考代码中nfft 2^nextpow2(nm)改动时需要同步检查频率轴计算。另外spectrogram窗口长度为128时时间分辨率约为0.128秒而cohen_tfd在winlen151时的等效窗长接近0.151秒两者分辨率相近正好可以互相印证调参结论。本文还有配套的精品资源点击获取