From: Eric Durant <edurant@engin.umich.edu>
Subject: Re: Vectorize Lagrange interpolator
Date: 20 Nov 1999 00:00:00 GMT
Message-ID: <3836305B.C6914591@engin.umich.edu>
Content-Transfer-Encoding: 8bit
References: <80vi9v$nml$1@nnrp1.deja.com>
To: cykhung@my-deja.com
X-Accept-Language: en
Content-Type: text/plain; charset=iso-8859-1
Organization: University of Michigan EECS
Mime-Version: 1.0
Newsgroups: comp.soft-sys.matlab

cykhung@my-deja.com wrote:
>Is there a way to vectorize the code to generate a Lagrange 
>interpolator. A Lagrange interpolator is actually a FIR filter 
>generating a fractional delay in the lowpass region. The impulse 
>response of this FIR filter is,
> 
>       N
>       --    D - k
>h(n) = ||   -------
>       ||    n - k
>      k = 0
>      k ~= n
> 
>where D = given fractional delay (constant)
>      N = filter order.
> 
>I can generate h(n) with a loop.
> 
>h = zeros(N+1,1);
>for n = 1:N+1
>   k = [1:n-1, n+1:N+1];
>   h(n) = cumprod(D - k) / cumprod(n - k);
>end
>
>Note that the index starts from 1 (Matlab thing).

Here's an algorithm for finding the Lagrange interpolation polynomial
(equivalently, filter) that isn't vectorized, but it does eliminate the
redundant multiplications of the straightforward implementation.  Also,
it doesn't require huge matrix storage.  It allows irregularly spaced
values of the independent variable; since your application doesn't
require this, you could simplify it.

I noticed the matrix solution given in this thread doesn't give the
answer I expect; perhaps it is a scaling issue.  In any case, I've
included a test case.  If the algorithm is working correctly,
all(polyval(p,x)==y) returns 1 (if we're lucky enough that roundoff
errors are strictly negligible) or, more generally,
sum([polyval(p,x)-y].^2) approaches 0 as the roundoff error approaches
0.

I can provide the derivation if you're interested.  It uses a neat,
Horner-inspired factorization of the Lagrange polynomial.

=================================================
function [p] = lagrange(x,y)
% Eric Durant <edurant@umich.edu>
% Wednesday 6 October 1999
%
% Find the interpolation polynomial for f(x) = y
% using the Lagrange interpolation method.

l = length(y);
for i=1:l,
    alpha(i) = y(i);
    for j=1:i-1,
        alpha(i) = (alpha(i) - alpha(j)) / (x(i) - x(j));
    end
end
ip = 1;
p = zeros(1,l);
for i=1:l-1,
    p(l-i+1:l) = p(l-i+1:l) + alpha(i) * ip;
    ip = conv(ip, [1 -x(i)]);
end
p = p + alpha(l) * ip;

%%% Test case:
% » x = 1:4;
% » y = [2 7 18 41];
% » p = lagrange(x,y)
% 
% p =
%      1    -3     7    -3
% 
% » polyval(p,x)
% 
% ans =
%      2     7    18    41
=================================================

Eric Durant
http://www.edurant.com/