I got a problem with scipy.signal.deconvolve() I tried the answer from Understanding scipy deconvolve but it didn't work. here's my code:
n = 3
s = 240
p = np.linspace(-s*n, s*n, s*n*2)
g = mlab.normpdf(p, 0, s)
f1 = g/max(g)
s = 1920
p = np.linspace(-s*n, s*n, s*n*2)
g = mlab.normpdf(p, 0, s)
g1 = g/max(g)
let the signal be box-like
signal = g1
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 everywhere
gauss = f1
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, _ = deconvolve( filtered, gauss )
the deconvolution has n = len(signal) - len(gauss) + 1 points
n = 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] = deconv
deconv = 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")
ax[1].plot(gauss, color="#68934e", label="gauss filter")
we need to divide by the sum of the filter window to get the convolution normalized to 1
ax[2].plot(filtered/np.sum(gauss), color="#325cab", label="convoluted")
ax[3].plot(deconv, color="#ab4232", label="deconvoluted")
# ax[3].set_yscale('log')
for i in range(len(ax)):
ax[i].set_xlim([0, len(signal)])
ax[i].legend(loc=1, fontsize=11)
if i != len(ax)-1 :
ax[i].set_xticklabels([])
plt.show()
The deconvoluted curve was suppose to be the same as the original gaussian.