This program fits an order-five polynomial to a set of 20000 experimental data points. It demonstrates the power of BBC BASIC's built-in array and matrix operations to perform high-speed calculations. It could be adapted to carry out other least-squares fitting tasks.
| Download POLYFIT.BBC | Run POLYFIT.EXE |
|---|
REM POLYFIT - fit a polynomial to a set of experimental data
REM demonstrating the use of BBC BASIC's matrix operations.
REM In this version the maximum order of polynomial is five.
REM Richard T. Russell, richard@rtrussell.co.uk, 12-Apr-2011
HIMEM = PAGE + 5000000
INSTALL @lib$+"ARRAYLIB"
INSTALL @lib$+"FNUSING"
INSTALL @lib$+"WINLIB5"
*FLOAT 64
file$ = @dir$+"polyfit.csv"
max% = 20000
REM Initialise screen:
SYS "SetWindowText", @hwnd%, "Least squares polynomial fitting"
VDU 23,22,720;480;8,16,16,128
VDU 5
ORIGIN 140,168
COLOUR 2,0,100,0
PROCaxes
REM Create controls:
button1% = FN_button("Load && plot file", 190, 12, 120, 20, FN_setproc(PROCload), 0)
button2% = FN_button("Fit polynomial", 390, 12, 120, 20, FN_setproc(PROCfit), 0)
SYS "EnableWindow", button2%, 0
REM Declare arrays:
DIM vector(5), matrix(5,5)
DIM x(max%), x2(max%), x3(max%), x4(max%), x5(max%)
DIM x6(max%), x7(max%), x8(max%), x9(max%), x10(max%)
DIM y(max%), xy(max%), x2y(max%), x3y(max%), x4y(max%), x5y(max%)
REM Idle loop:
REPEAT
WAIT 1
UNTIL FALSE
END
REM Draw axes and labels:
DEF PROCaxes
CLS
GCOL 0
FOR x = 0.0 TO 1.0 STEP 0.1
MOVE x*1200, 0 : PLOT 21, x*1200, 700
MOVE x*1200-24, -24 : PRINT FNusing("#.#", x);
NEXT
FOR y = 0 TO 6
MOVE 0, 700*y/6 : PLOT 21, 1198, 700*y/6
MOVE -80, 700*y/6+16 : PRINT FNusing("#.#", y);
NEXT
RECTANGLE 0, 0, 1198, 700
ENDPROC
REM Load and plot file:
DEF PROCload
file% = OPENIN(file$)
IF file% = 0 ERROR 0, "Could not open file " + file$
PROCaxes
SYS "EnableWindow", button2%, 0
index% = 0
x() = 0.0
y() = 0.0
GCOL 4
WHILE NOT EOF#file% AND index% <= max%
record$ = GET$#file%
comma% = INSTR(record$, ",")
IF comma% THEN
x(index%) = VAL(record$)
y(index%) = VAL(MID$(record$,comma%+1))
LINE x(index%)*1200, y(index%)*700/6, x(index%)*1200, y(index%)*700/6
index% += 1
ENDIF
ENDWHILE
CLOSE #file%
npts% = index%
SYS "EnableWindow", button2%, 1
ENDPROC
REM Fit and plot polynomial:
DEF PROCfit
TIME = 0
sum_x = SUM(x())
x2() = x() * x() : sum_x2 = SUM(x2())
x3() = x() * x2() : sum_x3 = SUM(x3())
x4() = x2() * x2() : sum_x4 = SUM(x4())
x5() = x2() * x3() : sum_x5 = SUM(x5())
x6() = x3() * x3() : sum_x6 = SUM(x6())
x7() = x3() * x4() : sum_x7 = SUM(x7())
x8() = x4() * x4() : sum_x8 = SUM(x8())
x9() = x4() * x5() : sum_x9 = SUM(x9())
x10() = x5() * x5() : sum_x10 = SUM(x10())
sum_y = SUM(y())
xy() = x() * y() : sum_xy = SUM(xy())
x2y() = x2() * y() : sum_x2y = SUM(x2y())
x3y() = x3() * y() : sum_x3y = SUM(x3y())
x4y() = x4() * y() : sum_x4y = SUM(x4y())
x5y() = x5() * y() : sum_x5y = SUM(x5y())
matrix() = \
\ npts%, sum_x, sum_x2, sum_x3, sum_x4, sum_x5, \
\ sum_x, sum_x2, sum_x3, sum_x4, sum_x5, sum_x6, \
\ sum_x2, sum_x3, sum_x4, sum_x5, sum_x6, sum_x7, \
\ sum_x3, sum_x4, sum_x5, sum_x6, sum_x7, sum_x8, \
\ sum_x4, sum_x5, sum_x6, sum_x7, sum_x8, sum_x9, \
\ sum_x5, sum_x6, sum_x7, sum_x8, sum_x9, sum_x10
vector() = \
\ sum_y, sum_xy, sum_x2y, sum_x3y, sum_x4y, sum_x5y
PROC_invert(matrix())
vector() = matrix().vector()
took = TIME / 100
GCOL 1
MOVE 80, -72
PRINT "y = " FNusing("##.###", vector(0)) FNusing("+##.###", vector(1)) " x" \
\ FNusing("+##.###", vector(2)) " x^2" FNusing("+##.###", vector(3)) " x^3" \
\ FNusing("+##.###", vector(4)) " x^4" FNusing("+##.###", vector(5)) " x^5"
GCOL 2
MOVE 80, -120
PRINT "Fitting order-five polynomial to " ;npts% " points took " ;
PRINT FNusing("#.##", took) " seconds"
GCOL 9
FOR x = 0.0 TO 1.0 STEP 0.001
y = 0
FOR power% = 0 TO 5
y += vector(power%) * x^power%
NEXT
PLOT x*1200, y*700/6
NEXT x
ENDPROC
