Human Activity Classification based on Smartphone Sensor Signals

Part1 - (Post) Data Visualization & feature extraction

Contents

Introduction

This example describes an analysis approach on accelerometer signals captured with a smartphone. The smartphone is worn by a subject during 6 different types of physical activity. The goal of the analysis is to build an algorithm that automatically identifies the activity type given the sensor measurements.

The example uses data from a recorded dataset, courtesy of the bat365 France employees (N=6 subjects to date + 1 dummy subject)

n_subjects = 7; % Automate later

Open full recording for a subject & look into the data

Retrieve a single acceleration component for a particular subject.

[acc, actid, actlabels, t, fs] = getRawAcceleration('SubjectID',3,...
    'Component','x','DataFolder','Data\Prepared_iPhone_32');

% Let's look at the same data, colored based on the activity type.
% Given this data:
%
% * We would like to be able to tell the difference between the different
%   activities, just based on the content of the signal.
% * Note in this case the coloring is based on existing knowledge (actid)
% * Labeled data can be used to "train" a classification algorithm so it
%   can it later predict the class of new (unlabeled) data.

% Visualize the same signal using a custom plotting function
plotAccelerationColouredByActivity(t, acc, actid, {'Vertical acceleration'},actlabels)

Amplitude only - Using a mean measure

Looking only at the raw values irrespective or their position in time can provide a first set of clues

% This can easily help separate e.g. Walking from Sitting
figure
plotCompareHistForActivities(acc, actid,'Sitting', 'Standing',actlabels)

Amplitude only - Using an RMS or standard deviation measure

This can help separate things like e.g. Walking and Standing

plotCompareHistForActivities(acc, actid,'Standing', 'Walking',actlabels)

Amplitude-only methods are often not enough... What next?

For example it would be very hard to distinguish between simply Walking and Running (close statistical moments)

Simple statistical analysis is often not sufficient. For signals one should also consider methods that measure signal variations over time

plotCompareHistForActivities(acc, actid,'Walking', 'Running',actlabels)

Time-domain analysis - preliminary considerations

There two main different types of causes behind our signals:

As we focus on classifying the physical activities, we should focus time-domain analysis on the effects of body movements. These are responsible for the most "rapid" (or frequent) variations in our signal.

In this specific case a simple average over a period of time would easily estimate the gravitational component, which could be then subtracted from the relevant samples to obtain the signal due to the physical movements.

For the sake of generality here we introduce an approach based on digital filters, which is much more general and can be reused in more challenging situations.

Digital filtering workflow

To isolate the rapid signal variations from the slower ones using digital filtering:

Filter out gravitational acceleration

As well as interactively, filters can be designed programmatically. In this case the function hpfilter was generated automatically from the Filter Design and Analysis Tool, but it could have just as well been created manually

% Design filter
fhp = hpfilter;

% "Decouple" acceleration due to body dynamics from gravity
ab = filter(fhp,acc);

% Compare the filtered signal with the original one
plotAccelerationColouredByActivity(t, [acc, ab], actid,...
    {'Original','High-pass filtered'},actlabels)

Focus on stationary activities

After filtering signals artefact were detected in stationnary signals (1 & 2)

% Assume we know the activity id for Walking is 3
stationnary = [1 2];
sel = ismember(actid,stationnary);

% Select only desired array segments for time vector
% and acceleration signal
ts = t(sel);
abs = ab(sel);

% Chop off 1st & last seconds
ts = ts(fs:end-fs);
abs = abs(fs:end-fs);

% Plot stationnary-only signal segment. Use interactive plot tools to zoom in
% and explore the signal
figure
subplot(211);
plot(ts, abs);title('Stationnary activities');
axis tight

% Moving-Average filtering (Low-Pass filtering)
abs_f = smooth(abs,fs);
subplot(212);
plot(ts, abs_f);title('Low-pass filtered');
axis tight

Plot Power Spetral Density & Validate activities differentiation

Use Welch method with its default options, using known sample frequency

% % When called without output arguments, the function pwelch plots the PSD
[~, f] = pwelch(abs,512,256,2048,fs);
f = f(f<5);

% E.g. are position and height in the two respsctive sets of peaks
% different for different activities?
P = zeros(1024,max(actid));
for k = 1:max(actid)
    p = getActivityPSD_New(ab,actid,k,fs);
    P(:,k) = p;
end

% % Re-Scale PSD to fit on same scale (visual trick)
% P_scaled = P./repmat(sqrt(sum(P.*P)),size(P,1),1);
ribbon(f,P); ylabel('Frequency (Hz)'); xlabel('Activity Type')
set(gca,'Xtick',1:max(actid),'XTickLabel',actlabels); axis tight

Validate consistency of PSD information across different subjects

Compare Power Spectral Density of signals for walking activity across all subjects in the dataset.

% This helper function uses a linear amplitude scale so PSD peaks
% visually stand out better

for i=3
    plotPSDForGivenActivity(i,n_subjects);
    axis tight
end

Automatic peak identification (Walking activity example)

The function findpeaks in Signal Processing Toolbox can be used to identify amplitude and location of spectral peaks

walking = 3;
abw = acc(actid == walking);

% Compute numerical values of PSD and frequency vector. When called with
% one or two output arguments, the function pwelch does not automatically
% plot the PSD
[p,f] = pwelch(abw,[],[],[],fs);

% Use findpeaks to identify values (amplitude) and index locations of peaks
[pks,locs] = findpeaks(p);

% Plot PSD manually and overlay peaks
figure
plot(f,db(p),'.-')
grid on
hold on
plot(f(locs),db(pks),'rs')
hold off
addActivityLegend(walking)
title('Power Spectral Density with Peaks Estimates')
xlabel('Frequency (Hz)')
ylabel('Power/Frequency (dB/Hz)')

Peak finding with more specific requirements

PSD with numerical output

[p,f] = pwelch(abw,[],[],[],fs);

% Find max 8 peaks, at least 0.25Hz apart from each other and with a given
% prominence value

fmindist = 0.25;                    % Minimum distance in Hz
N = 2*(length(f)-1);                % Number of FFT points
minpkdist = floor(fmindist/(fs/N)); % Minimum number of frequency bins

[pks,locs] = findpeaks(p,'npeaks',10,'minpeakdistance',minpkdist,...
    'minpeakprominence', 0.15);

% Plot PSD and overlay peaks
plot(f,db(p),'.-')
grid on
hold on
plot(f(locs),db(pks),'rs')
hold off
legend(actlabels(walking));
title('Power Spectral Density with Peaks Estimates')
xlabel('Frequency (Hz)')
ylabel('Power/Frequency (dB/Hz)')

Autocorrelation as a feature

Autocorrelation can also be powerful for frequency estimation. It is especially effective for estimating low-pitch fundamental frequencies

% xcorr with only one input will compute the autocorrelation
[c, lags] = xcorr(abw);

% Highlight the main t=0 peak (overal energy) and a few secondary peaks
% The position of the second highest peaks identifies the main period
tmindist = 0.3;
minpkdist = floor(tmindist/(1/fs));
[pks,locs] = findpeaks(c,'minpeakprominence',1e4,...
    'minpeakdistance',minpkdist);

% Plot autocorrelation and three key peaks
tc = (1/fs)*lags;
figure
plot(tc,c,'.-')
grid on
hold on
plot(tc(locs),c(locs),'rs')
hold off
xlim([-5,5])
legend(actlabels(walking));
title('Autocorrrelation with Peaks Estimates')
xlabel('Time lag (s)')
ylabel('Correlation')

Check position of first peak varies between different activities

Compare autocorrelation plots for same subject and different activity By zooming into the secong highest peaks, note how the respective second-peak positions are still spaced apart by more than a few samples.

plotCorrActivityComparisonForSubject(1, 'Walking', 'Running')

Feature extraction summary

After exploring interactively a few different techniques to extract descriptive features from this type of signals, we can collect all the analysis methods identified into a single function. The responsibility of this function is to extract a fixed number of features for each signal buffer provided as input.

type featuresFromBuffer
type featuresFromBuffer_codegen

function feat = featuresFromBuffer(at, fs)
% featuresFromBuffer Extract vector of features from raw data buffer
% 
% Copyright 2014-2015 The bat365, Inc.

% Initialize digital filter
persistent fhp
if(isempty(fhp))
    fhp = hpfilter;
    fhp.PersistentMemory = false;
end

% Initialize feature vector
feat = zeros(1,60);

% Remove gravitational contributions with digital filter
ab = filter(fhp,at);

% Average value in signal buffer for all three acceleration components (1 each)
feat(1:3) = mean(at,1);

% RMS value in signal buffer for all three acceleration components (1 each)
feat(4:6) = rms(ab,1);

% Autocorrelation features for all three acceleration components (3 each): 
% height of main peak; height and position of second peak
feat(7:15) = autocorrFeatures(ab, fs);

% Spectral peak features (12 each): height and position of first 6 peaks
feat(16:51) = spectralPeaksFeatures(ab, fs);

% Spectral power features (5 each): total power in 5 adjacent
% and pre-defined frequency bands
feat(52:60) = spectralPowerFeatures(ab, fs);

% --- Helper functions
function feats = autocorrFeatures(x, fs)
n_channels = size(x,2);
feats = zeros(1,3*n_channels);

[c,lags] =arrayfun(@(i) xcorr(x(:,i)),1:n_channels,'UniformOutput',false);

minprom = 0.0005;
mindist_xunits = 0.3;
minpkdist = floor(mindist_xunits/(1/fs));

% Separate peak analysis for all channels
for k = 1:n_channels
    [pks,locs] = findpeaks(c{k},...
        'minpeakprominence',minprom,...
        'minpeakdistance',minpkdist);
    
    tc = (1/fs)*lags{k};
    tcl = tc(locs);
    
    % Feature 1 - peak height at 0
    if(~isempty(tcl))   % else f1 already 0
        feats(n_channels*(k-1)+1) = pks((end+1)/2);
    end
    % Features 2 and 3 - position and height of first peak 
    if(length(tcl) >= 2)   % else f2,f3 already 0
        feats(n_channels*(k-1)+2) = tcl(2);
        feats(n_channels*(k-1)+3) = pks(2);
    end
end

function feats = spectralPeaksFeatures(x, fs)
n_channels = size(x,2);
mindist_xunits = 0.3;

feats = zeros(1,12*n_channels);

N = 4096;
minpkdist = floor(mindist_xunits/(fs/N));

% Cycle through number of channels
for k = 1:n_channels
    [p, f] = periodogram(x(:,k),rectwin(length(x)),4096,fs);
    [pks,locs] = findpeaks(p,'npeaks',20,'minpeakdistance',minpkdist);
    if(~isempty(pks))
        mx = min(6,length(pks));
        [spks, idx] = sort(pks,'descend');
        slocs = locs(idx);
        
        pks = spks(1:mx);
        locs = slocs(1:mx);
        
        [slocs, idx] = sort(locs,'ascend');
        spks = pks(idx);
        opks = spks;
        locs = slocs;
    end
    ofpk = f(locs);
    
    % Features 1-6 positions of highest 6 peaks
    feats(12*(k-1)+(1:length(opks))) = ofpk;
    
    % Features 7-12 power levels of highest 6 peaks
    feats(12*(k-1)+(7:7+length(opks)-1)) = opks;
end

function feats = spectralPowerFeatures(x, fs)
n_channels = size(x,2);

edges = [0.5, 1.5, 5, 10];
n_feats = length(edges)-1;
feats = zeros(1,n_feats*n_channels);

for k=1:n_channels
    [p, f] = periodogram(x(:,k),rectwin(length(x)),4096,fs);
    for kband = 1:n_feats
        feats(n_feats*(k-1)+kband) = sum(p( (f>=edges(kband)) & (f<edges(kband+1)) ));
    end
end

function feat = featuresFromBuffer_codegen(at, fs)
% featuresFromBuffer Extract vector of features from raw data buffer
% 
% Copyright 2014-2015 The bat365, Inc.

% Initialize digital filter
persistent dcblock corr spect f
if(isempty(dcblock))
    [s,g] = getFilterCoefficients(fs);
    dcblock = dsp.BiquadFilter('Structure','Direct form II transposed', ...
        'SOSMatrix',s,'ScaleValues',g);

    NFFT = 4096;
    spect = dsp.SpectrumEstimator('SpectralAverages',1,...
        'Window','Rectangular','FrequencyRange','onesided',...
        'SampleRate',fs,'SpectrumType','Power density',...
        'FFTLengthSource','Property','FFTLength',4096);
    f = (fs/NFFT)*(0:NFFT/2)';

    corr = dsp.Autocorrelator;
end

% Initialize feature vector
feat = zeros(1,60);

% Remove gravitational contributions with digital filter
ab = step(dcblock,at);

% Average value in signal buffer for all three acceleration components (1 each)
feat(1:3) = mean(at,1);

% RMS value in signal buffer for all three acceleration components (1 each)
feat(4:6) = rms(ab,1);

% Autocorrelation features for all three acceleration components (3 each): 
% height of main peak; height and position of second peak
feat(7:15) = autocorrFeatures(ab, corr, fs);

% Pre-compute spectra of 3 channels for frequency-domain features
af = step(spect,ab);

% Spectral peak features (12 per channel): value and freq of first 6 peaks
feat(16:51) = spectralPeaksFeatures(af, f);

% Spectral power features (3 each): total power in 3 adjacent
% and pre-defined frequency bands
feat(52:60) = spectralPowerFeatures(af, f);

% --- Helper functions
function feats = autocorrFeatures(x, corr, fs)
n_channels = size(x,2);
feats = zeros(1,3*n_channels);

c = step(corr, x);
lags = (0:length(x)-1)';

minprom = 0.0005;
mindist_xunits = 0.3;
minpkdist = floor(mindist_xunits/(1/fs));

% Separate peak analysis for all channels
for k = 1:n_channels
    [pks,locs] = findpeaks([0;c(:,k)],...
        'minpeakprominence',minprom,...
        'minpeakdistance',minpkdist);

    tc = (1/fs)*lags;
    tcl = tc(locs-1);

    % Feature 1 - peak height at 0
    feats(n_channels*(k-1)+1) = c(1,k);
    % Features 2 and 3 - position and height of first peak 
    if(length(tcl) >= 2)   % else f2,f3 already 0
        feats(n_channels*(k-1)+2) = tcl(2);
        feats(n_channels*(k-1)+3) = pks(2);
    end
end

function feats = spectralPeaksFeatures(xpsd, f)
n_channels = size(xpsd,2);
mindist_xunits = 0.3;

feats = zeros(1,12*n_channels);

minpkdist = floor(mindist_xunits/f(2));

% Cycle through number of channels
nfinalpeaks = 6;
for k = 1:n_channels
    [pks,locs] = findpeaks(xpsd(:,k),'npeaks',20,'minpeakdistance',minpkdist);
    opks = zeros(nfinalpeaks,1);
    if(~isempty(pks))
        mx = min(6,length(pks));
        [spks, idx] = sort(pks,'descend');
        slocs = locs(idx);

        pkssel = spks(1:mx);
        locssel = slocs(1:mx);

        [olocs, idx] = sort(locssel,'ascend');
        opks = pkssel(idx);
    end
    ofpk = f(olocs);

    % Features 1-6 positions of highest 6 peaks
    feats(12*(k-1)+(1:length(opks))) = ofpk;
    
    % Features 7-12 power levels of highest 6 peaks
    feats(12*(k-1)+(7:7+length(opks)-1)) = opks;
end

function feats = spectralPowerFeatures(xpsd, f)
n_channels = size(xpsd,2);

edges = [0.5, 1.5, 5, 10];
n_feats = length(edges)-1;

featstmp = zeros(n_feats,n_channels);
    
for kband = 1:length(edges)-1
    featstmp(kband,:) = sum(xpsd( (f>=edges(kband)) & (f<edges(kband+1)), : ),1);
end
feats = featstmp(:);

function [s,g] = getFilterCoefficients(fs)
coder.extrinsic('zp2sos')
[z,p,k] = ellip(7,0.1,60,0.4/(fs/2),'high');
[s,g] = coder.const(@zp2sos,z,p,k);

Parallel data preparation to train Machine-Learning algorithms

To train the model, assume we:

The buffered data is already available and stored in the file BufferedAccelerations.mat

Computing the features is a fairly efficient process, but it takes a while in this case because of the high number fo signal vectors available. The pre-computed set of feature vectors for all available signal buffers is available in the file BufferFeatures.mat To re-compute all features use the function extractAllFeatures, which will

extractAllFeatures can distribute the computations to a pool of workers if Parallel Computing Toolbox is installed

p = gcp('nocreate');
if isempty(p)
    p = parpool('local');
end

extractAllFeatures('Data\Prepared_iPhone_32','BufferFeatures60.mat');

Feature extraction starting for 24075 observations... 
...Done!
Total time elapsed: 187.166 seconds

Copyright 2014-2015 The bat365, Inc.
Published with MATLAB® R2015b