[SciPy-dev] reimplementation of lfilter

[Sturla Molden]

Ralf Gommers skrev:
> Sure, I know. The docs I linked have improvements over svn. If you 
> merge your changes, there'll be a conflict and those docs will not get 
> merged soon - I was hoping to avoid that.

Oh, the docs, not the code... Would something like this be better?

The Cython module doesn't need to have docs in it.



from linear_filter import lfilter as _lfilter

def lfilter(b, a, x, axis=-1, zi=None, direction=1, inplace=-1):
    """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. If axis is -1, the filter is applied to x
            flattened. (*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. If axis is None, zi should be a vector
          of length K = max(M,N), where M=len(b)-1 and N=len(a)-1. If
          an axis is not None, zi should either be a vector of length
          K, or an array with same shape as x, except for zi.shape[axis]
          which should be K. If zi is a vector, the same initial
          conditions are applied to each vector along the specified
          axis. Default: zi=None (initial rest).
    direction -- Filtered in the reverse direction if -1. (*Default* = 1)
    inplace -- If possible, modification to x is done implace. 
(*Default* = -1)

  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

    a[0]*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

    """
    if isscalar(a):
        a = [a]

    if axis < 0: axis = None
    if inplace < 0: inplace = None
   
    return _lfilter(b, a, x, axis=axis, zi=zi,
             inplace=inplace, direction=direction)



Regards,
S.M.