Deconvolution to Remove Gaussian Blur in 1D Signal (Wiener Filtering?)

Question

I've got a set of biology data that I'm trying to denoise (effectively, a population statistic can only be measured convolved with a gaussian of known width)

My problem is this: I can measure (f*g), and I can measure g(x), but I need to know f(x). Because this is real data with slight blips, I can't just use a deconvolve() function.

I've been told by an EE friend that I should be using 1-D Wiener Deconvolution. The problem is, I have no earthly idea how to implement this. Supposedly, there's a Matlab package for it, but I don't have access to Matlab; and all the Python examples I've been able to find (e.g. scikit-learn) are 2D only.

Is there something I'm missing? Is there an easy formula I can implement?

Please, don't just point me to the Wikipedia page/a textbook for 1D Wiener filtering. I've tried to read several of these, but they look like incomprehensible gibberish (I'm a biologist).

EDIT: Also what's up with point spread functions? I assume my PSF should be the known Gaussian g(x), and yet every example I can find of Gaussian deconvolution with Wiener filtering (all in 2d) use a uniform distribution as their PSF instead of a Gaussian distribution.

Thanks!

Answer 1

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 (Look at the SignalProcessing\Q2969 folder).

Resources

• Deconvolution of Synthetic 1D Signals - How To?
• Estimate the Filter Coefficients of 1D Filtration (Convolution).
• 1D Deconvolution with Gaussian Kernel (MATLAB).
• How Is the Formula for the Wiener Deconvolution Derived?

Answer 2

It's a system identification problem. As we know convolution is commutative: $f(n)*g(n)=g(n)*f(n)$. Now you have an unknown signal $f(n)$ filtered by a known system $g(n)$, and we say the output is $d(n)$. On the other hand, we can also say $d(n)$ is the output of a known signal $g(n)$, which is gaussian noise in your case, passed by an unknown system $f(n)$. You can identify the unknown system by Wiener filter or adaptive filter like least-mean-squares (LMS). The system block is as follows.

In your case we assume there is no interference and $y(n)=d(n)$. The error signal $e(n)=\hat{y}(n)-d(n)$, where $\hat{y}(n)$ is the output of the adaptive filter. Iterative methods such as gradient decsent are applied to minimize $e(n)$, and after the adaptive filter converges, $\hat{h}(n)$ is what you need.

Here I find a python package of LMS algorithm with some examples.

Answer 3

We have just released the SPOQ method code that does 1D deconvolution with a potentially time-varying kernel. It was initially meant for . The first release is only in Matlab (one should ckeck whether it works with Octave, a free clone), and we expect a release in Python soon. So, I am proposing it anyway.

There is one mode in the code to simulate spikes convolved with "finite Gaussian kernels" with known width(script Load_SPOQ_Data_Simulated.m). I will display the principles. We observe a 1D signal, made of a sum of Gaussians of different locations and amplitudes (top plot). A question is "How can we recover "locations and amplitudes" with noise and close Gaussians?

This problem is inspired by typical problems in analytical chemistry, especially in chromatography. A spike signal is convolved with an instrumental response or a lineshape, yielding the top picture. Then, can we recover the location and height of the original spikes?

The problem is often ill-posed. Yet, with t other assumptions, like sparsity, or the concept that a few peaks can be separated, there are some issues. Refernces:

• SPOQ: smooth, sparse ℓp-over-ℓq ratio regularization Matlab toolbox,
• SPOQ ℓp-Over-ℓq Regularization for Sparse Signal Recovery applied to Mass Spectrometry,
• SPOQ webpage.

Answer 4

Use a scipy.signal.deconvolve(signal, divisor) function: parameter signal is your recorded population growth data that you believe is "convolved with a gaussian of known width"; parameter divisor is an impulse response of your Gaussian filter.

Here is an example of using this function. A raw data (array original) is an optical density data time series of microbial population growth in microfluidic drops during unrestrained stage of growth. Let this data, in its raw form, is not accessible, as you claim it is inaccessible in your experiment. I simulate your measurement procedure by convolving this hidden data with an impulse response of the Gaussian filter -- for example, with a parameter $\sigma = 1$ (the "known width" of your Gaussian filter): $g(t) = \exp \left(-t^2/(2 \sigma ^2) \right)/\left(\sigma \sqrt(2\pi)\right)$. I cut off the impulse response tails (array impulse_response) at $\Delta t = 3$ before and after the maximum value, leaving in total 7 values in the impulse_response array. Now, the data you receive for analysis is in an array recorded.

In this simple scenario, the original data can be restored with the deconvolve function. In the example code, the restored data is returned in a recovered array. A remainder array can be used to estimate the error in the recovered data. In this example, the remainder values are very small in the interval size of the original data minus size of the impulse_response (as it should be) and much greater outside this interval, and I cut off a plot of the recovered data at the 23rd sample.

from scipy import signal
import matplotlib.pyplot as plt

original = [0.011,0.012,0.012,0.015,0.026,0.028,0.033,0.023,0.024,0.056,0.026,0.027,0.061,0.045,0.039,0.073,0.090,0.111,0.047,0.094,0.117,0.144,0.150,0.163,0.080,0.215,0.166]
plt.plot(original)
impulse_response = [0.011,0.135,0.607,1,0.607,0.135,0.011]
recorded = signal.convolve(impulse_response,original)
plt.plot(recorded)
recovered, remainder = signal.deconvolve(recorded, impulse_response)
plt.plot(recovered[0:23])
plt.show()

Output:

The original data color is blue, the filtered data color is orange, and the recovered data color is green. The original data graph is partially overlapped by the recovered data graph.

In this example, function deconvolve can be applied because the frequency response of the Gaussian filter is nonzero for all frequency values. For more complicated scenarios, it may become necessary to use the Wiener deconvolution. See an example of Python's implementation in https://gist.github.com/danstowell/f2d81a897df9e23cc1da.

In your scenario, you can recover unfiltered data using function deconvolve(). Remember that the subarray of the valid recovered data samples is shorter than the original array by the length of the impulse response array (seven samples in my example).