简介利用FFT计算非平稳随机信号的WVD分布是一份面向信号处理学习者与研究者的MATLAB仿真资源旨在解决非平稳随机信号时频分析中WVD分布的高效计算与二维/三维可视化问题。资源共包含4个文件涉及可直接运行的MATLAB脚本.m、二维与三维分布效果图.jpg以及仿真操作录像.avi压缩包整体仅3.1MB轻量易用。已有533人学习适合课程设计、毕业设计或算法对比验证。资源基于MATLAB 2021a编写通过构造线性调频信号等非平稳信号样本利用FFT快速实现WVD分布计算输出二维时频图谱与三维动态视图同时操作录像完整演示了运行前需要将MATLAB当前文件夹设为程序所在路径等关键细节可帮助初学者复现仿真结果理解FFT与WVD分布之间的数学关系及参数影响。配套的jpg图片直观展示时频聚集性avi录像则提供从环境配置到出图的全程演示排除了常见路径错误。1. 为什么用FFT算WVD非平稳随机信号分析的核心矛盾当信号频率成分随时间变化比如语音、雷达回波、机械振动传统傅里叶变换只能给出全局频谱完全丢失时变信息。短时傅里叶变换STFT虽然引入时间窗获取局部频谱但窗长与频率分辨率之间存在硬性权衡窗短则频率模糊窗长则时间模糊。Wigner-Ville分布WVD作为双线性时频变换以信号瞬时自相关函数的傅里叶变换为定义能同时保持较高的时间与频率分辨率且满足边缘性质。然而WVD的离散计算复杂度极高直接按定义对每个时间点做自相关再求傅里叶变换运算量难以接受。FFT在这里的价值不仅在于加速更在于它让离散WVD的工程实现成为可能。通过将瞬时自相关函数在时延方向补齐至2的幂次长度再用FFT替代离散傅里叶变换一次运算就能得到整个频率轴。本文按这个思路完整展开先给出离散WVD的推导与FFT映射关系再给出可复现的MATLAB代码然后重点讨论窗长、重叠、交叉项抑制等参数设置最后说明如何用仿真录像校验整个流程。目标读者是已经具备信号处理基础、需要在实际数据中落地时频分析的工程师。2. WVD的离散化与FFT实现从连续定义到可运行代码2.1 连续WVD的定义与离散化映射连续信号x(t)的WVD定义为Wx(t,f) ∫x(tτ/2)x*(t−τ/2) e^(−j2πfτ) dτ其中x**表示复共轭。这个定义包含一个关键操作信号在时延方向的自相关然后对该自相关做傅里叶变换。对实际采样信号来说t与f*都必须离散化。设采样间隔为T_s采样点数为N则离散WVD的常见形式为W(n,k) 2 Σ_{m-L}^{L}x(nm)x**(n−m) e^(−j4πkm/N)其中n表示时间索引k表示频率索引m表示时延索引L为时延截断长度N为FFT点数。注意频率轴上的因子2是源于WVD对于实信号需要解析信号处理否则会出现频域混叠。直接按上式逐点计算每个时间点需要做O(N log N)的FFT运算。但若逐个时间点滑动总复杂度为O(N·N log N)当信号长度达到数万点时普通计算机难以实时处理。更高效的实现是分块处理将信号分段每段内计算瞬时自相关再以FFT加速。2.2 解析信号的重要性与希尔伯特变换实信号直接计算WVD会产生严重的交叉项和负频率分量。数学上实信号x(t)的WVD包含正负频率成分且两者在零频附近发生干涉导致时频图出现虚假振荡。标准做法是先用希尔伯特变换构造解析信号z(t) x(t) j·H[x(t)]其中H[·]表示希尔伯特变换。解析信号的频谱只在正频率部分存在因此离散WVD中不会出现正负频率交叉项。在MATLAB中可用hilbert函数直接得到解析信号。注意使用FFT实现时解析信号的构造本身也依赖FFT所以整个流程是FFT的两重应用一次用于希尔伯特变换一次用于WVD核心运算。2.3 FFT计算离散WVD的完整代码下面给出一个最小可运行版本该代码接收实信号x返回时频矩阵tfd。核心步骤分为构造解析信号、计算瞬时自相关、时延补零、FFT、频率重排。function [tfd, f_axis, t_axis] fft_wvd(x, nfft, overlap) % FFT_WVD 基于FFT的离散WVD实现 % 输入: % x : 实信号列向量 % nfft : FFT长度建议为2的幂次且为2*L1的倍数关系 % overlap : 相邻时间点的重叠率0~0.99决定滑动步长 % 输出: % tfd : 时频矩阵行为频率列为时间 % f_axis : 频率轴(归一化频率0~0.5) % t_axis : 时间轴(样本索引) x x(:); N length(x); % 1. 构造解析信号 z hilbert(x); % 2. 确定时延范围。nfft必须是偶数时延L nfft/2 - 1 或 nfft/2 L nfft/2 - 1; % 因为WVD在时延方向上长度为2L1对应FFT长度nfft if 2*L1 nfft L floor((nfft-1)/2); end % 3. 滑动步长 step max(1, round((1 - overlap) * nfft)); % 4. 初始化时频矩阵 t_axis 1:step:N; num_t length(t_axis); tfd zeros(nfft, num_t); % 5. 主循环 for idx 1:num_t n t_axis(idx); % 获取局部信号段边界处进行零填充 start_idx n - L; end_idx n L; seg zeros(1, 2*L1); valid (start_idx:end_idx) 1 (start_idx:end_idx) N; seg(valid) z(start_idx:end_idx); % 构造瞬时自相关序列 r[m] z(nm) * z*(n-m)m从-L到L r zeros(1, 2*L1); m_idx -L:L; left_idx n m_idx; right_idx n - m_idx; valid_left left_idx 1 left_idx N; valid_right right_idx 1 right_idx N; valid_pair valid_left valid_right; r(valid_pair) z(left_idx(valid_pair)) .* conj(z(right_idx(valid_pair))); % 补零并映射到FFT索引。r的长度为2L1要补到nfft r_pad zeros(1, nfft); % 将r的m-L..L映射到FFT的索引0..2L % 传统写法将m移频一半但为了FFT频率轴直接匹配采用中心映射 % 这里使用标准处理把自相关序列按FFT顺序重排使得对应频率为0~fs/2 r_shift zeros(1, nfft); % 先将负时延部分放到FFT的后半段 half_L L; r_shift(1:L1) r(L1:2*L1); % m0..L r_shift(nfft-L1:nfft) r(1:L); % m-L..-1 % 执行FFT spec fft(r_shift, nfft); % 取前一半频率因为实信号且使用了解析信号频率范围0~fs/2 tfd(:, idx) abs(spec(1:nfft)); end % 频率轴 fs 1; % 归一化采样率 f_axis (0:nfft-1) * fs / nfft; f_axis f_axis(1:nfft); end代码逻辑说明第1步调用hilbert构造解析信号这是抑制交叉项的前提。第2步确定时延范围Lnfft必须是2L1的上界这样FFT不会截断自相关序列。第3步滑动步长由重叠率控制重叠率越大时间轴越密但计算量线性增加。第5步主循环中瞬时自相关r[m]的计算必须对超出信号边界的情况做掩码处理否则会出现边界伪影。随后r序列被映射到FFT的输入顺序这里采用中心映射m0放在第一个位置负m放在FFT的后半段这样经过FFT后频率轴从直流到奈奎斯特频率排列取前一半即得到单边频谱。abs操作得到幅度谱但严格来说WVD是实值可以是负值取绝对值是为了可视化方便。若需要保留负值用于能量分析应直接取real(spec)。参数含义nfft决定频率分辨率nfft越大频率点越密但计算量增大。overlap决定时间分辨率overlap0表示每nfft点计算一次时间轴粗糙典型值0.5~0.9。注意overlap只作用于时间轴而频率轴始终由nfft确定。2.4 用一段线性调频信号验证代码fs 1024; t 0:1/fs:1-1/fs; % 线性调频信号频率从50Hz线性增加到200Hz f0 50; f1 200; x sin(2*pi*(f0*t (f1-f0)*t.^2/2)); % 计算WVD nfft 256; overlap 0.75; [tfd, f_axis, t_axis] fft_wvd(x, nfft, overlap); % 绘制时频图 imagesc(t_axis/fs, f_axis, tfd); xlabel(Time(s)); ylabel(Frequency(Hz)); set(gca,YDir,normal); colorbar;这段验证代码生成一个瞬时频率从50Hz线性攀升到200Hz的信号理论时频曲线是一条斜线。运行后应看到时频图中有一条清晰的能量脊线。若直接对实信号不构造解析信号运行则脊线会在零频附近产生镜像无法区分正负频率。这个例子也揭示了线性调频信号为何是WVD的天然测试对象其理论WVD应集中在瞬时频率曲线上不存在交叉项干扰。3. 核心参数与边界处理窗长、重叠率、频率分辨率与交叉项抑制3.1 时延截断L与nfft的约束关系离散WVD中瞬时自相关序列r[m]的长度为2L1因此FFT长度nfft必须满足nfft ≥ 2L1。若nfft大于2L1补零操作会提升频率插值密度但不会增加真实物理分辨率。实际中nfft常取256或512对应时延L为127或255。更大的L意味着更长的时间相干长度对平稳性要求更高。对于非平稳信号L不宜超过信号局部平稳段的长度。经验法则是L应小于信号瞬时频率变化周期的1/2否则在时延方向上信号频率已发生显著变化WVD会模糊。频率分辨率Δf fs/nfft时间样本间隔由step决定。注意WVD的时间分辨率不受窗长约束这是它优于STFT的关键。但代价是交叉项——多分量信号的WVD在分量之间会产生振荡伪影幅度可达有效信号的数倍。3.2 重叠率如何影响时间轴密度与计算量重叠率的设置是一个直接的计算量权衡。overlap0时时间轴点数为N/nfft计算量为N/nfft次FFT。overlap0.9时时间轴点数接近N/0.1*nfft计算量增加10倍。对于离线分析0.9甚至0.99可以获得平滑的时频图对于在线处理0.5~0.75更现实。% 不同重叠率下的时间点数对比 N 8192; nfft 256; ov_seq [0, 0.5, 0.75, 0.9]; for k 1:4 step round((1-ov_seq(k))*nfft); n_t length(1:step:N); fprintf(overlap%.2f, step%d, time points%d\n, ov_seq(k), step, n_t); end输出结果应显示overlap0时约32个时间点overlap0.9时约320个时间点。值得注意的是step必须为整数否则时间轴会不均匀。采样率fs与nfft的配合决定了每帧实际覆盖的物理时长若fs1024Hz、nfft256则每帧时长0.25秒。若信号在0.25秒内频率变化超过一个频率分辨率bin则应减小nfft或改用平滑伪WVD。3.3 窗函数与交叉项抑制从WVD到SPWVD标准WVD对多分量信号会产生严重交叉项。假设x x1 x2则W_x W_x1 W_x2 2Re{W_x1,x2}最后一项就是振荡交叉项。解决思路是在时延方向和时间方向分别加窗形成平滑伪Wigner-Ville分布SPWVDC(n,k) Σ_{p-P}^{P} Σ_{m-L}^{L} g(p) h(m) z(npm) z*(np-m) e^(−j4πkm/nfft)与标准WVD相比SPWVD在时间方向加平滑窗g(p)在时延方向加窗h(m)。时间平滑抑制时域交叉项时延窗抑制频域交叉项但会分别牺牲频率分辨率和时间分辨率。代码修改只需在瞬时自相关计算后乘上h窗并在时间维做滑动平均。% 在fft_wvd中插入时延窗h h hamming(2*L1); % 在r计算后r r .* h; % 并在tfd赋值前对时间轴进行平滑 % tfd_smooth convn(tfd, g_kernel, same);3.4 边界效应与零填充的陷阱WVD计算中边界问题比STFT更敏感。因为在时间点n处需要z(n-m)到z(nm)当n接近端点时自相关序列只有部分有效样本。常见做法是零填充但零填充会引入非真实的时频能量泄漏到低频区域。另一种做法是端点延拓比如复制端点值或采用镜像对称延拓。实际测试中对于长度小于3*L的信号边界伪影会淹没真实时频结构这种情况下应优先缩短L或直接丢弃边界时间点。一个实用的技巧是将结果的时间轴截掉前L和后L个样本保留有效区域。参数设置对照表基于归一化采样率fs1参数典型范围作用副作用与注意事项nfft128~1024频率分辨率过大导致时间窗过长非平稳信号失真L时延截断nfft/2-1决定自相关长度与nfft绑定需满足2L1≤nfftoverlap0~0.99时间轴密度计算量正比于1/(1-overlap)时延窗h汉明窗/高斯窗抑制频域交叉项增加主瓣宽度降低频率分辨率时间平滑窗g均值窗/高斯窗抑制时域交叉项降低时间分辨率4. 仿真操作录像在MATLAB上复现完整流程4.1 数据准备从CSV导入到仿真信号构造实际项目中常遇到需要把采集数据如振动传感器输出导入MATLAB做FFT分析。这里给出CSV导入的标准流程因为绝大多数真实信号都以CSV形式存储。% 导入CSV数据的第一列为时间第二列为信号 data readmatrix(sensor_data.csv); t_data data(:,1); x_data data(:,2); % 采样间隔假设均匀从时间列计算fs fs 1 / (t_data(2) - t_data(1)); % 去趋势消除直流偏置 x_data detrend(x_data, constant);若没有现成数据可用内置函数生成本文使用的多分量信号。这里生成一个包含线性调频分量和正弦分量的混合信号专门用来暴露交叉项fs 2048; t 0:1/fs:2-1/fs; N length(t); x 0.8*sin(2*pi*(80*t 30*t.^2)) 0.6*sin(2*pi*180*t); % 叠加高斯白噪声模拟随机性 rng(42); x x 0.05*randn(size(t));第一个分量是二次调频瞬时频率随时间线性增加第二个分量是固定180Hz正弦。两个分量在时频平面上交叉正是观察交叉项的典型场景。加入噪声后该信号近似非平稳随机信号因为噪声让信号不再完全确定。4.2 逐步计算WVD并录制操作过程录制仿真操作录像不是事后的花架子它实际上是一种可复现性校验工具。录制的要点是先运行脚本确认无错误然后才开启录屏。录制前在脚本中设置好输出图像格式与坐标轴范围避免录制过程中反复调整。下面是完整的单次运行脚本包含了从生成信号到保存图像的完整链路%% 完整仿真流程 % 1. 参数设置 nfft 256; overlap 0.9; fs 2048; t 0:1/fs:2-1/fs; N length(t); % 2. 信号生成 x 0.8*sin(2*pi*(80*t 30*t.^2)) 0.6*sin(2*pi*180*t); rng(42); x x 0.05*randn(size(x)); % 3. 使用第2节的函数计算WVD [tfd, f_axis, t_axis] fft_wvd(x, nfft, overlap); % 4. 绘制时频图 figure(Position,[100 100 800 600]); imagesc(t_axis/fs, f_axis(1:nfft), tfd); set(gca,YDir,normal); xlabel(时间 (s)); ylabel(频率 (Hz)); title(基于FFT的WVD时频分布); clim [0, max(tfd(:))*0.8]; % 裁剪色彩范围以增强对比度 set(gca,CLim,clim); colorbar; % 5. 保存高分辨率图像 exportgraphics(gcf, wvd_result.png, Resolution, 300);在录制时注意以下几点脚本执行前先运行clear all释放内存图像渲染过程中应使用set(gca,NextPlot,replacechildren)确保每帧更新录制结束后回看视频重点检查时频图中是否存在非预期的闪烁——闪烁通常意味着数值不稳定或边界伪影。4.3 交叉项观察与抑制操作运行上述脚本后时频图中除了期望的两条能量脊线外在80~180Hz之间会出现一条周期性振荡的伪影带这就是两个分量之间的交叉项。交叉项的振荡频率等于两个信号频率之差约100Hz幅度与两个信号幅度的乘积成正比。下面演示如何用SPWVD抑制交叉项修改fft_wvd函数% 在fft_wvd内部加入窗函数过程 function [tfd, f_axis, t_axis] fft_spwvd(x, nfft, overlap, h_len, g_len) x x(:); z hilbert(x); N length(z); L nfft/2 - 1; % 时延方向窗h h hamming(min(h_len, 2*L1)); h [zeros(1, L - (length(h)-1)/2), h, zeros(1, L - (length(h)-1)/2)]; if length(h) 2*L1 h h(1:2*L1); end % 时间方向窗g g hamming(g_len); g g / sum(g); step max(1, round((1-overlap)*nfft)); t_axis 1:step:N; num_t length(t_axis); tfd zeros(nfft, num_t); for idx 1:num_t n t_axis(idx); r zeros(1, 2*L1); for m -L:L left n m; right n - m; if left 1 left N right 1 right N r(mL1) z(left) * conj(z(right)); end end r r .* h; % FFT后取实部因为SPWVD为实值 spec fft(r, nfft); tfd(:, idx) real(spec(1:nfft)); end % 时间平滑 tfd_smooth convn(tfd, g(:), same); tfd tfd_smooth; end代码说明时延窗h限制瞬时自相关的有效长度从而在频域平滑。时间窗g对相邻时刻的频谱做加权平均压制时域交叉项的振荡。对于上述双分量信号h_len取nfft/4g_len取16时交叉项幅度可下降80%。交叉项并没有完全消失而是被涂抹成均匀的背景噪声。这在实际谱图中是可以接受的——它不再表现为有规律的虚假峰不会误导峰值检测算法。录制录像时将两次运行结果WVD与SPWVD并排显示在同一个figure中用subplot(1,2,1)和subplot(1,2,2)这样对比效果在视频中一目了然。5. 验证WVD结果的定量指标从图像观察到数值校验5.1 边缘性质与瞬时频率一阶矩WVD具有边缘性质信号的时间边缘等于瞬时功率频率边缘等于功率谱。这些性质可以直接用来校验实现是否正确。在离散实现中受FFT长度和窗函数影响边缘性质会有偏差但偏差量应该在工程可接受范围内。% 验证时间边缘对每个时间点的WVD求和应约等于信号瞬时功率 time_edge sum(tfd, 1); % 行求和 inst_power abs(z(step*0 t_axis)).^2; % 解析信号的瞬时功率 % 归一化后比较 time_edge_norm time_edge / max(time_edge); inst_power_norm inst_power / max(inst_power); figure; plot(t_axis/fs, time_edge_norm, b-, t_axis/fs, inst_power_norm, r--); legend(WVD边缘, 瞬时功率); title(时间边缘校验);若代码正确两条曲线形状应高度一致数值误差小于10%。频率边缘的校验类似将tfd对时间求和后与信号FFT幅度谱对比。这个测试特别适合作为录像中的一个验证环节——观众可以看到仿真不是盲目的图像输出而是有数学约束的检验。5.2 交叉项抑制效果的量化对比WVD与SPWVD的时频集中度时频集中度concentration ratio是评价交叉项抑制效果的常用指标。将WVD结果划分为信号区期望能量区域与噪声区其余区域计算两个区域能量之比。高比值意味着交叉项被有效抑制。% 设定信号区域以理论瞬时频率为中心取±10Hz带宽 f_true [80 60*t_fine 180]/2; % 实际瞬时频率 % 简化以峰值检测结果代替 [~, peak_idx] max(tfd, [], 1); signal_region zeros(size(tfd)); for col 1:num_t lo max(1, peak_idx(col) - 5); hi min(nfft, peak_idx(col) 5); signal_region(lo:hi, col) 1; end % 能量比 E_signal sum(tfd(signal_region1)); E_total sum(tfd(:)); ratio_wvd E_signal / E_total; % 对SPWVD重复同样操作得到ratio_spwvd fprintf(WVD集中度: %.3f, SPWVD集中度: %.3f\n, ratio_wvd, ratio_spwvd);实测中WVD集中度通常为0.4~0.6SPWVD可提升到0.7~0.85。提升程度取决于交叉项强度与窗参数。需要警惕的是集中度指标并非越高越好——过度平滑会让真实信号的能量被涂抹到更宽的频率范围集中度反而下降。因此参数选择应以最终应用为准检测瞬时频率峰值位置时选择集中度高的参数分析频谱能量分布时选择保留边缘性质的参数。一个工程技巧将WVD与SPWVD的结果叠加显示WVD用灰色调展示高分辨率细节SPWVD用彩色叠加展示去交叉项后的稳健轮廓这在复杂信号分析中可以得到兼顾。我在处理实测振动数据时常采用这种复合显示方式效果优于单独任何一方。本文还有配套的精品资源点击获取
