The tail of scipy deconvolve

Question

I have a dataset for the function $g(t) = \int_{0}^{\infty}f(t-\tau)h(\tau)d\tau$ I would like to deconvolve. My assumptions are

• The signal $f(t) = 0$ before the start of observation, that is for $t < 0$
• The filter I use is causal, meaning that $h(t) = 0$ for $t < 0$ (for example, an exponential distribution function)

I have read this post, checked the proposed solution, and it works. My problem is that my filter function has a reasonably long tail, and the sampling frequency is rather coarse. My signal does not manage to recede to zero during the time of observation, as given in all the toy examples. The scipy.signal.deconvolve function seems to only deconvolve the part of the original signal where the convolved signal and the filter function fully overlap.

For my purposes, this is too wasteful.

Is it possible to use scipy.signal.deconvolve or another tool to deconvolve the entire dataset, not just the fully-overlapping part. I understand that, potentially, deconvolving the signal at the very end of the tail has larger error than the one in the middle. An ideal solution would propose a value of the deconvolved function for the entire duration of observation, and also report the expected deconvolution error as function of time. Any advice is welcome

Edit

Here is the minimal code to reconstruct undesired behaviour. The signal deconvolved using scipy.signal.deconvolve is shorter than the recorded convolved signal. I believe that, under suitable assumptions, the convolved signal contains enough information to reconstruct the original signal over the entire duration of recording. Note that

• The signal provided in the example below is an oversimplification of my actual signal. It is a continuously changing signal, which is recorded over a short period of time. It is never really zero
• The real signal I have is not zero before the observation either, but not much action happens there. I believe it can be approximated by a constant baseline, effects of which can be subtracted from the period of interest. This carries an intrinsic error, of course, and I'll try to estimate it at some point, but that is a topic for another question. For now, I would assume that the signal value before the start of observation is zero in order to isolate the main problem of this question

import numpy as np
import matplotlib.pyplot as plt
import scipy.signal

def exp_ker(t, tau):
    return np.exp(-t/tau)/tau

# Problem size
T_SIGNAL = 10
T_KERNEL = T_SIGNAL / 2

DT = 0.01
t_signal = np.arange(0, T_SIGNAL + DT, DT)
t_kernel = np.arange(0, T_KERNEL + DT, DT)
t_conv_full = np.arange(0, T_SIGNAL + T_KERNEL + DT, DT)
t_deconv_scipy = np.arange(0, T_SIGNAL - T_KERNEL + DT, DT)
NPOINT = len(t_signal)

# Box Signal
signal = np.zeros(NPOINT)
signal[np.sin(t_signal - 1.5)**2 > 0.5] = 1

# Kernel
TAU = 1.0
kernel = exp_ker(t_kernel, TAU)
kernel /= np.sum(kernel)

# Convolve
sig_conv_full = scipy.signal.convolve(signal, kernel, mode='full')
sig_conv_real = sig_conv_full[:NPOINT]

# Deconvolve
sig_deconv_full, sig_rem_full = scipy.signal.deconvolve(sig_conv_full, kernel)
sig_deconv_real, sig_rem_real = scipy.signal.deconvolve(sig_conv_real, kernel)

fig, ax = plt.subplots(nrows=2, ncols=3, figsize = (11,7), tight_layout=True)
ax[0][0].plot(t_signal, signal)
ax[1][0].plot(t_kernel, kernel)
ax[0][1].plot(t_conv_full, sig_conv_full)
ax[1][1].plot(t_signal, sig_conv_real)
ax[0][2].plot(t_signal, sig_deconv_full)
ax[1][2].plot(t_deconv_scipy, sig_deconv_real)

ax[0][0].set_title("Original signal")
ax[1][0].set_title("Filter")
ax[0][1].set_title("Full convolution")
ax[1][1].set_title("Realistic measured convolution")
ax[0][2].set_title("Deconvolution from full")
ax[1][2].set_title("Deconvolution from realistic")

ax[0][0].set_xlim([0, T_SIGNAL+T_KERNEL])
ax[1][0].set_xlim([0, T_SIGNAL+T_KERNEL])
ax[0][1].set_xlim([0, T_SIGNAL+T_KERNEL])
ax[1][1].set_xlim([0, T_SIGNAL+T_KERNEL])
ax[0][2].set_xlim([0, T_SIGNAL+T_KERNEL])
ax[1][2].set_xlim([0, T_SIGNAL+T_KERNEL])

plt.show()

Answer 1

I do not think there is a method to get a perfect reconstruction of the original signal by using a trimmed/cropped convolution. However, I am hoping there are approximations/assumptions you can make like the measured signal has finite energy.

If that is the case, what you can do instead is to find the inverse of the kernel, which you can then convolve with the measured signal.

You will, however, have some minor amount of error near the end, and this is unavoidable unless you have an infinitely long signal.

Following is what I did with your example:

inv_kernel = np.fft.ifft(1/np.fft.fft(kernel))

...

sig_deconv_full_2 = scipy.signal.convolve(sig_conv_real, inv_kernel, mode='full')
signal_estimate = sig_deconv_full_2[:NPOINT] #signal trimmed to the required length

---

This is the script run with your original signal. The squared error term is the squared error of signal_estimate and signal.

The impulse is expected as the measured signal suddenly ends. We can just trim the signal to get rid of the impulse.

---

I also tested it with a random signal I generated using NumPy's random library.

---