I'm getting weird units while double-integrating accelerometer data to extract displacement

I'm currently working on a research project, where we are trying to measure how much test subjects move their head while looking at pictures under different conditions. We are using an accelerometer attached to their head on a hat, and the head movements are meant to be very slight (only moving about an inch side-to-side). We recorded a few test trials just to see what the data looks like, and this is the script I wrote to get average displacement out of it:

clc, clearvars

hat_on_table = readtable('hat_on_table.csv');
looking_at_screen = readtable('looking_at_screen.csv');
talking = readtable('talking.csv');
walking = readtable('walking.csv');

avg_displacements = {func(hat_on_table), func(looking_at_screen),func(talking), func(walking)}
column_names = {'hat_on_table', 'looking_at_screen', 'talking', 'walking'};

data_table = cell2table(avg_displacements, 'VariableNames', column_names, 'RowNames', {'average displacement (m)'});

function avg_displacement = func(curr)
    time = curr.('epoc_ms_');
    time = time - time(1);
    % convert from ms to s
    time = time / 1000;
    
    acc_x = curr.('x_axis_g_');
    acc_y = curr.('y_axis_g_');
    acc_z = curr.('z_axis_g_');
    % convert from Gs to m/s^2
    acc_x = acc_x*9.80665;
    acc_y = acc_y*9.80665;
    acc_z = acc_z*9.80665;
    
    % find displacement for each axis
    
    % normalize the data
    acc_x = acc_x - mean(acc_x);
    time_and_acc_x = table(time, acc_x, 'VariableNames', {'time (s)', 'x-axis acceleration m/s^2'});
    % integrate first for velocity (m/s), then again for displacement (m)
    velocity_x = cumtrapz(time_and_acc_x.("time (s)"), time_and_acc_x.('x-axis acceleration m/s^2'));
    displacement_x = cumtrapz(time_and_acc_x.("time (s)"), velocity_x);

    % normalize the data
    acc_y = acc_y - mean(acc_y);
    time_and_acc_y = table(time, acc_y, 'VariableNames', {'time (s)', 'y-axis acceleration m/s^2'})
    % integrate first for velocity (m/s), then again for displacement (m)
    velocity_y = cumtrapz(time_and_acc_y.("time (s)"),time_and_acc_y.('y-axis acceleration m/s^2'));
    displacement_y = cumtrapz(time_and_acc_y.("time (s)"), velocity_y);

    % normalize the data
    acc_z = acc_z - mean(acc_z);
    time_and_acc_z = table(time, acc_z, 'VariableNames', {'time (s)', 'z-axis acceleration m/s^2'});
    % integrate first for velocity (m/s), then again for displacement (m)
    velocity_z = cumtrapz(time_and_acc_z.("time (s)"), time_and_acc_z.('z-axis acceleration m/s^2'));
    displacement_z = cumtrapz(time_and_acc_z.("time (s)"), velocity_z);
    

    % combine the three axes into one displacement value
    total_displacements = sqrt((displacement_x.^2)+(displacement_y.^2)+(displacement_z.^2));

    avg_displacement = mean(total_displacements)
end

The issue is that the numbers which we get for average displacements do not make sense for the conditions under which we recorded. Just leaving the hat completely still on the table translated to an average displacement of 3.9m, and looking at a laptop screen resulted in an average displacement of 87.7m. Are the units wrong? Or have I messed up somewhere while integrating/combining the data? Does it need to be filtered first? (Due to the small nature in the desired behavior, I thought filtering could lose the meaningful differences).

Any direction would be much appreciated!

[James Tursa]

If you leave the hat on a table, do the raw accel values make sense? I.e., do you measure pretty much constant reaction force to gravity and nothing else? And do the units make sense? Do you also have gyro data, or only accel data?

Answer 1

[Paul]

Hi @Nathaniel,

I only see one issue with the code, but I have some other thoughts for consideration.

Read in the data

warning('off')
hat_on_table = readtable('hat_on_table.csv');
looking_at_screen = readtable('looking_at_screen.csv');
talking = readtable('talking.csv');
walking = readtable('walking.csv');
warning('on')

Consider the hat on the table. Plot the measured acceleration.

temp = hat_on_table;
ameas = [temp.x_axis_g_ , temp.y_axis_g_ , temp.z_axis_g_];
time = (temp.epoc_ms_ - temp.epoc_ms_(1))/1000; 
figure
plot(time,ameas),xlabel('time');ylabel('g')

Recall that accelerometers don't measure acceleration relative to an inertial reeference frame. Rather, an accelerometer measures specific force, which is the difference between inertial acceleration and gravity. In other words, the ideal accelerometer would measur

ameas = atrue - g

where atrue is the component of the inertial acceleration along the sensing axis and g is the component of gravitational acceleration along the sensing axis.

Assume that the accelerometer axes are orthogonal. Assume that the z-axis of the accelerometer is orthogonal to the local level. The the gravity vector in the along the x-, y-, and z-axes of the accelerometer package is

g = [0 0 -1];

Plot the inertial acceleration, atrue, assuming ideal accelerometers

figure
plot(time,ameas + g),xlabel('Time');ylabel('g')

We see that the there is still some residual error left over. Accelerometers can have a variety of errors. We'll assume the simplest error model of a constant bias.

ameas = (atrue - g) + bias

With this assumption, we'll just estimate the bias as the mean of the residual errors, assuming acceleration due to earth rotation is negligible.

bias = mean(ameas - ([0 0 0] - g))

bias = 1×3

0.1148 -0.0689 -0.0019

Now, reconstruct atrue and plot it

atrue = ameas - bias + g;
figure
plot(time,atrue),xlabel('time');ylabel('g')
ylim([-0.01 0.01])

We see that there is quantization error that will also have an effect.

Compute the position via double integration

pos = cumtrapz(time,cumtrapz(time,9.80665*atrue,1),1);
figure
plot(time,pos),xlabel('Time');ylabel('m')

Even though that hat is not accelerating (not counting earth rotation effects), we see that our model still has residual errors (quantization, scale factor, nonorthogonality, etc.). The spec sheet for your accelerometers should give some insight.

In your code, you take the norm of the position as a function of time, and then take the mean of that result. I'm not sure what that quantity called "average displacement" is supposed to represent.

Look at walking.

temp = walking;
ameas = [temp.x_axis_g_ , temp.y_axis_g_ , temp.z_axis_g_];
time = (temp.epoc_ms_ - temp.epoc_ms_(1))/1000;
figure
plot(time,ameas),xlabel('Time');ylabel('g')

Here, it looks like the walking motion starts at around t = 4 seconds. But for t < 4, the measurements are not the same as for the hat on the table. Is that because the measurement errors are different, or because the z-axis of the hat is no longer orthogonal to the local level (either due to how it's placed on the head or maybe the subject's head is sligthly tilted or both)?

Let's assume that the bias errors are the same as for the hat on the table. Then we can use the t < 4 measurements to estimate the gravity vector in the "tilted" frame of reference, which will be the frame in which the position will be estimated

index = time <= 3.5;
g = mean(bias - ameas(index,:))

g = 1×3

0.5572 -0.1065 -0.9167

The norm of the estimated g vector is close to unity, which is good. We'll normalize it so it makes physical sense.

norm(g)

ans = 1.0780

g = g/norm(g);

Note that we could use the this gravity vectory to estimate the initial orientation of the accelerometer package relative the local level and then we could take the measurements back to the local level frame from the "tilted frame."

Compute true accelerration and plot

atrue = ameas - bias + g;
figure
plot(time,atrue),xlabel('Time');ylabel('g')

That all three measurements aren't zero mean for t < 4 suggest that the error model is probably not quite right.

Position

pos = cumtrapz(time,cumtrapz(time,9.80665*atrue,1),1);
figure
plot(time,pos),xlabel('Time');ylabel('m')

Of course in this case, the pos includes the effect of the subject walking, presumably along the x direction.. Not sure how you want to deal with that, but taking the time-average of the norm of pos probably isn't what's desired. Mabye a high pass filter on the accelerometer data would be appropriate? I assume the goal is to get the incremental position of the acclerometer relative to some reference point (like the center of the torso or something).

Also, while the subject is walking we're assuming thet the subject's head has zero angular velocity and that the accelerometer measurements are taken in a nonrotating frame (again, neglecting earth rotation). Absent gyro measuremements, not sure how to deal with that effect other than to make that assumption.

Consider all of the assumptions made just to get to this point.

@William Rose might have insights into this problem. IIRC, he's done work in the past that sounds similar to this.

avg_displacements = {func(hat_on_table), func(looking_at_screen),func(talking), func(walking)};
column_names = {'hat_on_table', 'looking_at_screen', 'talking', 'walking'};
data_table = cell2table(avg_displacements, 'VariableNames', column_names, 'RowNames', {'average displacement (m)'})

data_table = 1×4 table
                                hat_on_table    looking_at_screen    talking    walking
                                ____________    _________________    _______    _______

    average displacement (m)       3.9492            87.706          278.58     67.193 

function avg_displacement = func(curr)
    time = curr.('epoc_ms_');
    time = time - time(1);
    % convert from ms to s
    time = time / 1000;
    
    acc_x = curr.('x_axis_g_');
    acc_y = curr.('y_axis_g_');
    acc_z = curr.('z_axis_g_');
    % convert from Gs to m/s^2
    acc_x = acc_x*9.80665;
    acc_y = acc_y*9.80665;
    acc_z = acc_z*9.80665;
    
    % find displacement for each axis
    
    % normalize the data
    acc_x = acc_x - mean(acc_x);
    time_and_acc_x = table(time, acc_x, 'VariableNames', {'time (s)', 'x-axis acceleration m/s^2'});
    % integrate first for velocity (m/s), then again for displacement (m)
    velocity_x = cumtrapz(time_and_acc_x.("time (s)"), time_and_acc_x.('x-axis acceleration m/s^2'));
    displacement_x = cumtrapz(time_and_acc_x.("time (s)"), velocity_x);
    % normalize the data

    acc_y = acc_y - mean(acc_y);
    time_and_acc_y = table(time, acc_y, 'VariableNames', {'time (s)', 'y-axis acceleration m/s^2'});
    % integrate first for velocity (m/s), then again for displacement (m)
    velocity_y = cumtrapz(time_and_acc_y.("time (s)"),time_and_acc_y.('y-axis acceleration m/s^2'));
    displacement_y = cumtrapz(time_and_acc_y.("time (s)"), velocity_y);
    % normalize the data
    acc_z = acc_z - mean(acc_z);
    time_and_acc_z = table(time, acc_z, 'VariableNames', {'time (s)', 'z-axis acceleration m/s^2'});
    % integrate first for velocity (m/s), then again for displacement (m)
    velocity_z = cumtrapz(time_and_acc_z.("time (s)"), time_and_acc_z.('z-axis acceleration m/s^2'));
    displacement_z = cumtrapz(time_and_acc_z.("time (s)"), velocity_z);
    
    % combine the three axes into one displacement value
    total_displacements = sqrt((displacement_x.^2)+(displacement_y.^2)+(displacement_z.^2));
    avg_displacement = mean(total_displacements);
end