Problem: Given a six-sided die, fitted with accelerometers, we want to estimate which side is up by just looking at the accelerometer data.
Solution: We model each side by a multivariate Gaussian, and perform Maximum Likelihood (ML) estimation.
Concrete: Within the die, two dual-axis accelerometers are fitted, so we'll have 4 values, two of which should behave approximately the same. See the cubicles website for more information. The input vectors, Gaussians, and mean vectors are therefore four-dimensional.
[1: Gaussian Modelling] [2: Maximum Likelihood Estimation] [3: generating the PICC-code] [4: visualization] [5: adding gestures]
We start by reading the datasets that contain the sensordata from the cube, in this case for two different surfaces (table and carpet). By plotting, it becomes clear what parts of the dataset are the various states of the cube.
>> cubicle1 = textread('cubicle1.txt');
>> cubicle2 = textread('cubicle2.txt');
>> cubicle1 = cubicle1(:,1:4);
% just take the accelerometer values
>> cubicle2 = cubicle2(:,1:4);
% just take the accelerometer values
>> plot(cubicle1)
>> plot(cubicle2)
We get to the six different datasets (for each side being the one on top) by looking at the different states in the figures..
>> side1 = [cubicle1(1:150,:); cubicle2(1:150,:)];
>> side2 = [cubicle1(200:350,:); cubicle2(225:375,:)];
>> side3 = [cubicle1(423:573,:); cubicle2(417:567,:)];
>> side4 = [cubicle1(612:762,:); cubicle2(595:745,:)];
>> side5 = [cubicle1(810:940,:); cubicle2(800:950,:)];
>> side6 = [cubicle1(1000:1150,:); cubicle2(1000:1110,:)];
So we have more or less clear datasets to start from:
>> plot([side1; side2; side3; side4; side5; side6]);
And calculate the mean vector for each dataset, plus the covariance matrix:
>> mean1 = mean(side1);
>> mean2 = mean(side2);
>> mean3 = mean(side3);
>> mean4 = mean(side4);
>> mean5 = mean(side5);
>> mean6 = mean(side6);
>> cov1 = cov(side1);
>> cov2 = cov(side2);
>> cov3 = cov(side3);
>> cov4 = cov(side4);
>> cov5 = cov(side5);
>> cov6 = cov(side6);
Which gives:
>> disp(mean1); disp(mean2); disp(mean3); disp(mean4);disp(mean5);
disp(mean6);
504.0900 522.2233 462.0733 525.8733
564.2980 519.8775 522.8113 523.1358
502.8742 463.0762 518.4834 466.0728
446.5497 523.8907 516.3046 528.8411
507.7553 580.4716 520.1135 585.6738
506.6260 521.0802 577.2519 525.6069
>> disp(cov1); disp(cov2); disp(cov3); disp(cov4);disp(cov5); disp(cov6);
0.3430 -0.0302 -0.0200 0.0148
-0.0302 0.5687 -0.0465 0.2190
-0.0200 -0.0465 1.1585 -0.1010
0.0148 0.2190 -0.1010 0.5525
0.3959 -0.0830 -0.0133 0.0059
-0.0830 0.4268 0.1662 0.1164
-0.0133 0.1662 1.1835 0.0290
0.0059 0.1164 0.0290 0.4965
1.5256 0.1259 -0.1084 0.1886
0.1259 0.8546 0.0627 0.6456
-0.1084 0.0627 0.9815 0.0112
0.1886 0.6456 0.0112 1.0578
0.4278 -0.0693 0.0014 0.0013
-0.0693 0.4033 0.0168 0.1254
0.0014 0.0168 1.3487 -0.0577
0.0013 0.1254 -0.0577 0.5394
0.4061 -0.2045 0.0777 -0.1193
-0.2045 0.6415 -0.1320 0.2434
0.0777 -0.1320 1.1828 -0.2155
-0.1193 0.2434 -0.2155 0.6334
0.4189 -0.0963 -0.0165 -0.0403
-0.0963 0.3499 -0.0049 0.0393
-0.0165 -0.0049 0.9555 0.0075
-0.0403 0.0393 0.0075 0.3851
And with these values we can complete the constants we need for the Maximum Likelihood Estimation:
>> ln1 =
log(det(cov1)); ln2 = log(det(cov2)); ln3 = log(det(cov3)); ln4 =
log(det(cov4)); ln5 = log(det(cov5)); ln6 = log(det(cov6));
>> disp( [ln1 ln2 ln3 ln4 ln5 ln6] );
-2.2723 -2.4796 -0.3534 -2.1876 -2.0451 -3.0039
>> inv1 = inv(cov1); inv2 = inv(cov2); inv3 =
inv(cov3); inv4 = inv(cov4); inv5 = inv(cov5); inv6 = inv(cov6);
>> disp(inv1); disp(inv2); disp(inv3); disp(inv4); disp(inv5); disp(inv6);
2.9440 0.2209 0.0459 -0.1580
0.2209 2.0920 0.0152 -0.8325
0.0459 0.0152 0.8780 0.1533
-0.1580 -0.8325 0.1533 2.1724
2.6476 0.5779 -0.0473 -0.1642
0.5779 2.7713 -0.3670 -0.6350
-0.0473 -0.3670 0.8951 0.0342
-0.1642 -0.6350 0.0342 2.1630
0.6762 -0.0252 0.0775 -0.1060
-0.0252 2.1870 -0.1274 -1.3289
0.0775 -0.1274 1.0350 0.0530
-0.1060 -1.3289 0.0530 1.7748
2.4110 0.4496 -0.0129 -0.1116
0.4496 2.7615 -0.0627 -0.6499
-0.0129 -0.0627 0.7463 0.0945
-0.1116 -0.6499 0.0945 2.0153
2.9651 0.8529 -0.0614 0.2097
0.8529 2.0774 0.0636 -0.6160
-0.0614 0.0636 0.9062 0.2724
0.2097 -0.6160 0.2724 1.9476
2.5653 0.6848 0.0464 0.1979
0.6848 3.0744 0.0296 -0.2426
0.0464 0.0296 1.0477 -0.0185
0.1979 -0.2426 -0.0185 2.6427
After this 'training phase', we can hardcode the results (marked in yellow) in the code for the microcontroller, using the constant structures ln[], mean[], and inv[].
For estimating what class the incoming accelerometer values (x1,x2,x3,x4) belong to, we just use this formula (the marked items are the hardcoded bits we calculated above):
estimated side = minarg[ ln(determinant(covar))+ (x-mean)*inverse(covar)*(x-mean)]
Notice that the second term in the minarg function is the Mahalanobis metric between the current sensor values and the means.
In C sourcecode, we translate this to:
#define MAX_CLASS 6
const float mean[MAX_CLASS][4] = {{504.0900, 522.2233, 462.0733, 525.8733},
{564.2980, 519.8775, 522.8113, 523.1358},
{502.8742, 463.0762, 518.4834, 466.0728},
{446.5497, 523.8907, 516.3046, 528.8411},
{507.7553, 580.4716, 520.1135, 585.6738},
{506.6260, 521.0802, 577.2519, 525.6069}};
const float inv11[MAX_CLASS] = {2.9440, 2.6476, 0.6762, 2.4110, 2.9651, 2.5653};
const float inv12[MAX_CLASS] = {0.2209, 0.5779, -0.0252, 0.4496, 0.8529,
0.6848};
const float inv13[MAX_CLASS] = {0.0459, -0.0473, 0.0775, -0.0129, -0.0614,
0.0464};
const float inv14[MAX_CLASS] = {-0.1580, -0.1642, -0.1060, -0.1116, 0.2097,
0.1979};
const float inv22[MAX_CLASS] = {2.0920, 2.7713, 2.1870, 2.7615, 2.0774, 3.0744};
const float inv23[MAX_CLASS] = {0.0152, -0.3670, -0.1274, -0.0627, 0.0636,
0.0296};
const float inv24[MAX_CLASS] = {-0.8325, -0.6350, -1.3289, -0.6499, -0.6160,
-0.2426};
const float inv33[MAX_CLASS] = {0.8780, 0.8951, 1.0350, 0.7463, 0.9062, 1.0477};
const float inv34[MAX_CLASS] = {0.1533, 0.0342, 0.0530, 0.0945, 0.2724,
-0.0185};
const float inv44[MAX_CLASS] = {2.1724, 2.1630, 1.7748, 2.0153, 1.9476, 2.6427};
const float ln[MAX_CLASS] = {-2.2723, -2.4796, -0.3534, -2.1876, -2.0451,
-3.0039};
// estimate the most likely class
int get_ml(long int x1, long int x2, long int x3, long int x4 ){
int k;
signed long int mx1, mx2, mx3, mx4;
float mah1, mah2, mah3, mah4;
float res, prev_res;
int temp_min;
// go over all classes
prev_res = 10000000000000.0;
temp_min = MAX_CLASS+1;
for (k=0; k<MAX_CLASS; k++) {
// calculate temporary means
mx1 = x1 - mean[k][0]; mx2 = x2 -
mean[k][1];
mx3 = x3 - mean[k][2]; mx4 = x4 -
mean[k][3];
//printf("\nTemp means: %li %li %li %li\n\r",
mx1, mx2, mx3, mx4);
// calculate 4d Mahalanobis:
mah1 = (float)mx1*inv11[k];
mah1 += (float)mx2*inv12[k];
mah1 += (float)mx3*inv13[k];
mah1 += (float)mx4*inv14[k];
mah1 *= (float)mx1;
mah2 =
(float)mx1*inv12[k]+(float)mx2*inv22[k];
mah2 +=
(float)mx3*inv23[k]+(float)mx4*inv24[k];
mah2 *= (float)mx2;
mah3 =
(float)mx1*inv13[k]+(float)mx2*inv23[k];
mah3 +=
(float)mx3*inv33[k]+(float)mx4*inv34[k];
mah3 *= (float)mx3;
mah4 =
(float)mx1*inv14[k]+(float)mx2*inv24[k];
mah4 +=
(float)mx3*inv34[k]+(float)mx4*inv44[k];
mah4 *= (float)mx4;
res = mah1 + mah2 + mah3 + mah4 +
ln[k];
// get the minarg:
if (res < prev_res) {
prev_res =
res;
temp_min = k
+ 1;
}
}
return temp_min;
}
void main() {
int ret;
long int x1,x2,x3,x4;
setup_adc_ports(ALL_ANALOG);
setup_adc(ADC_CLOCK_INTERNAL);
while (1) {
set_adc_channel(0); delay_us(100); x1
= read_adc();
set_adc_channel(1); delay_us(100); x2
= read_adc();
set_adc_channel(2); delay_us(100); x3
= read_adc();
set_adc_channel(3); delay_us(100); x4
= read_adc();
ret = get_ml(x1,x2,x3,x4);
printf("result=%i\n\r",ret);
}
}
The variable 'temp_min' will then return the number of the most likely side that facing up.
Document created by Kristof Van Laerhoven, March 2003.