# Deconvolution of 1D Signals Blurred by a Gaussian Kernel

I have convolved a random signal with a a Gaussian and added noise (Poisson noise in this case) to generate a noisy signal. Now I would like to deconvolve this noisy signal to extract the original signal using the same Gaussian.

The problem is that I need a code which does the job of deconvolution in 1D. (I have already found some in 2D but my main aim is 1D).

Can you please suggest me some packages or programs that are able to do so? (Preferably in MATLAB)

Thanks in advance for the help.

• use the function deconv in MATLAB. – GOEKHAN GUEL Jul 28 '12 at 12:26
• doesn't work with added noise... – user1724 Jul 30 '12 at 16:19
• You can't deconvolve a signal. You can estimate an inverse convolution given two signals: system's impulse response and system output. Which one are you trying to do? – Phonon Aug 1 '12 at 13:57
• @Phonon: Pretty late with this comment, but there are blind deconvolution methods that don't require knowledge of the system impulse response. As you might imagine, you can do better if you do know the impulse response, though. – Jason R Oct 2 '12 at 12:37
• @JasonR Fair point. – Phonon Oct 2 '12 at 13:54

I've explained it once on StackOverflow.

Your signal can be represented as a vector, and convolution is multiplication with an N-diagonal matrix (where N is the length of the filter). I am assuming for the sake of the answer that the filter is much smaller than the signal

For example:

V1
V2
...
Vn


  [b1 b2 b3];


So the matrix is nxn: (Let it be called A):

[b2 b3 0  0  0  0.... 0]
[b1 b2 b3 0  0  0.... 0]
[0  b1 b2 b3 0  0.... 0]
.....
[0  0  0  0  0  0...b2 b3]


Convolution is:

A*v;


And de-convolution is

A^(-1) * ( A) * v;


Obviously, in some cases de-convolution is not possible. These are the cases when you have singular A. Even matrixes that are not singular, but close to being singular, can be problematic, as they will have large numeric error. You can estimate it by computing the condition number of the matrix.

If A has low condition, you can compute the inverse, and apply it on the result.

Now, let's see some examples in Matlab:

First, I've made a function that computes the convolution matrix.

function A = GetConvolutionMatrix(b,numA)
A = zeros(numA,numA);
vec = [b  zeros(1,numA-numel(b))];
for i=1:size(A,1)
A(i,:) = circshift(vec,[1 i]);
end
end


Now, let's try to see what happens with different kernels:

    b = [1 1 1];
A = GetConvolutionMatrix(b,10);
disp(cond(A));


The condition number is :

 7.8541


This one is problematic, as expected. After averaging, it is hard to get back the original signal.

Now, let's try some milder averaging:

b = [0.1 0.8 0.1];
A = GetConvolutionMatrix(b,10);
disp(cond(A));


The result is:

1.6667


That goes well with our intuition, mild averaging of original signal is much easier to reverse.

We can also see how the inverse matrix looks like:

 figure;imagesc(inv(A));


Here is one line from the matrix:

  0.0003   -0.0026    0.0208   -0.1640    1.2910   -0.1640    0.0208   -0.0026    0.0003   -0.0001


We can see that most of the energy in each line is concentrated in 3-5 coefficients around the center. Therefore, in order to deconvolve, we can simply convolve the signal again with this approximation:

   [0.0208   -0.1640    1.2910   -0.1640    0.0208]


This kernel looks interesting! It is a sharpening operator. Our intuition is correct, sharpening cancels out blur.

• This answer deserves more upvotes – dynamic Jan 25 '13 at 22:01
• Why do you think the matrix is tridiagonal? For Circular Convolution it will be circulant. In most cases it will be Toeplitz. Have a look on my solution. – Royi Nov 25 '19 at 15:23
• Read the answer - I am analyzing a case where the filter has 3 elements. In most cases in image processing, the filter is much smaller than the signal . So yes, it is a Toepliz matrix, but it is also N-diagonal, where N is the length of the filter. Circular convolution is also quite useless in image processing. – Andrey Rubshtein Nov 26 '19 at 9:37
• I've updated the answer to avoid any further confusion. – Andrey Rubshtein Nov 26 '19 at 9:39
• Have you seen Gaussian Kernel which is implemented in 3 samples? – Royi Nov 26 '19 at 11:53

If you have added random noise you cannot get the original signal... You can try to separate the signals in the frequency domain (if the noise and the signal are of different frequencies). But it seems that what you are searching for is a Wiener filter.

I think this is still an open problem.

There are numerous research papers that try to recover the original signal the best they can.

One classic approach is through Wavelet-based Methods.

There are also dictionary approaches like this one.

You can get a more in-depth view of the problem by following the research done by David L. Donho, Michael Elad, Alfred M. Bruckstein etc.

• A recent paper using complex Morlet wavelet by Nguyen, Farge & Schneider seems to yield good results. Google this bibliographic code: 2012PhyD..241..186N A friend of mine used this method with 2D wavelets on the interstellar medium with excellent results. I have yet to look into it in details. – PhilMacKay Oct 2 '12 at 14:11

If I understood the problem properly, we can formalize the problem as follows:

We have a signal model,

$y = Hx + \eta$

where $y$ is the observation, $H$ is the convolution operator, and $\eta$ is the noise. We want to estimate $x$ by using observation and the knowledge of characteristics of noise.

In this case, $\eta$ is simulated from a Poisson distribution. However, the above mentioned dictionary approaches underlie a Gaussian noise assumption. In this case, Gaussian is the convolution operator, not the noise.

I did not work on signal recovery under the Poisson noise, but I googled and found this paper may be useful. Similar approaches in that context can be useful for this problem.

Deconvolution of a noisy data is known to be an ill-posed problem, since the noise is arbitrarily magnified in the reconstructed signal. Therefore, a regularization method is required to stabilize the solution. Here, you can find a MATLAB package that addresses this issue by implementing the Tikhonov's regularization algorithm:

I will go to very beginning of the question. There are deconvolution functions in MATLAB which are used for image processing applications. However, you can also use these functions for 1D signals. For example,

% a random signal
sig_clean = zeros(1,200);
sig_clean(80:100)=100;

figure
subplot(1,3,1)
plot(sig_clean,'b-.','LineWidth',2)
legend('Clean Signal')

% convolve it with a gaussian
x=1:30;
h = exp(-(x-15).^2/20); h=h/sum(h);
sig_noisy = conv(sig_clean,h,'same');

sig_noisy = awgn(sig_noisy,0,'measured');

subplot(1,3,2)
plot(sig_noisy,'r')
hold on, plot(sig_clean,'b-.','LineWidth',3)
legend('Blurred and noise added signal','Clean Signal')


(sig_noisy = sig_clean * h + noise) Then why not deconvolve the output signal with the h function and obtain the (almost) input signal. I am using Wiener deconvolution here

sig_deconvolved=deconvwnr(sig_noisy,h,1);

subplot(1,3,3)
plot(sig_noisy,'r')
hold on, plot(sig_clean,'b-.','LineWidth',2)
hold on, plot(sig_deconvolved,'k--','LineWidth',2)
legend('Blurred and noise added signal','Clean Signal','Deconvolved Signal')


Alternatively, if you don't know the h function, but know the input and output, this time why not deconvolve the input signal with the output which will give the h^-1 function. Then you can use it as a filter to filter the noisy signal. (sig_clean = sig_noisy * h^-1)

h_inv=deconvwnr(sig_clean,sig_noisy,1);

figure;
subplot(1,2,1)
plot(h_inv)
legend('h^-^1')

sig_filtered=conv(sig_noisy,h_inv,'same');
subplot(1,2,2)
plot(sig_noisy,'r')
hold on, plot(sig_clean,'b-.','LineWidth',2)
hold on, plot(sig_filtered,'k--','LineWidth',2)
legend('Blurred and noise added signal','Clean Signal','Filtered Signal')


I hope it helps.

Convolution is the multiplication and summation of two signals one on to the other. I am talking about two deterministic signals. If you want to deconvolve one from the other then this corresponds to the solution of system of equations. As you might know system of equations are not always solvable. The system of equations can be overdetermined, underdetermined or exactly solvable.

In case you add some noise, then you loose some information and you can not get this information back. What you can do is again to solve the linear system of equations considering the fact that each coefficient is added a noise term. Or as you can see in another answer to your question, you might want to first estimate the original signal from the noisy signal and then try to solve the system of equations.

It is important to note that the noise is added to the multiplied and summed up coefficients. Therefore it might even be the case that your system of equation is eventually not uniquly solvable. To be sure that it is uniquely solvable your coefficient matrix should be square and of full rank.

This would be difficult to do. Convolution with a Gaussian is equivalent to multiplication with a Fourier Transform of the Gaussian in the frequency domain. This happens to be also a Gaussian of in essence this is a low pass filter and a really effective one at that. Once you add noise all the information that's in the "stop band" of the Gaussian is destroyed. There is no way to recover that.

De-convolution is essentially multiplying with the inverse of the frequency response. Here is the problem: The inverse of the frequency response gets really, really big where the original Gaussian is very small. At these frequencies you would basically amplify the noise by huge amounts. Even if everything would be completely noise free, you'd most likely run into numerical problems.

### Approaches

There are many methods for Deconvolution (Namely the degradation operator is linear and Time / Space Invariant) out there.
All of them try to deal with the fact the problem is Ill Poised in many cases.

Better methods are those which add some regularization to the model of the data to be restored.
It can be statistical models (Priors) or any knowledge.
For images, a good model is piece wise smooth or sparsity of the gradients.

But for the sake of the answer a simple parametric approach will be takes - -Minimizing the Least Squares Error between the restored data in the model to the measurements.

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

### Numerical Solution

Ensemble of Gaussian Blur Kernel was created.

The parameters are $n = 300$, $k = 31$ and $m = 270$.
The data is random and no noise were added.

In MATLAB the Linear System was solved using pinv() which uses SVD based Pseudo Inverse and the \ operator.

As one can see, using the SVD the solution is much less sensitive as expected.

Why is there an error?
Looking at a solution (For the highest STD):

As one could see the signal is restored very well except for the start and the end.
This is due to the use of Valid Convolution which tells us little on those samples.

### Noise

If we added noise, things would look differently!
The reason results were good before is due to the fact MATLAB could handle the DR of the data and solve the equations even though they had large condition number.

But large condition number means the inverse filter amplify strongly (To reverse the strong attenuation) some frequencies.
When those contain noise it means the noise will be amplified and the restoration will be bad.

As one could see above, now the reconstruction won't work.

### Summary

If one knows the Degradation Operator exactly and the SNR is very good, simple deconvolution methods will work.
The main issue of deconvolution is how hard the Degradation Operator attenuates frequencies.
The more it attenuates the more SNR is needed in order to restore (This is basically the idea behind Wiener Filter).
Frequencies which were set to zero can not be restored!

In practice, in order to have stable results on should add some priors.

The code is available at my StackExchange Signal Processing Q2969 GitHub Repository.

In general, one method to handle the issue that generalizes substantially to a problem of extracting two or more components is to take the spectra G¹, G² ⋯, Gⁿ of signals #1, #2, ..., #n, tabulate the total square Γ(ν) = |G¹(ν)|² + |G²(ν)|² + ⋯ + |Gⁿ(ν)|² at each frequency ν, and normalize G₁(ν) ≡ G¹(ν)* / Γ(ν), G₂(ν) ≡ G²(ν)* / Γ(ν), ..., G_n(ν) ≡ Gⁿ(ν)* / Γ(ν). The problem with ill-definedness and noise corresponds to the fact that Γ(ν) ~ 0 is possible for some frequencies ν. To handle this, add in another "signal" to extract G⁰(ν) = constant - the "noise" signal. Now Γ(ν) will be strictly bounded below. This is almost certainly connected with Tikhonov regularization, but I never found or established any equivalency result or other correspondence. It is simpler and more direct and intuitive.

Alternatively, you may treat the G's as vectors equipped with a suitable inner product e.g. «G,G'» ≡ ∫ G(ν)* G'(ν) dν, and take (G₀, G₁, ⋯, G_n) as a dual (e.g. the generalized inverse) of (G⁰, G¹, ⋯, Gⁿ) - assuming, of course, the component vectors are linearly independent.

For Gaussian deconvolution, one would set up n = 1, G⁰ = the "noise" signal and G¹ = the "Gaussian" signal.

The answer provided by Andrey Rubshtein will fail miserably in the presence of noise, as the described problem is very sensitive to noise and modeling errors. It is a good idea to construct a convolution matrix, but the use of regularization in the inversion is an absolute must in a problem like this. A very simple and straight forward regularization method (although computationally expensive) is the Truncated Singular Value Decomposition (TSVD). Methods like Tikhonov Regularization and Total Variation Regularization are worth checking out. Tikhonov regularization (and its general form) has a very elegant stacked form that is easy to implement in Matlab. Check out the book: Linear and nonlinear inverse problems with practical applications by Samuli Siltanen and Jennifer Mueller.

Actually, the question is not clear. But the answers carified what you've asked for. You can build a system of linear algebraic equations as some people advice, that is correct, but the matrix built on known signal is so-called poorly conditioned. That means when you try to invert it, the truncation errors kill solution and you receive random numbers in result. The common approach is constrained extremum. You minimize norm of solution ||x|| with constrain ||Ax - y|| < delta. So you are looking for such x with smallest norm that not allow difference between Ax and y to be large. It is very simple you need to add so-called parameter of regularization on the principal diagonal of matrix obtained in application of least squares. It is called Tikhonov regularization. I have coding samples that do that, but they are designed for different problems, so I can't provide them, but take my word it is really an elementary problem.