2

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()    

enter image description here

The deconvoluted curve was suppose to be the same as the original gaussian.

Community
  • 1
  • 1
  • look at [this](http://stackoverflow.com/questions/40615034/understanding-scipy-deconvolve) awww what should I write to have 30 symbols? – Dmitry Pukhov Apr 13 '17 at 21:05

0 Answers0