[SciPy-User] how to use signal.lfiltic?
Skipper Seabold
jsseabold@gmail....
Wed Nov 24 10:09:35 CST 2010
On Wed, Nov 24, 2010 at 10:39 AM, Skipper Seabold <jsseabold@gmail.com> wrote:
> This is mainly for my own understanding of going back and forth
> between signal processing language and time series econometrics.
>
> I don't see how to use lfiltic.
>
> Say I have a (known) output vector of errors and an input y vector. y
> follows a mean zero ARMA(p,q) process given by b and a below where p =
> q = 2. If I want to use lfilter to recreate the errors, forcing the
> first p outputs (errors) to be zero, then I need to solve the
> difference equations for the zi that do this which is given by zi
> below. But if I try to recreate this using lfiltic, it doesn't work.
> Am I missing the intention of lfiltic? Matlab's documentation, which
> a lot of the signal stuff seems to be taken from, suggests that y and
> x in lfiltic need to be reversed, but this also doesn't give the zi I
> want.
>
> from scipy import signal
> import numpy as np
>
> errors = np.array([ 0., 0., 0.00903417, 0.89064639, 1.51665674])
> y = np. array([-0.60177354, -1.60410646, -1.16619292, 0.44003132, 2.36214611])
>
> b = np.array([ 1. , -0.8622494 , 0.34549996])
> a = np.array([ 1. , 0.07918344, -0.81594865])
>
> # zi I want to produce errors = 0,0,...
>
> zi = np.zeros(2)
> zi[0] = -b[0] * y[0]
> zi[1] = -b[1] * y[0] - b[0] * y[1]
>
> zi
> # array([ 0.60177354, 1.08522758])
>
> e = signal.lfilter(b, a, y, zi=zi)
> e[0]
> # array([ 0. , 0. , 0.00903417, 0.89064639, 1.51665674])
>
> zi2 = signal.lfiltic(b,a, errors[:2], y[:2])
>
> zi2
> # array([-0.03533984, -0.20791273])
>
> e2 = signal.lfilter(b, a, y, zi=zi2)
>
> e2[0]
> # array([-0.63711338, -1.24269149, -0.41241704, -0.0899541 , 1.25042151])
>
> Skipper
>
Basically, I think this line in signal.lfiltic
for m in range(M):
zi[m] = sum(b[m+1:]*x[:M-m],axis=0)
should be
for m in range(M):
zi[m] = sum(-b[:m+1][::-1]*x[:m+1],axis=0)
I'm not sure about the next loop, since my output are zero, I didn't
have to solve for it.
Is this a bug or am I misunderstanding?
Skipper
More information about the SciPy-User
mailing list