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