9

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.

enter image description here

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;
Royi
  • 19,608
  • 4
  • 197
  • 238
Anand Mohan
  • 117
  • 9

3 Answers3

7

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.

Hilmar
  • 44,604
  • 1
  • 32
  • 63
  • 4
    In the case above, synthetic case, the SNR is approaching infinity. Since the Gaussian LPF is not zero the Deconvolution should work almost perfectly. To show it though, one needs the signals the OP used. – Royi Feb 04 '19 at 20:02
  • 1
    welcome to the 10K club, Hil. – robert bristow-johnson Feb 07 '19 at 21:16
  • 2
    The Fourier transform of a Gaussian is a Gaussian. Both have support $\mathbb R$. Hence, the high-frequency content is not quite "gone", it's simply attenuated to the point where it is non-recoverable. – Rodrigo de Azevedo Feb 10 '19 at 10:16
  • I think in perfect SNR world your statement is wrong. In reality it depends on the SNR on the frequencies with effective information. Though it is marked as "Accepted" I think it should be edited. – David Apr 18 '20 at 20:23
  • There is no such thing as perfect SNR world (unless you stay in the realm of pencil and paper) – Hilmar Apr 19 '20 at 16:05
  • @Hilmar, I think you're too conclusive with your statements which in this case are wrong. Every time you take a photo with your smartphone the ISPC applies deconvolution. Do you find the results unacceptable? Have you seen results of SOTA networks doing inverse problems? – Royi Apr 16 '22 at 14:32
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.

enter image description here

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:

enter image description here

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

Royi
  • 19,608
  • 4
  • 197
  • 238
  • 1
    Regarding the plots, are you sure about the output of the convolution? The kernel is wide, something like $5$ periods of the input signal. The output of the convolution looks too neat. – Rodrigo de Azevedo Feb 11 '19 at 16:28
  • I checked the code again and it seems valid. Do you see any issue in the code? I think the scaling is confusing. – Royi Feb 11 '19 at 16:44
  • 1
    I am not talking about the code. I am talking about the plots. The kernel is much wider than the input pulses. Is the issue in the scaling of the plots? – Rodrigo de Azevedo Feb 12 '19 at 06:41
  • But the plot is straight from the code. The plot is plot(vX, vC); where vC = conv2(vS, vH, 'same');. – Royi Feb 12 '19 at 12:51
  • 1
    Let's forget about the code for now. The problem is in the convolution. I would claim that a kernel that wide would produce an output in which there is a single, wiggly pulse. – Rodrigo de Azevedo Feb 12 '19 at 14:59
  • 1
    I think what's getting ion your way is the scaling of the axis relative to Sampling Rate. I will fix that later on. – Royi Feb 12 '19 at 15:25
1

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.