Reply by Steve Pope March 3, 20092009-03-03
<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
Reply by March 3, 20092009-03-03
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
Reply by Steve Pope March 3, 20092009-03-03
<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
Reply by March 3, 20092009-03-03
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
Reply by Martin Eisenberg March 3, 20092009-03-03
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.
Reply by Steve Pope March 3, 20092009-03-03
<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
Reply by Steve Pope March 3, 20092009-03-03
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 >Lagrangian interpolator solution? I have access to a pretty good >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
Reply by Martin Eisenberg March 3, 20092009-03-03
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.
Reply by March 3, 20092009-03-03
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
Reply by Dave March 3, 20092009-03-03
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 Lagrangian interpolator solution? I have access to a pretty good university library but have been unable to find an analytic solution that I can cite. Computationally, I realize this is probably not the most efficient way to resample a single data vector. The reason for wanting to find W is that the signal is corrupt with additive Gaussian noise. Unfortunately, noise in adjacent channels is correlated due to a poor receiver/amplifier design and I've been asked to evaluate the correlation after resampling, i.e. Cov' = W^T Cov W, and co-addition of the signals. I've done this numerically using a simulink model but I was asked if I could approach the problem with algebra and see what comes out.