简介这份资源是面向信号处理方向科研人员与工程师的时频分析MATLAB代码合集聚焦魏格纳分布及其相关时频变换的实现与验证适合已具备一定信号处理基础、希望深入理解非平稳信号时频特性的中高级学习者。压缩包共96个文件全部为.m脚本整体约133KB涵盖信号源生成、时频变换计算、魏格纳分布求解、交叉项处理与可视化绘图等模块可配合短时傅里叶变换、小波变换等方法对比分析。内容预览显示代码涉及多种经典时频分布函数与瞬时频率估计工具便于读者直接调用、修改并嵌入自己的实验流程。目前已有292人学习下载可作为课程设计、论文复现或工程算法原型的参考素材帮助读者在动手实践中理解时频分辨率的权衡与魏格纳分布交叉项问题的处理思路。1. 魏格纳分布到底在算什么从一段非平稳信号说起手里有一段轴承振动信号采样率 12.8 kHz转频 30 Hz 附近有故障冲击可你把它丢进 FFT得到的只是一条糊成一团的谱线——冲击发生的那一瞬间到底对应哪个频率完全看不出来。这就是非平稳信号的经典困境频率成分随时间变化而傅里叶变换把时间信息彻底积分掉了。魏格纳分布Wigner Distribution工程上常写 Wigner-Ville DistributionWVD就是为解决这类问题而生的时频分布工具它把一维时间信号映射到时间-频率二维平面让你能同时看到「什么时候」出现了「什么频率」。这份 part3.zip 里的时频分布代码核心就是魏格纳分布的实现与调用。它适合三类人做旋转机械故障诊断、要把冲击时刻和频率对上号的工程师做雷达/通信信号分析、需要高时频聚集度的算法同学以及正在学时频分析、被 STFT 和 WVD 的差异绕晕的学生。魏格纳分布最大的卖点是时频聚集度远高于短时傅里叶变换代价是会产生交叉项干扰这个矛盾贯穿整篇代码的每一个参数选择。下面从原理、代码结构、参数、避坑一路讲到怎么验证结果对不对。2. 魏格纳分布的数学骨架与代码映射为什么它比 STFT 聚集度高2.1 从定义式到离散实现连续魏格纳分布的定义是信号与其共轭在时间轴上做相关再对时延做傅里叶变换W(t, f) ∫ x(t τ/2) · x*(t − τ/2) · e^(−j2πfτ) dτ关键点在于那个 τ/2它让信号自己和自己做相关而不是像 STFT 那样拿一个窗函数去截。正因为没有窗频率分辨率不受窗长限制时频聚集度自然高。但代价也在这里——两个频率分量之间会产生交叉项出现在两者频率的中点位置这是 WVD 的固有属性不是 bug。离散实现时常见做法是对每个时间点 t取信号在 t 附近的一段构造解析信号去掉负频率避免混叠然后对时延轴做 FFT。代码里通常分三步希尔伯特变换得到解析信号、逐点构造相关序列、对相关序列做 FFT。下面是一段最小可跑的 Python 实现骨架import numpy as np from scipy.signal import hilbert def wigner_ville(x, fs): # x: 实信号, fs: 采样率 # 1. 解析信号抑制负频率带来的交叉项 z hilbert(x) N len(z) # 2. 时延轴长度取 N//2 保证对称 tau np.arange(-N//2, N//2) wvd np.zeros((N, N), dtypecomplex) # 3. 逐时间点构造相关序列并做 FFT for t in range(N): # 边界处理越界处补零避免索引报错 idx1 t tau // 2 idx2 t - tau // 2 valid (idx1 0) (idx1 N) (idx2 0) (idx2 N) corr np.zeros(N, dtypecomplex) corr[valid] z[idx1[valid]] * np.conj(z[idx2[valid]]) wvd[t, :] np.fft.fftshift(np.fft.fft(corr)) return wvd逻辑说明hilbert把实信号变成解析信号这是抑制交叉项的第一步不做的话正负频率会互相干扰。tau取-N//2到N//2保证时延轴关于零点对称。逐点循环里用valid掩码处理边界这是新手最容易翻车的地方——不处理边界直接索引会越界或者补零方式不对导致边缘出现虚假能量。最后fftshift把零频移到中心方便画图。参数说明fs只影响频率轴的刻度换算不参与计算本身N是信号长度直接决定计算量因为循环是 O(N²)N 上万时纯 Python 循环会慢到无法接受实际工程里要么用向量化要么用现成的tftb库。2.2 为什么工程上更常用平滑伪魏格纳分布原始 WVD 的交叉项在真实信号上往往严重到没法看。比如两个频率分量 f1 和 f2交叉项会出现在 (f1f2)/2 处而且幅度可能比真实分量还大。工程上的标准解法是加窗平滑得到平滑伪魏格纳分布SPWVDW_sp(t, f) ∫∫ g(u) · h(τ) · x(t − u τ/2) · x*(t − u − τ/2) · e^(−j2πfτ) du dτ多出来的 g(u) 是时间平滑窗h(τ) 是频率平滑窗。两个窗一加交叉项被压下去代价是时频聚集度下降——这就是时频分析里绕不开的「不确定性原理」权衡。代码里如果看到两个可调窗长参数基本就是 SPWVD 的实现。选型理由很直接如果你的信号是单分量或者分量间隔很远用原始 WVD聚集度最高如果是多分量、交叉项糊成一片老老实实上 SPWVD把窗长调到一个交叉项可接受、聚集度还能看的平衡点。我一般会先用原始 WVD 跑一遍看交叉项有多严重再决定平滑窗开多大。3. 把 part3.zip 跑起来环境、调用与参数怎么设3.1 环境准备与依赖确认拿到一个时频分布代码包第一步不是急着跑而是确认依赖。魏格纳分布的实现通常依赖numpy、scipy画图依赖matplotlib如果代码里用了tftbTime-Frequency Toolbox 的 Python 移植还得单独装。常见做法是先建虚拟环境避免和系统里的科学计算库版本打架python -m venv tfenv source tfenv/bin/activate # Windows 用 tfenv\Scripts\activate pip install numpy scipy matplotlib # 如果代码 import tftb再补一句 pip install tftb逻辑说明虚拟环境隔离是关键时频分析代码经常对numpy版本敏感尤其是np.fft的行为在不同版本间有细微差异。参数说明没有特殊版本要求时用最新稳定版即可但如果代码里写死了某个 API按报错回退版本。提示如果运行时报ModuleNotFoundError: No module named tftb说明代码用了这个库如果报的是cannot import name hilbert检查scipy是否装全。3.2 调用魏格纳分布函数的标准流程假设代码包里有一个wvd.py或者类似的主模块典型调用流程是读信号 → 预处理 → 算 WVD → 画时频图。下面是一段通用调用骨架你可以对照自己包里的函数名替换import numpy as np import matplotlib.pyplot as plt from wvd import wigner_ville # 按实际模块名替换 # 1. 读信号假设是单列文本或 npy sig np.loadtxt(vibration.txt) fs 12800 # 采样率必须和采集时一致 # 2. 去均值避免直流分量在零频处堆一大块能量 sig sig - np.mean(sig) # 3. 算 WVD wvd wigner_ville(sig, fs) # 4. 画图频率轴换算成 Hz freqs np.fft.fftshift(np.fft.fftfreq(len(sig), 1/fs)) plt.imshow(np.abs(wvd).T, aspectauto, originlower, extent[0, len(sig)/fs, freqs[0], freqs[-1]]) plt.xlabel(Time (s)); plt.ylabel(Frequency (Hz)) plt.colorbar(); plt.show()逻辑说明去均值这一步很多人省掉结果时频图零频处一条亮线把真实低频分量全盖住了。np.abs(wvd)取模是因为 WVD 是复值画图只看能量幅度。extent把像素坐标映射到真实时间和频率不设的话横纵轴是采样点索引没法读。参数说明fs必须和采集设备一致填错的话频率轴整体缩放故障特征频率就对不上了。aspectauto让图像自适应长宽比信号长的时候不加这个图会被压扁。3.3 三个必调参数窗长、时延范围、采样率参数作用调大后果调小后果平滑窗长 g/h压制交叉项交叉项弱但聚集度差聚集度高但交叉项明显时延范围 τ决定频率分辨率频率分辨率高计算量大分辨率低可能漏分量采样率 fs频率轴刻度频率范围宽可能混叠范围窄高频看不到时延范围这个参数最容易被忽略。理论上 τ 取满整个信号长度分辨率最高但计算量是 O(N²)实际会截断到一个合理值。截断太狠两个靠近的频率分量就分不开了。我一般先取信号长度的 1/4 到 1/2 试看目标分量能不能分开再定。采样率则是硬约束根据奈奎斯特能看到的最高频率是 fs/2。如果故障特征频率在 6 kHz 而 fs 只有 10 kHz那 5 kHz 以上全是混叠怎么调窗都没用只能重新采。4. 交叉项、边界效应与计算量魏格纳分布避坑清单4.1 现象时频图中间冒出一条不属于任何分量的亮线原因这是 WVD 的交叉项出现在两个真实分量频率的中点。比如信号里有 100 Hz 和 300 Hz交叉项就在 200 Hz。它不是噪声是双线性变换的固有产物。解决上 SPWVD加时间平滑窗 g(u) 和频率平滑窗 h(τ)。窗长从信号长度的 1/16 开始试逐步加大直到交叉项降到可接受。注意平滑窗会同时削弱真实分量别一味加大。4.2 现象时频图左右两端能量异常边缘发黑或发亮原因边界效应。计算 t 附近的 WVD 时t ± τ/2 会超出信号范围补零处理会让边缘的相关值偏小能量失真。解决常见做法是两端各丢弃 τ_max/2 长度的结果只保留中间可信区域。或者对信号做镜像延拓再算算完裁掉延拓部分。代码里如果没做边界处理边缘那几列直接不看。4.3 现象信号一长程序就卡死内存爆掉原因WVD 输出是 N×N 的复矩阵N10000 时就是 1 亿个复数内存直接几个 G。纯 Python 双重循环更是慢到离谱。解决分段计算把长信号切成若干段分别算 WVD 再拼接或者降采样只要目标频率在降采样后的奈奎斯特范围内就行。向量化实现能把速度提几十倍但内存问题还在分段是更稳的路子。4.4 现象换了台机器跑结果频率轴对不上原因fftfreq或fftshift的用法在不同代码里不一致有的用采样点数算有的用信号长度算差一个点频率刻度就偏。解决固定一套换算逻辑频率轴统一用np.fft.fftshift(np.fft.fftfreq(N, 1/fs))N 用实际参与 FFT 的长度。换机器前先跑一个已知频率的正弦信号验证比如 1 kHz 正弦看时频图亮线是不是正好在 1 kHz。4.5 现象解析信号做完时频图反而更乱了原因hilbert变换对非窄带信号或者含直流分量的信号效果不好解析信号构造失真负频率没压干净反而引入新干扰。解决先去掉直流分量再做hilbert如果信号带宽很宽考虑先带通滤波到目标频段再算 WVD。不是所有信号都适合直接上解析信号这一步要验证。5. 怎么验证你的魏格纳分布算对了三个可复现的检验技巧第一个技巧是单频正弦检验。造一个 1 kHz 的正弦fs 取 10 kHz算 WVD正确的时频图应该是一条水平亮线位置精确在 1 kHz时间方向均匀。如果亮线有倾斜或者位置偏移说明频率轴换算错了如果亮线很粗说明时延范围截得太短频率分辨率不够。这个检验五分钟就能做完但能排掉一大半低级错误。第二个技巧是双分量交叉项定位。造 1 kHz 加 3 kHz 两个正弦正确的 WVD 除了两条真实亮线还会在 2 kHz 处出现一条交叉项。这条交叉项的位置是理论可预测的正好是两个频率的中点。如果你看到的交叉项不在 2 kHz那要么是解析信号没做对要么是频率轴映射有问题。反过来如果你上了 SPWVD 而 2 kHz 那条线还在说明平滑窗没起作用检查窗函数是不是真的加进去了。第三个技巧是用已知的线性调频信号chirp验证时频聚集度。造一个从 1 kHz 线性扫到 5 kHz 的 chirpWVD 应该呈现一条清晰的斜线斜率对应调频速率。STFT 在同参数下这条线会明显更粗把两者并排画出来聚集度的差距一目了然。这也是判断代码有没有真正实现 WVD 而不是拿 STFT 糊弄的最直接方法。我自己的习惯是每换一个信号源或者改一次采样率先把这三个检验跑一遍确认代码没退化再上真实数据。时频分析这东西图一画出来看着都差不多但频率轴偏一点、交叉项位置错一点后面故障诊断的结论就全歪了。宁可前面花十分钟验证也别拿着错图去下结论。希望帮到你。本文还有配套的精品资源点击获取