Search code examples
pythonscipysignal-processingdeconvolution

Python deconvolution giving unexpected result

Below I have plotted the signal (Lifetime decay) I am trying to deconvolve from an impulse response function, i.e. the divider (IRF). So I should just get the decay a bit sharper. Here is an example of a topic I look at that gives what I need:

Understanding scipy deconvolve

enter image description here

Please not for my code, I am using only the peak of the divider (IRF), not the entire array sequence as shown on the image.

I am using the following code to do that:

IRF = IRF * (max(decay)/max(IRF))
# replace 0s 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 so familiar with the maths of deconvolution, so I am trying blindly different things, like trying different divider functions, 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 give me some idea if it's something in scipy.deconvolve I don't understand, what the reason could be for this strange behaviour, or even some high-level reading material that might help me out? Or if you think this problem is more physics-y than coding-related, a better place to ask such a question?

Solution

  • The problem is that deconvolve is a sort of polynomial division, it decomposes the output signal in $conv(h, x) + r$, if your signal is noisy it may give strange results. Also if the first sample in the inpulse response is small it tends to produce the exponentially growing output.

    What I would do for this problem is the division of FFTs.

    N = 2**(ceil(log2(len(IRF) + len(decay)))
filtered = ifft(fft(decay, N) / fft(IRF, N))