Deconvolution of a 1D Time Domain Wave Signal Convolved with Series of Rect Signals

Question

I have a synthesized signal (the bottom of the following figure), which is the convolution of the input signal (at the top) and the objective function (in the middle). The intention is to retrieve the objective function from the convoluted signal, when the input signal is known. From practical point of view, it seems to be an ill-conditioned problem, but I'm curious to know the expert opinions on how far one can get, and with what SP toolset. Specifically, I would like to know: 1). comparing the first wave packet and second, how to tell whether the first square in the objective function is one with shorter time span but the same amplitude, or with the same time span but smaller-amp; 2). for the 4th square, how to tell it's not two separate short squares rather than a single long one. I've also attached the Matlab codes I used to generate the figures. Thank you very much.

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

Answer 1

Solving a deconvolution isn't easy even in simulated environment not to mention in practice.

The main trick to solve it is using the proper model / prior for the problem and very good measurements (High SNR).

So basically, for deconvolution we're after:

$$ \hat{\boldsymbol{x}} = \arg \min_{\boldsymbol{x}} \frac{1}{2} {\left\| H \boldsymbol{x} - \boldsymbol{y} \right\|}_{2}^{2} $$

Where $ H $ is the matrix form of the convolution operator of the known signal and $ \boldsymbol{y} $ is our measurement samples.

Now, we need to talk about the Convolution form. The convolution has to deal with boundary conditions which might be crucial for good results.
We have few approaches (Basically extrapolation):

1. Extrapolate with Zeros - Assuming data outside the samples of the signals is zero.
2. Extrapolate by Nearest Neighbor (Also called Replicate) - The unknown value is extrapolated by the nearest known value.
3. Extrapolate by Periodic continuation - The data is assumed to be periodic. Hence any missing value is based on that.

The building of $ H $ must match that. Since you used, in your code conv() with no explicit mode it means you basically chose zeros and since your convolution output is full (The default) it means we see the transients as well in the output.

The solution to the above is given by:

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

The stability of this solution depends on the Condition Number of $ {H}^{T} H $.
Remark - One could solve this in Frequency Domain as well. Though it will require some touches as in Frequency Domain the model is the periodic model.

Let's have a look on the results:

First we can see the Condition Number is huge! You may think about the condition number as the amplification of the error. It means even the slightest noise will make things unstable. As can be seen, indeed even a white noise with a standard deviation of 1e-8 caused errors!

In practice, to deal with this instability we use some regularization.

$$ \hat{\boldsymbol{x}} = \arg \min_{\boldsymbol{x}} \frac{1}{2} {\left\| H \boldsymbol{x} - \boldsymbol{y} \right\|}_{2}^{2} + \lambda R \left( x \right) $$

Where $ R \left( \cdot \right) $ is the regularization function and $ \lambda $ is the regularization factor which balances between listening to the naïve deconvolution model or to the regularization model.

The regularization function must be chosen with respect to the prior knowledge we have about the signal of interest.
In your case, something clear about your signal is its piece wise smooth property. Its gradient is very sparse. Hence it is a perfect match to the Total Variation model:

$$ \hat{\boldsymbol{x}} = \arg \min_{\boldsymbol{x}} \frac{1}{2} {\left\| H \boldsymbol{x} - \boldsymbol{y} \right\|}_{2}^{2} + \lambda TV \left( x \right) = \arg \min_{\boldsymbol{x}} \frac{1}{2} {\left\| H \boldsymbol{x} - \boldsymbol{y} \right\|}_{2}^{2} + \lambda {\left\| D x \right\|}_{1} $$

Where $ D $ is the finite differences operator (Discrete Derivative).
This is a relatively easy task to solve. In my project Total Variation (TV) Regularized Least Squares - Solvers Analysis I implemented and compared few solvers.
For this example I used the ADMM based solver.

Here is the result for the TV regularization:

As can be seen, it totally overcome the (Very Low!) noise from above.
In real world (And higher noise levels) one needs to tweak the $ \lambda $ parameter. You will hardly recover perfect results, but they will be much better than doing the naive approach.

MATLAB Implementation

The full code is available on my StackExchange Signal Processing Q71822 GitHub Repository.
It includes the functions to build the convolution matrix from the samples and solve the TV problem.

Answer 2

Actually I have a similar problem like yours. But in mine, the objective function is not like rectangular pulse but just spikes as shown below. I work in ultrasoinc testing field. So, this example is kind of like an ultrasoinc signal.

But, first I tried with the code from Royi using the ADMM solver from the: answer. For me, I find his solution was painfully slow to find the objective function and sometimes there are random spikes atleast for the real data. Since I have large number of signals to deconvolve, I went with L1-MM algorithm. Royi once explained here the time taken between L1-MM vs ADMM solver where L1-MM converges faster.

For our problem where we know the input and output signal, Royi's ADMM solver code provided here worked. But to be honest, I had a hard time to understand the code to be able to relate with the theory. So, I stood no chance in modifying his L1-MM code (provided here) for my problem.

I had used the L1-MM code provided here by Ivan Selesnick. There the author had given explanation with examples and the derivation of the L1-MM solver equations. The code had variables named after the equations so it was easy to follow. The result of using L1-MM for the ultrasoinc signal is given below:

Here the spikes had not been identified like the real objective function but the location of these spikes was determined. For me atleast that is what I would like to achieve for now.Still, you can see that the convolved signal from the recovered objective function matched the original signal.

If your data was spikes instead of rectangles, here is the result from the L1-MM algorithm. (p.s., with rectangle objective function, the recovered objective function from L1-MM does not look correct. But you can give it a go.)

All the codes and data are available here: In my github

For my ultrasoinc guided wave data:

The measured data that I want to deconvolve:

The deconvolution performed with ADMM solver of Royi

The deconvolution performed with L1-MM solver from Ivan Selesnick:

L1-MM solver even with 10dB SNR showed good results: So, it can be used for denoising for my data:

Answer 3

This is relatively straight forward if you have decent signal to noise ratio at all frequencies. Since you have simulated this in double precision this can be easily done with something like

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)));

In practice that doesn't work that well since any real world measurement will have bad signal to noise ratio at some frequency and you need ALL frequencies to properly reconstruct the time domain signal.

In that case you may have to do bandpass reconstruction or try to fill in, interpolate or repair the bad frequency points.

If that's not good enough, you can set it up as a least square error problem in the time domain: Convolution can be expressed as a matrix multiplication and you just need to invert this equation.

I'm guessing that nothing will work particularly well with a wavelet input signal since they are very narrow band. Your input signal must have sufficient energy at all frequencies that are relevant to the signal you want to recover. The information at frequencies with insufficient energy is lost in the output signal.