[SciPy-dev] reimplementation of lfilter

David Cournapeau cournape@gmail....
Wed Sep 23 03:31:40 CDT 2009


On Wed, Sep 23, 2009 at 3:10 PM, Sturla Molden <sturla@molden.no> wrote:
> Sturla Molden skrev:
>>
>> Well, I changed all the code to C++, to use the std::complex type,
>> std::vector instead of malloc, and templates for specializing the filter
>> functions to all dtypes supported by lfilter. I'm quite happy with the way
>> the C++ looks :-)
>
> I forgot to add support for the object dtype.
>
> This also fixes a bug in the current signal.lfilter. This will crash the
> interpreter:
>
>  import numpy as np
>  import scipy
>  from scipy.signal import butter, lfilter
>  b,a = butter(4, .25)
>  x = np.random.randn(1000000).astype(object).reshape((1000,1000))
>  fx = lfilter(b,a,x, axis=1)
>
> The C++ version here does not crash on this. :-)
>
> I just ripped some code from current lfilter and cleaned it up slightly; I
> hope that is ok. (Cython was not happy about "cdef object *ptr", so I
> reluctantly had to put it all in C++.)  The reason my code does not crash
> like current lfilter, is copy-in copy-out for strided vectors. It seems
> lfilter messes up the strides and then segfaults.

I am all for improving the lfilter code, here a few comments:
 - it is much more likely that your improvements will be included if
you provide patches instead of rewrite of the full code - it makes
reviewing much easier.
 - I would also prefer having C instead of C++ as well - in this case,
C++ does not bring much since we have our "templating" system and you
don't use the STL much.
 - In any case, please do not use exception, it is not portable.
 - you cannot use Py_ssize_t, as it is python 2.5 >= feature - there
is nothing wrong with npy_intp, I don't understand your comment.
 - using cython is good - that could be a first patch, to replace the
existing manual wrapping by cython. You can declare pointer without
trouble in cython, you would have to be more explicit about the exact
problem you had.

cheers,

David


More information about the Scipy-dev mailing list