pythonscipysignal-processingconvolutiondeconvolution

A question about deconvolution of a signal using Python scipy


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


Solution

  • 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'])