Plotting and fitting

A raw sensor reading is not a measurement until you know what it means. A potentiometer on a robot joint returns ADC counts, not degrees. The fix is a calibration: record the counts at several known angles, fit a curve through them, and use that curve to convert every future reading. Octave gives you the three tools for this in a few lines: plot to look, polyfit to fit, polyval to apply.

The calibration data

You clamp the joint at seven angles measured with a protractor and write down the ADC count at each. This is the data we will use throughout.

The counts rise by about 9 per degree, but the steps are not quite even: 81, 87, 91, 88, 92, 97. A pot with a slightly non-linear track or a lever arm geometry would do this. A straight line may or may not be good enough, and we want a number to decide.

Fitting with polyfit

polyfit(x, y, n) finds the polynomial of degree n that minimises the sum of squared errors, and returns its coefficients with the highest power first. polyval(p, x) evaluates that polynomial. Under the hood polyfit builds a matrix of powers of x and solves it with the same backslash from the last lesson, which handles more equations than unknowns as a least-squares problem.

x = [0 10 20 30 40 50 60];            % reference angle, degrees
y = [102 183 270 361 449 541 638];    % ADC counts
p1 = polyfit(x, y, 1);
p2 = polyfit(x, y, 2);
printf("line:  slope %.3f, offset %.3f\n", p1);
printf("quad:  %.4f x^2 + %.3f x + %.3f\n", p2);

r1 = y - polyval(p1, x);
r2 = y - polyval(p2, x);
printf("line residuals: %s\n", mat2str(round(r1 * 100) / 100));
printf("rms residual: line %.3f, quad %.3f\n", sqrt(mean(r1.^2)), sqrt(mean(r2.^2)));

r = roots([p2(1), p2(2), p2(3) - 300]);
theta = r(r >= 0 & r <= 60);
printf("reading 300 -> angle %.2f deg\n", theta);

line:  slope 8.939, offset 95.250
quad:  0.0118 x^2 + 8.232 x + 101.143
line residuals: [6.75 -1.64 -4.04 -2.43 -3.82 -1.21 6.39]
rms residual: line 4.266, quad 1.237
reading 300 -> angle 23.37 deg

The line says 8.939 counts per degree with 95.25 counts at zero. The quadratic adds a small curvature term of 0.0118 and slightly changes the other two coefficients. The last two lines solve the calibration backwards: to turn a reading of 300 into an angle you set p2(x) = 300 and take the root. A quadratic has two roots, here about 23.37 and about -721.9, and only one lies inside the 0 to 60 degree range of the sensor, so the logical test picks it out.

Residuals tell the truth

A residual is measured minus fitted. The overall size is summarised by the root-mean-square (rms): 4.266 counts for the line and 1.237 for the quadratic. But the sign pattern is more informative than the total. Look at the line's residuals: positive at both ends (6.75 and 6.39) and negative through the middle. That is a bowl, a systematic shape the straight line cannot follow. Residuals from a good model look like random scatter, and the quadratic's do: they alternate between about -1.6 and +2.3 with no pattern.

Do not conclude that higher degrees are always better. With seven points, a degree 6 polynomial passes exactly through every one of them, residual zero, and then swings wildly between the points. It has memorised the noise. Prefer the lowest degree whose residuals look random, and when you can, test on points you did not fit.

Making the figure

Plots open in a window when you run these commands in Octave. Here is the script that draws the data, both fits, and the residuals.

xx = linspace(0, 60, 121);
plot(x, y, 'o', xx, polyval(p1, xx), '-', xx, polyval(p2, xx), '--');
xlabel('reference angle (deg)'); ylabel('ADC counts');
legend('measured', 'line', 'quadratic', 'location', 'northwest');
grid on;
print -dpng calibration.png

figure;
plot(x, r1, 'o-', x, r2, 's-');
xlabel('reference angle (deg)'); ylabel('residual (counts)');
legend('line', 'quadratic'); grid on;

Since a figure cannot be shown here, this is what you will see. The first figure has seven circles climbing from the lower left to the upper right. A solid line runs through them, passing slightly below the first and last circles and slightly above the middle ones. The dashed curve bends gently upward and passes within about two counts of every circle. At this scale the two curves look almost the same, which is exactly why you need the second figure. It shows the residuals: the line's curve starts near +7, dips to about -4 around 20 to 40 degrees and climbs back to +6, a clear bowl, while the quadratic's curve zigzags close to zero. The print command writes the picture to a PNG file, useful on a robot without a display.

The same maths in NumPy

The NumPy runner prints the same coefficients, so you can verify the numbers above and then change the data.

numpy_version.py

ResetRun

import numpy as np x = np.array([0, 10, 20, 30, 40, 50, 60], dtype=float) # reference angle, deg y = np.array([102, 183, 270, 361, 449, 541, 638], dtype=float) # ADC counts p1 = np.polyfit(x, y, 1) p2 = np.polyfit(x, y, 2) print("line: slope %.3f, offset %.3f" % tuple(p1)) print("quad: %.4f x^2 + %.3f x + %.3f" % tuple(p2)) r1 = y - np.polyval(p1, x) r2 = y - np.polyval(p2, x) print("line residuals: [%s]" % " ".join("%.2f" % v for v in r1)) print("rms residual: line %.3f, quad %.3f" % (np.sqrt(np.mean(r1 ** 2)), np.sqrt(np.mean(r2 ** 2)))) roots = np.roots([p2[0], p2[1], p2[2] - 300]) theta = roots[(roots.real >= 0) & (roots.real <= 60)].real[0] print("reading 300 -> angle %.2f deg" % theta) # The same quadratic through least squares, as in lesson 2 V = np.column_stack([x ** 2, x, np.ones(len(x))]) coef = np.linalg.lstsq(V, y, rcond=None)[0] print("least squares matches polyfit:", np.allclose(coef, p2))

Try np.polyfit(x, y, 6) and print the rms residual: it is essentially zero, and the fit is useless between the points. Also change one count by 20 and see how much the coefficients move; that is a feel for how noise in calibration data turns into error in every future reading.