去卷积去除一维信号中的高斯模糊(维纳滤波?)

信息处理 过滤 反卷积 1天
2022-01-19 09:27:44

我有一组试图去噪的生物学数据(实际上,人口统计数据只能用已知宽度的高斯卷积来测量)

我的问题是:我可以测量 (f*g),我可以测量 g(x),但我需要知道 f(x)。因为这是带有轻微斑点的真实数据,所以我不能只使用 deconvolve() 函数。

一位 EE 朋友告诉我,我应该使用 1-D Wiener Deconvolution。问题是,我不知道如何实现这一点。据说,它有一个 Matlab 包,但我无权访问 Matlab;并且我能找到的所有 Python 示例(例如 scikit-learn)都是 2D 的。

有什么我想念的吗?我可以实现一个简单的公式吗?

请不要只将我指向维基百科页面/一维维纳过滤的教科书。我试图阅读其中的几个,但它们看起来像是难以理解的胡言乱语(我是一名生物学家)。

编辑:点扩散函数又是怎么回事?我假设我的 PSF 应该是已知的高斯 g(x),但是我能找到的每个使用 Wiener 滤波的高斯反卷积示例(全部在 2d 中)都使用均匀分布作为其 PSF,而不是高斯分布。

谢谢!

4个回答

方法

有许多反卷积方法(即退化算子是线性的和时空不变的)。
在许多情况下,他们都试图处理问题是 Ill Poised 的事实。

更好的方法是那些为要恢复的数据模型添加一些正则化的方法。
它可以是统计模型(先验)或任何知识。
对于图像,一个好的模型是渐变的分段平滑或稀疏。

但是为了得到答案,将采用一种简单的参数化方法——最小化模型中恢复的数据与测量值之间的最小二乘误差。

模型

最小二乘模型很简单。
作为数据函数的目标函数由下式给出:

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

优化问题由下式给出:

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

其中$ x $是要恢复的数据,$ h $是模糊内核(在这种情况下为高斯),$ y $是给定测量值的集合。
该模型假设仅对卷积的有效部分给出测量值。即如果$ x \in \mathbb{R}^{n} $$ h \in \mathbb{R}^{k} $那么$ y \in \mathbb{R}^{m} $其中$ m = n - k + 1 $

这是有限空间中的线性运算,因此可以使用矩阵形式编写:

$$ \arg \min_{x} f \left( x \right) = \arg \min_{x} \frac{1}{2} \left\| H x - y \right\|_{2}^{2} $$

其中$H \in \mathbb{R}^{m \times n} $是卷积矩阵。

解决方案

最小二乘解由下式给出:

$$ \hat{x} = \left( {H}^{T} H \right)^{-1} {H}^{T} y $$

可以看出,它需要矩阵求逆。
充分解决这个问题的能力取决于运算符$ {H}^{T} H $的条件数,它服从$ \operatorname{cond} \left( H \right) = \sqrt{\operatorname{cond} \left ( {H}^{T} H \right)} $

条件数分析

这个条件编号的背后是什么?
可以使用线性代数来回答它。
但在我看来,更直观的方法是在频域中考虑它。

基本上,退化算子通常会衰减高频的能量。
现在,由于在频率上这基本上是元素乘法,有人会说反转它的简单方法是逆滤波器的元素除法。
嗯,就是上面做的。
问题出现在滤波器将能量实际上衰减到零的情况下。然后我们遇到了真正的问题......
这基本上是条件数告诉我们的,一些频率相对于其他频率衰减的程度。

在此处输入图像描述

上面可以看到条件数(使用 [dB] 单位)作为高斯滤波器 STD 参数的函数。
正如预期的那样,STD 越高,条件数越差,因为 STD 越高意味着 LPF 越强(最后下降的值是数值问题)。

数值解

创建了高斯模糊内核的集合。

在此处输入图像描述

参数为$ n = 300 $$ k = 31 $$ m = 270 $
数据是随机的,没有添加任何噪音。

在 MATLAB 中pinv(),使用基于 SVD 的伪逆和\算子来求解线性系统。

在此处输入图像描述

可以看出,使用 SVD 的解决方案不如预期的那么敏感。

为什么会出现错误?
寻找解决方案(对于最高 STD):

在此处输入图像描述

正如人们所看到的,除了开始和结束之外,信号恢复得非常好。
这是由于使用了有效卷积,它在这些样本上几乎没有告诉我们什么。

噪音

如果我们添加噪音,事情看起来会有所不同!
之前结果很好的原因是由于 MATLAB 可以处理数据的 DR 并求解方程,即使它们的条件数很大。

但是大的条件数意味着逆滤波器会强烈放大(以逆转强衰减)某些频率。
当那些包含噪音时,这意味着噪音会被放大并且恢复会很糟糕。

在此处输入图像描述

正如上面所见,现在重建将无法进行。

概括

如果一个人准确地知道降级算子并且 SNR 非常好,那么简单的反卷积方法将起作用。
反卷积的主要问题是退化算子衰减频率的难度。
它衰减得越多,恢复所需的 SNR 就越多(这基本上是Wiener Filter背后的想法)。
设置为零的频率无法恢复!

在实践中,为了得到稳定的结果,应该添加一些先验。

该代码可在我的StackExchange Signal Processing Q2969 GitHub 存储库中找到(查看SignalProcessing\Q2969文件夹)。

资源

这是一个系统识别问题。我们知道卷积是可交换的:$f(n)*g(n)=g(n)*f(n)$现在你有一个由已知系统$g(n )$过滤的未知信号$f(n )$ ,我们说输出是$d(n)$另一方面,我们也可以说$d(n)$是已知信号$g(n)$的输出,在您的情况下它是高斯噪声,由未知系统$f(n)$传递。您可以通过维纳滤波器或最小均方 (LMS) 等自适应滤波器来识别未知系统。系统块如下。

管理系统

在您的情况下,我们假设没有干扰并且$y(n)=d(n)$误差信号$e(n)=\hat{y}(n)-d(n)$,其中$\hat{y}(n)$是自适应滤波器的输出。应用梯度下降等迭代方法来最小化$e(n)$,在自适应滤波器收敛后,$\hat{h}(n)$就是你所需要的。

在这里,我找到了一个带有一些示例的 LMS 算法的 python 包。

我们刚刚发布了 SPOQ 方法代码,它使用潜在的时变内核进行 1D 反卷积。它最初是为 . 第一个版本仅在 Matlab 中(应该检查它是否与 Octave 一起工作,一个免费的克隆),我们预计很快会在 Python 中发布。所以,无论如何我都在提议。

代码中有一种模式可以模拟与已知宽度(脚本Load_SPOQ_Data_Simulated.m)的“有限高斯核”卷积的尖峰。我将展示这些原则。我们观察到一维信号,由不同位置和幅度的高斯总和组成(上图)。一个问题是“我们如何才能恢复带有噪声和接近高斯的“位置和幅度”?

SPOQ 高斯核稀疏信号恢复

这个问题的灵感来自分析化学中的典型问题,尤其是色谱中的问题。尖峰信号与仪器响应或线形卷积,产生顶部图像。那么,我们可以恢复原始尖峰的位置和高度吗?

分析化学中稀疏正峰的恢复

这个问题往往是不恰当的。然而,如果没有其他假设,例如稀疏性,或者可以分离几个峰值的概念,就会出现一些问题。参考文献:

使用scipy.signal.deconvolve(signal, divisor)函数:参数signal是您记录的人口增长数据,您认为“与已知宽度的高斯卷积”;参数divisor是高斯滤波器的脉冲响应。

这是使用此功能的示例。原始数据(阵列original)是在无限制生长阶段微流体液滴中微生物种群生长的光密度数据时间序列。让这些原始形式的数据无法访问,因为您声称在您的实验中无法访问它。我通过将这些隐藏数据与高斯滤波器的脉冲响应进行卷积来模拟您的测量过程——例如,使用参数$\sigma = 1$(高斯滤波器的“已知宽度”):$g(t) = \exp \left(-t^2/(2 \sigma ^2) \right)/\left(\sigma \sqrt(2\pi)\right)$我在最大值之前和之后在$\Delta t = 3$处切断了脉冲响应尾部(数组impulse_response),总共留下 7 个值impulse_response大批。现在,您收到的用于分析的数据位于一个数组recorded中。

在这个简单的场景中,可以使用该deconvolve功能恢复原始数据。在示例代码中,恢复的数据以recovered数组的形式返回。数组可remainder用于估计恢复数据中的误差。在这个例子中,剩余值在区间内非常小size of the original data minus size of the impulse_response(应该如此),而在区间外则大得多,我在第 23 个样本处截取了恢复数据的图。

from scipy import signal
import matplotlib.pyplot as plt

original = [0.011,0.012,0.012,0.015,0.026,0.028,0.033,0.023,0.024,0.056,0.026,0.027,0.061,0.045,0.039,0.073,0.090,0.111,0.047,0.094,0.117,0.144,0.150,0.163,0.080,0.215,0.166]
plt.plot(original)
impulse_response = [0.011,0.135,0.607,1,0.607,0.135,0.011]
recorded = signal.convolve(impulse_response,original)
plt.plot(recorded)
recovered, remainder = signal.deconvolve(recorded, impulse_response)
plt.plot(recovered[0:23])
plt.show()

输出:

图1

原始数据颜色为蓝色,过滤数据颜色为橙色,恢复数据颜色为绿色。原始数据图与恢复的数据图部分重叠。

在此示例中,deconvolve可以应用函数,因为高斯滤波器的频率响应对于所有频率值都是非零的。对于更复杂的场景,可能需要使用 Wiener 反卷积。请参阅https://gist.github.com/danstowell/f2d81a897df9e23cc1da中的 Python 实现示例

在您的场景中,您可以使用 function 恢复未过滤的数据deconvolve()请记住,有效恢复数据样本的子阵列比原始阵列短脉冲响应阵列的长度(在我的示例中为七个样本)。