I am trying to implement a deconvolution-based event detection algorithm in python, but scipy.signal.deconvolve doesn't seem to work in my case. Here is a basic example:

import numpy as np
from scipy.signal import convolve, deconvolve, unit_impulse

def event(x, start, a, t1, t2):
    """Event model based on two-exponent differential.

    Positional arguments:
    x - time array
    start - event start time
    a - event amplitude
    t1 - rise tau
    t2 - decay tau
    y = -np.exp(-(x - start) / t1) + np.exp(-(x - start) / t2)
    n = a / y.max()
    res = n * np.array([0 if t < start else i for t, i in zip(x, y)])

signal = unit_impulse(1000, idx=100)
ev = event(np.arange(300), 0, -10, 1, 5)
conv = convolve(signal, ev)
signal_recovered = deconvolve(conv, ev)[0]

Here, signal_recovered always ends up containing all NaNs, except for the first value, which is Inf. Why does it happen?

  • 1
    $\begingroup$ Have you tried looking at your variables after every step to identify where the NaNs are introduced? $\endgroup$ – MBaz May 24 '18 at 14:53
  • $\begingroup$ As you can see from the minimal example, NaNs don't appear before the deconvolution step. $\endgroup$ – Axon May 24 '18 at 16:04
  • $\begingroup$ Well, it's not obvious just from looking at your code. Anyway, now you may consider reading the docs for deconvolve to figure out if its arguments are valid, asking the author, or reading its source code. $\endgroup$ – MBaz May 24 '18 at 19:08
  • $\begingroup$ If for other reasons you have to deconvolve you have to, but deconvolution is always tricky, usually because of some division by zero. For a function such as you have and also data obtained in an experiment, iterative re-convolution (with a known function) with non-linear least squares works well. $\endgroup$ – porphyrin May 29 '18 at 7:53

Browse other questions tagged or ask your own question.