与一系列矩形信号卷积的一维时域波信号的反卷积

信息处理 matlab 反卷积 逆问题 全变
2022-01-04 06:26:29

我有一个合成信号(下图底部),它是输入信号(顶部)和目标函数(中间)的卷积。目的是在输入信号已知时从卷积信号中检索目标函数。从实际的角度来看,这似乎是一个病态问题,但我很想知道专家对一个人能走多远以及使用什么 SP 工具集的意见。具体来说,我想知道:1)。比较第一个波包和第二个波包,如何判断目标函数中的第一个方格是时间跨度较短但幅度相同的波包,还是时间跨度相同但幅度较小的波包;2)。对于第 4 个方格,如何判断它不是两个单独的短方格,而不是一个长方格。一世' 我还附上了我用来生成数字的 Matlab 代码。非常感谢。

在此处输入图像描述

clear; close all;

%% 10MHz incident signal;
fs = 1000e6;
f = 10e6;
wavelength = fs/f;
sig = wavemaker(3.5, f, fs);
figure; subplot(3,1,1);
plot(sig, 'LineWidth', 1); title(strcat('input signal, wavelength =',num2str(wavelength),' data pts'));
axis([0,3300,-2,2])

%% Sparse signal;
N = 3000;                    % N : length of signal
s = zeros(N,1);
k = [50:(50+wavelength*0.1) 500:(500+wavelength*0.6) 1200:(1200+wavelength*1.6) 2200:(2200+wavelength*4.6)];
s(k) = 1;
subplot(3,1,2);
plot(s, 'LineWidth', 1); title('distribution: objective function');
axis([0,3300,0,2])


%% convoluted signal
y = conv(sig,s);
subplot(3,1,3);plot(y, 'LineWidth', 1);hold on;
plot(abs(hilbert(y)), 'LineWidth', 1);
title('convoluted signal between input and districution');
xlim([0,3300])


function x = wavemaker(nCycles, fc, fs)
% function to generate wave packet;
nSample = round(fs / fc * nCycles);
ts      = 1 / fs;
T       = ts * nSample;
t_max   = ts * (nSample-1);
t       = 0: ts: t_max;
 
x       = sin( 2 * pi * fc .* t);
  
x = x.*hanning(nSample)';
end
3个回答

即使在模拟环境中解决反卷积也不容易,更不用说在实践中了。

解决它的主要技巧是使用正确的模型/先验来解决问题和非常好的测量(高 SNR)。

所以基本上,对于反卷积,我们追求的是:

x^=argminx12Hxy22

其中是已知信号的卷积算子的矩阵形式,是我们的测量样本。Hy

现在,我们需要谈谈卷积形式。卷积必须处理可能对获得良好结果至关重要的边界条件。
我们有几种方法(基本上是外推法):

  1. 用零推断 - 假设信号样本之外的数据为零。
  2. 通过最近邻推断(也称为复制)- 未知值由最近的已知值推断。
  3. 通过周期性延续推断 - 假设数据是周期性的。因此,任何缺失值都基于此。

的建筑必须与之相匹配。由于您使用了,在没有显式模式的代码中,这意味着您基本上选择了零,并且由于您的卷积输出是(默认),这意味着我们在输出中也看到了瞬态。Hconv()full

上述问题的解决方案如下:

x^=argminx12Hxy22=(HTH)1HTy

该解的稳定性取决于条件数备注- 也可以在频域中解决这个问题。尽管在频域中需要一些接触,但该模型是周期性模型。HTH

让我们看看结果:

在此处输入图像描述

首先我们可以看到条件数是巨大的!您可能会将条件数视为错误的放大。这意味着即使是最轻微的噪音也会使事情变得不稳定。可以看出,确实即使是带有标准偏差的白噪声也会1e-8导致错误!

在实践中,为了处理这种不稳定性,我们使用了一些正则化。

x^=argminx12Hxy22+λR(x)

其中是正则化函数,是正则化因子,它在聆听朴素反卷积模型或正则化模型之间取得平衡。R()λ

必须根据我们对感兴趣信号的先验知识来选择正则化函数。
在您的情况下,您的信号清楚的是它的分段平滑属性。它的梯度非常稀疏。因此,它与 Total Variation 模型完美匹配:

x^=argminx12Hxy22+λTV(x)=argminx12Hxy22+λDx1

其中有限差分算子(离散导数)。 这是一个相对容易解决的任务。在我的项目Total Variation (TV) Regularized Least Squares - Solvers Analysis我实施并比较了几个求解器。 对于这个例子,我使用了基于 ADMM 的求解器D

这是 TV 正则化的结果:

在此处输入图像描述

可以看出,它完全克服了上面的(非常低!)噪音。
在现实世界(和更高的噪声水平)中,需要调整参数。你很难恢复完美的结果,但它们会比采用幼稚的方法要好得多。λ

MATLAB 实现

完整代码可在我的StackExchange Signal Processing Q71822 GitHub 存储库中找到(查看SignalProcessing\Q71822文件夹)。
它包括从样本构建卷积矩阵和解决 TV 问题的功能。

其实我也有和你类似的问题。但在我的情况下,目标函数不像矩形脉冲,而只是如下所示的尖峰。我在超声波测试领域工作。所以,这个例子有点像超声波信号。

在此处输入图像描述

但是,首先我尝试使用来自Royi的代码,使用来自答案的 ADMM 求解器对我来说,我发现他的解决方案很难找到目标函数,有时至少对于真实数据会有随机峰值。由于我有大量信号要反卷积,所以我使用了 L1-MM 算法。Royi曾经在这里解释了 L1-MM 与 ADMM 求解器之间的时间,其中 L1-MM 收敛得更快。

对于我们知道输入和输出信号的问题,此处提供的 Royi 的 ADMM求解器代码有效。但老实说,我很难理解代码以便能够与理论联系起来。因此,我没有机会为我的问题修改他的 L1-MM 代码(在此处提供)。

我使用了Ivan Selesnick在此处提供的 L1-MM 代码。在那里,作者通过示例和 L1-MM 求解器方程的推导给出了解释。该代码具有以方程式命名的变量,因此很容易理解。使用 L1-MM 处理超声波信号的结果如下:

在此处输入图像描述

这里的尖峰没有像真正的目标函数那样被识别出来,但是这些尖峰的位置是确定的。至少对我来说,这就是我现在想要实现的。不过,您可以看到来自恢复的目标函数的卷积信号与原始信号相匹配。

如果您的数据是尖峰而不是矩形,这是 L1-MM 算法的结果。(ps,使用矩形目标函数,从 L1-MM 恢复的目标函数看起来不正确。但你可以试一试。)

在此处输入图像描述 在此处输入图像描述 在此处输入图像描述

所有代码和数据都可以在这里找到:在我的 github

对于我的超声波导波数据:

我要反卷积的测量数据: 在此处输入图像描述

使用 Royi 的 ADMM 求解器进行的反卷积 在此处输入图像描述

使用 Ivan Selesnick 的 L1-MM 求解器执行的反卷积: 在此处输入图像描述

即使 SNR 为 10dB 的 L1-MM 求解器也显示出良好的结果:因此,它可用于对我的数据进行去噪: 在此处输入图像描述 在此处输入图像描述 在此处输入图像描述

如果您在所有频率上都具有良好的信噪比,这相对简单。由于您已经以双精度对此进行了模拟,因此可以通过以下方式轻松完成

objectiveFunction  = real(ifft(fft(y,8192)./fft(sig',8192))); 
sDiff = s-objectiveFunction(1:length(s));
fprintf('Recovery error= %6.2fdB\n',10*log10(mean(sDiff.^2)./mean(s.^2)));

实际上,这并不能很好地工作,因为任何现实世界的测量在某些频率上都会有不良的信噪比,并且您需要所有频率来正确重建时域信号。

在这种情况下,您可能必须进行带通重建或尝试填充、内插或修复坏频点。

如果这还不够好,您可以将其设置为时域中的最小二乘误差问题:卷积可以表示为矩阵乘法,您只需要反转这个方程。

我猜小波输入信号不会特别好用,因为它们的频带非常窄。您的输入信号必须在与您要恢复的信号相关的所有频率上都具有足够的能量。能量不足频率的信息会在输出信号中丢失。