4
$\begingroup$

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()

Results of the code seen here

$\endgroup$
  • $\begingroup$ Are you saying that scipy.signal.deconvolve does not return a long enough signal? The example in the docs suggests it should be the right length... $\endgroup$ – Peter K. Jan 9 '19 at 20:48
  • $\begingroup$ Yes, I am. In the linked example the convolved signal has length 9, and deconvolved has length 8. In general, if the length of the convolved signal is $M$, and the length of the filter is $N$, it will produce a deconvolved signal of length $M-N+1$, which is the length of the full overlap between convolved signal and kernel. I want the deconvolved signal to have length $M$, not by adding extra zeros at the beginning or end, but by deconvolving the tail correctly $\endgroup$ – Aleksejs Fomins Jan 10 '19 at 8:39
  • $\begingroup$ Would it be possible to talk a little bit more about the problem you are faced with and then we can have a look at deconvolution as well. Are there multiple interferences in the signal from which you want to derive the excitation? Can the signal be considered zero from a point onwards? Where are these signals coming from? $\endgroup$ – A_A Jan 10 '19 at 9:44
  • $\begingroup$ I'll add a minimal example to the question shortly. $\endgroup$ – Aleksejs Fomins Jan 10 '19 at 10:03
  • $\begingroup$ Thank-you! That makes the question clearer. $\endgroup$ – Peter K. Jan 10 '19 at 12:49

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.