# [SciPy-dev] reimplementation of lfilter

Sturla Molden sturla@molden...
Tue Sep 22 01:21:44 CDT 2009

```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.scipy.org/pipermail/scipy-dev/attachments/20090922/f052db95/attachment-0001.pl
-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: _linear_filter.c
Url: http://mail.scipy.org/pipermail/scipy-dev/attachments/20090922/f052db95/attachment-0001.c
```