企业官网建设流程全解析

简介：SPWVD时频分析工具包是一份面向信号处理与故障诊断等场景的MATLAB源码资源，适合需要分析非平稳信号局部时频特征的研究人员与工程师。该工具包的核心是平滑伪维格纳-维尔分布实现代码，通过加窗平滑有效抑制传统维格纳-维尔分布的交叉项干扰，同时保持较好的时频聚集性。包内共1个m文件，压缩后仅1KB，属于轻量级算法函数；用户可直接调用或在此基础上按需修改窗函数、平滑长度等参数，用于语音识别、生物医学信号处理、旋转机械故障诊断等任务。目前已有1172人浏览学习，实用性获得一定认可。对希望避开复杂理论推导、快速上手SPWVD并借助MATLAB完成时频图绘制的读者而言，这份资源可节省大量编码时间，并可通过阅读源码理解预处理、PWVD计算、平滑及显示等完整流程，便于后续算法改进与应用开发。

1. SPWVD 时频分析：先看懂它凭什么压制交叉项

一个 10 秒的振动信号，里面有变转速的齿轮箱频率分量，还有间歇出现的轴承冲击。用短时傅里叶变换看，冲击会被窗长拉宽；直接用 WVD 看，分量之间又冒出一堆莫名的交叉项。这时你能想到的最实用的解法就是平滑伪 Wigner-Ville 分布（SPWVD）。SPWVD 时频分析是 MATLAB 信号处理里被反复提及却又很少有人讲透的工具：它在 WVD 的时频平面上同时做两个方向的平滑，把交叉项压下去，同时保留相对较高的时频聚集性。这篇内容适合正在做非平稳信号分析、故障监测或者雷达回波时频特征提取的工程师。我不会堆公式，但会把算法拆开，给出一个不依赖第三方工具箱也能跑的 MATLAB 实现。

2. SPWVD 时频分析原理：从 WVD 的双线性结构到可分离平滑

2.1 双线性时频分布与交叉项的来源

Wigner-Ville 分布的定义是单时间点双线性运算，对任意信号 (x(t)) 有

[ W_x(t, f) = \int_{-\infty}^{\infty} x\left(t+\frac{\tau}{2}\right) x^*\left(t-\frac{\tau}{2}\right) e^{-j2\pi f\tau}, d\tau ]

因为信号在其本身和自身的延时副本之间做相关，所以 WVD 没有窗口限制，时频聚集性理论上接近极限，脉冲、扫频都能在时频图上形成很细的脊线。但问题恰恰出在“双线性”这三个字上。如果信号由两个分量组成 (x(t) = x_1(t) + x_2(t))，展开时频分布会得到：

[ W_x(t,f) = W_{x_1}(t,f) + W_{x_2}(t,f) + 2\operatorname{Re}\left[ W_{x_1, x_2}(t,f) \right] ]

第三项就是交叉项。它出现在两个自项的坐标中点，幅度可以比自项还大，而且是振荡的，看起来像“幽灵分量”。在 MATLAB 里直接对一段多分量信号运行未平滑的 WVD，时频图中会出现明显的棋盘状干扰，和真实谱峰混在一起，导致瞬时频率提取完全不可用。因此单独使用 WVD 的场景极其有限，工程上几乎都会做某种平滑处理。

2.2 平滑核：为什么分成两个窗

平滑伪 Wigner-Ville 分布的基本思路是用一个二维低通滤波器对 WVD 结果做卷积：

[ \text{SPWV}x(t,f) = \int{-\infty}^{\infty}\int_{-\infty}^{\infty} G(t-t', f-f') , W_x(t',f'), dt', df' ]

如果二维核 (G(t,f)) 可以分解成两个一维窗 (g(t)) 和 (h(f)) 的乘积，那么平滑过程就能拆成先在时间轴上做一次卷积、再在频率轴上做一次卷积，计算复杂度低很多，参数也更直观。(g(t)) 负责抑制 WVD 在时间方向的振荡，平滑掉交叉项的规律性起伏；(h(f)) 负责在频率方向压低旁瓣和残留干扰。这个分解就是 SPWVD 最常被实现的形式。

在实际离散实现中，(g(t)) 用一段短时窗，中心点对应当前时刻，窗长越长，时间方向越平滑，但瞬时冲击会被拉宽；(h(f)) 同样用一段窗，窗长越长，频率方向越平滑，但两个靠得很近的频率分量会混在一起。时频聚集性和交叉项抑制程度本质上是矛盾的，SPWVD 能做的只是让用户在两个维度分别权衡，而不是像普通平滑伪分布那样只有一个全局参数。

2.3 WVD、PWVD、SPWVD 三者的对比与选型

一般如果我需要做故障诊断或雷达多普勒分析，会直接跳过 PWVD，因为多分量背景下它保留的交叉项依然明显。SPWVD 是性价比最高的起点，先固定一组中等窗长，比如频率窗 31 点、时间窗 31 点，再看结果决定是否加长或缩短。

3. 用 MATLAB 实现 SPWVD 时频分析：可运行的平滑伪 WVD 函数

3.1 实现思路：先用 WVD，再做可分离二维卷积

最常见的实现方式是用 for 循环逐时间点构造局部自相关，然后做 FFT 得到 WVD 频率切片。计算完 WVD 后，再把平滑核分别沿频率轴和时间轴做conv2卷积。这种思路的好处是把 WVD 和 SPWVD 拆成两个独立阶段，方便调试：先确认 WVD 有交叉项，再确认加卷积后交叉项被压掉。

另一个必须做的步骤是取解析信号。把实信号通过hilbert转成解析信号后，负频率成分为零，WVD 的交叉项会明显减少，而且自项幅度更集中。若不取解析信号，实信号的 WVD 会在负频率产生镜像项，平滑后依然会污染时频图。

3.2 myspwvd 函数完整代码

function [tfr, f] = myspwvd(x, fs, nfreq, flen, tlen) % myspwvd - MATLAB 实现平滑伪 Wigner-Ville 分布 % 输入: % x : 实信号行向量或列向量 % fs : 采样率 (Hz) % nfreq : 频率点数，建议取 256/512 等 2 的幂 % flen : 频率方向平滑窗长度，奇数 % tlen : 时间方向平滑窗长度，奇数 % 输出: % tfr : N x (nfreq/2) 的时频矩阵，N 为信号长度 % f : 频率轴向量 x = x(:).'; N = length(x); if mod(nfreq, 2) ~= 0 nfreq = nfreq + 1; end half = nfreq / 2; % 解析信号，抑制负频率镜像项 xa = hilbert(x); % ---- 第一步：计算 WVD 正频率部分 ---- wvd = zeros(N, half); for tidx = 1:N % 滞后 m 的范围受信号边界和 FFT 点数限制 m_max = min([tidx - 1, N - tidx, half - 1]); if m_max < 0 continue; end lag = (-m_max:m_max).'; % 列向量 pos = half + lag + 1; % 零滞后对应位置 half+1 % 双线性自相关 R = xa(tidx + lag) .* conj(xa(tidx - lag)); % 填充到长度为 nfreq 的序列，未填充处为零 Rw = zeros(nfreq, 1); Rw(pos) = R; % FFT 并移位，取正频率部分 spec = fftshift(fft(Rw)); wvd(tidx, :) = spec(half + 1:end).'; end wvd = real(wvd); % ---- 第二步：可分离二维平滑 = SPWVD ---- if mod(flen, 2) == 0, flen = flen + 1; end if mod(tlen, 2) == 0, tlen = tlen + 1; end g = hann(flen, 'periodic'); % 频率平滑窗 h = hann(tlen, 'periodic'); % 时间平滑窗 g = g(:)./sum(g); h = h(:)./sum(h); % 沿频率方向卷积：核必须是行向量 tfr = conv2(wvd, g.', 'same'); % 沿时间方向卷积：核必须是列向量 tfr = conv2(tfr, h, 'same'); % 频率轴 if nargout > 1 f = (0:half-1) * fs / nfreq; end end

这段代码把 SPWVD 的计算拆成清晰的三个阶段。hilbert(x)将实数信号转为解析信号，这样后续自相关只保留正频率；Rw的长度等于nfreq，滞后为零时放在half+1位置，未填充的滞后位置自动补零，这与标准 FFT 方法要求的一致。第二步的conv2(wvd, g.', 'same')用行向量g.'沿频率方向平滑，每列的频率谱线被窗口加权后叠加；接着用列向量h沿时间方向平滑，最终输出就是同时经过时间方向和频率方向平滑的 SPWVD。

参数方面，nfreq直接决定频率分辨率，也影响 FFT 计算量。flen太大会把两个相近的频率峰合并，太小平滑不足；tlen太大会把瞬态冲击抹平，太小则交叉项压不干净。我一般先取nfreq=512, flen=31, tlen=31，再观察结果微调。

3.3 参数怎么定：频率点数、平滑窗长与信号长度的关系

• 频率点数nfreq：通常设置为信号长度的 2 的幂，不超过信号长度。nfreq=256适合短数据，nfreq=1024适合 1 秒左右、采样率在 kHz 以上的信号。
• 频率平滑窗flen：需要在交叉项抑制与频率分辨率之间折中。对 60 Hz 和 120 Hz 两个分量，flen超过 65 就可能让两个峰粘连。
• 时间平滑窗tlen：影响瞬时变化能力。冲击、咳嗽、爆破音等非平稳事件，tlen应小于信号平稳段长度的三分之一。
• 信号过短时，边界处 WVD 本身误差大，平滑窗超出边界后的结果不可信，建议丢弃时频图两侧各tlen/2个时间点。

提示：conv2的'same'选项会让边缘处的卷积窗口自动截断，导致图中左右边缘的幅度略低。判断参数效果时，只看时频图中间 80% 的区域。

4. SPWVD 时频分析实战：线性调频和正弦叠加信号的参数对比

4.1 生成测试信号并调用 myspwvd

这里我构造一个多分量测试信号，由一段线性调频、一个固定频率正弦和少量白噪声组成。线性调频的瞬时频率随时间线性增加，正弦分量则保持恒频，两者叠加能同时检验 SPWVD 的跟踪能力和交叉项抑制效果。

fs = 1024; % 采样率 1024 Hz t = (0:511) / fs; % 0.5 秒信号 x = chirp(t, 50, t(end), 200, 'linear') ... + sin(2*pi*120*t) ... + 0.3 * randn(size(t)); % 直接调用上面的 myspwvd 函数 [tfr, f] = myspwvd(x, fs, 512, 31, 31); % 画图：时频平面 imagesc(t, f, abs(tfr)); axis xy; xlabel('时间 (s)'); ylabel('频率 (Hz)'); title('SPWVD 时频分布 (flen=31, tlen=31)'); colorbar;

chirp函数生成 50 Hz 到 200 Hz 的线性扫频信号，sin(2*pi*120*t)固定正弦频率是 120 Hz。myspwvd输出的矩阵行数是频率点数的一半，这里正好是 256 条频率线，从 0 到 512 Hz（采样率 1024 的一半）。imagesc纵轴顺序默认从低到高，所以需要axis xy让频率轴按常规方向显示。

运行后能看到时频图中有一条从左下向右上扫过的斜线，同时在中间高度有一条水平亮线。由于加入了 0.3 倍幅值的白噪声，时频图中会有细小颗粒噪点，但交叉项的条状干扰应该已经不明显。

4.2 不同平滑窗长对时频图的影响

要理解 SPWVD 的边界，可以直接对比不同flen和tlen下的输出。我用同一信号，分别跑三组参数，然后比较交叉项的可见程度和扫频脊线宽度。

figure; subplot(1,3,1); [tfr1, f1] = myspwvd(x, fs, 512, 11, 31); imagesc(t, f1, abs(tfr1)); axis xy; title('flen=11, tlen=31'); subplot(1,3,2); [tfr2, f2] = myspwvd(x, fs, 512, 31, 31); imagesc(t, f2, abs(tfr2)); axis xy; title('flen=31, tlen=31'); subplot(1,3,3); [tfr3, f3] = myspwvd(x, fs, 512, 61, 31); imagesc(t, f3, abs(tfr3)); axis xy; title('flen=61, tlen=31');

从表格可以得出，flen对交叉项影响最直接，因为频率方向的平滑相当于把 WVD 在频率域的低通能力放大。flen从 11 增到 61，交叉项确实消失，但线性调频的斜边也变毛糙了，看起来像有一层雾。使用时优先调小flen解决混叠，再用tlen处理时间方向的多余纹波。

4.3 为什么 SPWVD 会比 STFT 更“锐利”

短时傅里叶变换（STFT）等效于对信号加窗后做傅里叶变换，窗函数一旦选定，时间分辨率与频率分辨率的乘积固定受不确定性原理约束，而且窗内信号被强制看成平稳。SPWVD 的平滑是在双线性谱上进行的，它先计算理想聚集的 WVD，再做后置平滑，因此理论上其频率跟踪速度比 STFT 快得多。举个例子，对线性调频信号，STFT 需要把窗长控制在极短范围内才能避免扫频引起频谱展宽，但窗一短，频率分辨率又下降。SPWVD 使用后置平滑，能够在不直接牺牲瞬时频变响应的情况下压掉交叉项，这是它能在齿轮箱故障诊断、水下目标识别等场景中替代 STFT 的主要原因。

5. 验证 SPWVD 结果的三个技巧与常见坑

5.1 用瞬时频率脊线验证时频分布正确性

拿到时频矩阵后，第一件事不是美化颜色映射，而是验证峰值位置是否和理论瞬时频率一致。对上面生成的测试信号，理论瞬时频率为：

[ f_{\text{inst}}(t) = \begin{cases} 50 + (200-50)/0.5 \cdot t, & t \le 0.5s \ 120, & \text{恒定分量} \end{cases} ]

在 MATLAB 里可以把时频矩阵每一列的最大值索引转成频率，再和理论曲线对比：

[~, idx_max] = max(abs(tfr), [], 1); f_est = f(idx_max); % 理论瞬时频率 f_theory = 50 + (200-50) * t; % 线性调频分量 plot(t, f_est, 'o'); hold on; plot(t, f_theory, 'r-'); xlabel('时间 (s)'); ylabel('瞬时频率 (Hz)'); legend('SPWVD 峰值', '理论值');

如果时频分布正确，峰值轨迹应该贴着理论直线走。出现系统性偏移，通常是因为频率轴起点不是 0 Hz；出现明显分岔，则是交叉项残留或窗长过大导致两个分量融合。

5.2 频点宽度、边缘振荡和内存占用的处理

使用 SPWVD 时常见三个坑都与参数直接相关。第一个是频率轴宽度选得太小，导致时频图只有少量像素，两个分量在图上无法区分。把nfreq设置为信号长度的 2 的幂，一般能避免这个问。

第二个是边缘振荡。WVD 本身要求信号无限长，截断后首尾时段的滞后项估计不准，经过conv2后的边缘也会被窗截断。我一般只取时频图中央 70%-80% 的区域做脊线提取，边缘直接丢弃。

第三个是内存。信号长度 2 万点、nfreq=4096时，时频矩阵就达到 8000 万复数元素，占内存约 1.3 GB。如果内存吃紧，建议对信号分段处理，每段长度为 4096 点左右，段与段重叠 50%，再把时频结果拼起来。注意段边界处要去掉tlen个点，避免拼接痕迹。对长信号这种策略比一次算完全部要稳得多。

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