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;
  • 1
    $\begingroup$ Is it possible to provide the Gaussian curve you used for the convolution or its parameters (?). $\endgroup$
    – A_A
    Commented Feb 4, 2019 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$ Commented Feb 4, 2019 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$ Commented Feb 4, 2019 at 13:31
  • $\begingroup$ If you know the input is a square wave, why not quantize the reconstruction to obtain a square wave? $\endgroup$ Commented Feb 5, 2019 at 10:24
  • $\begingroup$ Could you also share the code to generate the other signal? $\endgroup$
    – Royi
    Commented Feb 6, 2019 at 19:33

3 Answers 3


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.

  • 4
    $\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
    Commented Feb 4, 2019 at 20:02
  • 1
    $\begingroup$ welcome to the 10K club, Hil. $\endgroup$ Commented Feb 7, 2019 at 21:16
  • 2
    $\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$ Commented Feb 10, 2019 at 10:16
  • $\begingroup$ 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. $\endgroup$
    – David
    Commented Apr 18, 2020 at 20:23
  • $\begingroup$ There is no such thing as perfect SNR world (unless you stay in the realm of pencil and paper) $\endgroup$
    – Hilmar
    Commented Apr 19, 2020 at 16:05

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


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.


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

  • 1
    $\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$ Commented Feb 11, 2019 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
    Commented Feb 11, 2019 at 16:44
  • 1
    $\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$ Commented Feb 12, 2019 at 6:41
  • 1
    $\begingroup$ 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. $\endgroup$ Commented Feb 12, 2019 at 14:59
  • 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
    Commented Feb 12, 2019 at 15:25

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.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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