高斯核模糊的一维信号的反卷积

#1 楼

我已经在StackOverflow上对其进行了解释。


您的信号可以表示为矢量，并且卷积与N对角矩阵相乘（其中N是滤波器的长度） 。为了回答这个问题，我假设滤波器比信号小得多。例如：

您的矢量/信号是：

V1
V2
...
Vn

您的过滤器（卷积元素）为： br />

  [b1 b2 b3];

卷积是：

[b2 b3 0  0  0  0.... 0]
[b1 b2 b3 0  0  0.... 0]
[0  b1 b2 b3 0  0.... 0]
.....
[0  0  0  0  0  0...b2 b3]

反卷积是

A*v;

显然，在某些情况下，无法进行反卷积。这是当您使用奇数A时的情况。即使不是奇数但接近奇数的矩阵也可能会出现问题，因为它们会产生较大的数值误差。您可以通过计算矩阵的条件数来估计它。

如果A具有较低的条件，则可以计算逆，并将其应用到结果中。 />
首先，我做了一个计算卷积矩阵的函数。

A^(-1) * ( A) * v;

现在，让我们尝试看看不同的内核会发生什么：

function A = GetConvolutionMatrix(b,numA)
    A = zeros(numA,numA);
    vec = [b  zeros(1,numA-numel(b))];
    for i=1:size(A,1)
        A(i,:) = circshift(vec,[1 i]);
    end
end

条件编号为：

    b = [1 1 1];
    A = GetConvolutionMatrix(b,10);
    disp(cond(A));

如预期的那样，这是有问题的。平均后，很难恢复到原始信号。

现在，让我们尝试更温和的平均：

 7.8541

结果是：

b = [0.1 0.8 0.1];
A = GetConvolutionMatrix(b,10);
disp(cond(A));

这与我们的直觉相吻合，对原始信号进行轻度平均更容易反转。

我们还可以看到反转的方式。矩阵看起来像：

1.6667

这里是矩阵的一行：

 figure;imagesc(inv(A));

我们可以看到，每条线中的大多数能量都集中在中心附近的3-5个系数中。因此，为了进行反卷积，我们可以简单地用这种近似方法再次对信号进行卷积：

它是一个锐化运算符。我们的直觉是正确的，锐化可以消除模糊。

#2 楼

如果添加了随机噪声，则无法获得原始信号...
您可以尝试在频域中分离信号（如果噪声和信号的频率不同）。但似乎您要搜索的是维纳滤波器。

#3 楼

我认为这仍然是一个未解决的问题。

有许多研究论文试图尽其所能恢复原始信号。

一种经典方法是通过基于Wavelet的方法方法。

也有类似的字典方法。 Michael Elad，Alfred M.Bruckstein等。

#4 楼

如果我对问题的理解正确，我们可以将问题形式化如下：

我们有一个信号模型，

$ y = Hx + \ eta $

其中$ y $是观测值，$ H $是卷积运算符，而$ \ eta $是噪声。我们要通过观察和噪声特征知识来估计$ x $。在这种情况下，根据泊松分布模拟$ \ eta $。但是，以上提到的字典方法是高斯噪声假设的基础。在这种情况下，高斯不是卷积运算符，而是噪声。在这种情况下，类似的方法对于解决此问题可能很有用。

#5 楼

众所周知，对噪声数据进行去卷积是一个不适当地的问题，因为在重构信号中噪声会被任意放大。因此，需要一种正规化方法来稳定溶液。在这里，您可以找到一个通过实现Tikhonov的正则化算法来解决此问题的MATLAB软件包：

https://github.com/soheil-soltani/TranKin。

#6 楼

我将开始讨论这个问题。 MATLAB中有用于图像处理应用程序的去卷积功能。
但是，您也可以将这些功能用于一维信号。例如，

% a random signal
sig_clean = zeros(1,200); 
sig_clean(80:100)=100;

figure
subplot(1,3,1)
plot(sig_clean,'b-.','LineWidth',2)
legend('Clean Signal')

% convolve it with a gaussian
x=1:30;
h = exp(-(x-15).^2/20); h=h/sum(h);
sig_noisy = conv(sig_clean,h,'same');

% and add noise
sig_noisy = awgn(sig_noisy,0,'measured');

subplot(1,3,2)
plot(sig_noisy,'r')
hold on, plot(sig_clean,'b-.','LineWidth',3)
legend('Blurred and noise added signal','Clean Signal')

（sig_noisy = sig_clean * h + noise）
然后为什么不使用h函数对输出信号进行反卷积并获得（几乎）输入信号。我在这里使用维纳反卷积

sig_deconvolved=deconvwnr(sig_noisy,h,1);

subplot(1,3,3)
plot(sig_noisy,'r')
hold on, plot(sig_clean,'b-.','LineWidth',2)
hold on, plot(sig_deconvolved,'k--','LineWidth',2)
legend('Blurred and noise added signal','Clean Signal','Deconvolved Signal')

或者，如果您不知道h函数，但是知道输入和输出，这次为什么不要将输入信号与输出解卷积，这将提供h^-1功能。然后，您可以将其用作过滤器来过滤噪声信号。帮助。

#7 楼

卷积是两个信号彼此相乘和相加。我说的是两个确定性信号。如果要彼此解卷积，则这对应于方程组的解。如您所知，方程组并不总是可解的。方程组可以是超定的，不确定的或完全可求解的。考虑到每个系数都加上一个噪声项这一事实，您可以再次解决线性方程组。或者，如您在问题的另一个答案中看到的那样，您可能想要首先从噪声信号中估计原始信号，然后尝试求解方程组。

请务必注意噪声被加到相加和相加的系数上。因此，甚至可能最终您的方程组无法唯一地求解。为确保它是唯一可解的，您的系数矩阵应为正方形且具有完整等级。

#8 楼

这将很难做到。与高斯的卷积等效于在频域中与高斯的傅立叶变换相乘。这恰好也是高斯，实际上这是一个低通滤波器，并且实际上是一个有效的滤波器。
添加噪声后，高斯“阻带”中的所有信息都会被破坏。无法恢复。

去卷积实际上是与频率响应的倒数相乘。这是问题所在：原始高斯很小的地方，频率响应的反函数变得非常大。在这些频率下，您基本上可以将噪声放大很多。即使一切都完全没有噪音，您也很可能会遇到数值问题。

#9 楼

通常，一种解决该问题的方法基本上可以概括为提取两个或多个分量的问题，即采用信号＃1，＃2，...，＃n的频谱G¹，G²，Gⁿ将总和制成表格平方Γ（ν）= |G¹（ν）|²+ |G²（ν）|²+⋯+ |Gⁿ（ν）|²在每个频率ν，并对G₁（ν））G¹（ν）* /Γ进行归一化（ν），G 2（ν）≡G2（ν）* /Γ（ν），...，G_n（ν）≡Gⁿ（ν）* /Γ（ν）。不确定性和噪声的问题与以下事实相对应：对于某些频率ν，Γ（ν）〜0是可能的。要处理此问题，请添加另一个“信号”以提取G⁰（ν）=常数-“噪声”信号。现在Γ（ν）将严格限制在下面。这几乎可以肯定与Tikhonov正则化有关，但是我从未发现或建立任何等效结果或其他对应关系。它更简单，更直接，更直观。

或者，您可以将G视为带有合适内积的向量，例如«G，G'»≡∫G（ν）* G'（ν）dν，并将（G₀，G₁，⋯，G_n）取为（G⁰，G¹，⋯，Gⁿ）的对偶（例如广义逆） -当然，假设分量矢量是线性独立的。

对于高斯反卷积，将设置n = 1，G⁰=“噪声”信号，G¹=“ Gaussian”信号。

#10 楼

方法

这里有许多用于解卷积的方法（即，退化算子是线性的并且时间/空间不变）。
所有这些方法都试图解决以下问题：很多情况下。

更好的方法是在要恢复的数据模型中添加一些正则化的方法。
可以是统计模型（先验）或任何知识。
对于图像，好的模型是逐段平滑的或稀疏的渐变。

但是为了回答这个问题，将采用一种简单的参数方法-最小化模型中恢复的数据与度量之间的最小二乘误差。

模型

最小二乘模型很简单。
目标函数作为数据的函数由下式给出： \ 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} $是卷积矩阵。 br />
解决方案

最小二乘解由下式给出：

$$ \ hat {x} = \ left（{H} ^ {T} H \ right）^ {-1} {H} ^ {T} y $$

可以看出，它需要矩阵求逆。
充分解决此问题的能力取决于遵守$ \ operatorname {cond} \ left（H \ right）= \ sqrt {\ operatorname {cond} \ left的运算符$ {H} ^ {T} H $的条件编号（{H} ^ {T} H \ right）} $。

条件数分析

该条件数的含义是什么？
可以使用线性代数来回答它。
但我认为更直观。方法将在频域中考虑。

降级算符通常会衰减高频能量。
现在，由于在频率上这基本上是元素明智的乘法，因此会说将其反转的简单方法是使用逆滤波器对元素进行除法。那么我们就有了实际的问题...
这基本上是条件编号告诉我们的，某些频率相对于其他频率衰减有多困难。



以上可以看到条件数（使用[dB]单位）是高斯滤波器STD参数的函数。最后是数值问题）。

数值解

高斯模糊核集合创建。



参数为$ n = 300 $，$ k = 31 $和$ m = 270 $。
数据是随机的，没有添加噪声。逆运算符和pinv()运算符。



可以看到，使用SVD，该解决方案不如预期那样敏感。有错误吗？
寻找解决方法（最高STD）：和结尾。
这是由于使用有效卷积对这些样本所知甚少。由于MATLAB可以处理数据的DR并即使方程具有很大的条件数，也可以求解方程，因此之前的结果很好。

但是条件数大意味着逆滤波器在某些频率上会强烈放大（以抵消强衰减）。
当包含噪声时，则意味着噪声会被放大，恢复效果很差。



如上所示，现在的重建将不起作用。

摘要


解卷积的主要问题是衰减算子衰减频率的难度。
它越衰减越多，则需要更多的SNR。进行还原（这基本上是维纳滤波器的原理）。
设置为零的频率无法还原！

在实践中，为了获得稳定的结果，应添加一些先验条件。

我的StackExchange Signal Processing Q2969 GitHub存储库中提供了代码。

#11 楼

实际上，这个问题尚不清楚。但是答案证实了您的要求。您可以按照某些人的建议构建线性代数方程组，这是正确的，但是基于已知信号的矩阵是所谓的条件差的。这意味着当您尝试将其求逆时，截断错误会杀死解决方案，并且您会收到随机数。常见的方法是约束极值。您可以最小化解范|| x ||约束|| Ax-y || <三角洲。因此，您正在寻找具有最小范数且不允许Ax和y之间的差异较大的x。非常简单，您需要在应用最小二乘法获得的矩阵的主对角线上添加所谓的正则化参数。这称为Tikhonov正则化。我有这样做的编码示例，但是它们是针对不同问题而设计的，因此我无法提供它们，但请注意这确实是一个基本问题。

#12 楼

在存在噪声的情况下，Andrey Rubshtein提供的答案将惨败，因为所描述的问题对噪声和建模错误非常敏感。构造卷积矩阵是一个好主意，但是在这样的问题中，绝对必须使用反演中的正则化。截断奇异值分解（TSVD）是一种非常简单直接的正则化方法（尽管计算量很大）。诸如Tikhonov正则化和总变化正则化之类的方法值得一试。 Tikhonov正则化（及其一般形式）具有非常优雅的堆叠形式，可以轻松在Matlab中实现。查阅本书：Samuli Siltanen和Jennifer Mueller撰写的《线性和非线性反问题及其实际应用》。