简介利用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时刻的频率切片。关键区别在这里短时傅里叶变换是先把信号乘窗函数加窗再做FFTWVD是对自相关核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[nm]·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]的长度是2M1FFT点数N_FFT必须不小于2M1否则会因为序列截断丢信息。更常见的做法是N_FFT取2的幂且大于2M1比如M128时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分布公式要点时频聚集性交叉项抑制适用场景WVD不加任何窗直接全滞后域FFT最高单分量接近理想无单分量信号、短观察窗内的特征分析伪WVDPWVD滞后域加窗h(τ)截断m略有下降时间方向分辨率受影响弱能压远距离交叉项工程常用默认值Chirp类信号平滑伪WVDSPWVD滞后窗h(τ) 频率平滑窗g(s)聚集性继续下降较强但时频弥散多分量非平稳随机信号重排SPWVD在SPWVD基础上做能量重排恢复部分聚集性较强且能量集中已知噪声较大、需要可视化的场景仅基于标题“利用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_maxFFT点数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_max1 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_max64~256决定滞后域窗长也就是时间方向的积分宽度频域更细但交叉项变多时间分辨降低N_FFT2的幂≥2*tau_max1控制频率采样点数频率轴更密运算量增加窗类型hamming/hanning滞后域平滑压制远距离交叉项旁瓣更低主瓣变宽是否取解析信号必须消除正负频率镜像不取时零频附近出现假条纹实战调试时先固定tau_max128N_FFT512观察时频图。如果斜线区域“糊成一片”说明tau_max偏大或窗太长适当减小到64如果交叉条纹明显说明窗旁瓣抑制不够改用kaiser窗并调beta参数。N_FFT通常不需要超过2048多余的零填充只改变插值密度不提升真实分辨率。一个容易被忽略的细节N_FFT取2的幂方便使用FFT但n的循环本身是串行for结构对N2048的信号要跑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_axis0), abs(tfr(:, f_axis0)).^2); axis xy; xlabel(时间/s); ylabel(频率/Hz); colorbar;绘图只保留f_axis0的正半轴因为解析信号在负频段能量趋近于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对零频附近的直流分量非常敏感残留的直流会在f0处形成亮带并掩盖低频信号。幅值归一化不是必须但它能让噪声标准差和后处理阈值在不同数据间保持一致录像讲述参数时也更有说服力。如果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式的说明而不是直接报参数值。重点讲清楚“这个参数变大时频图怎么变”让观看者能建立参数和图像之间的映射关系。录制像素至少1080pMATLAB窗口字号调大命令行窗口和图形窗口分别放左右两侧避免图形被遮挡。还有一个实务技巧先把脚本跑通一遍确定参数区间和图像输出稳定后再录。如果录制中发现数据异常导致程序报错不必重拍整段保留报错过程作为排错讲解反而更有价值但要在视频里明确标出这是预期展示的错误。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实现方式优点缺点适合场景MATLAB手写wvd_fft逻辑完全可控参数透明无额外依赖循环慢代码量大学习原理、定制算法MATLAB Time-Frequency Toolboxtfrwv等函数开箱即用实现优化过需额外安装工具包参数封装多快速试算、对比验证Python tftb开源免费基于NumPy便于集成到自动处理流程文档偏少版本间API有变动需要批量处理或工业部署验证手写函数正确性的方法是与工具箱输出对比。取同一段信号手写结果的能量峰值位置与工具箱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 endwin_f用元素和为1的归一化窗保证平滑不改变总能量。平滑宽度Nf取值5~11比较实用太小起不到抑制效果太大会把Chirp的斜线本身也抹平。SPWVD在时频图上的视觉特点是底色噪声更均匀真实分量边缘更“实”但频率方向分辨率与PWVD相比变粗。操作录像在展示交叉项抑制时应该采用“先PWVD后SPWVD”的对照方式同一段信号、同一坐标轴范围只切换平滑窗。录像里说“交叉项怎么看”不如直接展示“这个条纹是交叉项加平滑后它变淡了但主分量还在”。最后再补充一次边缘分布校验让观众确认平滑过程没有引入能量损失完整录像到此结束所有的脚本、参数和验证逻辑都在前述代码段中留下可复现路径。本文还有配套的精品资源点击获取