//Logo Image
AuthorsChe-Chang Yang (2010-09-25); recommended: Yeh-Liang Hsu (2011-02-11).
Note: This article is the Chapter 5 of Che-Chang Yang’s doctoral thesis “Development of a Home Telehealth System for Telemonitoring Physical Activity and Mobility of the Elderly”.

Chapter 5. Estimation of energy expenditure using a wearable accelerometry systems

This chapter presents a method of estimating the energy expenditure (EE) by using a wearable accelerometry system which includes a tri-axial accelerometer and a barometric pressure sensor. In structured ambulation courses, the accelerations (x-, y-, and z-axis) and the atmospheric pressure signals were recorded while the oxygen uptake, a variety of pulmonary parameters, heart rate, and walking speed were measured by means of a portable oxygen analyzer apparatus. A black-box estimation method was developed to estimate EE with a selection of features extracted from the accelerometry data, heart rate and walking speed. The preliminary estimation result is also shown.

5.1 Introduction

Estimation of energy expenditure (EE) is essential in the research of obesity, sports exercise, physical activity, and public health investigation. EE can be accurately assessed by means of the doubly-labeled water (DLW) method, or direct and indirect calorimetry in clinical-settings, which are however costly and cannot be implemented in a larger group study. Traditionally EE assessment tools commonly rely on the uses of dietary record or subjective investigation.

Motion sensors have been used to estimate EE. Mechanical pedometers are the simplest motion sensor for this purpose by detecting the steps taken during ambulation and locomotion. The numbers of steps are used to estimate the travelled distance, as well as the corresponding EE. Though being simple and low cost, those mechanical pedometers are unable to distinguish the intensity of motions and therefore such a major drawback leads to less accurate and compromised EE estimation performance. Accelerometers have been used to enhance the EE estimation accuracy and have become the ideal replacement of traditional pedometers. A vast amount of studies have been using accelerometers (or accelerometry) to measure physical activity, identify daily movement and adverse activities (i.e., the falls), as well as estimate EE.

In addition, the Biomedical System Laboratory (BSL) in the University of New South Wales (UNSW) has developed an wearable accelerometry measurement system for this purpose. Recently the attempt has been made to utilize a barometric pressure sensor in the existing wearable accelerometry measurement system to enhance the accuracy in classifying activities that produce position shift in altitude, such as walking on inclined surfaces or walking up/downstairs. Therefore, the resulting EE is expected to be better estimated due to the added altitude information. The trial study has shown an improvement and positive outcome in EE estimation.

This chapter presents a method of estimating the energy expenditure (EE) by using a wearable accelerometry system which includes a tri-axial accelerometer and a barometric pressure sensor. In structured ambulation courses, the accelerations (x-, y-, and z-axis) and the atmospheric pressure signals were recorded while the oxygen uptake, a variety of pulmonary parameters, heart rate, and walking speed were measured by means of a portable oxygen analyzer apparatus. A black-box estimation method was developed to estimate EE with a selection of features extracted from the accelerometry data, heart rate and walking speed. The preliminary estimation result is also shown.

5.2 Data collection

A “black-box” polynomial estimation model is used to estimate EE. To establish a database for training model parameters, test subjects (2 young healthy male subjects) were recruited in the data collection phase. This section describes the use of the instruments and protocol of the data collection.

5.2.1 Instruments

(1)  Triaxial accelerometer module

The triaxial accelerometry module (TAM. or Triax) is a small (70×54×14mm), light weight (57.5g) wearable device that measures body movements. The Triax uses a MEMS capacitive tri-axis accelerometer (MMA7260, Freescale) to measure accelerations in a 4g sensitivity range (default setting). A barometric pressure sensor (SCP1000, VTI Technologies) is also integrated to provide altitude information. A Class-I Bluetooth transceiver is used to enable real-time s data transmission in radio band. The Triax uses a microcontroller (MSP430F149, Texas Instrument) as a control and processing core. The signals from the tri-axis accelerometer and the barometric pressure sensor are sampled by the 12-bit analog-to-digital conversion (ADC) circuits of the microcontroller at the sampling rates of 80Hz, and 1.8Hz, respectively. A PC with Bluetooth connectivity is required and necessary to receive and log the real-time data transmitted from the Triax.

DSC04756.jpg

Figure 5.1 Triaxial accelerometry module

(2)  Cardio-Pulmonary Exercise Tester (K4b2)

The cardio-pulmonary exercise tester (CPET), (K4b2, COSMED, Italy) is a portable system that measures a variety of pulmonary parameters of human body during exercises and ambulation. The measured rate of oxygen uptake per mass (, in the unit of ) is used as a gold standard estimate of energy expenditure during exercise and metabolism. This device consists of a gas mask, a portable unit and a battery (Figure 5.2). The gas mask collects the breathing gas that is to be directed to the portable unit via a tube. A small turbine with an opto-reader (i.e., a flow metering sensor) is inserted in the mask to measure the volume of the breathing gas. The portable unit analyzes the gas volume (from the turbine signals) and gas components, O2 consumption and CO2 production, breath-by-breath in real-time. A series of pulmonary parameters are also derived. The data are logged in the portable unit. The recorded data can be downloaded to a PC via a RS-232 interconnect.

In addition, a GPS module is connected to the portable unit to measure walking speed. Heart rate (HR) was measured by the Polar T31 transmitter that was worn around the chest of the subjects. HR signals were wirelessly transmitted to a receiver probe connected to the portable unit. During the measurement, the subjects wore a harness to which all the components, except the Polar T31 transmitter, are attached to. For best measurement accuracy, a minimal warm-up time of 45 minutes for the K4b2 is required before system calibration. The K4b2 weighs 2.5kg in total and it does not causes significant burdens to the subjects during data collection. The calibration procedure for K4b2 consists of: (1) Room air calibration, (2) Reference gas calibration, (3) Gas delay calibration, and (4) Turbine calibration. The calibration configuration may vary when the components are replaced, disinfected and even re-connected. Ill-calibrated or non-calibrated system will result in less accurate measurement and analysis. Therefore, calibrations must be carried out before every measurement of single subject, or whenever any activated component is replaced.

DSC04625-6.jpg

Figure 5.2 CPET K4b2

5.2.2 Data collection protocol

The data collection consists of 3 phases for each subject:

(1)  Measuring  at rest state

In this data collection phase, the subjects’  at a 3-min resting stance state was measured individually in the laboratory. The data collected in this phase is used to establish the base line of the metabolism rate at resting state. Note that in this phase the Triax is not used because the subjects always kept still stance state and therefore the accelerometry data and its corresponding features are deemed identical for all the subjects.

(2)  Determining walking cadences

The subjects’ walking cadences (steps/min.) were measured individually before the 3rd data collection phase. The subjects’ cadences were measured by asking them to walk at their won normal paces along a level and 92m straight pathway without obstacles while their steps taken and the time spent were recorded. The cadences at normal walking paces can be obtained, and the waking paces 20% reduced and 20% increased from their normal walking paces were assumed as their slow and fast cadences.

(3)  Collecting  in structured ambulation

This data collection task was conducted in an outdoor environment (walkway between Scientia Building and Red Centre in the UNSW campus), as depicted in Figure 5.3. Wearing both the K4b2 and Triax, the subjects were directed to walk on a level walkway (near Scientia Building) and then to perform the sequential ambulation states across different terrains as listed in the test protocol (Table 5.1). Note that the ambulation in the terrain states 7 to 12 (reverse, uphill) are actually the reverse sequence of the states 1 to 6 (forward, downhill). Between the uphill and downhill ambulation, a minimum break of 2 minutes is given for them to stabilize their metabolism to their resting state. During data collection, the subjects are asked to walk at regulated (slow, normal, and fast) cadences. For each subject, data during their ambulation at slow, normal, and fast walking cadences were recorded.

Before data collection, the locations of terrain transits are shown to the subjects so that they have clear information about the entire test terrain. During the data collection in this phase, an experiment instructor carried a laptop PC to receive Triax signals, and the experiment instructor followed behind the subjects to avoid interfering the subjects. In the meanwhile, the subjects were asked to press the “marker” button on the portable unit to annotate the terrain transits along the entire walkway while the instructor annotates the transit on the Triax data logger on the laptop PC.

  

Figure 5.3 Terrain illustration

Table 5.1 The test protocol in the 3rd phase

No.

Terrain

1

Walk on level (1)

Forward, downhill

2

Walk downstairs: stair (1), and (2)

3

Walk on inclined ramp: slope (1)

4

Walk on inclined ramp: slope (2)

5

Walk on inclined ramp: slope (3)

6

Walk on level (2)

Break

7

Walk on level ground (2)

Reverse, uphill

8

Walk on inclined ramp: slope (3)

9

Walk on inclined ramp: slope (2)

10

Walk on inclined ramp: slope (1)

11

Walk upstairs: stair(2), and (1)

12

Walk on level (1)

5.3 Developing the estimation model

This section describes the development of the EE estimation model. A black-box model is used to provide the EE estimate from the inputs of features. The development tasks include data preparation, feature selection and extraction, and model parameter selection.

5.3.1 Data preparation

For the obtained in the 1st data collection phase from all the subjects, the average rate of oxygen uptake per unit mass,  is calculated and is used as the baseline at resting state. Because the subjects kept standing still in this phase, accelerometry data are assumed identical for all the subjects.

Although the data from the Triax and K4b2 were synchronously and continuously recorded, each ambulation state can be identified and separated according to the annotated markers in the data. As a result, the data recorded between each adjoining terrain state (as the segments of dot lines shown along the terrain path in Figure 5.3) is excluded and then the remaining data (referring to Table 5.1) can be concatenated sequentially.

5.3.2 Feature selection and extraction

The x-, y-, and z-axis output signals were media-filtered (window length n=3) to remove artifact signal spikes before the analog-to-digital conversion (ADC) that converts the signals into accelerations (ax, ay, and az in the unit of g=9.81 ) and tilts (tx, ty, and tz in the unit of degree). Because the MMA7260QT tri-axis accelerometer is gravity-responsive, low pass filtering (3rd-order elliptical IIR, fs=0.25Hz) is applied to the accelerations to separate gravity component (, , and ). Body accelerations (,, and ) can be estimated by subtracting the gravity components from the acceleration outputs (ax, ay, az). Acceleration features were extracted using sliding window with a constant window length l=200 samples (equal to 2.5s at sampling rate fs=80Hz) and 50% overlapped in the window length.

Several signal features from the collected data were identified as follows:

(1)  Signal magnitude area (SMA)

SMA is defined as the sum of integrals of the signal magnitudes of the three body acceleration.

                                                       (5.1)

(2)  Vertical acceleration energy (VE)

VE is defined as the sum of the square of vertical acceleration (x-axis)

                                                                                            (5.2)

(3)  Peak frequency (PF)

Peak frequency is derived by applying Fast Fourier transfer to the vertical acceleration.

                                                                            (5.3)

(4)  Peak to peak magnitudes (P2Px, P2Py, P2Pz)

    

                                                     (5.4)

(5)  Trunk tilt angles (TAx, TAy, TAz)

                                                                                  (5.5)

(6)  Acceleration pair correlations (Cxy, Cyz, Cxz)

The coefficients of correlation between x-, and y-axis, y- and z-axis, and x- and z-axis were calculated

                                                                                      

                                                                                       

                                                                                (5.6)

(7)  Differential pressure (DP)

Although the barometric pressure data is recorded synchronously with accelerations that are sampled at 80Hz, the actual sampling rate of barometric pressure is approximately 1.8Hz, which is far slower than that of accelerations. Therefore, the barometric pressure data was firstly median-filtered (n=29) and then up-sampled to match the sampling rate of 80Hz by linear interpolation.

Because the response of the barometric pressure sensor is not fast and its output is liable to external atmospheric interference, the pressure differences is obtained from the interpolated pressure data by calculating the differences between two neighboring 8-sec averages before and after each pressure sample.  

                                                     (5.7)

The differential pressure DP is calculated from the smoothened barometric pressure:

                                                                               (5.8)

5.3.3 Oxygen uptake data and the additional features

During the data collection, the additional data, heart rate (HR, bpm) and walking speed (SP, m/s2) is also recorded simultaneously with the rate of oxygen uptake. Note that the sampling intervals of rate of oxygen uptake, heart rate and walking speed measured by K4b2 are not uniform and depends on the varying breath-by-breath intervals of the subjects. The sampling intervals of K4b2 are longer than that of the Triax. Therefore linear interpolation (upsampling to approximate fs=80Hz) is applied to these data so as to match the sampling rate of the Triax signals.

5.3.4 Estimation model

All the extracted feature sequences compiled as in the previous section were used to estimate the rate of volume of oxygen uptake by using the output-error model (OE). The OE model structure can be expressed as Equation (5.9) and Figure 5.3 shows the model structure.

                                                                    (5.9)

Figure 5.3 The OE model structure

The polynomial B(q) and F(q) are:

                             (5.10)

                                                                (5.11)

The parameters nb and nf are the orders of the OE model and nk is the delay. When [nb nf nk] is specified, the coefficients [] and [] of the polynomial B(q) and F(q) can be determined in by using prediction error minimization method (PEM) conducted in MATLAB System Identification ToolBox. Each feature was separately estimated with a selection of order range: nb=1 to 20; nf=1 to 20; and nk=0, 5, 9. As a result, 21×21×3=1323 iterations was calculated for each feature vector. The order set that has minimal mean square error (MSE, in the unit of ) was selected.

                                                                   (5.12)

 

                                                               (5.13)

5.3.5 Data preparation for the training phase

After completing feature extraction, data preparation for each feature vectors is required for Equation (8) and (9) show the data arrangements for K4b2 and Triax data ( and , respectively) for signal estimation processing. In Equation (x),  is the average rate of oxygen uptake measured during 3-minute still standing of all the subjects. ,  and  are  the rate of oxygen uptakes during slow, normal and fast walking speed of the sth subject, respectively. The Triax data of all the subjects was arranged as the Equation (x), in which ,  and  are the features generated from Triax signals during slow, normal and fast walking speed of the sth subject. is the feature generated from Triax signals during standing still. Note that  is not the average of the features of all the subjects because there is no apparent difference between Triax signals measured from subjects during still standing. Therefore  is a assumed identical here. Also the feature of peak frequency (PF) and differential pressure (DP) are assumed zero.

In each walking cadence, the data collected in forward walking (A to B) and reverse walking (B to A) is then concatenated. Finally the concatenated data at slow, normal, and fast walking cadences is arranged sequentially. Triax data is automatically processed in MATLAB and the gas, heart rate, and speed data from K4b2 is manually classified.

                                        (5.14)

                                          (5.15)

Because the sliding windowing technique of the window length l=200 samples with 50% overlap is applied to the acceleration and barometric pressure sampled/up-sampled at 80Hz, the sampling rate of the features therefore result in approximately 0.8Hz (equals to sampling interval Ts=1.25s) . However, the sampling rate of  in K4b2 is not uniform and depends on breath-by-breath intervals of the subjects that vary individually. Similarly, the heart rate, speed data and other pulmonary parameters are recorded simultaneously.

5.4 Estimation results

As shown in the Table 5.2, there were 18 features extracted from the Triax (and the heart rate and walking speed).  All the features were used to estimate the  separately. In order to find the order (nb, nf, nk) that best approximates the , an iterative procedure was conducted to approximate  by using a selected order combination of nb, nf, nk (giving the ranges of nb=1 to 20, nf=1 to 20, and nk=1, 5, 9). Therefore, there are 20×20×3=1200 iterations for each feature approximation. The mean squared error (MSE) of  and was used as the performance metric. The orders that generated minimal were selected as the best order. Table 5.2 lists the best orders for all the features, and the corresponding MSE and the coefficient of correlations.

Table 5.2 Feature estimation results

Feature

Order [nb nf nk]

MSE

r

SMA

[ 10 17 5]

30.5022

0.5131

PF

[14 10 1]

27.6993

0.6219

DP

[20 11 5]

16.0892

0.7813

VE

[3 8 9]

28.3137

0.5557

P2Px

[11 19 5]

26.9749

0.5936

P2Py

[18 6 9]

30.7391

0.5410

P2Pz

[19 10 9]

27.6973

0.5819

HR

[4 11 1]

11.829

0.8445

SP

[15 11 5]

24.4691

0.6600

TAx

[16 2 1]

28.8679

0.5489

TAy

[2 12 5]

23.5378

0.6634

TAz

[5 10 1]

22.2780

0.6767

Cxy

[18 18 1]

42.4292

0.3347

Cxz

[3 17 9]

44.3567

0.4132

Cyz

[2 13 1]

28.6351

0.5576

AVx

[10 8 1]

28.3839

0.5581

AVy

[20 9 9]

31.1053

0.4928

AVz

[4 19 5]

30.8676

0.5823

All the features provide different performance in estimating . As a result, the estimation performance of the model may vary with different selection of the feature combinations. Figure 5.4 shows the overall estimation performance of the model by using the features of peak frequency (PF), differential pressure (DP), heart rate (HR), walking speed (SP), y-axis tilt angle (TAy), z-axis tilt angle (TAz), and x-axis peak to peak magnitude (P2Px). The over best order for the 7 selected features is [8,11,1] and it corresponds to the MSE of 7.937 and the coefficient of correlation of 0.9024.

best fit

Figure 5.4 The estimated total (blue) and the measured (red)

Appendix 1. Quick guide to MMA7260QT tri-axis accelerometer

A 5.1 Key features

The MMA7260QT (Freescale Semiconductor) is a micro-machined, capacitive low-g triaxial accelerometer. It uses a small package (6×6×1.45mm 16-lead QFN) and can be operated with low power voltage (2.2V-3.6V) and low current consumption(500uA; Sleep mode: 3uA). It offers selectable sensing sensitivity of ±1.5g, 2g, 4g, and 6g with integral low-pass filter. These superior features make MMA7260QT applicable to portable electronics, navigation or robotics. Figure A5.1 show the QFN package and its pin connections are shown in Figure A5.2. Although a 16-lead package is used, only 8 connections are required, including VDD, VSS (GND), x-, y-, z-axis outputs, Sleep Mode and g-Select1 and 2. The other 8 pins that noted as “N/C” can be left float or unconnected.

Figure A5.1 MMA 7260QT tri-axis accelerometer in QFN package

 

Figure A5.2 Pin connections of MMA7260QT

A 5.2 Converting accelerations

Figure A5.3 shows the axis orientation of the MMA7260QT. The positive signs along x-, y-, and z-axis (with arrows indicated) define the direction that the sensor is accelerated to. The outputs from the MMA7260QT are analog signals with maximal bandwidth response of 350Hz (x- and y-axis) and 150Hz (z-axis). For any axis with no applied acceleration, its output is equal to half the supply voltage (VDD). The output voltage increases from the half VDD level when the sensor is accelerated in the positive direction along its sensitive axis. On the contrary, the signal output is below the half VDD level when the sensor is accelerated in negative direction (or decelerated) along its sensitive axis.

Figure A5.3 Axis orientation of the MMA7260QT

For a typical VDD=3.3V application, the zero-acceleration output is 0.5×3.3=1.65V. When the sensor is accelerated, the outputs of the sensitive axes deviate from 1.65V and the variation is according to the selected sensitivity S (mV/g, voltage per gravity) as shown in Table A5.1. For example, if 2g sensitivity is selected, its sensitivity is 600mV/g (g is gravity in the amount of 9.81m/s2) and the voltage within the sensitivity range changes linearly with the measured acceleration (Acc). Equation (A5.1) explains the relationship. If an axis output is 2.35V, the actual measured acceleration along that axis is Acc=. The Equation (A5.1) is valid for x-, y-, and z-axis outputs.

Sensitivity can be selected with 2 logic inputs connected to pin g-Select1 and g-Select2 (Referring to Figure A5.2). The sensitivity can be changed at anytime during operation. The g-Select pins of the MMA7260QT can be configured with high (1) or low (0) status by microcontroller outputs, as shown in the Table 1. The g-Select pins can be left unconnected for applications only requiring a 1.5g sensitivity.

                                                                 (A5.1)

Table A5.1 Sensitivity selector

g-Range

Sensitivity (mV/g)

g-Select1

g-Select2

1.5g

800

0

0

2g

600

1

0

4g

300

0

1

6g

200

1

1

The Sleep Mode pin can be connected to a logic inputs for mode switch. Set this pin low to enable MMA7260QT in Sleep Mode that will only consumed trickle current. A high logic input at this pin will switch the sensor to normal operation mode.

A 5.3 Tilt sensing

The MMA7260QT can respond to gravity or constant acceleration due to its capacitive detection principle and mechanism. When gravity is perpendicular to an axis, its axis output is zero-acceleration and therefore is half the VDD (i.e., 1.65V for typical 3.3V application). When gravity is parallel to an axis and the gravity direction is toward the positive direction of that axis, its axis output is half the VDD plus the selected sensitivity (Table A5.1). For 2g sensitivity application if the sensor is placed horizontally on a surface with top side upward (referring to Figure A5.1), the x- and y-axis outputs are 1.65V but the z-axis output is 2.25V due to applied gravity. Table A5.2 further illustrates the MMA7260QT outputs when it is placed at different orientations,

Table A5.2 The MMA7260QT outputs with respect to different sensor orientations (VDD=3.3V, 2g sensitivity)

Orientation

Outputs (Vout)

x-axis

1.65

2.25

1.65

1.05

1.65

1.65

y-axis

2.25

1.65

1.05

1.65

1.65

1.65

z-axis

1.65

1.65

1.65

1.65

2.25

1.05

The gravity response capability of the MMA7260 is useful for accurate tilt sensing with respect to any orthogonal planes. Assume the ,  and  are the tilt angles of x-, y-, and z-axis with respect to horizon, respectively with known accelerations obtained from Equation (A5.1), all the three tilt angles follow sinusoidal relationship that is expressed in Equation (A5.2).

                                                                                   (A5.2)

Note that Equation (A5.2) stays most accurate when the MMA7260QT is hold static or in the condition without varying acceleration. The resolution (acceleration changed per degree, i.e., the slope the sinusoidal curve) for any axis also varies with tilt angles due to the sinusoidal relationship. Take the x-axis for example, the maximal resolution can be obtained when its tilt  increases from 0° or 180°, and the minimal resolution occurs at  approaches 90° or 270°. Therefore a modified tilt calculation as Equation (A5.3) is suggested and is valid and applicable because it combines other axis outputs and therefore a maximal resolution of tilt sensing can be retained across any rotation and orientation with respect to any axis.

                                                               (A5.3)