MATLAB Examples

# Human Activity Classification based on Smartphone Sensor Signals

Part1 - (Post) Data Visualization & feature extraction

## 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 Mathworks 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:

• One to do with "fast" variations over time, due to body dynamics (physical movements of the subject)
• The other, responsible for "slow" variations over time, due to the position of the body with respect to the vertical gravitational field

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:

• Design a high-pass filter (e.g. using the Filter Design and Analysis Tool, fdatool, in MATLAB)
• Apply the filter to the original signal

## 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 MathWorks, 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 MathWorks, 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:

• Re-organise the acceleration signals into shorter buffers of fixed length L, each labeled with a single activity ID
• Extract a vector of features for each Lx3 signal buffer [ax, ay, az] using the function featuresFromBuffer
• Provide the model with two sets of feature vectors and corresponding activity ID

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

• Read the buffered signals from BufferedAccelerations.mat
• Compute a feature vector for each buffer using featuresFromBuffer
• Save all feature vectors into the file BufferFeatures.mat

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 ```