[SciPy-user] Solving "difference" equations

[Travis Oliphant]

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.