信息处理 - 估计一维过滤的过滤系数（卷积） - 吾爱随笔录 - 问答

估计一维过滤的过滤系数（卷积）

信息处理 matlab 过滤器冲动反应频域反卷积

2022-01-13 20:22:56

我有一个输出信号这是一个输入信号卷积与脉冲响应函数与一些添加噪声： $y$ $x$ $\star$ $h$ $n$

y (t) = h (t) ⋆ x (t) + n (t)

$y(t) = h(t) \star x(t) + n(t)$

我知道输入信号和输出信号并且想计算脉冲响应函数。我发现反卷积不像卷积那么简单，因为输入信号包含零，然后不会定义频域中的除法。如果找到两种方法，请在互联网上寻找“反卷积”的方法：维纳反卷积和正则化反卷积。维纳反卷积似乎更容易理解，所以我想尝试在 Matlab 中实现它（Matlab 函数给我关于输入信号在第一个条目处为零的错误，如果我阅读帮助文件，它似乎只适用于多项式？）。 $x$ $y$ $h$ deconv

因此，根据Wikipedia - Wiener Deconvolution解释，您想要找到以便： $g$

\hat{x} (t) = g (t) ⋆ y (t)

$\hat{x}(t) = g(t) \star y(t)$

但是在如何计算的定义中，他们使用了原始方程中的所有变量。他们也只展示了如何找到但可能和可以交换，因为卷积是可交换的，但我不确定两个向量的正确长度。目前它们的长度相同。 $G$ $x$ $x$ $h$

G (f) = \frac{H^{*} (f) S (f)}{| H (f) |^{2} S (f) + N (f)}

$G(f) = \frac{H^*(f) S(f)}{|H(f)|^2 S(f) + N(f)}$

在哪里：

$H = {\tt fft}(h)$
$G = {\tt fft}(g)$
$S =$ 的功率谱密度？这是吗？ $x$ ${\tt fft}(x)$
$N =$ 的平均功率谱密度，不太明白这是什么 $n$

我的问题是如何在不知道脉冲响应函数的情况下获得它（因为它在的定义和原始方程中都是如此）？由于我知道输入和输出，因此与找到具有已知脉冲响应函数的输入应该没有太大区别。我想在 Matlab 中做到这一点。 $h$ $G$

1个回答

这是一个很好的问题。
我将尝试使用 2 种方法（基本相同）来解决它。

解决方案是最小二乘解决方案：

\hat{h} = \arg min_{h} \frac{1}{2} {‖ h * x - y ‖}_{2}^{2}

$\hat{h} = \arg \min_{h} \frac{1}{2} \left\| h \ast x - y \right\|_{2}^{2}$

我们假设数据以有限离散形式给出（实际上是这样）。
卷积在有效模式下完成（如在 MATLABvalid属性中）。

直接解决方案

可以将上面的内容写成：

\hat{h} = \arg min_{h} \frac{1}{2} {‖ X h - y ‖}_{2}^{2}

$\hat{h} = \arg \min_{h} \frac{1}{2} \left\| X h - y \right\|_{2}^{2}$

基本上使用卷积的交换属性可以将上述作为构建卷积矩阵。 $x$ $x$

现在，解决方案是通常的最小二乘：

\hat{h} = \arg min_{h} \frac{1}{2} {‖ X h - y ‖}_{2}^{2} = {(X^{T} X)}^{- 1} X^{T} y

$\hat{h} = \arg \min_{h} \frac{1}{2} \left\| X h - y \right\|_{2}^{2} = \left( {X}^{T} X \right)^{-1} {X}^{T} y$

请注意，对于长信号，矩阵是巨大的。
因此，这种方法仅适用于小信号。

迭代解决方案

好吧，可以使用梯度下降法来最小化成本函数：

f (h) = \frac{1}{2} {‖ h * x - y ‖}_{2}^{2}

$f \left( h \right) = \frac{1}{2} \left\| h \ast x - y \right\|_{2}^{2}$

梯度由（相关是卷积的伴随算子）给出：

\frac{d}{d h} f (h) = (h * x - y) ⋆ h

$\frac{d}{d h} f \left( h \right) = \left( h \ast x - y \right) \star h$

其中代表互相关。 $\star$

在 MATLAB 代码中，它由下式给出：

hObjFun = @(vH) 0.5 * mean( (conv2(vX, vH, 'valid') - vY) .^ 2 );
vG      = conv2(flip(vX, 1), (conv2(vX, vHEst, 'valid') - vY), 'valid');

所有数据都在列向量中（conv2()使用它是因为它更快，因为conv()它只是它的包装器）。

这种方法速度快，不受大小的限制（或者更小，准确地说）。

完整的 MATLAB 代码在我的StackExchange Signal Processing Q18993 GitHub 存储库中（查看SignalProcessing\Q18993文件夹）。

其它你可能感兴趣的问题

上一篇使用python实时在线过滤数字信号下一篇卡尔曼滤波器协方差