5
$\begingroup$

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;
|improve this question
$\endgroup$
  • 1
    $\begingroup$ Is it possible to provide the Gaussian curve you used for the convolution or its parameters (?). $\endgroup$ – A_A Feb 4 '19 at 12:29
  • $\begingroup$ also, considering the result looks exactly like Gibb's phenomenon in its purest textbook shape, the exact code you used to convolve and deconvolve might be relevant. $\endgroup$ – Marcus Müller Feb 4 '19 at 13:28
  • $\begingroup$ oh, and, especially: how you plotted this. Wild guess: you've plotted the real part of the IFFT, and due to ambiguities that's not what you think it should be. $\endgroup$ – Marcus Müller Feb 4 '19 at 13:31
  • $\begingroup$ If you know the input is a square wave, why not quantize the reconstruction to obtain a square wave? $\endgroup$ – Rodrigo de Azevedo Feb 5 '19 at 10:24
  • $\begingroup$ Could you also share the code to generate the other signal? $\endgroup$ – Royi Feb 6 '19 at 19:33
7
$\begingroup$

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.

|improve this answer
$\endgroup$
  • 2
    $\begingroup$ 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. $\endgroup$ – Royi Feb 4 '19 at 20:02
  • 1
    $\begingroup$ welcome to the 10K club, Hil. $\endgroup$ – robert bristow-johnson Feb 7 '19 at 21:16
  • 1
    $\begingroup$ 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. $\endgroup$ – Rodrigo de Azevedo Feb 10 '19 at 10:16
2
$\begingroup$

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.

|improve this answer
$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – Rodrigo de Azevedo Feb 11 '19 at 16:28
  • $\begingroup$ I checked the code again and it seems valid. Do you see any issue in the code? I think the scaling is confusing. $\endgroup$ – Royi Feb 11 '19 at 16:44
  • $\begingroup$ 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? $\endgroup$ – Rodrigo de Azevedo Feb 12 '19 at 6:41
  • $\begingroup$ But the plot is straight from the code. The plot is plot(vX, vC); where vC = conv2(vS, vH, 'same');. $\endgroup$ – Royi Feb 12 '19 at 12:51
  • 1
    $\begingroup$ I think what's getting ion your way is the scaling of the axis relative to Sampling Rate. I will fix that later on. $\endgroup$ – Royi Feb 12 '19 at 15:25
1
$\begingroup$

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.

|improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.