title: Inverse filtering using Python
body: <p>Given an impulse response <code>h</code> and output <code>y</code> (both one-dimensional arrays), I'm trying to find a way to compute the inverse filter <code>x</code> such that <code>h * x = y</code>, where <code>*</code> denotes the convolution product.</p>
  
  <p>For example, suppose that the impulse response <code>h</code> is <code>[1,0.5]</code>, and the output is a step function (that is, consists of all <code>1</code>s). One can figure out that the first coefficients should be <code>[1, 0.5, 0.75]</code>, which yields an output of <code>[1, 1, 1, 0.375]</code>. The last term contains an error, but this is not really an issue as I'm only concerned with the output up to a certain maximum time.</p>
  
  <p>I would like to automate and 'scale up' this inverse filtering to a longer, more complex impulse response function. So far, the only way I've figured out to get the coefficients is by generating the series expansion of the Z-transform using sympy. (Note that the Z-transform of a step function is 1/(1-z)).</p>
  
  <p>However, I've noticed that sympy is quite slow in computing the coefficients: it takes 0.8 seconds even for a simple, short example, as demonstrated by the script below:</p>
  
  <pre><code>import numpy as np
  from scipy import signal
  from sympy import Symbol, series
  import time
  
  h = np.array([1,0.5])           # Impulse response function
  y = np.array([1,1,1,1,1])       # Ideal output is a step function
  
  my_x = np.array([1,0.5,0.75])   # Inverse filter response (calculated manually)
  
  my_y = signal.convolve(h,my_x)  # Actual output (should be close to ideal output)
  
  print(my_y)
  
  start = time.time()
  z = Symbol('z')
  print(series(1/((1-z)*(1+z/2)),z))
  end = time.time()
  print("The power series expansion took "+str(end - start)+" seconds to complete.")
  </code></pre>
  
  <p>This produces the output:</p>
  
  <pre><code>[ 1.     1.     1.     0.375]
  1 + z/2 + 3*z**2/4 + 5*z**3/8 + 11*z**4/16 + 21*z**5/32 + O(z**6)
  The power series expansion took 0.798881053925 seconds to complete.
  </code></pre>
  
  <p>In short, the coefficients of the power expansion match the desired filter response, but it seems cumbersome to me to use sympy to do this. Is there a better way to compute inverse filter coefficients in Python?</p>