[SciPy-Dev] Splines

[Charles R Harris]

On Mon, Mar 18, 2013 at 11:17 AM, denis <denis-bz-py at t-online.de> wrote:

> Charles R Harris <charlesr.harris <at> gmail.com> writes:
>
>
> > On Wed, Mar 13, 2013 at 12:54 PM, denis <denis-bz-py <at> t-online.de>
> wrote:
>
> > Chuck,
> >   is "prefilter=True will filter out high frequencies in the data"
> > not correct ?
> >
> > No, I suspect it increases the high frequencies ;) For cubic splines
> convolving the coefficients with [1/6, 2/3, 1/6] will reproduce the sample
> values, so to get the spline coefficients you need to deconvolve what is
> essentially a low pass filter.
>
> Chuck,
>   you're right; could spline_filter1d be buggy ?
>
>     f .5  x [ 100 -100  100 -100  100 -100  100 -100  100 -100]
>     np.convolve: [ 50 -33  33 -33  33 -33  33 -33  33 -50]
>     ndimage.convolve: [ 50 -33  33 -33  33 -33  33 -33  33 -50]
>     ndimage.spline order 2: [ 200 -200  200 -200  200 -200  200 -200  200
> -200]
>     ndimage.spline order 3: [ 300 -300  300 -300  300 -300  300 -300  300
> -300]
>     ndimage.spline order 4: [ 480 -480  480 -480  480 -480  480 -480  480
> -480]
>
>
Unlikely, probably there is a scaling factor floating about somewhere. You
need to be more precise about what you did for me to say more.

NI_SplineFilter1D in $scipysrc/ndimage/src/ni_interpolation.c has
>
>    case 2:
>         npoles = 1;
>         pole[0] = sqrt(8.0) - 3.0;
>         break;
>     case 3:
>         npoles = 1;
>         pole[0] = sqrt(3.0) - 2.0;
>
>
The roots of 1/6 + 2/3*x + 1/6*x**2 are -2 +/- sqrt(3), which is where the
pole in the (factored) inverse IIR filter for cubic splines, case 3, comes
from. Offhand, I don't see where spline_filter1d handles the boundary
conditions, which I assume are only the mirrored case. The usual way to
handle that, sometimes obscured in the literature, is to treat the boundary
conditions as a perturbation and use the
Sherman-Morrison<http://en.wikipedia.org/wiki/Sherman%E2%80%93Morrison_formula>formula
to deal with them. That is somewhat clearer if one sets up the
circulant matrix for periodic boundary conditions and uses its exact
factorization into circulant matrices as a starting point.

This brings up the whole question of boundary conditions for uniform
splines. The version in signal processing seems to be mirrored (like the
DCT II), and I assume the same for spline_filter1d. The reference for
spline_filter1d looks to be Unser's early
paper<http://www.google.com/url?sa=t&rct=j&q=&esrc=s&source=web&cd=1&cad=rja&ved=0CDMQFjAA&url=http%3A%2F%2Fusers.fmrib.ox.ac.uk%2F~jesper%2Fpapers%2Ffuture_readgroups%2Funser9302.pdf&ei=mlNHUZzMMcHwigKajYH4DQ&usg=AFQjCNF7RqY0hwNrZgAzod-1tnHffE3otw&bvm=bv.43828540,d.cGE>.
I think we need at least periodic conditions for closed loops, perhaps
reflected (DCT-I like), and among the remaining options would be
not-a-knot, and endpoint derivatives. Some of these get a bit tricky for
higher order splines, so perhaps the spline_filter1d choice of degree 5 as
the maximum is reasonable. Note that high degree uniform splines go over to
sinc interpolation methods and show the Gibbs phenomenon.

Chuck
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/scipy-dev/attachments/20130318/7836ab5b/attachment.html>