This is the function which computes the parabolic approximations. The function receives an input vector, x and should fill three vectors, a0, a1 and a2 which represent the three coefficients of the parabolic approximation.
The computation of the approximation is done with the
aid of a Maple program
which will compute the value of Compute_as.
As can be seen from the formula manipulation, to compute the
parabolic least squares approximation, three sums are needed
(sx, sjx and sj2x).
These can be computed for the index i based on the
ones for i-1.
Procedure( MovingParabApprox,
Params( [n,m] = Integer, [x,a0,a1,a2] = Array(n,Float)),
Block(
Declaration(
i = [Local,Integer],
[sx,sjx,sj2x] = [Local,Float] ),
If( m < 3, Error("Parabolic order not large enough") ),
If( n >= 1, Block( Assign(a0[1],x[1]), Assign(a1[1],0),
Assign(a2[1],0))),
If( n >= 2, Block( Assign(a0[2],x[2]), Assign(a1[2],x[2]-x[1]),
Assign(a2[2],0) )),
Comment( "The definitions of sx, sjx, and sj2x are:
m - 1 m - 1 m - 1
----- ----- -----
\ \ \ 2
) x[i - j] = sx, ) j x[i - j] = sjx, ) j x[i - j] = sj2x
/ / /
----- ----- -----
j = 0 j = 0 j = 0 "),
If( n >= 3, Block(
Assign( sx, x[1]+x[2] ),
Assign( sjx, x[1] ),
Assign( sj2x, x[1] ) )
),
ForLoop( i, 3, 1, min(n,m), Block(
Assign( sj2x, sj2x+2*sjx+sx ),
Assign( sjx, sjx+sx ),
Assign( sx, sx+x[i] ),
subs( m=i, Compute_as))
),
ForLoop( i, m+1, 1, n, Block(
Assign( sj2x, sj2x-x[i-m]*(m-1)^2 ),
Assign( sjx, sjx-x[i-m]*(m-1) ),
Assign( sx, sx-x[i-m] ),
Assign( sj2x, sj2x+2*sjx+sx ),
Assign( sjx, sjx+sx ),
Assign( sx, sx+x[i] ),
Compute_as)
),
Return
)
);
Notes:
For i=1 and i=2, the general formulas do not work, and a direct assignment is best.
The comment has been cut-and-pasted from the output of the symbolic algebra system.
The computation of the partial sums is different for i<=m and i>m. For the latter, the last term has to be subtracted off.
Two values of Compute_as are needed. For the first the limit of the sum is i and for the second is m. The values computed are based on m, hence the substitution.