简介：这份资源围绕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·t，f0 是起始频率，k 是调频斜率。把时频图上的每个能量点 (t, f) 映射到 (f0, k) 空间，每个点对应 (f0, k) 平面上的一条直线，真实信号的 (f0, k) 处会形成累积峰值。这种参数化比极坐标更直观，因为 f0 和 k 就是我们要估计的物理量，不需要再做坐标反变换。
实际实现时，f0 和 k 的取值范围要预先划定。f0 从 0 到采样率的一半，k 从 −k_max 到 +k_max，k_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), dtype=complex) # 对每个时间点 n，计算滞后相关 for n in range(N): # 滞后 m 的范围受边界限制 m_max = min(n, N-1-n) m = np.arange(-m_max, m_max+1) # 对称相关 r(m) = z(n+m) * conj(z(n-m)) r = x[n+m] * np.conj(x[n-m]) # 补零到 N 点后 FFT 得到该时刻的频率分布 r_padded = np.zeros(N, dtype=complex) 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)。对 N=1000 的信号，单次运行在普通笔记本上大约几秒，可以接受。如果信号更长，建议用 GPU 或者分帧处理。
3.3 在时频图上做 Hough 变换
拿到 WVD 矩阵后，先取模值或者模值的平方作为能量图，再做 Hough 累积。这里用 (f0, k) 参数化。
def hough_lfm(W, fs, k_max, n_k=200): """ 在时频图上做 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 的计算量随信号长度平方增长。N=1000 时几秒能跑完，N=4000 时可能要几分钟。但信号截短了，时频图上的直线变短，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 轴的分辨率在时频图上对应的斜率变化量匹配。
举个例子，T=1 s，fs=1000 Hz，N=1000，频率分辨率 1 Hz。如果 k_max=500 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/T，k_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_iter=3): """迭代检测多个 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 的每个参数都理解到位，再上多分量。希望帮到你。
本文还有配套的精品资源，点击获取
