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!
Answers:
Thank you for visiting the Q&A section on Magenaut. Please note that all the answers may not help you solve the issue immediately. So please treat them as advisements. If you found the post helpful (or not), leave a comment & I’ll get back to you as soon as possible.
Method 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:
All methods was sourced from stackoverflow.com or stackexchange.com, is licensed under cc by-sa 2.5, cc by-sa 3.0 and cc by-sa 4.0