简介这份资源围绕Wigner-Hough变换展开面向从事非平稳信号处理、时频分析与故障诊断的工程师及科研人员帮助理解如何将Wigner分布与Hough变换结合以抑制交叉项干扰并检测时频域中的显著特征。压缩包共4个文件以3个m脚本文件和1个txt说明文件为主m文件用于实现Wigner分布计算、Hough变换投票与峰值检测等核心流程txt文件提供许可或使用说明整体约3KB体量轻便便于快速阅读与调试。目前已有205人学习下载适合作为入门与验证的参考。通过运行脚本读者可掌握从时频图生成、参数空间投票到特征回溯的完整思路理解周期性不明显或瞬时变化信号的识别方法并在此基础上迁移到自身项目中进行二次开发与实验验证。1. Wigner-Hough 变换从时频图里把斜线“捞”出来的工程做法雷达、声呐、通信侦察里经常遇到一类信号持续时间不长频率随时间线性变化业内叫线性调频信号简称 LFM。单看时域波形它跟噪声几乎没区别单看频谱它是一条被展宽的包络也看不出什么门道。真正能把这类信号从低信噪比背景里拎出来的是把 Wigner-Ville 分布和 Hough 变换串起来用也就是标题里的 Wigner-Hough 变换。Wigner-Ville 分布负责把一维时间信号变成二维时频图LFM 在这张图上表现为一条斜直线Hough 变换负责在二维图里做直线积分把这条斜线累积成一个尖峰。峰值位置对应起始频率和调频斜率峰值高度对应信号能量。这套组合在信噪比低到 -10 dB 量级时仍然能出结果是很多从业者处理非平稳信号时的默认选项。这篇笔记面向已经会写 FFT、但还没把时频检测跑通的工程师从原理到代码到参数到踩坑一步步把 Wigner-Hough 落地。2. Wigner-Ville 分布与 Hough 变换为什么这两个要拼在一起2.1 Wigner-Ville 分布到底给了我们什么Wigner-Ville 分布WVD是 Cohen 类时频分布里分辨率最高的一种它对信号做的是中心对称的双线性变换。离散形式下对长度为 N 的实信号 x(n)解析信号 z(n) 的 WVD 定义为W(n, k) Σ_m z(n m/2) · z*(n − m/2) · e^(−j2πkm/N)这个式子里 n 是时间索引k 是频率索引m 是滞后变量。跟短时傅里叶变换STFT比WVD 不需要窗函数所以不存在时间分辨率和频率分辨率互相牵制的问题理论上可以同时做到最好。代价是双线性带来的交叉项当信号里有多个分量时两两之间会产生虚假的时频能量位置在两个真实分量中间幅度还可能比真实分量高。对单个 LFM 信号来说WVD 的结果是一条清晰的斜线能量高度集中。这就是后面 Hough 变换能起作用的前提——如果时频图上的能量是散的直线积分就积不出尖峰。工程上常见的做法是先做解析信号变换用 Hilbert 变换去掉负频率再算 WVD避免正负频率之间的交叉项污染整个图。2.2 Hough 变换在时频图上的参数化Hough 变换本来是图像处理里检测直线的工具。标准形式用极坐标参数 (ρ, θ) 表示一条直线ρ x·cosθ y·sinθ。放到时频图上x 是时间y 是频率一条 LFM 对应的直线在 (ρ, θ) 空间里就是一个点反过来时频图上的所有点都在 (ρ, θ) 空间里投票真实直线对应的那个 (ρ, θ) 会累积出最大值。但时频图上的 LFM 检测有个更自然的选择直接用频率和时间的关系参数化。LFM 的瞬时频率 f f0 k·tf0 是起始频率k 是调频斜率。把时频图上的每个能量点 (t, f) 映射到 (f0, k) 空间每个点对应 (f0, k) 平面上的一条直线真实信号的 (f0, k) 处会形成累积峰值。这种参数化比极坐标更直观因为 f0 和 k 就是我们要估计的物理量不需要再做坐标反变换。实际实现时f0 和 k 的取值范围要预先划定。f0 从 0 到采样率的一半k 从 −k_max 到 k_maxk_max 由信号可能的最大调频斜率决定。这个范围划得越准累积矩阵越小计算越快但划得太窄会漏掉真实信号。我一般会先用 FFT 粗估一下信号占用的频段再定 f0 的范围。2.3 为什么不能只用 WVD 或只用 Hough只用 WVD低信噪比下时频图被噪声淹没人眼都看不清斜线更别说自动检测。只用 Hough输入是原始时域信号的话根本没有直线可检测。两者结合的逻辑是WVD 把信号能量从一维搬到二维并且沿直线集中Hough 沿直线做积分等效于对信号做相参积累噪声是非相参的积分后信噪比增益大致正比于直线长度。这就是为什么这套方法在低信噪比下仍然有效。代价是计算量。WVD 的复杂度是 O(N²)Hough 的复杂度取决于参数空间的大小通常是 O(N_t · N_f · N_k)N_k 是斜率的分辨率格点数。对长信号这个计算量不小。工程上常见的优化是先对 WVD 做阈值处理只保留超过门限的点再送 Hough能砍掉大部分无效投票。3. 用 Python 把 Wigner-Hough 跑通从解析信号到峰值检测3.1 生成测试用的 LFM 信号先构造一个干净的 LFM 信号加上高斯白噪声方便后面验证算法效果。import numpy as np import matplotlib.pyplot as plt def generate_lfm(fs, duration, f0, k, snr_db): 生成带噪声的 LFM 信号 fs: 采样率 Hz duration: 信号时长 s f0: 起始频率 Hz k: 调频斜率 Hz/s snr_db: 信噪比 dB t np.arange(0, duration, 1/fs) # 瞬时相位是频率的积分 phase 2 * np.pi * (f0 * t 0.5 * k * t**2) s np.exp(1j * phase) # 解析形式直接构造复信号 # 按信噪比加复高斯白噪声 signal_power np.mean(np.abs(s)**2) noise_power signal_power / (10**(snr_db / 10)) noise np.sqrt(noise_power/2) * (np.random.randn(len(t)) 1j*np.random.randn(len(t))) return t, s noise fs 1000 # 采样率 1 kHz duration 1.0 # 1 秒 f0 100 # 起始频率 100 Hz k 200 # 调频斜率 200 Hz/s snr_db -5 # 信噪比 -5 dB t, x generate_lfm(fs, duration, f0, k, snr_db)这里直接构造复解析信号省掉了 Hilbert 变换那一步。实际处理实信号时需要先做 Hilbert 变换得到解析信号否则 WVD 会出现正负频率的交叉项。参数上f0 和 k 的选择要保证瞬时频率 f0 k·t 始终在 0 到 fs/2 之间否则会产生混叠。上面这组参数下频率从 100 Hz 线性变到 300 Hz在 500 Hz 的奈奎斯特频率以内没问题。3.2 计算 Wigner-Ville 分布WVD 的离散实现有几种写法核心是对每个时间点做滞后方向的相关再 FFT。下面这个版本用矩阵化写法比双重循环快很多。def wigner_ville(x): 计算离散 Wigner-Ville 分布 x: 解析信号长度 N 返回: WVD 矩阵形状 (N, N)行是时间列是频率 N len(x) wvd np.zeros((N, N), dtypecomplex) # 对每个时间点 n计算滞后相关 for n in range(N): # 滞后 m 的范围受边界限制 m_max min(n, N-1-n) m np.arange(-m_max, m_max1) # 对称相关 r(m) z(nm) * conj(z(n-m)) r x[nm] * np.conj(x[n-m]) # 补零到 N 点后 FFT 得到该时刻的频率分布 r_padded np.zeros(N, dtypecomplex) r_padded[:len(r)] r wvd[n, :] np.fft.fftshift(np.fft.fft(r_padded)) return wvd W wigner_ville(x)这段代码里 m_max 的处理是关键在时间轴两端可用的滞后范围会缩小如果不做边界处理直接取全范围会引入虚假能量。补零到 N 点是为了让每个时间点的频率分辨率一致。fftshift 把零频移到中心方便后面画图和做 Hough 时频率轴的映射。计算量上外层循环 N 次每次 FFT 是 O(N log N)总体 O(N² log N)。对 N1000 的信号单次运行在普通笔记本上大约几秒可以接受。如果信号更长建议用 GPU 或者分帧处理。3.3 在时频图上做 Hough 变换拿到 WVD 矩阵后先取模值或者模值的平方作为能量图再做 Hough 累积。这里用 (f0, k) 参数化。def hough_lfm(W, fs, k_max, n_k200): 在时频图上做 Hough 变换检测 LFM W: WVD 矩阵 (N, N) fs: 采样率 k_max: 最大调频斜率绝对值 Hz/s n_k: 斜率分辨率格点数 返回: 累积矩阵 acc, f0 轴, k 轴 N W.shape[0] mag np.abs(W) # 只保留超过门限的点减少投票量 threshold np.mean(mag) 2 * np.std(mag) t_idx, f_idx np.where(mag threshold) # 时间轴和频率轴的实际值 t_vals t_idx / fs f_vals (f_idx - N//2) * fs / N # 参数空间 f0_axis np.linspace(0, fs/2, N) k_axis np.linspace(-k_max, k_max, n_k) acc np.zeros((len(f0_axis), len(k_axis))) # 对每个超过门限的点投票 for t, f in zip(t_vals, f_vals): # f f0 k*t f0 f - k*t f0_candidates f - k_axis * t # 找到落在 f0 轴范围内的索引 valid (f0_candidates 0) (f0_candidates fs/2) idx np.round(f0_candidates[valid] / (fs/2) * (len(f0_axis)-1)).astype(int) k_idx np.where(valid)[0] acc[idx, k_idx] mag[t_idx[0], f_idx[0]] # 用能量加权 return acc, f0_axis, k_axis门限那一步是工程上的关键优化。WVD 矩阵里大部分点的能量接近噪声水平如果全部拿去投票计算量翻几倍不说噪声的随机投票还会抬高累积矩阵的底噪让真实峰值不那么突出。用均值加两倍标准差做门限能滤掉大部分噪声点同时保留信号能量集中的区域。投票时用能量加权而不是简单计数是为了让强信号点的贡献更大弱噪声点的贡献更小。参数 n_k 控制斜率分辨率。n_k 太小相邻斜率分不开峰值会展宽n_k 太大每个格点分到的投票少峰值幅度下降而且计算量增加。经验上 n_k 取 100 到 500 之间比较合适具体看信号时长和调频斜率范围。k_max 要根据先验知识定如果完全不知道信号可能的最大调频斜率可以先设一个较宽的范围跑一遍看峰值落在哪里再缩小范围精跑。3.4 峰值检测与参数估计累积矩阵出来后找最大值的位置就是估计的 (f0, k)。k_max 500 acc, f0_axis, k_axis hough_lfm(W, fs, k_max) # 找全局最大值 peak_idx np.unravel_index(np.argmax(acc), acc.shape) f0_est f0_axis[peak_idx[0]] k_est k_axis[peak_idx[1]] print(f真实值: f0{f0} Hz, k{k} Hz/s) print(f估计值: f0{f0_est:.1f} Hz, k{k_est:.1f} Hz/s)如果信号里只有一个 LFM 分量全局最大值就够了。多个分量时需要做峰值提取找到最大值后把该峰值附近的一个邻域置零再找下一个最大值直到峰值低于某个门限。邻域大小取决于参数分辨率一般取累积矩阵中峰值半高宽的两倍。这一步没有后悔药邻域取小了会把同一个峰值的旁瓣当成第二个信号取大了会漏掉靠得近的两个真实信号。4. 参数怎么设WVD 长度、Hough 分辨率与门限的取舍4.1 信号长度与 WVD 计算量的平衡WVD 的计算量随信号长度平方增长。N1000 时几秒能跑完N4000 时可能要几分钟。但信号截短了时频图上的直线变短Hough 积分增益下降低信噪比下检测概率会掉。我一般会先估计信号的大致持续时间如果 LFM 只占整个采样时长的一小段就先做粗检测定位到大致时间段再截取那一段做精细 WVD。粗检测可以用 STFT虽然分辨率差但计算快能快速找到信号存在的区间。另一个思路是分帧做 WVD 再拼接但帧与帧之间的交叉项会污染拼接结果需要加窗和重叠处理实现起来比直接做全长 WVD 麻烦。除非信号特别长N 超过 10000否则不建议分帧。4.2 Hough 参数空间的分辨率选择f0 轴的分辨率通常跟 WVD 的频率分辨率对齐也就是 fs/N。k 轴的分辨率没有固定公式取决于信号时长 T 和允许的估计误差。如果要求调频斜率估计误差小于 Δk那么 k 轴的格点间距应该小于 Δk。但格点太密会导致每个格点累积的投票数减少峰值幅度下降。一个经验公式是k 轴格点数 n_k 取 T · k_max / (fs/N) 的量级也就是让 k 轴的分辨率跟 f0 轴的分辨率在时频图上对应的斜率变化量匹配。举个例子T1 sfs1000 HzN1000频率分辨率 1 Hz。如果 k_max500 Hz/s那么在 1 秒内频率变化 500 Hz对应 500 个频率格点。k 轴如果取 200 个格点每个格点对应 2.5 Hz/s 的斜率变化在 1 秒内对应 2.5 Hz 的频率变化比频率分辨率略粗可以接受。如果取 1000 个格点每个格点对应 0.5 Hz/s比频率分辨率还细但投票数会分散峰值幅度下降。我一般取 200 到 500 之间根据实际信噪比微调。4.3 门限设置对检测概率的影响WVD 门限设得太高弱信号点被滤掉Hough 累积的峰值幅度不够检测不到。设得太低噪声点大量参与投票累积矩阵底噪抬高真实峰值被淹没。均值加两倍标准差是一个保守的起点实际使用时可以画出门限后的时频图看看信号斜线是否还完整。如果斜线断成几截说明门限偏高降到均值加一倍标准差试试。如果时频图上全是散点说明门限偏低升到均值加三倍标准差。注意门限应该基于 WVD 模值的统计特性来定而不是固定值。不同信噪比下噪声的模值分布不同固定门限在信噪比变化时会失效。5. 避坑与排查Wigner-Hough 落地时最容易翻车的五个地方5.1 交叉项把真实峰值压下去了现象时频图上除了信号斜线还出现多条平行的虚假斜线Hough 累积后最大值对应的参数跟真实值对不上。原因WVD 是双线性变换当信号里有多个分量或者信号本身有幅度调制时分量之间会产生交叉项。交叉项的位置在两个真实分量中间幅度可能比真实分量还高。如果交叉项恰好也形成一条斜线Hough 会把它当成真实信号。解决先做解析信号变换去掉负频率分量能消掉正负频率之间的交叉项。如果信号本身有多个 LFM 分量考虑用平滑伪 Wigner-Ville 分布SPWVD或者 Choi-Williams 分布这些改进形式通过加核函数抑制交叉项代价是时频分辨率略有下降。另一个办法是先估计信号分量个数对消掉已知分量后再做 WVD。5.2 频率轴映射搞反了现象估计出的 f0 和 k 跟真实值符号相反或者 f0 落在负频率区域。原因WVD 做完 fftshift 后频率轴的中心是零频左边是负频率右边是正频率。如果 Hough 变换里频率轴的映射没有减去 N//2或者减的方向反了就会把正频率当成负频率。解决在计算 f_vals 时确认公式是 (f_idx - N//2) * fs / N。画一张时频图标出频率轴的实际值用已知频率的正弦信号验证一下。这个坑很隐蔽因为时频图看起来是对的只是坐标轴标错了但 Hough 累积时所有投票都偏了。5.3 信号时长估计错误导致 k 轴范围不够现象Hough 累积矩阵的最大值出现在 k 轴的边缘估计出的调频斜率刚好等于 k_max 或 −k_max。原因k_max 设小了真实信号的调频斜率超出了搜索范围峰值被截断在边缘。解决先设一个较大的 k_max 跑一遍看峰值是否落在边缘。如果是扩大 k_max 重跑。如果峰值在中间说明 k_max 合适。另一个办法是用 FFT 粗估信号带宽 B 和时长 T调频斜率的量级大约是 B/Tk_max 取这个值的两到三倍。5.4 噪声功率估计不准导致门限失效现象低信噪比下检测不到信号高信噪比下检测到一堆虚假峰值。原因门限用的是 WVD 模值的均值和标准差但 WVD 模值的分布不是高斯的尤其在信号存在时信号区域的模值会拉高均值和标准差导致门限被抬高弱信号点被滤掉。解决用 WVD 矩阵的中间区域没有信号的时间段估计噪声的均值和标准差而不是用全图。或者用中位数代替均值中位数对异常值更鲁棒。另一个办法是自适应门限先设一个低门限跑一遍 Hough找到候选峰值再根据候选峰值周围的能量分布调整门限重跑。5.5 峰值邻域抑制把真实信号吃掉了现象多分量场景下检测到第一个信号后第二个信号检测不到。原因峰值抑制的邻域取大了把第二个信号的峰值也置零了。或者两个信号的参数在累积矩阵里靠得太近邻域重叠。解决邻域大小应该根据累积矩阵中峰值的半高宽来定而不是固定值。先测量单个孤立峰值的半高宽取两倍作为抑制半径。如果两个信号参数确实很近考虑用 CLEAN 算法估计出第一个信号的参数后在时域重构该信号并减去再对残余信号重新做 WVD 和 Hough。这样能避免在累积矩阵里做邻域抑制带来的相互影响。6. 进阶技巧用重构对消提升多分量检测能力单分量检测跑通后多分量场景是下一个坎。两个 LFM 信号的参数如果相差不大WVD 上的两条斜线会交叉交叉点附近的交叉项能量很高Hough 累积时两个峰值会互相干扰。我试过几种做法最稳的是重构对消。思路很直接检测到第一个峰值后用估计的 (f0, k) 重构一个理想 LFM 信号幅度用峰值处的累积值反推相位用 f0 和 k 积分得到。然后在原始时域信号里减去这个重构信号对残余信号重新做 WVD 和 Hough。如果残余信号里还有第二个 LFM它的峰值就会干净地露出来。def reconstruct_lfm(t, f0, k, amp): 根据估计参数重构 LFM 信号 phase 2 * np.pi * (f0 * t 0.5 * k * t**2) return amp * np.exp(1j * phase) def iterative_detection(x, t, fs, k_max, n_iter3): 迭代检测多个 LFM 分量 x_res x.copy() results [] for i in range(n_iter): W wigner_ville(x_res) acc, f0_axis, k_axis hough_lfm(W, fs, k_max) peak_idx np.unravel_index(np.argmax(acc), acc.shape) f0_est f0_axis[peak_idx[0]] k_est k_axis[peak_idx[1]] amp_est np.abs(acc[peak_idx]) / len(t) # 粗略幅度估计 results.append((f0_est, k_est, amp_est)) # 重构并减去 x_recon reconstruct_lfm(t, f0_est, k_est, amp_est) x_res x_res - x_recon # 如果残余能量接近噪声水平提前退出 if np.mean(np.abs(x_res)**2) 1.5 * np.mean(np.abs(x)**2) * 10**(-snr_db/10): break return results幅度估计那一步是近似的用累积峰值除以信号长度。更准的做法是用最小二乘拟合固定 f0 和 k在时域上求最优幅度和相位。但迭代检测里用粗略估计就够了因为减不干净的部分会在下一轮被重新检测只要不把真实信号减过头就行。减过头的情况发生在幅度估计偏大时残余信号里会出现负的 LFM 分量WVD 上表现为一条反相的斜线Hough 累积后峰值位置不变但符号相反。判断方法是看残余信号的 WVD 上有没有负能量区域。迭代次数一般取 3 到 5 次。每次迭代后检查残余信号的能量如果降到噪声水平以下就提前退出。如果迭代到最大次数还有明显峰值说明要么信号分量超过预期要么某次幅度估计偏差太大导致对消不干净。验证检测结果是否可靠我习惯做两件事一是把估计的 (f0, k) 对应的直线画回时频图上看是否跟信号斜线重合二是在时域上重构所有检测到的分量跟原始信号做差看残余信号的频谱是否平坦。如果残余频谱还有明显尖峰说明漏检了分量。这两个验证步骤花不了几分钟但能避免把虚假峰值当成真实信号报出去。这套方法我从单分量调到多分量前后踩了大概两周的坑大部分时间花在门限和峰值抑制的调参上。后来发现与其在累积矩阵上做复杂的峰值处理不如回到时域做重构对消逻辑更干净参数更少。如果你也在做时频检测建议先把单分量跑通确认 WVD 和 Hough 的每个参数都理解到位再上多分量。希望帮到你。本文还有配套的精品资源点击获取