4. 离散反卷积的矩阵形式与直接求解

在连续时间域里，卷积写成积分；到了数字信号处理中，我们面对的是一串采样点，卷积就变成了求和。把这段求和按行展开、按列对齐，就得到一个非常清晰的线性方程组。本节的核心任务，就是把这个方程组写成矩阵形式，理解直接逆滤波在什么条件下会失效，并用矩阵求逆与伪逆完成一次基础反卷积实验。

一、从卷积求和到矩阵方程

设输入序列为 x[n]（长度 N），系统冲激响应为 h[n]（长度 M），输出为 y[n]。离散线性卷积定义为：

y[n] = Σ_{k=0}^{M-1} h[k] · x[n-k],  n = 0,1,...,N+M-2

把每个 y[n] 写成一行，把 x[0..N-1] 作为未知列向量，则 h 的各个移位版本就构成了系数矩阵的行。以 M=3、N=4 为例，输出长度为 6，矩阵形式为：

⎡ y0 ⎤   ⎡ h0  0   0   0 ⎤   ⎡ x0 ⎤
⎢ y1 ⎥   ⎢ h1  h0  0   0 ⎥   ⎢ x1 ⎥
⎢ y2 ⎥ = ⎢ h2  h1  h0  0 ⎥ · ⎢ x2 ⎥
⎢ y3 ⎥   ⎢ 0   h2  h1  h0⎥   ⎢ x3 ⎥
⎢ y4 ⎥   ⎢ 0   0   h2  h1⎥
⎣ y5 ⎦   ⎣ 0   0   0   h2⎦

这个矩阵称为卷积矩阵（Toeplitz 结构：每条对角线上的元素相同）。它的行数大于列数（6×4），是一个“高瘦”的长方阵。反卷积就是已知 y 和 H，反求 x，即解方程 Hx = y。

离散卷积的矩阵形式：Hx = y 卷积矩阵 H（6×4，Toeplitz） h0000 h1h000 h2h1h00 0h2h1h0 00h2h1 000h2 × 未知 x（4×1） x0 x1 x2 x3 = 观测 y（6×1） y0 y1 y2 y3 y4 y5 求解 路径 A：方阵求逆 x = H⁻¹y（仅当 H 为方阵且可逆） 条件：det(H) ≠ 0，无噪声 路径 B：伪逆（最小二乘） x̂ = (HᵀH)⁻¹Hᵀ y 适用：H 为长方阵，y 含噪声 风险：HᵀH 病态时误差放大 直接逆滤波失效条件： ① H 奇异或接近奇异（det≈0）→ 逆不存在或数值爆炸；② 输出含噪声 y = Hx + n → H⁻¹n 被高频放大；③ 卷积矩阵条件数大 → 微小扰动导致解严重偏离。

二、直接逆滤波为什么会失效

最朴素的想法是：既然 y = Hx，那 x = H⁻¹y 不就行了？问题在于，这个“求逆”在真实场景中往往不成立，原因可以归为三类。

第一，矩阵可能根本不可逆。 卷积矩阵是长方阵（行数 N+M−1 大于列数 N），严格来说没有通常意义的逆矩阵。只有在特定边界条件下（例如循环卷积、方阵形式）才可能构造出方阵。若冲激响应 h 的频谱存在零点（例如 h = [1, -1]，其频响在直流处为零），对应的卷积矩阵就是奇异的，逆不存在。

第二，噪声被无限放大。 实际观测总是 y = Hx + n。直接求逆得到 x̂ = H⁻¹y = x + H⁻¹n。卷积矩阵的特征值对应 h 的频响幅度；当某些频段 |H(f)| 很小，1/|H(f)| 就很大，该频段的噪声被成倍放大。这就是逆滤波的经典病态问题：信号被恢复的同时，噪声也淹没了结果。

第三，条件数决定敏感度。 矩阵条件数 κ(H) = σ_max / σ_min（奇异值之比）衡量解对扰动的敏感程度。κ 越大，y 的微小误差会导致 x̂ 的巨大偏差。对于低通特性的 h，高频段奇异值接近零，条件数极大，直接求逆几乎不可用。

三、用伪逆做基础反卷积实验

当 H 是长方阵时，方程 Hx = y 通常没有精确解（方程数多于未知数）。此时的标准做法是最小二乘：寻找使 ‖Hx − y‖² 最小的 x。对目标函数求导并令其为零，得到正规方程 HᵀHx = Hᵀy，其解为：

x̂ = (HᵀH)⁻¹ Hᵀ y  =  H⁺ y

其中 H⁺ 称为 Moore-Penrose 伪逆。下面给出一个可直接运行的教学示例（参数为合理假设，非实测数据）。设真实输入 x = [1, 2, 1, 0]，系统 h = [1, 0.5, 0.25]，先正演得到 y，再叠加少量噪声，最后用伪逆恢复。

import numpy as np

h = np.array([1.0, 0.5, 0.25])
x_true = np.array([1.0, 2.0, 1.0, 0.0])

# 构造卷积矩阵 H（6×4）
N, M = len(x_true), len(h)
H = np.zeros((N + M - 1, N))
for i in range(N + M - 1):
    for j in range(N):
        k = i - j
        if 0 <= k < M:
            H[i, j] = h[k]

y_clean = H @ x_true
rng = np.random.default_rng(0)
y_noisy = y_clean + rng.normal(0, 0.01, y_clean.shape)

# 伪逆求解
x_hat = np.linalg.pinv(H) @ y_noisy
print("条件数 κ(H) =", np.linalg.cond(H))
print("恢复 x̂ =", np.round(x_hat, 4))

运行后可以观察到两个关键现象。其一，κ(H) 的数值直接反映问题难度：若 h 的低频增益远大于高频增益，条件数会偏大，恢复误差随之上升。其二，即使噪声标准差只有 0.01，恢复出的 x̂ 也可能在尾部出现明显偏差——这正是伪逆在病态矩阵上仍会放大噪声的体现。判断实验是否“成功”，不能只看 x̂ 与 x_true 是否接近，还要结合条件数、残差 ‖Hx̂ − y‖ 以及噪声水平综合评估。

伪逆给出的是最小二乘意义下的最优解，但它并不“抗噪”。当 HᵀH 接近奇异时，(HᵀH)⁻¹ 的计算本身就会引入数值不稳定。工程上通常改用 SVD 截断、Tikhonov 正则化或迭代方法（如共轭梯度、Landweber 迭代）来抑制噪声放大，这些将在后续章节展开。

四、实验中的实用检查点

为了让基础反卷积实验可复现、可判断，建议在代码中固定以下检查项：构造 H 后先打印其形状与秩，确认是否为满列秩；计算条件数，若 κ 超过约 1e4，就要预期直接求逆结果不可靠；对比无噪声与含噪声两种输入下的恢复误差，观察误差随噪声方差的增长斜率；最后用残差范数验证解是否满足最小二乘意义下的最优性。把这些量记录下来，才能把“跑通了”变成“看懂了”。

矩阵形式的价值不仅在于能写出方程，更在于它把反卷积的困难量化了：秩、奇异值、条件数这些指标，直接告诉我们直接逆滤波在什么条件下会失效，以及伪逆能在多大程度上兜底。理解这一点，是后续引入正则化与迭代反演方法的必要前提。

← 上一章 📖 返回目录 下一章 →