Saturday, January 06, 2007

Interpolation

For those that aren't familiar with the term, and since I'm not too good at creating my own definitions (go on, ask Andy Kramek!), I thought I'd borrow a definition from wikipedia to get started:

"In the mathematical subfield of numerical analysis, interpolation is a method of constructing new data points from a discrete set of known data points." (Source)


If you follow that source link to the wikipedia article, you can read lots of nice things about interpolation and how it can be used in the math and science fields. About 8 years ago, I needed such a function: I needed the ability to find a set of values (x) that fell between two known values (f(x)). I wrote the following linear interpolation algorithm in FoxPro to get my answers:


*-- this sets up a test cursor with some data points
CREATE CURSOR crDataSet (fld_x n(4,2), fld_fx n(4,2))
INDEX ON fld_x TAG fld_x
INSERT INTO crDataSet (fld_x , fld_fx) VALUES (1,1.75)
INSERT INTO crDataSet (fld_x , fld_fx) VALUES (3,2.00)
INSERT INTO crDataSet (fld_x , fld_fx) VALUES (4,3.00)
INSERT INTO crDataSet (fld_x , fld_fx) VALUES (5,4.00)
INSERT INTO crDataSet (fld_x , fld_fx) VALUES (7,4.25)

? interpolate(2,"crDataSet") && returns 1.875
? interpolate(3,"crDataSet") && returns 2.000
? interpolate(3.1,"crDataSet") && returns 2.100
? interpolate(3.9,"crDataSet") && returns 2.900
? interpolate(10,"crDataSet") && returns 4.625

FUNCTION interpolate
PARAMETERS tnX, tcDataSet

LOCAL lcfld_x , lcfld_fx AS String
LOCAL lcOldNear , lcOldExac AS String
LOCAL X1, Y1, X2, Y2 AS Number
LOCAL lFound AS Logical

*-- name of fields in tcDataSet cursor
lcfld_x = FIELD(1)
lcfld_fx = FIELD(2)

*-- important to adjust these settings
lcOldNear = SET('NEAR')
lcOldExac = SET('EXACT')
SET NEAR ON
SET EXACT OFF

SELECT (tcDataSet)

lFound = SEEK(tnX,tcDataSet)
SET NEAR &lcOldNear
SET EXACT &lcOldExac

IF lFound
*-- if found, then just return the point
RETURN EVAL(tcDataSet + '.' + lcfld_fx)
ELSE
*- not found, so need to interpolate
cXStr = tcDataSet + '.' + lcfld_x
cYStr = tcDataSet + '.' + lcfld_fx

*-- if the seek put us at eof, back up 1
SELE (tcDataSet)
IF EOF(tcDataSet)
SKIP -1
ENDIF

*-- evaluate the first point
X2 = EVAL(cXStr)
Y2 = EVAL(cYStr)

*-- evaluate the second point
IF RECNO(tcDataSet) = 1
SKIP +1
ELSE
SKIP -1
ENDIF
X1 = EVAL(cXStr)
Y1 = EVAL(cYStr)

*-- draw a line and return the result!
RETURN ymxb(x1,y1,x2,y2,tnX)
ENDIF

ENDFUNC

FUNCTION ymxb
PARAMETERS X1,Y1,X2,Y2,X3

*-- y = mx + b
LOCAL nSlope , B1

nSlope = (Y1-Y2)/(X1-X2)
B1 = -((nSlope) * X1 - Y1)
RETURN (nSlope) * X3 + B1

ENDFUNC


There you have it!

A couple of notes: For space reasons, I've omitted some key error checking and bullet-proofing (I started doing it, but realized that an extra 20 lines of code for this example was overkill). But it would be a good idea, for example to verify that the columns in the cursor passed were numbers. Come to think of it, you should also very that the correct parameters were passed!

In addition, I remember from math class that such a linear method is not the most accurate method. I am using a straight line between points and not a curve. I suppose a complimentary function could be added along with ymxb to get a curve between the two points, but until I need to do that -- I'm not game!

No comments: