简介这份资源围绕Wigner-Hough变换展开面向从事非平稳信号处理、时频分析与故障诊断的工程师及科研人员帮助理解如何将Wigner分布与Hough变换结合以抑制交叉项干扰并检测时频域中的显著特征。压缩包共4个文件以3个m脚本文件和1个txt说明文件为主m文件用于实现Wigner分布计算、Hough变换投票与峰值检测等核心流程txt文件提供许可或使用说明整体约3KB体量轻便便于快速阅读与调试。目前已有205人学习下载。通过其中的MATLAB代码读者可掌握从二维自相关运算生成时频图、参数空间投票到回溯提取时频特征的完整思路并理解该方法在周期性不明显或瞬时变化信号识别中的应用价值适合作为时频分析入门与算法验证的参考素材。1. 从一段被噪声淹没的雷达回波说起Wigner-Hough 到底在解决什么一段实测雷达回波信噪比掉到 -5 dB 以下目标做匀加速运动回波在时频面上是一条斜着弯曲的脊线。你盯着频谱看什么都看不出来你换成短时傅里叶变换时间窗一短频率分辨率就崩窗一长脊线又被抹平。这时候 Wigner-Hough 变换就该上场了它先用 Wigner-Ville 分布把信号铺到时频面上再用 Hough 变换把时频面里那条直线或近似直线的能量脊累加成一个尖峰。噪声是散的脊线是连续的累加之后信噪比能拉出十几 dB 的增益。这个方案适合谁做雷达/声呐/通信信号检测的工程师做机械故障诊断里处理非平稳振动信号的还有做生物医学信号比如心电、肌电里 chirp 成分提取的。它解决的核心问题是在低信噪比下检测和参数估计线性调频LFM或多项式相位信号。代价也很直接——Wigner-Ville 有交叉项Hough 有计算量两个坑叠在一起参数没调好就是一张糊图。下面按我实际落地的顺序把这条路走一遍。2. Wigner-Ville 分布为什么它比短时傅里叶更适合做 Hough 的前置2.1 时频聚集性决定了 Hough 能不能累出峰Hough 变换的本质是把时频面上一条曲线上的点做积分。积分要出峰前提是这条曲线上的能量足够集中。短时傅里叶变换受海森堡不确定性原理约束时间分辨率和频率分辨率此消彼长一条 LFM 脊线在 STFT 面上是有一条宽度的带Hough 积分时带宽内的能量被摊薄峰值高度上不去。Wigner-Ville 分布WVD是双线性时频分布对单分量 LFM 信号它的时频聚集性理论上可以达到任意好——一条 LFM 在 WVD 面上就是一条理想的直线冲激脊。这就是 Wigner-Hough 的立论基础WVD 负责把信号能量压到最窄的脊上Hough 负责把这条脊上的能量累加成一个可检测的峰。但双线性带来一个绕不开的代价交叉项。两个分量之间会产生一个虚假的时频脊位置在两个真实分量中间振荡剧烈。如果直接对多分量信号做 WVD 再 Hough交叉项也会被累加出现虚假峰。所以工程上常见做法是先做 Cohen 类核函数平滑比如伪 Wigner-Ville 分布 PWVD、平滑伪 Wigner-Ville 分布 SPWVD牺牲一点聚集性换交叉项抑制。2.2 离散 WVD 的实现与解析信号预处理直接对实信号做 WVD正负频率之间会产生交叉项。标准做法是先构造解析信号把负频率干掉。下面是我常用的最小实现不依赖任何工具箱纯 NumPy。import numpy as np from scipy.signal import hilbert def wvd(sig, fs): 离散 Wigner-Ville 分布 sig: 实信号或解析信号一维 fs : 采样率 返回: t (时间轴), f (频率轴), W (时频矩阵, 行频率, 列时间) # 1. 构造解析信号抑制正负频率交叉项 z hilbert(sig) if np.isrealobj(sig) else sig N len(z) # 2. 时延轴从 -(N-1) 到 (N-1) tau np.arange(-(N-1), N) # 3. 初始化时频矩阵 W np.zeros((N, N), dtypecomplex) # 4. 逐时刻计算瞬时自相关 r(t, tau) z(ttau/2) * conj(z(t-tau/2)) for ti in range(N): for tau_i, tau_val in enumerate(tau): idx1 ti tau_val // 2 idx2 ti - tau_val // 2 if 0 idx1 N and 0 idx2 N: W[ti, tau_i] z[idx1] * np.conj(z[idx2]) # 5. 对 tau 维做 FFT得到频率维 W np.fft.fftshift(np.fft.fft(W, axis1), axes1) f np.fft.fftshift(np.fft.fftfreq(N, 1/fs)) t np.arange(N) / fs return t, f, W.T # 转置成 行频率, 列时间逻辑说明第 1 步用 Hilbert 变换构造解析信号这是抑制交叉项最省事的一步不做的话正负频率会在零频附近打架。第 4 步是 WVD 的定义式核心瞬时自相关函数在时延维上采样。第 5 步对时延维做 FFT 就得到频率维。参数说明fs只影响频率轴刻度不影响矩阵数值。N是信号长度注意这个双重循环是 O(N²)N 超过 2048 就会明显卡顿实际工程里要么分段处理要么用 FFT 基的快速算法。我一般会把 N 控制在 512 到 1024 之间配合滑窗。提示如果信号是多分量且交叉项严重把z换成先做 SPWVD 平滑的结果或者直接在 WVD 矩阵上做二维低通滤波再送 Hough。3. Hough 变换把时频面上的脊线累加成峰3.1 参数空间的选择为什么用极坐标而不是斜率截距时频面上的直线可以写成 f k·t b但斜率 k 在直线接近垂直时会趋于无穷参数空间发散。标准做法是转成极坐标ρ t·cosθ f·sinθ。每条时频面上的直线对应参数空间 (ρ, θ) 上的一个点而时频面上的一个点对应参数空间里的一条正弦曲线。所有正弦曲线的交点就是直线参数。对 LFM 信号时频脊线是 f f0 μ·t其中 μ 是调频率。映射到极坐标后θ 和 μ 有对应关系ρ 和 f0 有对应关系。所以 Hough 参数空间里的峰值位置直接给出起始频率和调频率估计。3.2 从 WVD 矩阵到 Hough 参数空间的完整代码def wigner_hough(sig, fs, theta_res180, rho_res200): Wigner-Hough 变换主流程 sig : 输入实信号 fs : 采样率 theta_res: 角度分辨率参数空间列数 rho_res : 距离分辨率参数空间行数 返回: H (Hough 参数空间), theta, rho, 以及峰值对应的 (f0, mu) t, f, W wvd(sig, fs) Wmag np.abs(W) # 1. 归一化并做软阈值压制噪声底 Wmag Wmag / (Wmag.max() 1e-12) Wmag[Wmag 0.1] 0 # 2. 构建 Hough 参数空间 Nt, Nf Wmag.shape theta np.linspace(-np.pi/2, np.pi/2, theta_res) # rho 范围由时频面尺寸决定 rho_max np.sqrt(Nt**2 Nf**2) rho np.linspace(-rho_max, rho_max, rho_res) H np.zeros((rho_res, theta_res)) # 3. 遍历时频面所有非零点投票 t_idx, f_idx np.nonzero(Wmag) for ti, fi in zip(t_idx, f_idx): val Wmag[ti, fi] for thi, th in enumerate(theta): r ti * np.cos(th) fi * np.sin(th) ri int((r rho_max) / (2 * rho_max) * (rho_res - 1)) if 0 ri rho_res: H[ri, thi] val # 4. 找峰值 ri, thi np.unravel_index(np.argmax(H), H.shape) rho_peak rho[ri] theta_peak theta[thi] # 5. 反算调频率和起始频率需根据时频面索引与实际物理量映射 # 时频面索引 ti 对应时间 t[ti]fi 对应频率 f[fi] # 直线在索引域: fi k * ti b, 极坐标 rho ti*cos fi*sin # 由 theta_peak 可解出 k -cos(theta)/sin(theta)sin 不为 0 时 if abs(np.sin(theta_peak)) 1e-6: k_idx -np.cos(theta_peak) / np.sin(theta_peak) # 索引斜率转物理调频率 dt t[1] - t[0] df f[1] - f[0] mu k_idx * df / dt else: mu np.inf return H, theta, rho, mu, theta_peak, rho_peak逻辑说明第 1 步的软阈值是关键WVD 面里噪声底虽然低但数量多不压制的话投票会被噪声淹没。第 3 步是标准 Hough 投票每个非零点对参数空间里所有 θ 投一票票重是该点的幅度。第 5 步的反算需要小心索引和物理量的映射这是最容易翻车的地方。参数说明theta_res决定调频率分辨率180 够用要更细就上 360。rho_res决定起始频率分辨率200 是折中。阈值 0.1 是经验值信噪比低时降到 0.05信噪比高时提到 0.2。注意第 3 步的双重循环在 N1024、theta_res180 时大约 180 万次投票Python 纯循环要几秒。生产环境用向量化或 C 扩展或者先对时频面做稀疏化只保留局部极大值点。4. 参数怎么设采样率、窗长、阈值与分辨率的联动4.1 采样率和信号长度的约束关系WVD 的频率分辨率是 fs/NHough 的调频率分辨率受限于时频面的时间和频率分辨率。要检测的最小调频率 μ_min 满足在信号时长 T 内频率变化量 Δf μ_min·T 必须大于频率分辨率 fs/N否则脊线在时频面上根本看不出斜率。反过来最大可检测调频率受 Nyquist 约束μ_max·T fs/2。我一般这样定先确定要检测的调频率范围 [μ_min, μ_max]然后反推需要的信号时长 T fs/(N·μ_min)以及采样率 fs 2·μ_max·T。这两个不等式一联立采样率和长度的可行域就出来了。4.2 窗长和重叠的选择WVD 本身不需要窗但实际信号往往分段处理。段长 N 决定频率分辨率重叠率决定时间轴平滑度。我的经验是段长取信号中一个 chirp 持续时间的 1.5 到 2 倍重叠 50%。段太长交叉项严重段太短频率分辨率不够。参数典型值影响调整方向段长 N512~1024频率分辨率 fs/N信噪比低时加大重叠率50%时间轴连续性调频率变化快时加大theta_res180~360调频率分辨率精度要求高时加大rho_res200~500起始频率分辨率同上幅度阈值0.05~0.2投票稀疏度信噪比低时降低4.3 交叉项抑制的三种工程手段第一种是解析信号预处理前面代码已经做了。第二种是核函数平滑用 PWVD 或 SPWVD 替代 WVD代价是脊线变宽Hough 峰值降低。第三种是时频面后处理对 WVD 矩阵做二维中值滤波或形态学开运算把细碎的交叉项抹掉再送 Hough。我一般先用第一种不够再上第二种。第三种会改变能量分布用的时候要重新标定阈值。三种手段不要同时上否则脊线被抹得连自己都认不出来。5. 避坑与排查Wigner-Hough 落地时最容易翻车的五个地方5.1 现象Hough 参数空间出现多个峰分不清哪个是真目标原因WVD 交叉项在时频面上形成了虚假脊线Hough 把它们也累加了。多分量信号或者信号本身有谐波时特别明显。解决先确认信号分量数。单分量直接上 WVD多分量先做 SPWVD 平滑。如果平滑后还有多峰检查是不是信号本身就有多个 chirp 分量那就不是交叉项而是真实目标需要按峰值高度和脊线连续性做二次筛选。5.2 现象峰值位置对但反算出的调频率偏差大原因索引域到物理域的映射搞错了。时频面的行索引对应频率列索引对应时间但np.fft.fftshift之后频率轴有偏移np.nonzero返回的索引和物理量之间差一个偏移量。解决反算前先把索引转成物理量再算斜率不要直接在索引域算完再乘系数。我习惯在 Hough 投票时直接用物理坐标 (t, f) 而不是索引 (ti, fi)虽然慢一点但不会错。5.3 现象低信噪比下 Hough 面一片模糊找不到峰原因WVD 面对噪声没有处理噪声点在参数空间里均匀投票把真峰的对比度拉低了。解决投票前做软阈值把幅度低于噪声底 2 到 3 倍的点置零。噪声底可以用时频面幅度中位数估计。另外可以只对时频面的局部极大值点投票点数减少一个数量级真峰反而更突出。5.4 现象程序跑得极慢N2048 时几分钟出不来原因WVD 的 O(N²) 双重循环加上 Hough 的 O(N²·theta_res) 投票两个瓶颈叠在一起。解决WVD 用 FFT 基的快速算法或者用scipy.signal.spectrogram的变体近似。Hough 投票用 NumPy 向量化把所有非零点一次性对 theta 数组做外积。再不行就降 N分段处理。5.5 现象单分量信号 Hough 峰值也不尖锐原因WVD 面里脊线本身有宽度尤其是信号有幅度调制或者频率不是严格线性时脊线在时频面上是弯的Hough 的直线模型匹配不上。解决如果调频率是常数但脊线弯检查是不是信号有幅度包络导致 WVD 脊线展宽。如果调频率本身在变多项式相位信号Hough 要换成广义 Hough 或者多项式 Hough参数空间维度要加。6. 进阶技巧用脊线提取替代全平面投票把速度拉起来全平面 Hough 投票的浪费在于时频面里 90% 以上的点都是噪声它们对真峰的贡献是负的。一个更聪明的做法是先做脊线提取只对脊线上的点投票。脊线提取可以用时频面的局部极大值追踪或者用图像处理里的细化算法。我常用的流程是WVD 面 → 二维高斯平滑 → 局部极大值检测 → 沿时间轴做脊线追踪相邻时刻极大值位置差不超过阈值→ 只对追踪到的脊线点做 Hough 投票。这样投票点数从 N² 降到 N 量级速度提升两个数量级而且因为噪声点被排除峰值对比度反而更高。验证方法拿一个已知调频率和起始频率的仿真 LFM 信号加不同信噪比的高斯白噪声跑 100 次蒙特卡洛统计调频率估计的均方根误差。信噪比从 -10 dB 到 10 dB看误差曲线在哪个点开始急剧上升那个点就是这个参数配置下的检测门限。我一般要求门限比理论值差不超过 3 dB超了就回去调阈值和分辨率。还有一个技巧是参数空间的非相干积累把多段短信号的 Hough 面直接相加因为真峰位置一致噪声峰位置随机相加之后真峰相对高度提升。这相当于用时间换信噪比适合信号持续时间长但单段信噪比低的情况。最后说个我自己的习惯每次调完参数先把 WVD 面和 Hough 面都画出来看一眼不要只看最终估计值。很多时候估计值看着对但时频面上脊线已经断了换个信号就翻车。图上看一眼比跑十次蒙特卡洛都管用。希望帮到你。本文还有配套的精品资源点击获取
