Best way to vectorize operation having input and o

Question

My goal is to vectorize the following operation in numpy,

y[n] = c1*x[n] + c2*x[n-1] + c3*y[n-1]

If n is time, I essentially need the outputs depending on previous inputs as well as previous outputs. I'm given the values of x[-1] and y[-1]. Also, this is a generalized version of my actual problem where c1 = 1.001, c2 = -1 and c3 = 1.

I could figure out the procedure to add the first two operands, simply by adding c1*x and c2*np.concatenate([x[-1], x[0:-1]), but I can't seem to figure out the best way to deal with y[n-1].

Answer 1

One may use an IIR filter to do this. scipy.signal.lfilter is the correct choice in this case.

For my specific constants, the following code snippet would do -

  from scipy import signal
  inital = signal.lfiltic([1.001,-1], [1, -1], [y_0], [x_0])
  output, _ = signal.lfilter([1.001,-1], [1, -1], input, zi=inital)

Here, signal.lfiltic is used to to specify the initial conditions.

Answer 2

Just by playing around with cumsum:

First a little function to produce your expression iteratively:

def foo1(x,C):
    x=x.copy()
    for i in range(1,x.shape[0]-1):
        x[i]=np.dot(x[i-1:i+2],C)
    return x[1:-1]

Make a small test array (I first worked with np.arange(10))

In [227]: y=np.arange(1,11); np.random.shuffle(y)
# array([ 4,  9,  7,  8,  2,  6,  1,  5, 10,  3])

In [229]: foo1(y,[1,2,1])
Out[229]: array([ 29,  51,  69,  79,  92,  99, 119, 142])
In [230]: y[0] + np.cumsum(2*y[1:-1] + 1*y[2:])
Out[230]: array([ 29,  51,  69,  79,  92,  99, 119, 142], dtype=int32)

and with a different C:

In [231]: foo1(y,[1,3,2])
Out[231]: array([ 45,  82, 110, 128, 148, 161, 196, 232])
In [232]: y[0]+np.cumsum(3*y[1:-1]+2*y[2:])
Out[232]: array([ 45,  82, 110, 128, 148, 161, 196, 232], dtype=int32)

I first tried:

In [238]: x=np.arange(10)
In [239]: foo1(x,[1,2,1])
Out[239]: array([  4,  11,  21,  34,  50,  69,  91, 116])
In [240]: np.cumsum(x[:-2]+2*x[1:-1]+x[2:])
Out[240]: array([  4,  12,  24,  40,  60,  84, 112, 144], dtype=int32)

and then realized that the x[:-2] term wasn't needed:

In [241]: np.cumsum(2*x[1:-1]+x[2:])
Out[241]: array([  4,  11,  21,  34,  50,  69,  91, 116], dtype=int32)

If I was back in school I probably would have discovered this sort of pattern with algebra, rather than a numpy trial-n-error. It may not be general enough, but hopefully it's a start.