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):
- Extrapolate with Zeros - Assuming data outside the samples of the signals is zero.
- Extrapolate by Nearest Neighbor (Also called Replicate) - The unknown value is extrapolated by the nearest known value.
- 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 (Look at the SignalProcessing\Q71822 folder).
It includes the functions to build the convolution matrix from the samples and solve the TV problem.
