BorlandTalk.com

Ronaldo Souza · Guest

Hi,

I need to interpolate a curve through a noisy signal. This curve is just
for trend visualization, I don't need the equation for the curve itself
and it might be many curves. The curve does not have to pass through
each point, just "in between" is what I'm looking for.

x = data points
- = interpolated curve

x x
x +----+
/ x \ x
/ +--+
x \ x

At first I thought a simple running average would be ok, but it turns
out I can have as few as 15 data points and, in these cases, I couldn't
get it to look right. Then I thought a polynomial least squares fit
would do the job, and it almost did. I wanted to find the degree of the
best fit automatically, using the Correlation Coefficient or the Average
Percent Absolute Error, but my "best fit" criteria isn't related to
minimum error and I can't define it mathematically, it's more a "look &
feel" thing. I ran fits from 2nd to 8th degree and the best fit using CC
or APAE would always be of the 8th order but the curve would have
undesired bumps and/or valleys. Visually I could assert that a 4th - 6th
degree fit would be exactly what I wanted, but could not come up any
algorithm to simulate whatever made me choose these lower order fits.

Please take a look at 2 pictures I uploaded to the attachments group.
The first is a 6th degree polynomial which "looks best". The second is
an 8th degree polynomial with a CC = 1.0 and the lowest APAE, but notice
the spurious "valley" on the left side.

What I need is either some algorithm to choose the best fit in a "visual
sense" or another interpolating algorithm altogether.

Any help will be greatly appreciated.

TIA,
Ronaldo

David Ninnes · Guest

Ronaldo Souza wrote:

Ronaldo Souza · Guest

David Ninnes wrote:

Nils Haeck · Guest

Yes, this is a common problem with interpolation of datapoints, when using
polynomial fits. The higher the degree, the more valleys and hills you get,
and you only interpolated.. you didn't even look yet at extrapolating the
data (you would be scared!). This is caused by the polynomial being
influenced by all the points in every point in the curve, so in theory the
curve on the complete left side is influenced by a variation of points on
the complete right side, with very unexpected results at times.

I think the best solution for this problem is not a polynomial fit, but a
b-spline curve. Here you can choose a certain degree, and each part of the
curve is only influenced by a number of "close" nodes (the nodes are your
points). You could also use a piece-wise fit using bezier curves (which is
at some point the same). It all depends a bit on how smooth your curve must
be.. up till which order the curve must be continuous (0-th order: just a
curve that is continuous, 1-th order: the slope is also continuous.. etc).

If you have many points, you could also look into simplifying your data
first.. see e.g. this piece of code I wrote:
http://www.simdesign.nl/Components/DouglasPeucker.html

I also have code for bspline and bezier curve matching, but that's not free.
For b-splines I'd advise you to pick up a good graphics book, as it is not
that straight-forward.

Nils

Rimvydas Paulavicius · Guest

This from my (very) old 1992 year collection:

procedure ChoBand( N, MP1 : integer; A, L : floating_array_ref_type );
{-------------------------------------------------------------------|
| Triangualizes banded symmetrical matrix A. |
|-------------------------------------------------------------------|
| Reference : |
| Wilkinson J. H., Reinsch C. |
| Handbook for automatic computations. Linear algebra. |
| New York Springer-Verlag Berlin, 1976, Algorithm I.4 |
|-------------------------------------------------------------------|
| Data : |
| A - lower triangle elements of banded symmetrical matrix where |
| A[*,MP1] are the diagonal elements. |
| For example if MP1=4 then A is processed as follows : |
| |
| * * * A[1,1] |
| * * A[2,1] A[2,2] |
| * A[3,1] A[3,2] A[3,3] |
| A[4,1] A[4,2] A[4,3] A[4,4] |
| ----------------------------------------------------- |
| A[N,N-3] A[N,N-2] A[N,N-1] A[N,N] |
|-------------------------------------------------------------------|
| Rezult : , |
| L - triangualization of matrix A = L * L . Has the same structu-|
| re as A. |
|-------------------------------------------------------------------}
var
P, i, r, j, q, k : integer;
Y : double;
begin
for i := 1 to N do begin
if i > MP1
then P := 1
else P := MP1-i+1;
r := i-MP1+P;
for j := P to MP1 do begin
q := MP1-j+p;
Y := A^[(i-1)*MP1+j];
for k := P to j-1 do begin
Y := Y - L^[(i-1)*MP1+k] * L^[(r-1)*MP1+q];
Inc(q);
end;{ k }
if j = MP1
then L^[(i-1)*MP1+j] := 1 / Sqrt( Y )
else L^[(i-1)*MP1+j] := Y*L^[r*MP1];
Inc(r);
end;{ j }
end;{ i }
end;{ ChoBand }

procedure ChoSol( N, MP1, R : integer; L, B, X : floating_array_ref_type );
{-------------------------------------------------------------------|
| , |
| Solves the system L * L * X = B . |
|-------------------------------------------------------------------|
| Reference : |
| Wilkinson J. H., Reinsch C. |
| Handbook for automatic computations. Linear algebra. |
| New York Springer-Verlag Berlin, 1976, Algorithm I.4 . |
|-------------------------------------------------------------------|
| Data : , |
| L * L - triangualization of banded symmetrical matrix. |
| B[1..R,1..N] - right sides |
|-------------------------------------------------------------------|
| Rezult : |
| X[1..R,1..N] - the solution |
|-------------------------------------------------------------------}
var
M , j, i, P, q, k : integer;
Y : double;
begin
M := Pred( MP1 );
for j := 1 to R do begin
{ Solve the system L * Y = B }
for i := 1 to N do begin
if i > M
then P := 1
else P := MP1-i+1;
Y := B^[ Pred(i)*R+j ];
q := i;
for k := M downto P do begin
Dec(q);
Y := Y - L^[ Pred(i)*MP1+k ] * X^[ Pred(j)*N+q ];
end;{ k }
X^[ Pred(j)*N+i ] := Y * L^[ i*MP1 ];
end;{ i }
{ Solve the system LT * X = Y }
for i := N downto 1 do begin
if (N-i)> M
then P := 1
else P := MP1-N+i;
Y := X^[ Pred(j)*N+i ];
q := i;
for k := M downto P do begin
Inc( q );
Y := Y - L^[ Pred(q)*MP1+k ] * X^[ Pred(j)*N+q ];
end;{ k }
X^[ Pred(j)*N+i ] := Y * L^[ i*MP1 ];
end;{ i }
end;{ j }
end;{ ChoSol }

procedure SmoothCubicSpline( N : integer ;
X, F : floating_array_ref_type;
Alpha : floating;
S : floating_array_ref_type );
{--------------------------------------------------------------------
| Smooths function values by using cubic spline approximations |
|-------------------------------------------------------------------|
| References : |
| 1) Vasilenko V.A. Teorija Splain-funkcij. Novosibirsk, NGU, 1978.|
| 2) Makarov B.L. Splain-aproksimacija funkcij. |
| Moskva, Vysshaja shkola, 1983 |
|-------------------------------------------------------------------|
| Data : |
| X[1..N] - argument values |
| F[1..N] - function to be smoothed values |
| Alpha - dumping factor |
|-------------------------------------------------------------------|
| Results : |
| S[1..N] - smoothed function values, might be the same as F |
--------------------------------------------------------------------}
var
H, A, B, L, W : floating_array_ref_type;
Beta : floating;
i : integer;
{ start }
begin
if N <= 0 then Exit;
S^[1] := F^[1];
if N <= 1 then Exit;
S^[2] := F^[2];
if N <= 2 then Exit;
{-------------------------------------------------------------------|
| form a linear system : |
| , |
| (A + Beta*G*G ) * W = G*F , where |
| |
| N * Alpha |
| Beta = ----------------- , |
| 6*Average(1/H**3) |
| |
| A = (H +H ) / 3 , A = H / 3 , |
| i,i i i+1 i,i+1 i+1 |
| |
| G = 1/H , G = -(1/H + 1/H ) , G = 1/H . |
| i,i i i,i+1 i i+1 i,i+2 i+1 |
| |
|-------------------------------------------------------------------}
GetMem( H, (N-1)*SizeOf(floating) );
GetMem( A, (N-2)*3*SizeOf(floating) );
GetMem( B, (N-2)*SizeOf(floating) );
GetMem( L, (N-2)*3*SizeOf(floating) );
GetMem( W, (N-2)*SizeOf(floating) );
Beta := 0;
for i := 1 to N-1 do begin
H^[i] := X^[i+1]-X^[i];
if H^[i] <= 0 then begin
Writeln('SmoothCubicSpline -> argument values aren''t in
ascending order');
Halt;
end;
Beta := Beta+Exp(-3*Ln(H^[i]));
end;{ i }
Beta := N*Alpha/(6*Beta/(N-1));
for i := 3 to N-2 do begin
A^[1+Pred(i)*3] := Beta/(H^[i-1]*H^[i]);
end;{ i }
for i := 2 to N-2 do begin
A^[2+Pred(i)*3] :=
H^[i]/6-Beta*(1/H^[i-1]+2/H^[i]+1/H^[i+1])/H^[i];
end;{ i }
for i := 1 to N-2 do begin
A^[3+Pred(i)*3] := (H^[i]+H^[i+1])/3 +

2*Beta*(1/Sqr(H^[i])+1/(H^[i]*H^[i+1])+1/Sqr(H^[i+1]));
end;{ i }
for i := 1 to N-2 do begin
B^[i] := (F^[i+2]-F^[i+1])/H^[i+1]-(F^[i+1]-F^[i])/H^[i];
end;{ i }
{ solve the system }
ChoBand( N-2, 3, A, L );
ChoSol( N-2, 3, 1, L, B, W );
{ evaluate smoothed F values }
S^[1] := F^[1]-Beta* W^[1]/H^[1];
S^[2] := F^[2]-Beta*(-W^[1]*(1/H^[1]+1/H^[2])+W^[2]/H^[2]);
for i := 3 to N-2 do begin
S^[i] := F^[i] - Beta*

(W^[i-2]/H^[i-1]-W^[i-1]*(1/H^[i-1]+1/H^[i])+W^[i]/H^[i]);
end;{ i }
S^[N-1] :=
F^[N-1]-Beta*(W^[N-3]/H^[N-2]-W^[N-2]*(1/H^[N-2]+1/H^[N-1]));
S^[N] := F^[N]-Beta*W^[N-2]/H^[N-1];
{ finish }
FreeMem( W, (N-2)*SizeOf(floating) );
FreeMem( L, 3*(N-2)*SizeOf(floating) );
FreeMem( B, (N-2)*SizeOf(floating) );
FreeMem( A, 3*(N-2)*SizeOf(floating) );
FreeMem( H, (N-1)*SizeOf(floating) );
{ exit }
end;{ SmoothCubicSpline }

Rimvydas

Ronaldo Souza wrote: