Below I have plotted the signal (Lifetime decay) I am trying to deconvolve from a known impulse response function (IRF), as well as the IRF itself. I'm using scipy.signal.deconvolve.

enter image description here

Please note for my deconvolution code, I am using only the peak of the IRF, not the entire array sequence as shown on the plot above.

Here is the code I am trying to use:

IRF = IRF * (max(decay)/max(IRF))
# replace zeros to avoid error message 
IRF = np.where(IRF == 0, 0.1, IRF)    
decay = np.where(decay == 0, 0.1, decay)  
# take only the quotient part of the result 
deconv = scipy.signal.deconvolve(decay, IRF)[0]
# "padding" the deconvolved signal so it has the same size as the original signal 
s = int((len(decay)-(len(deconv)))/2)  ## difference on each side 
deconv_res = np.zeros(len(decay))   
end = int(len(decay)-s-1)  # final index
deconv_res[s:end] = deconv
deconv = deconv_res
# convolved normalized to decay height for plotting 
deconv_n = deconv * (max(decay)/max(deconv))

The IRF is an array of float64, the signal is an array of uint16.

I admit I'm not familiar with the maths of deconvolution in detail (I'm just writing a software for microscopy image analysis), so I am blindly changing different aspects of the code, like trying different divider functions, normalization, number of points of the IRF, etc, but nothing is producing anywhere near as expected. The last result I got looks like this (see plot of the original signal and what the signal it tried to deconvolve..)

enter image description here

Could anyone just point me in the right direction as to what could be the issue or what I can read to solve my problem? I see scipy uses inverse filter for the deconvolution, is it not the right method for time domain data? Could you recommend any other way to do this impulse response deconvolution in Python? (I couldn't find any other package for this) ..

Thanks in advance!

  • $\begingroup$ Note that you can use np.pad to do the padding more simply numpy.org/doc/stable/reference/generated/numpy.pad.html $\endgroup$
    – endolith
    Oct 3, 2021 at 13:53
  • $\begingroup$ I am unclear as to what you're trying to do. You have an original unknown signal that has passed through the IRF, and this convolution produces the decay signal? What do you mean by "I am using only the peak of the IRF, not the entire array sequence"? Both inputs should be sequences. Have you tried converting the uint16 array to float? Unsigned int arrays can produce strange results if they wrap around, etc. Convert it to float for testing and only remove the conversion for efficiency after you've gotten it working as expected. $\endgroup$
    – endolith
    Oct 3, 2021 at 13:56
  • $\begingroup$ Thanks, sorry I wasn't clear enough. The known "decay" here is the final signal, result of the convolution of the known IRF and an unknown inherent signal. I am trying to deconvolve the known IRF from this known decay to get this unknown inherent signal (the "true" decay that will be fit later). Many thanks! I'll try the conversion now $\endgroup$
    – ISquared
    Oct 4, 2021 at 7:28


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.