python recursive vectorization with timeseries

Question

I have a Timeseries (s) which need to be processed recursively to get a timeseries result (res). Here is my sample code:

res=s.copy()*0  
res[1]=k # k is a constant  
for i in range(2,len(s)):  
    res[i]=c1*(s[i]+s[i-1])/2 +c2*res[i-1]+c3*res[i-2]

where c1,c2,c3 are constants. It works properly but I'd like to use vectorization and I tried with:

res[2:]=c1*(s[2:]+s[1:-1])/2+c2*res[1:-1]+c3*res[0:-2]  

but I get "ValueError: operands could not be broadcast together with shapes (1016) (1018) "
if I try with

res=c1*(s[2:]+s[1:-1])/2+c2*res[1:-1]+c3*res[0:-2]  

doesn't give any error, but I don't get a correct result, because res[0] and res[1] have to be initialized before the calculation will take place. Is there a way to process it with vectorization?
Any help will be appreciated, thanks!

Answer 1

This expression

    res[i] = c1*(s[i] + s[i-1])/2 + c2*res[i-1] + c3*res[i-2]

says that res is the output of a linear filter (or ARMA process) with input s. Several libraries have functions for computing this. Here's how you can use the scipy function scipy.signal.lfilter.

From inspection of the recurrence relation, we get the coefficients of the numerator (b) and denominator (a) of the filter's transfer function:

b = c1 * np.array([0.5, 0.5])
a = np.array([1, -c2, -c3])

We'll also need an appropriate initial condition for lfilter to handle res[:2] == [0, k]. For this, we use scipy.signal.lfiltic:

zi = lfiltic(b, a, [k, 0], x=s[1::-1])

In the simplest case, one would call lfilter like this:

y = lfilter(b, a, s)

With an initial condition zi, we use:

y, zo = lfilter(b, a, s, zi=zi)

However, to exactly match the calculation provided in the question, we need the output y to start with [0, k]. So we'll allocate an array y, initialize the first two elements with [0, k], and assign the output of lfilter to y[2:]:

y = np.empty_like(s)
y[:2] = [0, k]
y[2:], zo = lfilter(b, a, s[2:], zi=zi)

Here's a complete script with the original loop and with lfilter:

import numpy as np
from scipy.signal import lfilter, lfiltic


c1 = 0.125
c2 = 0.5
c3 = 0.25

np.random.seed(123)
s = np.random.rand(8)
k = 3.0

# Original version (edited lightly)

res = np.zeros_like(s)
res[1] = k  # k is a constant  
for i in range(2, len(s)):  
    res[i] = c1*(s[i] + s[i-1])/2 + c2*res[i-1] + c3*res[i-2]


# Using scipy.signal.lfilter

# Coefficients of the filter's transfer function.
b = c1 * np.array([0.5, 0.5])
a = np.array([1, -c2, -c3])

# Create the initial condition of the filter such that
#     y[:2] == [0, k]
zi = lfiltic(b, a, [k, 0], x=s[1::-1])

y = np.empty_like(s)
y[:2] = [0, k]
y[2:], zo = lfilter(b, a, s[2:], zi=zi)

np.set_printoptions(precision=5)
print "res:", res
print "y:  ", y

The output is:

res: [ 0.       3.       1.53206  1.56467  1.24477  1.08496  0.94142  0.84605]
y:   [ 0.       3.       1.53206  1.56467  1.24477  1.08496  0.94142  0.84605]

---

lfilter accepts an axis argument, so you can filter an array of signals with a single call. lfiltic does not have an axis argument, so setting up the initial conditions requires a loop. The following script shows an example.

import numpy as np
from scipy.signal import lfilter, lfiltic
import matplotlib.pyplot as plt


# Parameters
c1 = 0.2
c2 = 1.1
c3 = -0.5
k = 1

# Create an array of signals for the demonstration.
np.random.seed(123)
nsamples = 50
nsignals = 4
s = np.random.randn(nsamples, nsignals)

# Coefficients of the filter's transfer function.
b = c1 * np.array([0.5, 0.5])
a = np.array([1, -c2, -c3])

# Create the initial condition of the filter for each signal
# such that
#     y[:2] == [0, k]
# We need a loop here, because lfiltic is not vectorized.
zi = np.empty((2, nsignals))
for i in range(nsignals):
    zi[:, i] = lfiltic(b, a, [k, 0], x=s[1::-1, i])

# Create the filtered signals.
y = np.empty_like(s)
y[:2, :] = np.array([0, k]).reshape(-1, 1)
y[2:, :], zo = lfilter(b, a, s[2:], zi=zi, axis=0)

# Plot the filtered signals.
plt.plot(y, linewidth=2, alpha=0.6)
ptp = y.ptp()
plt.ylim(y.min() - 0.05*ptp, y.max() + 0.05*ptp)
plt.grid(True)
plt.show()

Plot: