[SciPy-dev] reimplementation of lfilter

[Sturla Molden]

Looking at this atrocious piece of C:

http://svn.scipy.org/svn/scipy/trunk/scipy/signal/lfilter.c.src

I got some inspiration for reusing some of my other C code to
reimplement scipy.signal.lfilter. I basically just had to write
this small C function + a Cython wrapper:


static void lfilter_1d(double *restrict b, double *restrict a,
    double *x, double *y, double *zi, npy_intp K, npy_intp N)
{
    double Z[K];
    double *restrict pz = zi, *restrict z = Z;
    register npy_intp n, k;
    for (n=0; n<N; n++) {
        register double yn, xn = x[n];
        yn = y[n] = b[0] * xn + pz[0];
        for (k=0; k<K-1; k++)
            z[k] = b[k+1] * xn + pz[k+1] - a[k+1]*yn;
        z[K-1] = b[K]*xn - a[K]*yn;
        double *tmp = pz;
        pz = z;
        z = tmp;
    }
    if (pz != zi)     
        memcpy(zi, pz, K*sizeof(double));       
}


Relative speed compared to scipy.signal.lfilter:

One processor:
1D array, 1000000 floats:              150%
2D array, 1000 x 1000 floats, axis=0:   80%
2D array, 1000 x 1000 floats, axis=1:  150%

It is faster than lfilter except for 2D with axis=0. For
2D with axis=1 it is probably faster than lfilter due to
copy-in copy-out optimization. For arrays with more than
one dimension, multithreading helps as well:

Four processors:
1D array, 1000000 floats:              150%
2D array, 1000 x 1000 floats, axis=0:  160%
2D array, 1000 x 1000 floats, axis=1:  250%

Multithreading obviously does not help for filtering single
vector, as there is no work to share.

Some other improvements over signal.lfilter, apart from readable
Cython code:

* GIL is released. The C code is not littered with Python C
  API calls that prevents this.

* Filtering can be done in reverse, i.e. no need for fliplr or
  flipud for filtering backwards.

* Arrays can be filtered in-place.

* Cleaned up docstring.

Still it only work as a "proof of concept" with double precision
floats. Adding support for other dtypes would be easy, just
replace double with any of:

float
float _Complex
double _Complex
long double
long double _Complex



Regards,
Sturla Molden







-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: linear_filter.pyx
URL: <http://mail.python.org/pipermail/scipy-dev/attachments/20090922/f052db95/attachment.ksh>
-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: _linear_filter.c
URL: <http://mail.python.org/pipermail/scipy-dev/attachments/20090922/f052db95/attachment.c>