Deconvolution of Synthetic 1D Signals - How To?

Question

I convolved a square wave with a Gaussian wave using linear convolution. Can I get the original square wave back by deconvolving my output with the Gaussian function?

I took the FFT of both signals, divided and then took the IFFT to get back the square wave. The output I got looked like random noise. I thought this might be due to the fact that I am doing a division and some denominator values might be very low, hence causing this error. I tried to set up a threshold and obtained the result depicted below.

How can I improve the result?

Edit 1 : Adding the code used to generate Gaussian signal. Was using matlab.

fs      =  200;           % Sampling frequency

t1 = 1/fs : 1/fs : n1/fs; % Where n1 is the length of the square wave signal

s2   =  4*gaussmf(t1, [ 0.4 4.5 ]);   % Generating Gaussian signal

Edit 2 : Adding code used to generate square wave. I am afraid, it's a clumsy code!

fs      =  200;                           % Sampling Frequency
t       =  1/fs : 1/fs : 3;
n       =  length(t);

sq(1 : round(n/3) )  = eps;                               
sq( (round(n/3) + 1) : 2*round(n/3) ) = 3;
sq( (2*round(n/3) + 1) : n ) = eps;

sq( (3*round(n/3) + 1) : 4*round(n/3) ) = 3;
sq( (4*round(n/3) + 1) : 5*round(n/3) ) = eps;
sq( (5*round(n/3) + 1) : 6*round(n/3) ) = 3;

sq( (6*round(n/3) + 1) : 7*round(n/3) ) = eps;
sq( (7*round(n/3) + 1) : 8*round(n/3) ) = 3;
sq( (8*round(n/3) + 1) : 9*round(n/3) ) = eps;

Answer 1

You cant't recover the original signal through deconvolution.

A Gaussian kernel is in essence a lowpass filter, i.e. it will remove information at higher frequencies from the signal. Once it's gone, it's gone and you can't recover it.

This problem shows up as "divide by zero" or "divide by a very small number", which then amplifies numerically noise of the original convolution to a massive amount.

In order to recover the original signal the system must have "sufficient signal to noise ratio in the bandwidth of interest". What exactly that is, depends heavily on your specific application.

You can de-convolve from a low-shelf or high-shelf filter, but not from a low-pass or anything with a stop band or lots of attenuation at specific frequencies.

You'd also lose some information at the beginning and end of the signal.

Answer 2

Most of the information is given in my answer to 1D Deconvolution with Gaussian Kernel (MATLAB) (Which is related to Deconvolution of 1D Signals Blurred by Gaussian Kernel).

Model

The least squares model is simple.
The objective function as a function of the data is given by:

$$ f \left( x \right) = \frac{1}{2} \left\| h \ast x - y \right\|_{2}^{2} $$

The optimization problem is given by:

$$ \arg \min_{x} f \left( x \right) = \arg \min_{x} \frac{1}{2} \left\| h \ast x - y \right\|_{2}^{2} $$

Where $ x $ is the data to be restored, $ h $ is the Blurring Kernel (Gaussian in this case) and $ y $ is the set of given measurements.
The model assumes the measurements are given only for the valid part of the convolution. Namely if $ x \in \mathbb{R}^{n} $ and $ h \in \mathbb{R}^{k} $ then $ y \in \mathbb{R}^{m} $ where $ m = n - k + 1 $.

This is a linear operation in finite space hence can be written using a Matrix Form:

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

Where $ H \in \mathbb{R}^{m \times n} $ is the convolution matrix.

Solution

The Least Squares solution is given by:

$$ \hat{x} = \left( {H}^{T} H \right)^{-1} {H}^{T} y $$

As can be seen it requires a matrix inversion.
The ability to solve this adequately depends on the condition number of the operator $ {H}^{T} H $ which obeys $ \operatorname{cond} \left( H \right) = \sqrt{\operatorname{cond} \left( {H}^{T} H \right)} $.

Condition Number Analysis

What's behind this condition number?
One could answer it using Linear Algebra.
But a more intuitive, in my opinion, approach would be thinking of it in the Frequency Domain.

Basically the degradation operator attenuates energy of, generally, high frequency.
Now, since in frequency this is basically an element wise multiplication, one would say the easy way to invert it is element wise division by the inverse filter.
Well, it is what's done above.
The problem arises with cases the filter attenuates the energy practically into zero. Then we have real problems...
This is basically what's the Condition Number tells us, how hard some frequencies were attenuated relative to others.

Above one could see the Condition Number (Using [dB] units) as a function of the Gaussian Filter STD parameter.
As expected, the higher the STD the worse the condition number as higher STD means stronger LPF (Values going down at the end are numerical issues).

Since you problem is basically a linear combination of Gaussian Kernels the answer stays the same.

Here is a simple example:

In the example above the reconstruction is perfect.
This is a result of the Convolution Kernel having Low Condition Number.
If you dial up the Gaussian Kernel STD things will deteriorate.

The code is available at my StackExchange Signal Processing Q55284 GitHub Repository (Look at the SignalProcessing\Q55284 folder).

Answer 3

The keyword for deconvolution from signals and filters with known statistics is "Wiener Filter". Basically you want to minimize the expected (square-norm here) error of the reconstructed signal taking into account the available signal-to-noise ratios. Your "noise" here will typically not be floating point errors but rather systematic errors due to using FFT windowing etc.