# Deconvolution of 1D signals

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 alreday 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 a tridiagonal matrix.

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

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.

-

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.

-

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:

-

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.

-

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.

-