Understanding scipy deconvolve

After some trial and error I found out how to interprete the results of scipy.signal.deconvolve() and I post my findings as an answer.

Let's start with a working example code

import numpy as npimport scipy.signalimport matplotlib.pyplot as plt# let the signal be box-likesignal = np.repeat([0., 1., 0.], 100)# and use a gaussian filter# the filter should be shorter than the signal# the filter should be such that it's much bigger then zero everywheregauss = np.exp(-( (np.linspace(0,50)-25.)/float(12))**2 )print gauss.min()  # = 0.013 >> 0# calculate the convolution (np.convolve and scipy.signal.convolve identical)# the keywordargument mode="same" ensures that the convolution spans the same#   shape as the input array.#filtered = scipy.signal.convolve(signal, gauss, mode='same') filtered = np.convolve(signal, gauss, mode='same') deconv,  _ = scipy.signal.deconvolve( filtered, gauss )#the deconvolution has n = len(signal) - len(gauss) + 1 pointsn = len(signal)-len(gauss)+1# so we need to expand it by s = (len(signal)-n)/2#on both sides.deconv_res = np.zeros(len(signal))deconv_res[s:len(signal)-s-1] = deconvdeconv = deconv_res# now deconv contains the deconvolution # expanded to the original shape (filled with zeros) #### Plot #### fig , ax = plt.subplots(nrows=4, figsize=(6,7))ax[0].plot(signal,            color="#907700", label="original",     lw=3 ) ax[1].plot(gauss,          color="#68934e", label="gauss filter", lw=3 )# we need to divide by the sum of the filter window to get the convolution normalized to 1ax[2].plot(filtered/np.sum(gauss), color="#325cab", label="convoluted" ,  lw=3 )ax[3].plot(deconv,         color="#ab4232", label="deconvoluted", lw=3 ) for i in range(len(ax)):    ax[i].set_xlim([0, len(signal)])    ax[i].set_ylim([-0.07, 1.2])    ax[i].legend(loc=1, fontsize=11)    if i != len(ax)-1 :        ax[i].set_xticklabels([])plt.savefig(__file__ + ".png")plt.show()    

This code produces the following image, showing exactly what we want (Deconvolve(Convolve(signal,gauss) , gauss) == signal)

Some important findings are:

• The filter should be shorter than the signal
• The filter should be much bigger than zero everywhere (here > 0.013 is good enough)
• Using the keyword argument mode = 'same' to the convolution ensures that it lives on the same array shape as the signal.
• The deconvolution has n = len(signal) - len(gauss) + 1 points. So in order to let it also reside on the same original array shape we need to expand it by s = (len(signal)-n)/2 on both sides.

Of course, further findings, comments and suggestion to this question are still welcome.

As written in the comments, I cannot help with the example you posted originally. As @Stelios has pointed out, the deconvolution might not work out due to numerical issues.

I can, however, reproduce the example you posted in your Edit:

That is the code which is a direct translation from the matlab source code:

import numpy as npimport scipy.signalimport matplotlib.pyplot as pltx = np.arange(0., 20.01, 0.01)y = np.zeros(len(x))y[900:1100] = 1.y += 0.01 * np.random.randn(len(y))c = np.exp(-(np.arange(len(y))) / 30.)yc = scipy.signal.convolve(y, c, mode='full') / c.sum()ydc, remainder = scipy.signal.deconvolve(yc, c)ydc *= c.sum()fig, ax = plt.subplots(nrows=2, ncols=2, figsize=(4, 4))ax[0][0].plot(x, y, label="original y", lw=3)ax[0][1].plot(x, c, label="c", lw=3)ax[1][0].plot(x[0:2000], yc[0:2000], label="yc", lw=3)ax[1][1].plot(x, ydc, label="recovered y", lw=3)plt.show()