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

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.