[SciPy-user] Solving "difference" equations
Travis Oliphant
oliphant at ee.byu.edu
Sun Aug 21 15:02:43 CDT 2005
Robert Kern wrote:
>Kumar Appaiah wrote:
>
>
>>Dear Scipy users,
>>I would like to know whether SciPy has a facility allowing one to
>>solve discrete-time difference equations. For example, equations like
>>the following:
>>
>>y[n] = y[n - 1] + x[n], n is an integer.
>>
>>I tried using odeint, but it solves it as a continuous time equation
>>(or so I think...). Setting hmax to 1.0 also doesn't help.
>>
>>What is the right way?
>>
>>
>
>Difference equations are rather trivial to do yourself via a for-loop. I
>don't think there are any special algorithms that need to be applied.
>
>
>
Robert is right, but there are some filtering tools that may be useful
because they do the loop in C.
The tools are in scipy.signal called lfilter (linear filter).
info(scipy.signal.lfilter)
lfilter(b, a, x, axis=-1, zi=None)
Filter data along one-dimension with an IIR or FIR filter.
Description
Filter a data sequence, x, using a digital filter. This works for many
fundamental data types (including Object type). The filter is a direct
form II transposed implementation of the standard difference equation
(see "Algorithm").
Inputs:
b -- The numerator coefficient vector in a 1-D sequence.
a -- The denominator coefficient vector in a 1-D sequence. If a[0]
is not 1, then both a and b are normalized by a[0].
x -- An N-dimensional input array.
axis -- The axis of the input data array along which to apply the
linear filter. The filter is applied to each subarray along
this axis (*Default* = -1)
zi -- Initial conditions for the filter delays. It is a vector
(or array of vectors for an N-dimensional input) of length
max(len(a),len(b)). If zi=None or is not given then initial
rest is assumed. SEE signal.lfiltic for more information.
Outputs: (y, {zf})
y -- The output of the digital filter.
zf -- If zi is None, this is not returned, otherwise, zf holds the
final filter delay values.
Algorithm:
The filter function is implemented as a direct II transposed structure.
This means that the filter implements
y[n] = b[0]*x[n] + b[1]*x[n-1] + ... + b[nb]*x[n-nb]
- a[1]*y[n-1] + ... + a[na]*y[n-na]
using the following difference equations:
y[m] = b[0]*x[m] + z[0,m-1]
z[0,m] = b[1]*x[m] + z[1,m-1] - a[1]*y[m]
...
z[n-3,m] = b[n-2]*x[m] + z[n-2,m-1] - a[n-2]*y[m]
z[n-2,m] = b[n-1]*x[m] - a[n-1]*y[m]
where m is the output sample number and n=max(len(a),len(b)) is the
model order.
The rational transfer function describing this filter in the
z-transform domain is
-1 -nb
b[0] + b[1]z + ... + b[nb] z
Y(z) = ---------------------------------- X(z)
-1 -na
a[0] + a[1]z + ... + a[na] z
Also, to change the initial conditions you may find
signal.lfiltic helpful
info(scipy.signal.lfiltic)
lfiltic(b, a, y, x=None)
Given a linear filter (b,a) and initial conditions on the output y
and the input x, return the inital conditions on the state vector zi
which is used by lfilter to generate the output given the input.
If M=len(b)-1 and N=len(a)-1. Then, the initial conditions are given
in the vectors x and y as
x = {x[-1],x[-2],...,x[-M]}
y = {y[-1],y[-2],...,y[-N]}
If x is not given, its inital conditions are assumed zero.
If either vector is too short, then zeros are added
to achieve the proper length.
The output vector zi contains
zi = {z_0[-1], z_1[-1], ..., z_K-1[-1]} where K=max(M,N).
-Travis O.
More information about the SciPy-user
mailing list