```<clay@claysturner.com> wrote:

>On Mar 3, 2:03&#2013266080;pm, spop...@speedymail.org (Steve Pope) wrote:

>> I do not actually think the standard formula is mis-behaved
>> if you compute it in the correct order.

>To get an idea of the issue is look at using equally spaced data and
>when you get up to around 50 points, the dynamic range of the weights
>will be near 15 orders of magnitude. At this point double precision
>floating point will start to have serious round off issues.

Thanks.

I've never needed to use a Lagrangian interpolator with anywhere
near that many points.  Often just four points results in
a signal distortion of -60 dBc.  I would randomly guess that
50 points would give you a residual that is hundreds of dB down.

Depending on scenario, etc.

Steve
```
```On Mar 3, 2:03&#2013266080;pm, spop...@speedymail.org (Steve Pope) wrote:
> <c...@claysturner.com> wrote:
> >On Mar 3, 12:27&#2013266080;pm, spop...@speedymail.org (Steve Pope) wrote:
> >> I have a question: is the "barycentric" form a
> >> better-behaved numerical formula for the same interpolator
> >> as a standard Lagrangian, or is it a different interpolator function?
> >Yes the barycentric form is numerically superior to the standard
> >Lagrangian form. Ideally they will both give the same results.
>
> Thanks.
>
> I do not actually think the standard formula is mis-behaved
> if you compute it in the correct order.
>
> The formula for a Lagrange coefficient is similar to that of
> a binomial coefficient in that there are a large number of
> factors in both numerator and denominator; if you were to compute
> either of these by itself you might run into a range problem,
> but you don't have to.
>
> For example one of the formulas I gave is a constant times:
>
> &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080;(prod over t=1..n)(p + n/2 - t)
> &#2013266080; &#2013266080; &#2013266080; &#2013266080; __________________________________________________
>
> &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; ((n-2)/2 + k)! (n/2 - k)!
>
> It is clear there are n terms in the numerator, and n-1
> in the denominator, so by pairing these up you should
> avoid any range problem even for large n.
>
> This is not to say the barycentric formula isn't numerically
> superior, just that the original formula is also reasonable.
>
> Steve

To get an idea of the issue is look at using equally spaced data and
when you get up to around 50 points, the dynamic range of the weights
will be near 15 orders of magnitude. At this point double precision
floating point will start to have serious round off issues. The
barycentric form (form 2) is neat in that the numerator and the
denominator's errors tend to cancel each other out. The "standard"
Lagrange form is quite problematic.

Clay

```
```<clay@claysturner.com> wrote:

>On Mar 3, 12:27&#2013266080;pm, spop...@speedymail.org (Steve Pope) wrote:

>> I have a question: is the "barycentric" form a
>> better-behaved numerical formula for the same interpolator
>> as a standard Lagrangian, or is it a different interpolator function?

>Yes the barycentric form is numerically superior to the standard
>Lagrangian form. Ideally they will both give the same results.

Thanks.

I do not actually think the standard formula is mis-behaved
if you compute it in the correct order.

The formula for a Lagrange coefficient is similar to that of
a binomial coefficient in that there are a large number of
factors in both numerator and denominator; if you were to compute
either of these by itself you might run into a range problem,
but you don't have to.

For example one of the formulas I gave is a constant times:

(prod over t=1..n)(p + n/2 - t)
__________________________________________________

((n-2)/2 + k)! (n/2 - k)!

It is clear there are n terms in the numerator, and n-1
in the denominator, so by pairing these up you should
avoid any range problem even for large n.

This is not to say the barycentric formula isn't numerically
superior, just that the original formula is also reasonable.

Steve
```
```On Mar 3, 12:27&#2013266080;pm, spop...@speedymail.org (Steve Pope) wrote:
> <c...@claysturner.com> wrote:
> >Here is an excerp from a document I wrote
> >a while back that includes this process. I hope this answers your
> >question.
> >http://www.claysturner.com/polyresample.pdf
>
> Thanks. &#2013266080;I have a question: is the "barycentric" form a
> better-behaved numerical formula for the same interpolator
> as a standard Lagrangian, or is it a different interpolator function?
>
> Steve

Hello Steve,

Yes the barycentric form is numerically superior to the standard
Lagrangian form. Ideally they will both give the same results.

Clay

```
```Steve Pope wrote:

> I have a question: is the "barycentric" form a better-behaved
> numerical formula for the same interpolator as a standard
> Lagrangian, or is it a different interpolator function?

It's the same interpolator, but tamer and cheaper to compute.

Martin

--
Quidquid latine scriptum est, altum videtur.
```
```<clay@claysturner.com> wrote:

>Here is an excerp from a document I wrote
>a while back that includes this process. I hope this answers your
>question.

>http://www.claysturner.com/polyresample.pdf

Thanks.  I have a question: is the "barycentric" form a
better-behaved numerical formula for the same interpolator
as a standard Lagrangian, or is it a different interpolator function?

Steve
```
```Dave  <Confused.Scientist@gmail.com> wrote:

>On Mar 2, 10:51&#2013266080;pm, spop...@speedymail.org (Steve Pope) wrote:

>> There is a closed form for the required elements of W for
>> a Lagrangian interpolator of arbitrary order. &#2013266080;(It is a
>> little long to type in and post.)

>Indeed, the idea is to precompute W.  Do you have a reference for the
>university library but have been unable to find an analytic solution
>that I can cite.

I'll type it in.  Abscissas spaced by one, point of interest p
lying in (0,1), n-point interpolation:

f(p) = (sum over k) A(n,k,p) * f[k]

if n even, k =  -(n-2)/2 ... n/2

if n odd, k = -(n-1)/2 ... (n-2)/2

if n even, A(n,k,p) =

(-1)^(n/2 + k) * (prod over t=1..n)(p + n/2 - t)
__________________________________________________

((n-2)/2 + k)! (n/2 - k)! (p-k)

if n odd, A(n,k,p) =

(-1)^((n-1)/2 + k) * (prod over t=0..n-1)(p + (n-1)/2 - t)
__________________________________________________

((n-1)/2 + k)! ((n-1)/2 - k)! (p-k)

From:

Abrahamowitz and Stegun, _Handbook of Mathematical Functions_,
National Bureau of Standards, 4th printing December 1965, Ch. 25.

This should work, but let me know if it doesn't seem to be correct.

Steve
```
```Here's a more extensive treatment of barycentric interpolation:

http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.15.5097

Martin

--
Quidquid latine scriptum est, altum videtur.
```
```On Mar 2, 9:40&#2013266080;pm, Dave <Confused.Scient...@gmail.com> wrote:
> I'm trying to define a couple of common signal interpolation/
> resampling methods as a matrix operation, i.e.
>
> &#2013266080; y = W z
>
> where,
>
> &#2013266080; y is the interpolated data vector of dimension m x 1
> &#2013266080; W is the interpolation matrix of dimension m x l
> &#2013266080; z is the original data vector of dimension l x 1
>
> I've worked out W for linear and nearest neighbour interpolation
> (easy) and I was wondering if more complicated approaches (cubic,
> quintic, Shannon, Lanczos) can be defined analytically? &#2013266080;I don't
> expect they are trivial to type out but I'd be grateful if anyone can
> point me to a suitable text.
>
> Thanks.

Hello Dave,

Yes this can be easily done. Here is an excerp from a document I wrote
a while back that includes this process. I hope this answers your
question.

http://www.claysturner.com/polyresample.pdf

IHTH,
Clay S. Turner

```
```On Mar 2, 10:51&#2013266080;pm, spop...@speedymail.org (Steve Pope) wrote:

> There is a closed form for the required elements of W for
> a Lagrangian interpolator of arbitrary order. &#2013266080;(It is a
> little long to type in and post.)

Indeed, the idea is to precompute W.  Do you have a reference for the