A question about deconvolution of a signal using Python scipy

713 views Asked by At

I am trying to learn a bit of signal processing , specifically using Python. Here's a sample code I wrote.

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import deconvolve

a = np.linspace(-1,1,50)
b = np.linspace(-1,1,50)**2
c = np.convolve(a,b,mode='same')
quotient,remainder = deconvolve(c,b);
plt.plot(a/max(a),"g")
plt.plot(b/max(b),"r")
plt.plot(c/max(c),"b")
plt.plot(remainder/max(remainder),"k")
#plt.plot(quotient/max(quotient),"k")

plt.legend(['a_original','b_original','convolution_a_b','deconvolution_a_b'])

In my understanding, the deconvolution of the convoluted array should return exactly the same array 'a' since I am using 'b' as the filter. This is clearly not the case as seen from the plots below.

I am not really sure if my mathematical understanding of deconvolution is wrong or if there is something wrong with the code. Any help is greatly appreciated!

enter image description here

1

There are 1 answers

2
pygri On BEST ANSWER

You are using mode='same', and this seems not to be compatible with scipy deconvolve. Try with mode='full', it should work much better.

Here the corrected code:

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import deconvolve

a = np.linspace(-1,1,50)
b = np.linspace(-1,1,50)**2
c = np.convolve(a,b,mode='full')
quotient,remainder = deconvolve(c,b)
plt.plot(a,"g")
plt.plot(b,"r")
plt.plot(c,"b")
plt.plot(quotient,"k")
plt.xlim(0,50)
plt.ylim(-6,2)

plt.legend(['a_original','b_original','convolution_a_b','deconvolution_c_b'])