信息处理 - 高斯滤波器不稳定差分的近似逆 - 吾爱随笔录

高斯滤波器不稳定差分的近似逆

信息处理过滤器设计卷积反卷积高斯差

2022-02-25 00:03:06

我正在尝试反转高斯（DoG）滤波器的差异。逆是不稳定的，所以我试图找到一个应用于特定输入的近似值。DoG 过滤器增加了边缘的对比度，我试图找到一个降低边缘对比度的逆。

近似的倒数不是很清楚，我的尝试让我接近但并不完全在那里。

输入

我的输入是一维的一系列上升步骤。它是这样创建的：

luminances = np.repeat(np.linspace(0, 1, 5), 100)

看起来像这样：

狗过滤器

我的 DoG 过滤器是用以下代码制作的：

filter_len = 81
excitatory_gaussian = scipy.signal.gaussian(filter_len, 1.4)
inhibitory_gaussian = scipy.signal.gaussian(filter_len, 16) * -0.043
filter = excitatory_gaussian + inhibitory_gaussian
filter /= sum(filter)
plt.plot(filter)

它看起来像这样：

卷积

当我将输入与 DoG 滤波器进行卷积时，我看到边缘的对比度更高，正如预期的那样：

我可以使用convolveor创建相同的输出lfilter：

deconvolved_luminances = signal.convolve(luminances, filter)
deconvolved_luminances = signal.lfilter(filter, [1], luminances)

反转尝试

我想找到 DoG 滤波器的倒数，这样，当将此输入与反向滤波器进行卷积时，输出是递增的，边缘对比度降低，并且可以与我的原始 DoG 滤波器进行卷积以输出原始输入的近似值.

信号.反卷积

似乎是因为 DoG 滤波器的逆滤波器不稳定，所以使用 signal.deconvolve 会导致输出错误。

deconvolved_luminances, remainder = scipy.signal.deconvolve(luminances, filter)

lfilter正如预期的那样，使用而不是时的输出是相同的deconvolve：

deconvolved_luminances = scipy.signal.lfilter([1], filter, luminances)

反相过滤器的一半

看起来过滤器的右半部分是稳定的。如果我将这个右半部分 ( half_filter_r = filter[len(filter)//2:]) 与我的信号 ( ) 反卷积s_deconvolved_r = signal.lfilter([1], half_filter_r, s)，我会得到一个有意义的输出，如下所示：

然而，取过滤器的左半部分 ( half_filter_l = filter[:len(filter)//2 + 1]) 并执行相同操作会产生一个输出，表明过滤器的左半部分不稳定：

我可以假设左半部分的适当输出只是右半部分输出的倒数s_deconvolved_l = max(s_deconvolved_r) - s_deconvolved_r[::-1]（s_deconvolved_r + s_deconvolved_l

但是，假设左半部分的输出基于右半部分似乎是可疑的，并且似乎应该有更好的方法。为什么左半边不稳定而右半边稳定？此外，当我将过滤器与上述输出进行卷积时，我得到的结果看起来与原始过滤器相似但不太正确（达到了我允许发布的图像的限制，但它在 colab 笔记本中）。

谢谢你的帮助！

3个回答

如果你有一个过滤器$\sum a_n z^n$，反卷积过滤器$1/ \sum a_n z^n$尝试通过添加过滤器的缩放版本来取消输入的最早样本，然后继续下一个样本（如高斯消除过程）。

如果滤波器中的第一个样本很小，则应用于滤波器的增益必须很大，这将导致后续样本变得更高。

是的，原件filter在单位圆之外有根：

所以倒置使用时会不稳定。

一种方法可能是形成最小等效相位：

minphase_filter = signal.minimum_phase(filter)

这将保证单位圆内的根。

产生的脉冲响应为：

这会产生滤波后的信号：

这允许s稳定地去卷积：

这是一个棘手的问题。粗略地说，这里的根本原因是原始过滤器“破坏”了无法恢复的信息，因此您必须做出一些取舍

设计反演的最佳方法是将反演公式化为最小二乘误差问题，其中误差标准用于为您的特定应用进行权衡。

假设您的传递函数是$H(z)$。假设我们要设计“某种逆”过滤器$G(z)$。所以理想情况下我们会有

$$H(z) \cdot G(z) = 1$$ 或 $$ G(z) = \frac{1}{H(z)}$$

但是，这仅在$H(z)$表现得相当好的情况下才有效，而您的则不然。特别是它不是最小相位，它在高频处具有“密集梳”，如果反转，将导致增益超过 100 dB 的窄带峰值。

您可以尝试估计您认为可以合理反转滤波器的频率范围。这是您的过滤器的函数，也是预期的本底噪声（物理或数字），然后定义目标函数$T(z)$。$T(z)$的一个不错的选择可能是一个低通滤波器，它的截止值与您的原始滤波器相似，它还应该包括足够的体延迟以允许反相成为因果关系。我可能会从两倍的过滤器长度开始。然后你可以通过最小化来解决$G(z)$

$$|G(z) \cdot H(z) - T(z)|^2 = min$$

这对于 FIR 滤波器来说非常简单，但对于需要一些智能搜索算法的 IIR 来说却很困难。

您可能还必须为最大系数值和/或增益添加一些加权函数和潜在的一些非线性约束。

我将从 FIR 开始，然后开始摆弄目标函数、权重、约束、滤波器长度、批量延迟，看看这是否“足够好”以供您做。

其它你可能感兴趣的问题

上一篇估计散粒噪声 - 起源是什么？下一篇音频，仅根据内容估计编解码器的压缩比