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
1
+50
$\begingroup$

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.

enter image description here


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

enter image description here


| improve this answer | |
$\endgroup$
  • $\begingroup$ Thanks for your work :). So you mentioned that the calculation of the inverse kernel has certain assumptions. Can you comment on what assumptions these are. In particular, can the stability of this procedure be guaranteed (or at least estimated) if the kernel is guaranteed to be finite and smooth? $\endgroup$ – Aleksejs Fomins Jun 3 at 12:59
  • $\begingroup$ The main assumption is that the Fourier Transform of the kernel exists and that it is non-zero everywhere in the frequency domain, which I think is the case for most smooth causal functions. Essentially, we want to be able to get back the original signal, so in the filtering process, no frequency information must be lost. This will not be true if the FT of the kernel is zero for certain frequencies (eg: an ideal LPF). $\endgroup$ – 2vrk1504 Jun 4 at 6:21
  • $\begingroup$ I will also like to add that longer your kernel is (more number of samples), more accurately you compute the inverse FT and hence, a lesser error between the original signal and the deconv output is achieved. $\endgroup$ – 2vrk1504 Jun 4 at 6:23
  • $\begingroup$ Ok, this really helps. Thanks again $\endgroup$ – Aleksejs Fomins Jun 4 at 12:53

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.