Given an impulse response h and output y (both one-dimensional arrays), I'm trying to find a way to compute the inverse filter x such that h * x = y, where * denotes the convolution product.

For example, suppose that the impulse response h is [1,0.5], and the output is a step function (that is, consists of all 1s). One can figure out that the first coefficients should be [1, 0.5, 0.75], which yields an output of [1, 1, 1, 0.375]. 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.

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

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:

import numpy as np
from scipy import signal
from sympy import Symbol, series
import timeh = np.array([1,0.5])           # Impulse response function
y = np.array([1,1,1,1,1])       # Ideal output is a step functionmy_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.")

This produces the output:

[ 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.

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?

This is called deconvolution: scipy.signal.deconvolve will do this for you. Example where you know the original input signal x:

import numpy as np
import scipy.signal as signal# We want to try and recover x from y, where y is made like this:
x = np.array([0.5, 2.2, -1.8, -0.1])
h = np.array([1.0, 0.5])
y = signal.convolve(x, h)# This is called deconvolution:
xRecovered, xRemainder = signal.deconvolve(y, h)
# >>> xRecovered
# array([ 0.5,  2.2, -1.8, -0.1])
# >>> xRemainder
# array([  0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
#          0.00000000e+00,  -1.38777878e-17])# Check the answer
assert(np.allclose(xRecovered, x))
assert(np.allclose(xRemainder, 0))

It sounds like you don’t know the original signal, so your xRemainder will not be 0 to machine precision—it’ll represent noise in the recorded signal y.