Deconvolution to Remove Gaussian Blur in 1D Signal (Wiener Filtering?) 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!
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.
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?