# Deconvolution with Python in real life

I have measured a signal which is convolved with the profile of the measuring apparatus. Now I want to remove this contribution to get the "real" signal. I am trying to do this with Python. For this I am using the scipy.signal.deconvolve function. To understand how to use it I started with this question in Stack Overflow, in particular the answer from Cleb. So I wrote a toy script to play with.

# Toy example

import numpy as np
import scipy.signal
import plotly.graph_objects as go
from plotly.subplots import make_subplots

x = np.arange(0., 20.01, 0.01)
original = np.zeros(len(x))
original[900:1100] = 1
filter_sigma = .2
filter_from_measuring_apparatus = np.exp(-(x-x.mean())**2/2/filter_sigma**2)
filter_from_measuring_apparatus = filter_from_measuring_apparatus[filter_from_measuring_apparatus>.1]
filter_for_deconv = filter_from_measuring_apparatus

measured = scipy.signal.convolve(original, filter_from_measuring_apparatus, mode='full') / filter_from_measuring_apparatus.sum()
deconvoluted, remainder = scipy.signal.deconvolve(measured, filter_for_deconv)
deconvoluted *= filter_for_deconv.sum()

signals = [original, filter_from_measuring_apparatus, measured, filter_for_deconv, deconvoluted]
labels = ['original', 'filter_from_measuring_apparatus', 'measured', 'filter_for_deconv', 'deconvoluted']
fig = make_subplots(rows = len(signals), shared_xaxes = True, vertical_spacing = 0.02)
for fig_row,(signal, label) in enumerate(zip(signals, labels)):
go.Scatter(
x = [i for i in range(len(signal))],
y = signal,
name = label,
),
row = fig_row+1,
col = 1,
)

fig.show()


The output from this toy example looks like this:

which is what we all expect, so fine. I am cutting the Gaussian filter_from_measuring_apparatus with >0.1 because in the answer from ImportanceOfBeingErnest he says that "the filter should be such that it's much bigger then zero everywhere", which is understandable as at some point it will be dividing.

# Real life measuring apparatus but still no noise

So now I want to move to real life. In real life the filter by the measuring apparatus can be anything, in particular a Gaussian that is not "much bigger than zero everywhere". So I modify only the filter_from_measuring_apparatus:

import numpy as np
import scipy.signal
import plotly.graph_objects as go
from plotly.subplots import make_subplots

x = np.arange(0., 20.01, 0.01)
original = np.zeros(len(x))
original[900:1100] = 1
filter_sigma = .2
filter_from_measuring_apparatus = np.exp(-(x-x.mean())**2/2/filter_sigma**2)
filter_for_deconv = filter_from_measuring_apparatus[filter_from_measuring_apparatus>.1] # The filter I will use for deconv is still "much bigger than zero".

measured = scipy.signal.convolve(original, filter_from_measuring_apparatus, mode='full') / filter_from_measuring_apparatus.sum()
deconvoluted, remainder = scipy.signal.deconvolve(measured, filter_for_deconv)
deconvoluted *= filter_for_deconv.sum()
deconvoluted[(deconvoluted<-2)|(deconvoluted>2)] = float('NaN') # Remove so we can see in the plot.

signals = [original, filter_from_measuring_apparatus, measured, filter_for_deconv, deconvoluted]
labels = ['original', 'filter_from_measuring_apparatus', 'measured', 'filter_for_deconv', 'deconvoluted']
fig = make_subplots(rows = len(signals), shared_xaxes = True, vertical_spacing = 0.02)
for fig_row,(signal, label) in enumerate(zip(signals, labels)):
go.Scatter(
x = [i for i in range(len(signal))],
y = signal,
name = label,
),
row = fig_row+1,
col = 1,
)

fig.show()


and now the output looks like this:

Why is this happening? I would expect an output very similar to the previous case as the measured signal is almost the same (only longer).

# Adding noise, keeping toy filter for measuring apparatus

Now let me go back to the toy example with a fake apparatus filter, but add a (very small) noise to the measurement:

import numpy as np
import scipy.signal
import plotly.graph_objects as go
from plotly.subplots import make_subplots

x = np.arange(0., 20.01, 0.01)
original = np.zeros(len(x))
original[900:1100] = 1
filter_sigma = .2
filter_from_measuring_apparatus = np.exp(-(x-x.mean())**2/2/filter_sigma**2)
filter_from_measuring_apparatus = filter_from_measuring_apparatus[filter_from_measuring_apparatus>.1]
filter_for_deconv = filter_from_measuring_apparatus

measured = scipy.signal.convolve(original, filter_from_measuring_apparatus, mode='full') / filter_from_measuring_apparatus.sum()
measured += 0.001*(np.random.randn(len(measured))**2)**.5 # Add some noise to the measurement, as I am measuring power the noise is only positive.
deconvoluted, remainder = scipy.signal.deconvolve(measured, filter_for_deconv)
deconvoluted *= filter_for_deconv.sum()
deconvoluted[(deconvoluted<-2)|(deconvoluted>2)] = float('NaN') # Remove so we can see in the plot.

signals = [original, filter_from_measuring_apparatus, measured, filter_for_deconv, deconvoluted]
labels = ['original', 'filter_from_measuring_apparatus', 'measured', 'filter_for_deconv', 'deconvoluted']
fig = make_subplots(rows = len(signals), shared_xaxes = True, vertical_spacing = 0.02)
for fig_row,(signal, label) in enumerate(zip(signals, labels)):
go.Scatter(
x = [i for i in range(len(signal))],
y = signal,
name = label,
),
row = fig_row+1,
col = 1,
)

fig.show()


which produces as output:

WTF?!?!?! The noise is 1/1000 as strong as the signal, almost invisible in the measured signal plot, how is it possible that produces such a drastic crap after the deconvolution?

Of course combining a real-life-filter for the measuring apparatus and noise simultaneously only worsens the results.

So my question is, how is the deconvolve function supposed to be used in real life cases? Or how should I do what I want to do in Python? I would appreciate an answer working out my toy example.

• Usually with real data iterative convolution with an assumed/known function to describe the process is used for the reasons you have discovered. This iterative approach is widely used, for example, in time correlated single photon counting. Aug 17 '21 at 10:06

Deconvolving in the time domain is equivalent to dividing in the frequency domain, i.e. we try to recover the input as

$$X(\omega) = \frac{Y(\omega)}{H(\omega)}$$

If there are frequencies where $$|H(\omega)|$$ is very small, this becomes an ill defined problem since you are approaching "zero divided by zero". That's what you are seeing in your second example. At high frequency the transfer function of your filter is so small that even the numerical noise of double precision arithmetic throws of the result.

Once you add noise, it becomes worse. Let's assume $$Y(\omega) = H(\omega)X(\omega)+N(\omega)$$. We get

$$X(\omega) = \frac{Y(\omega)}{H(\omega)} = \frac{X(\omega)+N(\omega)}{H(\omega)} = X(\omega) + \frac{N(\omega)}{H(\omega)}$$

The deconvolution amplifies the noise by the inverse of the filter transfer function. If your transfer function is down by 60 dB, your noise will be amplified by 60 dB, which is exactly what you are seeing. A Gaussian is a very effective low pass filter !

Deconvolution in real life

.. is tricky. If you transfer function is reasonably flat over a wide enough frequency range than time domain deconvolution can work. But it always requires careful inspection of the signal to noise ratio (SNR) as a function of frequency and might need "customization" to work around frequency bands with bad SNR.

The underlying problem is fundamental. If the filter attenuates bands of the input signal so they are below the noise floor, then the information is lost and not recoverable.

• Thanks for your answer. In your notation $H$ is the FT of my filter_for_deconv, $Y$ is the FT of my measured signal and $X$ is the FT of my deconvolved signal. The problems of "close to 0 division" should arise when $H$ (or filter_for_deconv) is close to 0. If $H$ is the same for the 3 of my examples, why in the first one the problem does not appear? Aug 12 '21 at 13:41
• I think I now see the reason: $$Y_{\text{estimated}}=X\frac{H_{\text{instrument}}}{H_{\text{model}}}+\frac{N}{H_{\text{model}}}$$ if $H_{\text{instrument}}\equiv H_{\text{model}}$ they end up canceling each other... Now I am intrigued, is it too complicated to deconvolve a Gaussian from such a simple and with really really high SNR signal? This has to be a very common problem, and solved decades ago. How would you remove the contribution of the Gaussian from the measured signal in my example? I refuse to believe that it is not possible... Aug 12 '21 at 13:48
• Convolve a Gaussian with a high frequency sine wave and see what you get: after the initial transient the output is essentially zero. There is nothing left to deconvolve. Aug 12 '21 at 19:31