简介:这份资源是面向信号处理方向科研人员与工程师的时频分析MATLAB代码合集,聚焦魏格纳分布及其相关时频变换的实现与验证,适合已具备一定信号处理基础、希望深入理解非平稳信号时频特性的中高级学习者。压缩包共96个文件,全部为.m脚本,整体约133KB,涵盖信号源生成、时频变换计算、魏格纳分布求解、交叉项处理与可视化绘图等模块,可配合短时傅里叶变换、小波变换等方法对比分析。内容预览显示代码涉及多种经典时频分布函数与瞬时频率估计工具,便于读者直接调用、修改并嵌入自己的实验流程。目前已有292人学习下载,可作为课程设计、论文复现或工程算法原型的参考素材,帮助读者在动手实践中理解时频分辨率的权衡与魏格纳分布交叉项问题的处理思路。
1. 魏格纳分布到底在算什么:从一段非平稳信号说起
手里有一段轴承振动信号,采样率 12.8 kHz,转频 30 Hz 附近有故障冲击,可你把它丢进 FFT,得到的只是一条糊成一团的谱线——冲击发生的那一瞬间到底对应哪个频率,完全看不出来。这就是非平稳信号的经典困境:频率成分随时间变化,而傅里叶变换把时间信息彻底积分掉了。魏格纳分布(Wigner Distribution,工程上常写 Wigner-Ville Distribution,WVD)就是为解决这类问题而生的时频分布工具,它把一维时间信号映射到时间-频率二维平面,让你能同时看到「什么时候」出现了「什么频率」。
这份 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), dtype=complex) # 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, dtype=complex) 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,交叉项会出现在 (f1+f2)/2 处,而且幅度可能比真实分量还大。工程上的标准解法是加窗平滑,得到平滑伪魏格纳分布(SPWVD):
W_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,如果代码里用了tftb(Time-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, aspect='auto', origin='lower', 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必须和采集设备一致,填错的话频率轴整体缩放,故障特征频率就对不上了。aspect='auto'让图像自适应长宽比,信号长的时候不加这个图会被压扁。
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 的复矩阵,N=10000 时就是 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 的 chirp,WVD 应该呈现一条清晰的斜线,斜率对应调频速率。STFT 在同参数下这条线会明显更粗,把两者并排画出来,聚集度的差距一目了然。这也是判断代码有没有真正实现 WVD 而不是拿 STFT 糊弄的最直接方法。
我自己的习惯是:每换一个信号源或者改一次采样率,先把这三个检验跑一遍,确认代码没退化再上真实数据。时频分析这东西,图一画出来看着都差不多,但频率轴偏一点、交叉项位置错一点,后面故障诊断的结论就全歪了。宁可前面花十分钟验证,也别拿着错图去下结论。希望帮到你。
本文还有配套的精品资源,点击获取