0

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?

ISquared
  • 364
  • 4
  • 22

1 Answers1

1

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))
Bob
  • 13,867
  • 1
  • 5
  • 27
  • 1
    Still have to handle a division by zero issue, though. For example, take an impulse convolved with a sinc. In the frequency domain, that's a (periodic) complex exponential multiplied by a rectangular window. Deconvolving by a sinc to get the impulse would be equivalent to dividing by the rectangular window, which has zeros. It makes sense why this is an error, there's no unique answer to reversing a windowing of a periodic function. – BatWannaBe Oct 04 '21 at 07:51
  • 1
    That's right for the general case. For this case, since the filter has a short response with sharp transitions there should be no zeros in the frequency domain :) let's see what the author says. – Bob Oct 04 '21 at 08:03
  • It worked very well (given the noisy data) by transforming to frequency domain and back, many thanks.I've selected your answer to be the solution. If you don't mind, could you explain to me why it doesn't as this here example with the gaussian filter, is it because the shape of the signal is not sharp? https://stackoverflow.com/questions/40615034/understanding-scipy-deconvolve – ISquared Oct 04 '21 at 21:41
  • Thank you. The explanation of the difference has more todo with noise, if the leading coefficient is small and there is noise in the filter or the data the noise excite unstable modes. https://dsp.stackexchange.com/a/78325/55891. The deconvolve does is not optimizing the least square error, it decomposing your signal in a filtered part and a constant part. – Bob Oct 05 '21 at 10:35