Deconvolution of 1D signals

Question

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.

Answer 1

I've explained it once on StackOverflow, but I decided to expand the answer a little bit.

---

Your signal can be represented as a vector, and convolution is multiplication with a tridiagonal matrix.

For example:

Your vector/signal is:

V1
V2
...
Vn

Your filter (convolving element) is:

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

Answer 2

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.

Answer 3

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.

Answer 4

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.

Answer 5

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.

Answer 6

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.