Recognizing Peaks in Accelerometer Data (Matlab)

Problem: Given accelerometer data, when and what size of peaks do we get in the signal(s)?.

Solution: We parse the accelerometer traces by first recognizing activity by thresholding the size and length of a running variance. Then, the available peak's are measured (by integration) and classified.

Using the method on this page, a toolkit (available upon request) provides visual and numerical feedback on the rough nature of the peaks in the signals. The figure on the right gives a sneak preview of what the toolkit looks like at the moment.

1. Activity Detection

The running variance is calculated over a sliding window n. As the variance is proportional to Σx²-(Σx)²/n, all we need to keep is a running sum and a running sum of squares. Try for instance:

A = textread('data.txt');

n = 10;
run_var = zeros(size(A,1),size(A,2));
var = zeros(size(A,1),size(A,2));

for i=1:length(A)-n,

sum = zeros(1,size(A,2));
sumsq = zeros(1,size(A,2));

for j=1:n,

sum = sum + A(i+j-1,:);
sumsq = sumsq + A(i+j-1,:).^2;

end

run_var(i,:) = sumsq - (sum.^2)/n;
var(i,:) = diag(cov(A(i:i+n-1,:)))';

end

subplot(3,1,1); plot(A);
subplot(3,1,3); plot(run_var);
subplot(3,1,2); plot(var);

A threshold is then empirically found to switch between activity detection:

t = 500;

for i=1:length(A)-n,

sum = zeros(1,size(A,2));
sumsq = zeros(1,size(A,2));
test = 0;
for j=1:n,
     sum = sum + A(i+j-1,:);
     sumsq = sumsq + A(i+j-1,:).^2;
end

for j=1:size(A,2),
     test = test + (sumsq(j) - (sum(j)^2)/n);
end
if test > 500,
     run_var(i,:) = 240;
end

end

plot([ A(1:length(run_var),:) run_var]);

Finally, the length of the found activity is the definite parameter to ignore most false positives (see graph).

2. Getting the Peak Parameters

For every piece of activity we found so far, we try to abstract the peaks that occur by:

subtracting a baseline from the data
finding the absolute area under all peaks
finding peaks by the net area the have under them (using a low threshold and the ratio to the absolute area)

There is a small delay (because of the running variance calculation), but this doesn't really matter too much if the window-size is small enough: (zoom from the previous plot)

These are the values we get from this particular gesture, for instance:

A = A(1400:1500,1:2);
n = 10; t = 500;
run_var = zeros(size(A,1));
last_action = [];

for i=1:length(A)-n,

sum = zeros(1,size(A,2));
sumsq = zeros(1,size(A,2));
test = 0;

for j=1:n,
    sum = sum + A(i+j-1,:);
    sumsq = sumsq + A(i+j-1,:).^2;
end

for j=1:size(A,2),
    test = test + (sumsq(j) - (sum(j)^2)/n);
end
if test > 500,
    run_var(i) = 240;
end

if run_var(i) > 0,
    last_action = [last_action; A(i,:)];
else
    % parse last action:

    %baselines:
     last_baseline = mean(last_action);
    %global area:
     last_globalsum = 0;
     for j=1:length(last_action),
          last_globalsum = last_globalsum + abs( last_action(j,:) - last_baseline );
     end
    %peaks:
     last_peak1 = 0;
     last_peak2 = 0;
     peaks1 = [];
     peaks2 = [];

     for j=1:length(last_action),
          if sign( last_action(j,1) - last_baseline(1) ) ~= sign(last_peak1) ,
           % new peak, so display last peak's info:
             peaks1 = [ peaks1 last_peak1 ];
             last_peak1 = last_action(j,1) - last_baseline(1);
          else
             last_peak1 = last_peak1 + (last_action(j,1) - last_baseline(1));
          end
          if sign( last_action(j,2) - last_baseline(2) ) ~= sign(last_peak2) ,
           % new peak, so display last peak's info:
             peaks2 = [ peaks2 last_peak2 ];
             last_peak2 = last_action(j,2) - last_baseline(1);
          else
             last_peak2 = last_peak2 + (last_action(j,2) - last_baseline(2));
          end
     end

     peaks1 = [ peaks1 last_peak1 ];
     peaks2 = [ peaks2 last_peak2 ];
     %last_action = [];
     %last_baseline = [];
end

end

figure(1);
plot([ A(1:length(run_var),:) run_var]);
figure(2);
plot([last_action ones(size(last_action))*diag(last_baseline)]);
disp(['last global sum: ' num2str(last_globalsum) ]);
disp(['peaks:']);
disp(peaks1');
disp(peaks2');

last global sum: 564.0645 402.3871
peaks:
0
-191.5484
282.0323
-90.4839

0
7.4194
6.4194
8.4194
7.4194
136.9032
2.4194
-177.4839
80.7419

We then get the parameters by applying a ratio threshold T / lambda where lambda is for instance 10:

global area in X axis: 564.0645

global area in Y axis: 402.3871

peaks in X axis:

negative, 191.5484
positive, 282.0323
negative, 90.4839

peaks in Y axis:

positive, 136.9032
negative, 177.4839
positive, 80.7419

Document created by Kristof Van Laerhoven, May 2003.