方法
有许多反卷积方法(即退化算子是线性的和时空不变的)。
在许多情况下,他们都试图处理问题是 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
文件夹)。
资源