1

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

CasualPythoner
  • 69
  • 1
  • 9

1 Answers1

1

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'])
pygri
  • 647
  • 6
  • 17
  • Your code works perfectly for this case. However, the actual data I will receive and will have to work with won't really work with mode='full'. Any ideas why mode ='same' is an issue? I am really confused by this. – CasualPythoner Oct 16 '20 at 16:51
  • is not that it is an issue per se, is just that it is not supported by scipy deconvolve function - if this is such a big problem for you you should roll out your own convolution/deconvolution function (is not that hard, look up how to do that, I think most implementation are Fast Fourier Transformation (FFT) based) – pygri Oct 17 '20 at 17:08