简介这份资源是面向信号处理方向科研人员与工程师的时频分析MATLAB代码合集聚焦魏格纳分布及其相关时频变换的实现与验证适合已具备一定信号处理基础、希望深入理解非平稳信号时频特性的学习者。压缩包共96个文件全部为.m脚本整体约133KB涵盖信号源生成、时频变换计算、魏格纳分布求解、交叉项处理以及时频图可视化等模块并配有示例脚本便于对照运行。资源围绕短时傅里叶变换、小波变换与魏格纳分布展开代码中涉及瞬时频率估计、模糊函数、平滑伪魏格纳分布等具体算法可帮助读者从理论公式过渡到可执行的数值实验。目前已有292人学习下载适合作为课程设计、论文复现或工程预研的参考素材通过阅读与调试这些脚本能够加深对时频分辨率、交叉项抑制等关键问题的理解并积累可复用的分析工具。1. 从 part3.zip 说起魏格纳分布到底能解决什么信号问题如果你手头有一个part3.zip里面塞着几段非平稳信号和一份时频分布代码大概率你正卡在同一个问题上傅里叶变换告诉你信号里有哪些频率却不告诉你这些频率在什么时刻出现。风机轴承的冲击、雷达回波的调频、脑电的瞬态爆发全是频率随时间变的信号普通频谱图一平均就把时间信息抹平了。魏格纳分布Wigner Distribution工程上常写 WVDWigner-Ville Distribution就是冲着这个痛点来的它属于 Cohen 类二次时频分布在时频平面上同时给出时间和频率的联合能量分布而且对线性调频信号的理论时频聚集性接近最优这是短时傅里叶变换STFT很难做到的。代价也很明确——它是双线性的多分量信号之间会产生交叉项cross-term这也是所有 WVD 代码里必须处理的核心矛盾。这篇笔记就围绕part3.zip里那套时频分布代码把魏格纳分布从原理、实现、参数到踩坑讲透适合已经会写 Python 或 MATLAB、想真正把时频分布跑出可用图的人。2. 魏格纳分布的原理与选型为什么不用 STFT 而选 WVD2.1 从瞬时自相关到 WVD 的定义式魏格纳分布的出发点不是加窗分段而是瞬时自相关函数。对解析信号 $x(t)$其 WVD 定义为$$W_x(t,f)\int_{-\infty}^{\infty} x\left(t\frac{\tau}{2}\right)x^*\left(t-\frac{\tau}{2}\right)e^{-j2\pi f\tau},d\tau$$把它和 STFT 对比一下就明白差异在哪。STFT 是 $x(\tau)$ 乘一个窗 $h(\tau-t)$ 再做傅里叶变换窗长一旦定死时间分辨率和频率分辨率就被海森堡不确定性原理锁死你只能二选一。WVD 没有窗它把信号在时间轴上对称地掰开——一半往前 $\tau/2$一半往后 $\tau/2$做自相关再变换。这个对称结构让它在理论上对线性调频LFM信号能形成一条理想的直线脊聚集性远好于任何固定窗的 STFT。代价来自二次两个字。因为 $W_x$ 是信号的双线性泛函当信号是多个分量之和 $xx_1x_2$ 时$$W_x W_{x_1}W_{x_2}2,\text{Re}{W_{x_1,x_2}}$$前两项是各自的自项我们想要的第三项就是交叉项它会出现在两个真实分量时频位置的几何中点幅度可能比自项还大振荡剧烈。这就是为什么直接对多分量信号跑 WVD图上经常糊成一片。2.2 离散化别直接对实信号做 WVD落到代码第一个坑是解析信号。实信号的 WVD 在正负频率上对称会产生无意义的交叉项所以标准做法是先做希尔伯特变换取解析信号。第二个坑是离散 WVD 的采样率问题定义式里 $x(t\tau/2)$ 涉及半整数索引直接离散会引入混叠。工程上通用做法是先把信号在时域插值 2 倍上采样再计算离散 WVD最后频率轴只取一半。下面是一段可直接抄的 Python 最小实现用scipy.signal.hilbert取解析信号用numpy做离散 WVDimport numpy as np from scipy.signal import hilbert, resample def wvd(x, fs, n_freqNone): 计算离散魏格纳分布 x : 输入实信号 (1D array) fs : 采样率 (Hz) n_freq : 频率轴点数, 默认取信号长度 返回 : (t, f, W) t为时间轴, f为频率轴, W为时频矩阵 (freq x time) x np.asarray(x, dtypefloat) N len(x) # 1. 取解析信号, 消除负频交叉项 xa hilbert(x) # 2. 时域2倍上采样, 避免半整数索引混叠 xa2 resample(xa, 2 * N) N2 len(xa2) if n_freq is None: n_freq N # 3. 逐时间点构造瞬时自相关, 再做FFT W np.zeros((n_freq, N), dtypecomplex) lags np.arange(-N, N) # 滞后 tau 范围 for ti in range(N): # 对称取 x(ttau/2) 与 x*(t-tau/2) idx_p ti lags # 上采样后索引步进为1 idx_m ti - lags valid (idx_p 0) (idx_p N2) (idx_m 0) (idx_m N2) r np.zeros(len(lags), dtypecomplex) r[valid] xa2[idx_p[valid]] * np.conj(xa2[idx_m[valid]]) # 4. 对滞后维做FFT得到频率维 spec np.fft.fftshift(np.fft.fft(r, n2 * n_freq)) W[:, ti] spec[n_freq // 2: n_freq // 2 n_freq] f np.linspace(0, fs / 2, n_freq) t np.arange(N) / fs return t, f, W逻辑说明第 1 步hilbert把实信号变成解析信号这是所有 WVD 实现的前提不做的话正负频率交叉项会污染整个图。第 2 步resample做 2 倍上采样是为了让t±τ/2落在整数索引上否则你要处理半整数插值代码复杂且容易错。第 3 步对每个时间点构造滞后序列r注意idx_p和idx_m是镜像关系这正是 WVD 对称结构的体现。第 4 步对滞后维做 FFT 得到频率维fftshift后只取正频率一半。参数说明fs必须和采集时一致否则频率轴全错n_freq决定频率分辨率取信号长度是常见默认但信号很长时比如 10 万点内存会爆需要分段或降采样。lags范围取[-N, N]是完整定义实际中常截断到[-N/2, N/2]来抑制交叉项这就是后面要讲的加窗 WVD。2.3 平滑伪魏格纳分布SPWVD抑制交叉项的主力原始 WVD 对多分量信号几乎不可用工程上真正落地的是它的加窗版本。最常用的是平滑伪魏格纳分布SPWVD在时间方向和频率方向各加一个窗$$SPW_x(t,f)\int\int h(\tau)g(u),x\left(tu\frac{\tau}{2}\right)x^*\left(tu-\frac{\tau}{2}\right)e^{-j2\pi f\tau},du,d\tau$$g(u)是时间平滑窗h(τ)是频率平滑窗。两个窗一加交叉项被大幅压制代价是自项也被展宽时频聚集性下降。这就是时频分析里绕不开的取舍交叉项和分辨率不可兼得只能按你的信号特性调窗长。def spwvd(x, fs, win_t31, win_f31): 平滑伪魏格纳分布, 用两个汉宁窗分别平滑时间和频率 from scipy.signal import hilbert, resample x np.asarray(x, dtypefloat) N len(x) xa hilbert(x) xa2 resample(xa, 2 * N) N2 len(xa2) gt np.hanning(win_t) # 时间平滑窗 hf np.hanning(win_f) # 频率平滑窗 W np.zeros((N, N), dtypecomplex) for ti in range(N): acc np.zeros(2 * N, dtypecomplex) for ui in range(-win_t // 2, win_t // 2 1): t_shift ti ui if t_shift 0 or t_shift N: continue lags np.arange(-N, N) idx_p t_shift lags idx_m t_shift - lags valid (idx_p 0) (idx_p N2) (idx_m 0) (idx_m N2) r np.zeros(len(lags), dtypecomplex) r[valid] xa2[idx_p[valid]] * np.conj(xa2[idx_m[valid]]) acc gt[ui win_t // 2] * r spec np.fft.fftshift(np.fft.fft(acc, n2 * N)) W[:, ti] spec[N // 2: N // 2 N] f np.linspace(0, fs / 2, N) t np.arange(N) / fs return t, f, W逻辑说明外层循环遍历时间点内层循环在时间方向用gt加权累加多个时间偏移的瞬时自相关这就是g(u)的作用。频率方向的平滑h(τ)在实现里通常通过对滞后序列加窗实现这里为了代码清晰先省略实际工程中会在r上再乘一个频率窗。参数说明win_t和win_f是核心调参对象。win_t越大时间方向平滑越强交叉项压得越狠但时间分辨率越差win_f同理作用于频率方向。经验值信号分量间隔近时窗要小分量间隔远、交叉项严重时窗可以大。一般从 31 或 51 起步看图上交叉项是否还在再微调。3. 用 part3.zip 里的时频分布代码跑通第一个例子3.1 构造一个能看出问题的测试信号不要一上来就拿真实数据跑先用合成信号验证代码对不对。构造一个两分量 LFM 信号两个分量在时频平面上交叉这样交叉项一定会出现方便你判断抑制效果。import numpy as np import matplotlib.pyplot as plt fs 1000.0 # 采样率 1kHz T 1.0 # 时长 1s t np.arange(0, T, 1 / fs) N len(t) # 分量1: 频率从 50Hz 线性升到 200Hz f1 50 150 * t / T s1 np.cos(2 * np.pi * np.cumsum(f1) / fs) # 分量2: 频率从 200Hz 线性降到 50Hz f2 200 - 150 * t / T s2 np.cos(2 * np.pi * np.cumsum(f2) / fs) x s1 s2逻辑说明用cumsum而不是直接f*t来生成相位是因为 LFM 的瞬时频率是相位的导数直接乘会得到错误的调频斜率这是新手最常翻车的地方。两个分量一个升一个降在时频图上会形成一个 X 形交叉交叉点就是交叉项最严重的位置。3.2 对比 WVD 与 SPWVD 的时频图t_axis, f_axis, W wvd(x, fs, n_freq256) _, _, Ws spwvd(x, fs, win_t51, win_f51) fig, axes plt.subplots(1, 2, figsize(12, 5)) for ax, M, title in zip(axes, [W, Ws], [WVD, SPWVD]): ax.pcolormesh(t_axis, f_axis, np.abs(M), shadingauto, cmapjet) ax.set_xlabel(Time (s)) ax.set_ylabel(Frequency (Hz)) ax.set_title(title) plt.tight_layout() plt.show()逻辑说明np.abs(M)取模得到能量分布pcolormesh画时频图。左图 WVD 你会看到两条清晰的斜线但中间交叉区域有一团剧烈振荡的伪影那就是交叉项。右图 SPWVD 交叉项被压下去但两条线变粗了这就是分辨率的代价。参数说明n_freq256控制频率轴点数越大图越细但计算越慢。win_t51是时间平滑窗长你可以改成 11、101 分别跑一遍直观感受窗长对交叉项和分辨率的影响。这一步建议手动改参数多跑几次比看任何理论都管用。3.3 从 part3.zip 迁移到自己的数据part3.zip里的时频分布代码大概率是 MATLAB 或 Python 的脚本集合迁移时注意三件事。第一采样率必须从你的采集配置里读不能沿用示例里的 1000。第二如果原代码直接对实信号做 WVD务必在前面补上希尔伯特变换否则你会看到正负频率对称的鬼影。第三数据长度如果不是 2 的幂FFT 会补零频率轴刻度要按实际n_freq重算别直接用fs/2除以点数。# 迁移模板: 读自己的数据 import scipy.io as sio data sio.loadmat(your_signal.mat)[signal].ravel() fs_real 25600.0 # 改成你的真实采样率 t_r, f_r, W_r spwvd(data, fs_real, win_t101, win_f101)逻辑说明ravel()把列向量拉平MATLAB 存的数据经常是 (N,1) 形状不拉平后面索引会出错。win_t和win_f要根据你的信号带宽和时长重新定采样率越高、信号越长窗可以适当加大。4. 魏格纳分布代码的避坑与排查清单4.1 图上出现对称鬼影现象时频图在正频率和负频率各有一条一模一样的脊或者频率轴出现镜像。原因直接对实信号做了 WVD没有取解析信号。实信号的频谱共轭对称WVD 会把负频分量也映射进来。解决在计算前统一加xa hilbert(x)所有后续操作基于xa。这是最高频的翻车点没有之一。4.2 交叉项压不下去现象调大win_t后交叉项还在甚至自项也糊了。原因交叉项的位置和强度取决于两个分量的时频距离如果两个分量靠得很近任何平滑窗都会同时伤到自项。解决先确认交叉项是不是真的来自双线性如果信号本身是单分量那图上的杂散可能是数值问题比如上采样没做。多分量且靠得近时考虑改用其他分布比如 Choi-Williams 分布或重排谱图不要死磕 WVD。4.3 频率轴刻度对不上现象明明信号主频在 200Hz图上脊线却出现在别的位置。原因离散 WVD 做了 2 倍上采样和 FFT 补零频率轴不能简单用fs/2均分。解决频率轴按np.linspace(0, fs/2, n_freq)生成n_freq是你 FFT 后实际取的频率点数不是原始信号长度。如果用了resample确认fs传的是原始采样率而不是上采样后的。4.4 长信号内存爆炸现象信号几万点跑 WVD 直接 MemoryError。原因朴素实现是O(N^2)甚至O(N^3)的复杂度每个时间点都要构造长度2N的滞后序列。解决分段处理每段加重叠窗或者用 FFT 加速的快速 WVD 算法。工程上更实际的做法是先降采样到你能接受的长度或者直接用 SPWVD 的分块实现。别硬扛几万点以上的信号用朴素循环跑 WVD 是不现实的。4.5 窗长调参没有方向现象win_t从 11 试到 201不知道哪个对。原因没有量化指标全靠肉眼。解决用两个客观指标辅助——自项脊线的能量集中度脊线峰值/旁瓣均值和交叉项区域的能量占比。先固定一个指标扫一遍窗长选指标拐点。经验上窗长取信号长度的 1/20 到 1/10 是个合理起点再根据交叉项残留微调。5. 进阶用重排和自适应核把魏格纳分布做到可用原始 WVD 和固定窗 SPWVD 都有明显短板真正在工程里落地时我一般会往两个方向走。第一个是重排reassignment把每个时频点的能量重新分配到其瞬时频率和群延迟的质心位置能在不牺牲太多分辨率的前提下锐化脊线。第二个是自适应核设计让核函数根据信号本身的模糊函数自动调整而不是用固定的汉宁窗。重排的核心思路是对 SPWVD 结果额外计算每个点的时间重心和频率重心然后把能量搬过去。def reassign_spwvd(x, fs, win_t51, win_f51): 对SPWVD结果做重排, 锐化时频脊线 t, f, W spwvd(x, fs, win_t, win_f) mag np.abs(W) # 时间重心: 沿时间轴的加权质心 t_idx np.arange(len(t)) t_centroid (mag * t_idx[None, :]).sum(axis1) / (mag.sum(axis1) 1e-12) # 频率重心: 沿频率轴的加权质心 f_idx np.arange(len(f)) f_centroid (mag * f_idx[:, None]).sum(axis0) / (mag.sum(axis0) 1e-12) # 重排: 把能量搬到重心位置 Wr np.zeros_like(W) for ti in range(len(t)): for fi in range(len(f)): if mag[fi, ti] 1e-6: continue ti_new int(np.clip(t_centroid[ti], 0, len(t) - 1)) fi_new int(np.clip(f_centroid[fi], 0, len(f) - 1)) Wr[fi_new, ti_new] W[fi, ti] return t, f, Wr逻辑说明t_centroid和f_centroid分别计算每个频率切片的时间重心和每个时间切片的频率重心这是重排的标准做法。1e-12是防止除零。重排后能量集中到脊线上图会明显变锐。注意重排对噪声敏感信噪比低时可能把噪声也锐化成假脊所以重排前最好先做一次轻度平滑。参数说明重排本身没有额外窗参数但它继承 SPWVD 的win_t和win_f。实践中我会先用较大的窗做 SPWVD 压交叉项再重排恢复分辨率这样两个目标都能兼顾一部分。验证重排效果的方法很简单看同一时刻的频率切片重排后主峰应该更窄、旁瓣更低。如果重排后出现孤立亮点多半是噪声被锐化了回去把窗调大或先降噪。一个我踩过的坑重排后的时频图不能直接用来做能量积分因为重排改变了能量的空间分布总能量虽然守恒但局部积分会失真。如果你的下游任务是提取瞬时频率或做分量分离重排图很好用如果要做能量统计还是用原始 SPWVD。这个区别我在一个轴承故障诊断项目里吃过亏当时用重排图算能量比结果和实际故障程度对不上排查半天才发现是重排的锅。最后说个习惯每次改完窗参数或换分布我都会固定跑那组两分量 LFM 合成信号看交叉项和脊线宽度有没有退化。合成信号是黑匣子测试真实数据是验收测试两个都不能省。时频分布这东西参数玄学多没有回归测试改着改着就不知道哪版是对的了。希望帮到你。本文还有配套的精品资源点击获取