企业官网建设流程全解析

简介：利用FFT计算非平稳随机信号的WVD分布，是一份适合信号处理与时频分析学习者、MATLAB使用者的实用仿真资料。资源提供完整MATLAB程序fft_WD.m，演示基于FFT的Wigner-Ville分布计算流程，运行后可得到二维与三维WVD分布图像，直观展示非平稳信号的时频聚集特征。压缩包内共4个文件，包括m脚本、操作录像avi以及两张结果预览jpg，整体大小约3.1MB，文件结构简洁，便于快速部署与验证。配套录像展示了在MATLAB 2021a中的完整操作与当前文件夹路径设置等注意事项，可有效降低入门门槛。该资源已有533人学习，适合需要理解WVD原理、完成课程实验或进行算法复现的读者参考。

1. 从非平稳随机信号到WVD：直接做FFT为什么不够

处理振动台采集的结构响应、语音或雷达回波时，频率成分会随时间移动。直接用FFT做一次频谱分析，只能得到一个时间段的平均能量分布，突发冲击、频率爬升这类细节会被抹平。非平稳随机信号需要的是“频率随时间怎么变”，也就是时频分布。

Wigner-Ville分布（WVD）在这个问题里是绕不开的候选。它对单分量信号有极高的时频聚集性，时频谱线的分辨率远超短时傅里叶变换；但代价是二次型分布特有的交叉项，一组频率分量之间会出现“幽灵能量”。很多人第一次用FFT实现WVD，看到时频图里出现原本不存在的条纹，第一反应是程序写错了，其实更可能是没做解析信号预处理，或者没有加平滑窗。这篇文章按“公式离散化 → MATLAB实现 → 参数调优 → 录像记录”的顺序，讲清楚怎么用FFT把WVD跑起来，以及哪些现象是真实算法特性，哪些是数值错误。

2. 用FFT实现WVD的数学原理与离散化路径

2.1 WVD的定义：从瞬时自相关到频谱

连续信号x(t)的WVD定义是：

W_x(t,f) = ∫ x(t+τ/2)·x*(t-τ/2)·e^{-j2πfτ} dτ

这个式子可以拆成两步来理解。先把x(t)在时间t处做“对称配对”：取t前面和后面各τ/2距离的两个样点，相乘得到R(t,τ) = x(t+τ/2)·x*(t-τ/2)，这相当于对t时刻的“局部相似性”做了一次瞬时自相关；然后对这个延时变量τ做傅里叶变换，得到的就是t时刻的频率切片。

关键区别在这里：短时傅里叶变换是先把信号乘窗函数加窗，再做FFT；WVD是对自相关核R(t,τ)做FFT。自相关核已经显式包含过去了和未来来的信号，所以WVD天然拥有时间方向的“干涉”能力，也正因为这种干涉，多分量信号会产生交叉项。视频里看到的时频图，横轴是时间t，纵轴是频率f，颜色深浅对应W_x(t,f)的幅度。

对非平稳随机信号来说，x(t)本身是随机过程的一次样本实现，WVD依然可以逐时刻计算，不要求信号平稳。这个特性正是它能替代FFT做非平稳分析的根本原因。

2.2 离散化：为什么能变成逐时刻的FFT

实际信号是采样得到的离散序列x[n]，采样周期T_s = 1/f_s。把连续定义离散化，延时τ用整数m表示，t对应整数n，那么延时域的自相关核变成：

R_n[m] = x[n+m]·x*[n-m]

对这串序列按m做DFT，就得到n时刻的频率切片。DFT用FFT实现，复杂度是O(N_FFT log N_FFT)，对每个时刻重复计算，总复杂度为O(N·N_FFT log N_FFT)，在PC上处理几十万点完全可行。

离散化时有一个隐藏约束：m的取值范围如果取-M到M，那么R_n[m]的长度是2M+1，FFT点数N_FFT必须不小于2M+1，否则会因为序列截断丢信息。更常见的做法是N_FFT取2的幂且大于2M+1，比如M=128时N_FFT选512或1024，相当于做了零填充，频域插值更密，频谱看起来更平滑。

FFT输出的顺序是0到N_FFT-1，对应频率从0到f_s到(N_FFT-1)/N_FFT。WVD的核函数对称，真实的频谱以f_s/2为中心折叠，所以大部分实现会在FFT之后做fftshift，把零频移到中间，再映射到频率轴(-f_s/2, f_s/2]。这一段映射关系是仿真录像参数区最常出错的地方，后面代码里会专门处理。

2.3 解析信号为什么必须提前处理

直接用实信号x[n]算WVD，交叉项会出现在正负频率之间。实信号的频谱关于零频对称，WVD里(t, f)处的能量和(t, -f)处相互干涉，时频图在零频附近会出现强烈的虚假分量。解决办法是先用希尔伯特变换构造解析信号：

z[n] = x[n] + j·H{x[n]}

H表示希尔伯特变换，z[n]的频谱只有正频率分量，负频率被清零，这样WVD只剩单边谱，交叉项也被压制了一部分。MATLAB的hilbert函数返回值本身就是复数解析信号，不是实数包络，直接用即可。如果信号本来就是复数基带信号，则不需要再做这一步。

这里要区分两个概念：取解析信号是为了避免±f交叉项，不是为了让信号“更平滑”。很多初学者看到hilbert输出感觉有点意外，就直接取实部用，反而把关键一步弄丢了。

2.4 不同变体怎么选：WVD、伪WVD、平滑伪WVD

仅基于标题“利用FFT计算非平稳随机信号的WVD分布”，实战中第一版先用PWVD，把时间、频率分辨率调到肉眼可接受，再决定要不要上SPWVD。下面两章的MATLAB实现即以PWVD为主线，末尾给SPWVD的平滑扩展。

3. 用MATLAB脚本落地：非平稳随机信号的WVD仿真

3.1 构造非平稳随机信号：chirp叠加白噪声

先造一个带频率爬升和随机成分的测试信号，验证时频图里能否同时看到“斜线”和“噪声底噪”。采样率设为1024 Hz，时长2秒，频率从50 Hz线性扫到300 Hz，另加高斯白噪声。

fs = 1024; % 采样率 t = 0:1/fs:2-1/fs; % 时间轴 N = length(t); % 总采样点数 f0 = 50; % 起始频率 f1 = 300; % 结束频率 x = chirp(t, f0, t(end), f1, 'linear'); x = x + 0.3 * randn(size(x)); % 叠加白噪声，噪声标准差0.3 z = hilbert(x); % 解析信号，后续WVD用复信号

chirp生成了线性调频信号，频率随时间单调上升；randn加的是平稳高斯白噪声，两者叠加后就是典型的非平稳信号模型。hilbert返回的复数序列，实部是原信号，虚部是希尔伯特变换结果，z^2的幅度近似原信号包络的平方。用z做WVD时，时频图也不会出现负频率镜像。

如果手里是实测CSV或采集仪数据，这段代码的x换成对应的数据列即可，后续只依赖采样率和信号本身，不关心来源。

3.2 核心函数：用FFT逐时刻计算WVD切片

下面这个函数是整篇文章的核心。输入解析信号x，滞后窗半宽度tau_max，FFT点数N_FFT，输出时频矩阵tfr和归一化频率轴。循环遍历每个时刻n，对自相关核做FFT。

function [tfr, f_axis] = wvd_fft(x, tau_max, N_FFT) % 用FFT逐时刻计算伪WVD分布 % x: 解析信号，行向量 % tau_max: 滞后量m的最大值，决定时间窗宽度 % N_FFT: FFT点数，建议 >= 2*tau_max+1 N = length(x); len_lag = 2 * tau_max + 1; % 滞后轴总长度 win = hamming(len_lag).'; % 滞后域加窗，抑制远处交叉项 tfr = zeros(N, N_FFT); % 时频矩阵初始化 m = -tau_max : tau_max; % 滞后索引 for n = 1:N idx_plus = n + m; % t + tau/2 对应索引 idx_minus = n - m; % t - tau/2 对应索引 valid = (idx_plus >= 1) & (idx_plus <= N) & ... (idx_minus >= 1) & (idx_minus <= N); R = zeros(1, len_lag); R(valid) = x(idx_plus(valid)) .* conj(x(idx_minus(valid))) .* win(valid); spec = fftshift(fft(R, N_FFT)); % 零频移到中心，对应频率轴 tfr(n, :) = spec; end f_axis = (0:N_FFT-1) / N_FFT - 0.5; % 归一化频率，-0.5对应-fs/2 end

逻辑拆开看：idx_plus和idx_minus分别取出n时刻左右各m个样点的索引，valid把越界的索引位置标记为无效，自相关核R在无效处保持0，这就是边界处理。FFT前乘窗win，等效于在滞后域做截断平滑，实际上得到的就是伪WVD。fftshift将零频从索引1搬到中心位置，使f_axis能正确对应负半轴频率。

调用时注意，tfr是复数，绘图用real(tfr)或abs(tfr)^2。WVD理论上应为实数，但数值计算截断会产生很小的虚部，直接取实部即可，abs则会把正负值的差异抹掉。多数论文图显示的是aes(tfr)的平方或实部，具体取决于展示意图。

3.3 参数怎么定：tau_max、N_FFT与窗的取舍

实战调试时先固定tau_max=128，N_FFT=512，观察时频图。如果斜线区域“糊成一片”，说明tau_max偏大或窗太长，适当减小到64；如果交叉条纹明显，说明窗旁瓣抑制不够，改用kaiser窗并调beta参数。N_FFT通常不需要超过2048，多余的零填充只改变插值密度，不提升真实分辨率。

一个容易被忽略的细节：N_FFT取2的幂方便使用FFT，但n的循环本身是串行for结构，对N=2048的信号要跑2048次FFT，总耗时在百毫秒到秒级。若信号超过十几万点，需要分段处理或改用Time-Frequency Toolbox里基于矩阵运算的实现，否则录像时会看到明显的卡顿。

3.4 绘制时频图与排查常见错误

tau_max = 128; N_FFT = 512; [tfr, f_axis] = wvd_fft(z, tau_max, N_FFT); figure; imagesc(t, f_axis(f_axis>=0), abs(tfr(:, f_axis>=0)).^2); axis xy; xlabel('时间/s'); ylabel('频率/Hz'); colorbar;

绘图只保留f_axis>=0的正半轴，因为解析信号在负频段能量趋近于0，显示出来也是噪声，影响视觉判断。imagesc前两个参数是坐标轴，第三个是颜色矩阵，注意方向要配合axis xy，否则图形上下翻转。

常见错误有几种。若整幅图在零频附近出现对称条带，说明用了实信号x而不是解析信号z。若斜线周围出现规律性“排骨纹”，是滞后域窗太长或未加窗导致的交叉项。若图像里有垂直的亮线，则是信号首尾的边界效应，valid置零保护已经有了，只有当tau_max过大导致有效数据占比太低时才会明显，减少tau_max即可缓解。若颜色只有噪点没有斜线，先检查chirp信号的幅度是否被噪声淹没，把噪声系数从0.3降到0.05试跑一次。

4. 仿真操作录像的关键步骤：数据准备与过程记录

4.1 把CSV导入到MATLAB中做FFT仿真

实测信号通常以CSV形式保存，第一列是时间戳，第二列是幅值。录像前先把数据读进来并按规范格式对齐时间轴，避免在录像过程中因为数据格式问题中断操作。

data = readmatrix('sensor_signal.csv'); t_raw = data(:, 1); x_raw = data(:, 2); fs = 1 / mean(diff(t_raw)); % 由时间差估算实际采样率 x_raw = x_raw - mean(x_raw); % 去直流 x_raw = x_raw / max(abs(x_raw)); % 幅值归一化到[-1,1]

readmatrix能自动识别表头和数据区，比csvread更稳健。去直流很重要，WVD对零频附近的直流分量非常敏感，残留的直流会在f=0处形成亮带并掩盖低频信号。幅值归一化不是必须，但它能让噪声标准差和后处理阈值在不同数据间保持一致，录像讲述参数时也更有说服力。

如果CSV的时间戳不是均匀间隔，直接插值到均匀时间轴再计算fs，否则diff(t_raw)的均值会产生偏差，FFT结果同样会失真。这一步在一个可复现的脚本里写完，录像时只执行不修改，能避免多次拍摄的口径不一。

4.2 仿真操作录像里应该录什么内容

仿真操作录像的常见做法是把执行过程分为三段：数据加载与预处理、WVD计算、参数交互调优。录制工具用MATLAB自带的上方Record按钮或第三方录屏，但重点不在工具，而在录什么：

第一段：执行readmatrix，输入CSV路径，展示数据规模和fs输出 第二段：运行wvd_fft，展示tau_max和N_FFT的赋值过程 第三段：逐步调整tau_max从256降到64，观察时频图交叉项变化

录像的解说词应配合参数变化讲，建议用paragrah式的说明，而不是直接报参数值。重点讲清楚“这个参数变大，时频图怎么变”，让观看者能建立参数和图像之间的映射关系。录制像素至少1080p，MATLAB窗口字号调大，命令行窗口和图形窗口分别放左右两侧，避免图形被遮挡。

还有一个实务技巧：先把脚本跑通一遍，确定参数区间和图像输出稳定后再录。如果录制中发现数据异常导致程序报错，不必重拍整段，保留报错过程作为排错讲解反而更有价值，但要在视频里明确标出这是预期展示的错误。

4.3 仿真发散与数值不稳定的处理

热词里提到“仿真发散”，WVD的数值环境里也有类似现象：tfr矩阵出现NaN或Inf，时频图整体变白，或能量随时间指数增长。

第一类原因是信号中包含极端幅值或零点。x(t)幅值出现NaN时，自相关核会扩散，需要在计算前用isfinite检查数据。第二类原因是某段信号幅度突然冲高，比如开关脉冲，使tfr局部能量骤增，显示范围被拉宽后其他区域变得不可见。处理方法是计算后做能量归一化或直接用中位数截断显示上限。

if any(~isfinite(x)) error('输入信号包含NaN或Inf，请先清洗数据'); end figure; surf_abs = abs(tfr).^2; median_val = median(surf_abs(:)); imagesc(t, f_axis, min(surf_abs, 20*median_val));

用中位数的20倍作为显示上限，属于鲁棒可视化策略，能压制尖峰脉冲对色标的影响，同时保留正常的时频结构。这个做法比简单设置caxis上限更稳定，因为不同信号的能量尺度差异很大。真正要修的数据问题，则要靠前置的异常值剔除，WVD是一种二次型变换，对离群点会做平方放大，任何input侧的噪声毛刺在时频图里都会被放大成亮斑。

4.4 工具对比：手写MATLAB、Time-Frequency Toolbox与Python tftb

验证手写函数正确性的方法是与工具箱输出对比。取同一段信号，手写结果的能量峰值位置与工具箱tfrwv一致即可，幅度有细微差异正常，因为窗函数和归一化定义可能不同。工具对比不需要引入额外依赖，只要在开发环境里跑一个简单的difference检查就能确认。

如果后续准备在Simulink里做在线或硬件协同仿真，MATLAB脚本计算WVD只是算法原型，可以把wvd_fft转成MATLAB Function块，把信号源替换成Simulink时间序列模块。这种迁移到嵌入式或仿真的路径，比在Simulink里直接写嵌套for循环要顺得多。

5. 验证算法与抑制交叉项：边缘分布校验和SPWVD平滑

时频图肉眼看起来“像那么回事”还不能证明实现正确。WVD有两条边缘分布性质，可以作为校验手段：对频率积分得到瞬时功率|x(t)|^2，对时间积分得到功率谱|X(f)|^2。用这两条性质检查数值实现，能快速定位代码里频移或窗函数出错的位置。

marginal_t = sum(abs(tfr).^2, 2) / N_FFT; marginal_f = sum(abs(tfr).^2, 1); figure; subplot(2,1,1); plot(t, abs(z).^2 / max(abs(z).^2)); hold on; plot(t, marginal_t / max(marginal_t), 'r--'); legend('解析信号瞬时功率', 'WVD时间边缘');

如果两条曲线形状偏差过大，优先检查时间边缘对应的频率求和范围是否被截断到正半轴。绘图时只显示了正频段，但求和时要全频率，否则能量对不上。频率边缘的验证需要做整段FFT对比，判断WVD在时间方向上的累积能量是否与原始信号频谱一致。这两条边缘性质都通过，就可以判定核心FFT实现基本可靠。

确认基础实现无误后，再看交叉项抑制。非平稳随机信号场景里交叉项往往会掩盖真实分量的边界，SPWVD在PWVD的基础上增加频率方向的平滑窗g(s)，压制速度更快的变化：

function [tfr, f_axis] = spwvd_fft(x, tau_max, N_FFT, win_f) len_lag = 2 * tau_max + 1; win_t = hamming(len_lag).'; Nf = length(win_f); % 频率平滑窗长度，通常为奇数 win_f = win_f / sum(win_f); % 归一化 [tfr_raw, f_axis] = wvd_fft(x, tau_max, N_FFT); tfr = zeros(size(tfr_raw)); for n = 1:size(tfr_raw, 1) tmp = conv(tfr_raw(n, :), win_f, 'same'); % 频率方向卷积平滑 tfr(n, :) = tmp; end end

win_f用元素和为1的归一化窗，保证平滑不改变总能量。平滑宽度Nf取值5~11比较实用，太小起不到抑制效果，太大会把Chirp的斜线本身也抹平。SPWVD在时频图上的视觉特点是：底色噪声更均匀，真实分量边缘更“实”，但频率方向分辨率与PWVD相比变粗。

操作录像在展示交叉项抑制时，应该采用“先PWVD后SPWVD”的对照方式：同一段信号、同一坐标轴范围，只切换平滑窗。录像里说“交叉项怎么看”不如直接展示“这个条纹是交叉项，加平滑后它变淡了但主分量还在”。最后再补充一次边缘分布校验，让观众确认平滑过程没有引入能量损失，完整录像到此结束，所有的脚本、参数和验证逻辑都在前述代码段中留下可复现路径。

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