Citation |

- Permanent Link:
- http://ufdc.ufl.edu/UFE0047295/00001
## Material Information- Title:
- Positioning of Accelerometers in Inertial Inclinometers and Its Pre-Processing
- Creator:
- Zhao, Linghan
- Place of Publication:
- [Gainesville, Fla.]
Florida - Publisher:
- University of Florida
- Publication Date:
- 2014
- Language:
- english
- Physical Description:
- 1 online resource (12 p.)
## Thesis/Dissertation Information- Degree:
- Master's ( M.S.)
- Degree Grantor:
- University of Florida
- Degree Disciplines:
- Mechanical Engineering
Mechanical and Aerospace Engineering - Committee Chair:
- CRANE,CARL D,III
- Committee Co-Chair:
- SCHUELLER,JOHN KENNETH
- Graduation Date:
- 8/9/2014
## Subjects- Subjects / Keywords:
- Acceleration ( jstor )
Accelerometers ( jstor ) Calibration ( jstor ) Coordinate systems ( jstor ) Gyroscopes ( jstor ) Kalman filters ( jstor ) Low pass filters ( jstor ) Mathematical vectors ( jstor ) Sensors ( jstor ) Signals ( jstor ) Mechanical and Aerospace Engineering -- Dissertations, Academic -- UF accelerometer -- calibration -- filter -- gyroscope -- positioning - Genre:
- bibliography ( marcgt )
theses ( marcgt ) government publication (state, provincial, terriorial, dependent) ( marcgt ) born-digital ( sobekcm ) Electronic Thesis or Dissertation Mechanical Engineering thesis, M.S.
## Notes- Abstract:
- Inclinometers based on inertial sensors are widely applied in modern industries. For better performance in accuracy and stability, usually the fusion of information provided by accelerometers and gyroscope is conducted. The algorithm in which the readings from two types of sensors are combined requires previously determined information about the position of accelerometers on object to be measured. However, the determination of positions of accelerometers relative to rotation center is not always available in practice due to geometry issue or flexibility of the surface where sensors are to be attached. To deal with this problem, a positioning algorithm for determining the positions of accelerometers using information provided by sensors themselves is presented here. The positioning process includes a complementary filter consists of Kalman filter and Savitzky-Golay filter to reduce noise the output of gyroscope, an analytical solution of five parameters to be determined in positioning, and an Extended Kalman filter to converge to the value of parameters directly. The positioning process can be conducted on-line once certain parts of sensor calibration are completed. To test the performance of complementary filter, Allan Variance is introduced. A number of experiments were conducted using a pair of MEMS dual-axis accelerometers and one MEMS single-axis gyroscope attached to a link with a rotating joint to validate the performance of positioning method proposed here. The results show that the algorithm is really applicable in practice to obtain the distances from accelerometers to the joint, and reach the degree of positioning error normally less than 10%. Moreover, due to concerns over accuracy of the method presented here, the calibration process for gyroscope and accelerometers is carefully conducted. To avoid the possible problem of commonly used least square method while applied in calibration of accelerometer, the determinations of sensitivities and offset are conducted independently. A two-location method for calibrating sensitivities is presented and applied, providing a good performance in calibration. The novelty of the positioning process is that it reduces the reliance of inertial measurement on external information. Once the sensors are soundly calibrated, they can be attached to any place on object directly, and see where they are placed by themselves. The application of the positioning method is foreseen in human body tracking, gait analysis, robot and manipulator working in aggressive condition, remote control, etc. ( en )
- General Note:
- In the series University of Florida Digital Collections.
- General Note:
- Includes vita.
- Bibliography:
- Includes bibliographical references.
- Source of Description:
- Description based on online resource; title from PDF title page.
- Source of Description:
- This bibliographic record is available under the Creative Commons CC0 public domain dedication. The University of Florida Libraries, as creator of this bibliographic record, has waived all rights to it worldwide under copyright law, including all related and neighboring rights, to the extent allowed by law.
- Thesis:
- Thesis (M.S.)--University of Florida, 2014.
- Local:
- Adviser: CRANE,CARL D,III.
- Local:
- Co-adviser: SCHUELLER,JOHN KENNETH.
- Electronic Access:
- RESTRICTED TO UF STUDENTS, STAFF, FACULTY, AND ON-CAMPUS USE UNTIL 2015-02-28
- Statement of Responsibility:
- by Linghan Zhao.
## Record Information- Source Institution:
- UFRGP
- Rights Management:
- Applicable rights reserved.
- Embargo Date:
- 2/28/2015
- Resource Identifier:
- 968786325 ( OCLC )
- Classification:
- LD1780 2014 ( lcc )
## UFDC Membership |

Downloads |

## This item has the following downloads:
ZHAO_L.pdf
ZHAO_L_pdf.txt FileType2_A_Triaxial_Accelerometer_and_Portable_Data_Processing_Unit_for_the_Assessment_of_Daily_Physical_Activity_pdf.txt |

Full Text |

PAGE 1 POSITIONING OF ACCELEROMETERS IN INERTIAL INCLINOMETERS AND IT S PRE PROCESSING By LINGHAN ZHAO A THESIS PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE UNIVERSITY OF FLORIDA 2014 PAGE 2 Â© 2014 Linghan Zhao PAGE 3 To my Parents Th e i r support, encouragement, and constant love have sustained me throughout my life. PAGE 4 4 ACKNOWLEDGMENTS I would like to express profound gratitude to my advisor, Prof essor Carl Crane, for his invaluable support, encouragement, supervision and useful suggestions throughout this research work. His moral support and continuous guidance enabled me to complete my work successfully. I would like to thank Prof essor John Schueller for serving on my committee. I am especially indebted to my parents, D octor Diansheng Zhao and Mrs Xiaowei Yu, for their love and support ever since my childhood. My family has been my ultimate support and I could not have done this without them. I would like to thank wonderful and inspiring teachers who have taught me over the ye ars, during my high school, during my undergraduate studies at Zhejiang University , Zhejiang University of Technology and my graduate studies here at University of Florida . I thank the fellow student Darsan Patel at the Center of Intelligent Machines and R obotics. I learned a great deal from him. I would also like to thank my classmates in Department of Mechanical and Aerospace Engineering, Joseph Yuan and Kecheng Tang. Special thanks to my dearest, oldest friend Zheng Ge. PAGE 5 5 TABLE OF CONTENTS page ACKNOWLEDGMENTS ................................ ................................ ................................ .. 4 LIST OF TABLES ................................ ................................ ................................ ............ 7 LIST OF FIGURES ................................ ................................ ................................ .......... 8 ABSTRACT ................................ ................................ ................................ ................... 10 CHAPTER 1 INTRODUCTION ................................ ................................ ................................ ....... 12 2 BASIC PROCEDURE OF ACCELEROMETER POSITIONING ................................ 25 Basic Principle ................................ ................................ ................................ ........ 25 Calibration of Sensors ................................ ................................ ............................. 26 Data P reprocessing ................................ ................................ ................................ 28 Positioning Algorithm ................................ ................................ .............................. 30 3 CALIBRATION OF GYROSCOPE AND DRIFT REDUCTION ................................ .. 35 Offset and Long term Drift Compensation ................................ .............................. 36 Sensitivity ................................ ................................ ................................ ................ 41 4 CALIBRATIONS OF ACCELEROMETERS ................................ .............................. 43 Least Square Method for Calibration ................................ ................................ ...... 43 Zero acceleration Offset ................................ ................................ ......................... 47 Dual position Calibration Method ................................ ................................ ............ 50 Validation of the Calibrated Parameters ................................ ................................ . 54 5 ONLINE ERROR REDUCTION OF GYROSCOPE ................................ ................... 56 Allan Variance and Allan Deviation ................................ ................................ ......... 56 Procedure Outline of Online Error Reduction ................................ .......................... 60 Determination of cutoff frequency for complementary filter ................................ ..... 64 Low Pass Filter for Obtaining LF and HF Components ................................ .......... 68 Original HF Signal Preconditioning ................................ ................................ ......... 74 ARMA Model for Preconditioned HF Signal ................................ ............................ 78 Kalman Filter for Preconditioned HF Signal ................................ ............................ 8 1 Result of Kalman Filter and Complementary Filter ................................ ................. 83 6 POSITIONING OF ACCELEROMETERS ................................ ................................ . 88 Analytical Result ................................ ................................ ................................ ..... 88 PAGE 6 6 Statistic al Filters ................................ ................................ ................................ ...... 95 Position ing using Extended Kalman Filter ................................ ............................... 99 7 CONCLUSION AND FU TURE WORK ................................ ................................ .... 104 LIST OF REFERENCES ................................ ................................ ............................. 106 BIOGRAPH ICAL SKETCH ................................ ................................ .......................... 116 PAGE 7 7 LIST OF TABLES Table page 4 1 Result of accelerometers offset calibration . ................................ ........................... 50 4 2 Results of calibrated sensitivities of accelerometers ................................ .............. 54 5 1 L ook up table for runs test given N = 3600. ................................ ........................... 77 5 2 N umber of runs of gyroscope signals . ................................ ................................ .... 77 6 1 Results of the analytical solutions, probability density method, Case 1 . ................ 94 6 2 Results of the analytical sol utions, probability density method, Case 2 . ................ 94 PAGE 8 8 LIST OF FIGURES Figure page 1 1 Sensor s configuration of VDI ................................ ................................ ................. 21 2 1 Configuration and coordinate systems for positioning of two accelerometers . ....... 25 2 2 Calibration pro cedure for gyroscope and accelerometer. ................................ ....... 28 2 3 Positioning program in cascade form for links chain starts from the base link. ...... 34 3 1 Gyroscope output sequence recorded in 1 hour . ................................ ................... 37 3 2 Long term drift in gyroscope output sequence. ................................ ...................... 37 3 3 Outputs of gyroscope and encoder during a pair of rotations ................................ . 38 3 4 Integral of gyro scope signal using original offset vs . the one using offset compensated by integral correction method . ................................ ...................... 40 3 5 G yroscope sensitivity results in 34 tests. ................................ ............................... 42 4 1 Eight accelerometer o rientations in z ero g bias c alibration . ................................ ... 47 4 2 O utput of accelerometer A1 for bias calibration. ................................ .................... 48 4 3 O utput of accelerometer A 2 for bias calibration ................................ ..................... 49 4 4 Sensor Axes and its equivalent orthogonal coordinate system . ............................. 51 4 5 Outputs of four accelerometer axes in a pair of static situations . ........................... 52 4 6 Residual of acceler ometer based on the calibrated parameters ............................ 55 5 1 Error percentage of Allan Variance vs. ratio of N and M . ................................ ....... 59 5 2 Schematic power density frequency distribution of information signal and noise/drift signal in output of gyroscope. ................................ ............................ 61 5 3 Complementary filter structure for gyroscope signal processing. ........................... 63 5 4 ADEV curve s of gyroscope signals in 20 tests . ................................ ...................... 65 5 5 Periodic component in gyroscope signal ................................ ................................ 66 5 6 PSD of original gyroscope signal vs . PSD of Savitzky Golay filtered gyroscope signal . ................................ ................................ ................................ ................. 73 PAGE 9 9 5 7 Result of normplot of output of gyroscope . ................................ ............................. 78 5 8 Autocorrelation function and P artial autocorrelation function of gyroscope noise .. 79 5 9 AIC values of 14 candidates of ARMA, AR, MA models . ................................ ....... 81 5 10 High frequency part of o riginal gyroscope signal and its Kalman filtered s equence ................................ ................................ ................................ ............ 84 5 11 ADEV of original static gyroscope sign al vs . ADEV of Kalman filtered static signal . ................................ ................................ ................................ ................. 84 5 12 Result of Kalman filter to process the gyroscope signal, Case 1 .......................... 85 5 13 Result of Kalman filter to process the gyroscope signal, Case 2. ......................... 86 6 1 Architecture and data flow in experiment of accelerometers positioning. ............... 88 6 2 Result of EKF of kinematic and distance parameters ................................ ............. 90 6 3 Calculated value sequences of four distance parameters r x,y and l x,y ..................... 91 6 4 Probability distribution of four distance parameters r x,y and l x,y ............................... 92 6 5 Outline of Kalman filter ................................ ................................ ........................... 97 6 6 EKF filtered angular velocity signal ................................ ................................ ...... 101 6 7 EKF filtered angular acceleration signal ................................ ............................... 101 6 8 EKF filtered values of four distance parameters r x,y and l x,y . ................................ 102 PAGE 10 10 Abstract of Thesis Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Master of Science POSITIONING OF ACCELEROMETERS IN INERTIAL INCLINOMETERS AND IT S PRE PROCESSING By Linghan Zhao Aug ust 2014 Chair: Carl Crane Major: Mechanical Engineering Inclinometers based on fusion of inform ation provided by accelerometer s and gyroscope s are widely applied in industries. The algorithm in which the signals were combined requires previously determined position of accelerometers on su bject. However, the determination of positions wa s not always available in practice due to geometry issue or flexibility of mounting surface . To deal with this, a positioning program for positioning accelerometers using information from sensors themselves is presented. The program includes a complementary filter to preprocess the output from gyroscope, analytical solution s of the position s to be determined and an Extended Kalman filter to improve performance of solution s . T he calibration s for sensors w ere carefully conducted for accuracy of positioning . S ensitivities and offset of accelerometer were calibrated i ndependently t o avoid problem s of regularly used least square method . A dual position method for calibrating sensitivities wa s applied for better performance. To test performance of complementary filter, Allan Variance wa s in troduced . E xperiments were conducted using two MEMS dual axis accelerometers and one MEMS single axis gyroscope to validate the method . The results show ed that the PAGE 11 11 algorithm wa s really applicable in practice to obtain the distances from accelerometers to the joint, and the relative positioning error was less than 10%. The novelty of the positioning process is that it eliminates the reliance of inertial measurement on external information. Once calibrated, the sensors can be attached to sub ject directly, and see where they are by themselves. The application is foreseen in human body tracking, gait analysis, robot and manipulator working in aggressive condition, remote control, etc. PAGE 12 12 CHAPTER 1 IN TR ODUCTION The motion characteristics of a three dimensional object can be described by six independent variables, three linear motions along the three perpendicular coordinate axes in the space and three rotational movements with respect to the axes. Measurement of the three orientations generates interest among researchers. I n order to accurately measure the rotational motion characteristics of an object, an angular motion sensor is required, which operates in essentially different principle from the mechanism of li near motion sensors. The angular motion sensor s can be typically categorized into two parts, inclinometers and goniometer s. T he former ones measure the angle of object relative to a natural vector, gravity or earth magnetic field, while the latter ones usu ally give the relative angle between two objects. H owever, since inclinometers could obtain relative angle using a pair of sensors, while a goniometer can be an inclinometer if given a predefined reference vector aligns with a natural vector, the two types of sensors can be deployed equivalently in practice . W e discuss inclinometers primarily in this thesis, but do not distinguish strictly between inclinometers and goniometer s if it is not necessary . I nclinometer is a type of a ngular motion sensor that can be applied for measuring angle of inclination or orientation . C urrently inclinometers based on various working principles [1 6] are widely used in ambulatory system to obtain the orientation of an object which is usually an important parameter in practice, for such tasks as gait analysis [ 7 ], human j oint a ngle t racking and activities monitoring [ 8 10 ], aircraft attitude estimation [ 11 ], a ttitude c ontrol of a s pherical r olling r obot [ 12 ] , load estimation [ 13 ] and level of acti v ity [ 14 ].I nclinometer s measure angles and convert them into electrical PAGE 13 13 signals, from which measurement and control system can derive the orientation of objects. They greatly vary with reference to accuracy, cost, reliability, and lifetime [15] . T o sense an angle , one direc t way is usin g an encoder , magnetic rotary encoder or optical rotary encoder . T he encoder meas ures the absolute or relative value of orientation on its shaft. T he advantage of encoder is primarily its significant accuracy. H owe ver, t he m ain limitation of an enco der that it is not always applicable in practice, because an encoder is a type of touch sensor. The shaft and the body of an encoder should be fixed to both side s of a joint respectively . The req uired mechanical contact interface makes the sensor , and subject s to wear and failures if deployed in harsh situation. Thus , the shaft should be fixed with relatively strict mounting requirements . F or some cases, it is even impossible to use an enco der, e.g., a human joint , or a joint that has been occupied by other equipment . T he characteristics of high price and high precision also limit application of encoder in low cost low accuracy cases. Apart from this, some types of relative inclination senso r based on other mechanical principles are proposed [16,17]. However, they also suffer the problem of mechanical contact and invasive character mentioned above. A nother type of sensors which can be the alternative of encoders is wearable sensor. I t avoids the problem of maintaining and precision suffered by encoder. Wearable sensors for measuring movement of joints have been studied for several applications [ 18,19 ]. In most studies, sensors are sewn to a piece of fabric and then 20,21 ]. However, uncertainties and drawbacks due to the fiber tensions and resistance alterations were reported [ 22 ] , and the problems of wear ing and localized space occupations still exist . PAGE 14 14 T o solve th ese problem s , some types of non touch sensing equipment were developed and applied. Typical non touch inclinometers are magnet based sensors, vision based sensors or inertial sensors, such as magnetic compass , optic marks system detected by camera or combination of gyroscopes and acceleromete rs . The magnetic sensors operate on the earth magnetic vector [6] . They find the orientation of sensor from the projection of earth magnetic vector on the sensitive axis. For vision based system s, computer vision techniques can be used to estimate motion p arameters directly from video footage without special markers, but these approaches are less accurate than optical systems [ 23 ]. In optical systems, m otion is captured by us ing a group of marks attached to a sub ject as reflectors and cameras placed around the su bject capturing light ray s reflected from the marks to record positions [24] . O nce a proper coordinate system is established, the change of position of marks detected by camera can be translated into the change of angle. Optical systems are common ly used and accurate in tracking movement [25, 4] . Unlike the mechanical sensors, the magnet based sensors and vision based sensors allow the angle measuring process to be absolutely decoupled from the physical joint axis or actuator axis by simply making sen sors [26], and consequently, the invasive character and maintaining requirement are avoided. I n addition, wires and tethers are not required for the common optical system. However, t he disadvantages of the two types of non touch s ystems above are that, for the magnetic one, the sensors should not be used in circumstances where ferromagnetic materials exists in the vicinity of magnetometer ; for the optic al one, the light signal may be jammed due to the environment or the motion, lea ding to the los s of the frames captured by camera . Thus the optical system require s PAGE 15 15 a clear line of sight between the source and the sensor . I t is expensive and can be used only in structured environments [27] . In order to provide an alternative means fo r motion tracking, new angular sens ing systems using different technologies have been develop ed . T he non touch inclinometers based on the inertial principle do not suffer the problems of magnetic disturbance or light jamming and laboratory setups , making a wider application field for them. Inertial sensors are designed for measuring acceleration and angular rate . The two basic types of inertial sensors corresponding to the two parameters to be sensed are accelerometers and gyrosco pes respectively . I nertial sensors have been popular in guidance, navigation, motion tracking and control applications for years [28 38] . In early days, due to the cost, size and complexity of inertial equipment, it was used primarily in fields of aviation and space industries or military applications [39,40] . I n recent years, however, new technologies for inertial sensors have been developed to achieve equivalent performance at lower cost [41] . In particular, rapid development of micro electromechanical sy stems (MEMS) with high accuracy, high reliability and multiple functionalities has provided a powerful tool set for inertial sensing [42 44]. Since the widespread availability of MEMS sensors today , inertial sensing technology provides low life cycle cost, small size, low production cost, and large volume manufacturing [ 41 ], making the sensors dominating inclinometer field in civil industries to a greater extent. In addition, the advance of the MEMS industry driven by consumer products has greatly improved the cost performance ratio of miniature accelerometers and gyroscopes. Over the past two decades, research on miniaturized inertia l motion sensors has received extensive attention, and continues to be an active domain [45, 46] . PAGE 16 16 O nce inertial sensors , especially MEMS inertial sensors are used to implement inclinometers, several typical configurations of gyroscope and accelerometers are available for obtaining proper information to find the angle, including accelerometer only configuration, gyroscope onl y configuration as well as the complementary configuration using accelerometers and gyroscope together. F or the accelerometers only configuration , a tri axial accelerometer can be used alone as an inclinometer in space [ 47 50 ], or a dual axial acceleromet er as well as two single axial sensors for angle in plane [ 22, 51,52 ]. F or some cases they were also used together with other sensors to find space orientation [53], or grouped to find attitudes integrally [54]. O nce the object to be tested is in static si tuation, there is only gravitational force applied to the sensor. The measured acceleration component on each sensitive axis represents the angle between the corresponding axis and the line of gravity vector. H owever, the problem of this method is that the result may be significantly deviated if there is relatively large extra acceleration compared to gravity applied to the object, making it less accurate in daily application since it is not easy or even impossible to eliminate all non gravity acceleration components. T o avoid the influence of extra acceleration on sensors in measuring angle, two optional ways are usually applied, the common mode rejection (CMR) method and the gyroscope only method . F or the CMR method, two dual axis accelerometers are mounted on adjacent links, both ideally attached to the joint center . Obviously the joint angle can be obtained from the difference in the directions of two acceleration vectors sensed by the two accelerometers respectively [55,56]. The CMR method ha s sign ificant errors if applied in rapid rotat ing cases, because usually it is not applicable to place two accelerometers at the center of joint connecting two links and to attach the PAGE 17 17 two sensors to adjacent links respectively [ 8 ] . T o eliminate the error introdu ced in CMR method, three alternatives are introduced based on the CMR, th e CMR with gyro integration (CMRGI) method [57] , the CMR with gyro d iffe rentiation (CMRGD) method [58] and the distributed CMR (DCMR) method [59,60] . H owever, the CMRGI and CMRGD methods subject to the problem of reduced accuracy over time due to drift accumulation and significant noise introduced by differentiation respectively. Theoretically , the DCMR method can calculate all joint angles without suffe ring the error due to non ideal placement faced by the CMR method [15]. However, the DCMR method was designed to find the relative angle between two coordinate systems established by configuration of accelerometers attached to two adjacent links. Therefore DCMR is a kinematics approach which depends on coordinate systems only, but does not provide any dynamics information which is necessary for dynamic control. F or the gyroscope only configuration, it can be used to obtain the angular velocity directly [61] as well as to calculate the angle. S ince the output of gyroscope is usually angular rate, clearly the current inclination angle can be determined by integral of gyroscope signal and the initial angle previously obtained in some other way s . F or some modern gyroscopes, the output of the sensor is directly the angle, making it more convenient in application as inclinometer. T he advantage of gyroscope only configuration is that the orientation can be derived without being affected by gravity or any inertial ef fect which causes problem in accelerometer only configuration [45]. O ne remarkable application of the gyroscope only configuration is the s trapdown integration algorithm used in inertial navigation [ 62 ] . I n this algorithm, the reading of gyroscope is integrated to find the inclination of computed virtual platform. H owever, one outstand drawback of the gyroscope only configuration is the accumulation of error over time PAGE 18 18 since integral is needed here. In pract ice , the integrals involved should be delicately handled so as to reduce the calculation errors as much as possible [63]. O nce the gyroscope is applied in angle measurement, usually some extra information sources such as magnetic compass, external optical marks are required for correction of the result from gyroscope only. Moreover , since inclinometers are commonly used to find absolute angle, the initial angle determined by other method is always required . T o reduce the problem of accelerometer only , gyros cope only configuration as well as the CMR, CMRGI and CMRGD methods i n calculating angles , one promising alternative is to estimate the orientation by fusion of accelerometer and gyroscope output , which is also the basic principle in which the vestibule sy stem of human works [64 68] . B y using this technology, the accelerometer and gyroscope acts as external information provider mutually to help reducing drift or noise in their outputs [69,70]. T he architecture combining the two types of inertial sensors above is usually called Inertial Measurement Units (IMUs) . Most likely contain a triaxial accelerometer and a triaxial gyroscope, IMUs are currently popular commercial devices used to sense movement and orientat ion of the moving body that can help calculate joint angles or inclination angles [12, 22, 71] . The design of commercial IMUs module s are commonly based on the principles of low cost, miniaturization, low power consumption, and modularity [27] . Aside from commercial IMU modules, currently combinations of inertial sensors in various fusion principles have been applied to various fields in which the parameters of orientation or inclination are needed, such as estimation of orientation and inertial effect to h uman head [5], measuring angle and angular rate of human knee joint [7,72], att itude e stimation for r obotic and h uman m otion t racking [73 75], estimation of s tate for f eedback c ontrol of h ydraulic m anipulator [76,77], t hree d imensional human g ait PAGE 19 19 a nalysis [78] and shape identification of bionic robot [79] . T he information provided by IMUs or combinations of sensors could be kinematic parameters like what obtained by the DCMR [5,22] , or dynamic parameters [12, 71 ]. It may be worthwhile to remark that a gyroscope which senses angular rate can be replaced by a group of accelerometers, making the combination an architecture consisting p u r e ly accelerometers [80,81]. A ctually, IMU is not a concept determining certain way in which the signals from two type s of sensors are processed in the sensor s architecture . F or angle measurement applications, i nformation obtained from the sensors in IMU is processed by an angle estimation method which may either be proposed by the developer or depend on the manufactory. Many algorithms for the angle estimation were developed [26, 82, 83]. Especially, Kalman filter based estimation has become common practice in IM Us application , either fusing information from types of inertial sensors within IMUs [84 86], or fusing informa tion from IMUs and other types of sources [12, 87 89] . S ince Kalman filtering is a recursive estimation procedure that uses sequential sets of measurements , it can improve the estimate of angle online at each step by taking the prior estimates and new sensor outputs [90] . There are multiple implementations and variations based on the Kalman filter or E , U nscented (UKF) using different estimations and optimization approaches [ 14, 39, 9 1 95 ]. T o avoid the s ingulariti es introduced by Euler representation which is regularly used in links configuration, quaternion is also introduced with Kalman filter in some researches [27, 82]. B ecause of the advantage of Kalman filter when applied in inertial based orientation estimat ion, the Kalman filter as well as the EKF is used in the program of this thesis. I n the conventional way using combination of gyroscope and accelerometer as an IMU , although error was compensated or corrected by information from PAGE 20 20 accelerometer, the long term integral of gyroscope signal wa s not avoided [ 93 95 ] . T he same problem also exists in some angle estimators based on Kalman filter, since the state update equations were established in an integral manner [93, 94]. However, similar to DCMR pr inciple, the integral process in IMUs can be cancelled by applying direct algebraic algorithm using outputs of accelerometers. One program called Vestibular Dynamic Inclinometer (VDI) was proposed in 2011, in which the inclination angle relative to the acc eleration vector of joint which was defined as Dynamic Equilibrium Axis (DAE) was obtained from measured accelerations [ 96 ] . I n this method, the problem of applying integral to gy roscope was absolutely avoided . T he advantage of VDI over the DCMR is that VD I obtains the inclination angle which can be further used in dynamics control strategy, making a b etter performance of sensors in inclinometer . An important question related to inertial sensor based motion recognition is sensor placement and displacement. C omparing the method of fusion of sensors and the ones based on one type of sensor only shows that, the fusion process as well as DCMR method requires previously determined positio n information of accelerometers, while it is not the case in most of one typ e of sensor methods. I f either one accelerometer o r one gyroscope is attached to the su bject whose inclination is to be measured, the angle can be obtained without information of location . On the contrary, the cost of fusion method being more accurate is t he requirement of exact information of sensors where they are attached to the su bject. Especially for some method, including VDI, a special geometry configuration of the sensors set is required to implement the algorithm. To deal with this, the conventiona l way is to measure the distance between sensors and joint prior to the application, and put the accelerometers at correct locations to satisfy the configuration if required. PAGE 21 21 T o demonstrate the requirement of locations of sensors fusion methods, we can co nsider the basic principle of VDI as a typical example of inertial sensor based inclinometer . T he configuration is shown in Figure 1 1 [ 96 ]. G iven the measured acceleration vectors a L and a R by both accelerometers, we have the difference and mean of the vectors as, where is the sum vector of gravitational acceleration g and the acceleration of joint O. is the angle of joint from vector, d is distance between the two accelerometers, is the angular rate obtained by gyroscope, l is the leng th of link shown in Figure 1 1 , and Ãª r , Ãª form the coordinate system fixed to the link. Figure 1 1 . Sensor s configuration of VDI , including one gyroscope denoted by G and two accelerometers denoted by L and R [96] . PAGE 22 22 B y mathematical manipulation, we have the following results, F rom the five equations above, obviously the angle and its 2 order rate as well as the absolute value of can be obtained b y 1,2 from accelerometers and from gyroscope. T his is the principle of inclinometer based on inertial sensors. Clearly, the inclinometer shown above requires two previously obtained distances l and r , which is considered as known parameters in measurement. The model used for VDI also implies that the two accelerometers should be placed at required and fixed locations on the link . Although there are differences among specified location requirements in various inertial inclinometers, the positions of sensors should be fixed is the basic rule. The vast majority of research in this area assumes well However, t he r equirement or assumption above may not be satisfied in real application, since it is not always available to place the sensors at required or fixed locations, especially in field of human body attitude tracking, the sensors could not be well fixed to body surface of people due to the flexibility of skin or clothes. Also in some cases, the locations are not easil y obtained previously using traditional measurement such as ruler. The tool for positioning may be costly or even PAGE 23 23 unavailable for the subject due to geometry issue, or size issue, or safety issue where people have to use remote operation [97] . T hus, eliminating the previous positioning procedure by using some other ways to deal with the locations of sensors online could be a promising way to improve the robustness of inertial system. However, there are few discussions on the sensors localization or positioning issue reported so far. The current way to improve the robustness of inertial system without given sensors position information is, with the help of gyroscope output, to determine if the measured acceleration is dominated by rotational effec t or translational effect [98]. I f rotational effect dominates the acceleration, this signal is simply ignored at the moment, instead, use the gyroscope signal only to monitor movement. The priori knowledge of the subject kinematic characteristics was also required. H owever, the problem of ignoring the requirement for is the loss of exact information which can be provid ed by sensors themselves. Moreover, the required assumption that the subject does not perform fast rotations and simultaneously also limits the feasibility of the method above in a wide range of cases . In order to reduce the reliance of inertial inclinometer on the previous auxiliary measurement in every application case and avoid the require ment of attaching the sensors at specified and fixed locations, an on line positioning method which depends on the information provided by gyroscope and accelerometers themselves to determine the locations of sensors is proposed and validated in following s ections. Primarily, this research involves the development of a prototype accelerometer positioning system together with careful calibration of inertial sensors and a complementary filter for drift compensation to gyroscope. T he accelerometer positioning algorithm is based on analytical solution and EKF. Since one underlying problem of the PAGE 24 24 commonly used least square method for accelerometer calibration was noticed. T he calibrations of offset and sensitivities were implemented separately . D ue to the drift p roblem which is usual in gyroscope, the output of it was carefully checked and processed, giving an error model which was then applied in Kalman filter to reduce the drift. T he performance of filter was validated using experimental data by Allan Variance. Physical experiments have been conducted to demonstrate the utility of the proposed positioning algorithm system. These tests did the position ing of two accelerometers fixed to a dual dimensional link with a rotation joint using the proposed method . Sensor data w ere provided to both analytical solution and EKF to obtain and compare the results by the two ways. Compared to the ignore signals if way mentioned previously, no signal neglect as well as kinematic assumption was made in this research, making it a wider field of applicatio n. By this positioning process, the combination of inertial sensors could be applied directly to the subject to be measured without previous procedure obtaining the configuration and gives the position parameters onlin e . Consequently, a foreseeable advantage of this method is improvement for the robustness of inertial inclinometer once applied in long term deployment where sensor shifts cannot be avoided. Moreover, better usability and user acceptance can be obtained if the PAGE 25 25 CHAPTER 2 B ASIC PROCEDURE OF ACCELEROMETER POSITIONING Basic Principle Figure 2 1 . Configuration and coordinate systems for positioning of two accelerometers . O nce the sensors were attached to an object which rotates about the rotation center , or the joint, we have the configuration shown in Figure 2 1 . Note that the position of gyroscope is not required in inertial measurement , thus it is not shown in the figure. T wo accelerometers A1 and A2 are arbitrarily attached . T he X and Y sensitive axes of both accelerometers are also arbitrarily selected as shown. The rotation angle between the two accelerometer axes system is denote d as . Note that t he non orthogonality of real sensitive axes of each accelerometer is not taken in consideration so far, that is, X and Y axes of one sensor shown in Figure 2 1 are mutually perpendicular. A ctually the X Y coordinate systems shown here depend on their PAGE 26 26 corresponding real sensitive axes. T he relationship will be discussed later in Dual position Calibration Method . The advantage of introducing this relationship is obtaining approximately identical cross sensitivities between X and Y axes of an accelerometer . V ec tor r A1A2 is the distance vector between A1 and A2 whose length is l . r A1 and r A2 are radius vectors of A1 and A2 respectively with respect to rotation center. r is the distance from rotation center to the midpoint of r A1A2 . T o determine the positions of A1 and A2 is to find th e two vectors r A1 and r A2 . It is obvious that once a unified global coordinate system is defined, determining the two radius vectors is equivalent to determining the relative position r A1A2 from the accelerometer to the other as well as the vector from joint center to midpoint of r A1A2 . I n other words, once the projections l x,y and r x,y of distance l and r onto the global coordinate system are obtained, each of the components of r A1 and r A2 in the global coordinate system is simply either sum or difference of l x,y / 2 and r x,y . N ote that the global system could be arbitrarily selected. F or the sake of simplicity, the global system was s elected to be aligning with the X Y coordinate system of one accelerometer in this research. A s shown in Figure 2 1 , the selected accelerometer to be base of global coordinate system is A1. B y est ablishing this global coordinate system attached to A1 , the task for positioning is to find l x,y and r x,y in this system using information from accelerometers A1, A2 and a gyroscope which is not shown in Figure 2 1 . Calibration of Sensors B efore the accelerometers and gyroscope being attached to the joint to be measured, the both types of sensors should be calibrated for bias voltage by zero input and sensitivities of sensitive axes. PAGE 27 27 N ormally external information provided by auxiliary sensor is required as reference signal for the calibration process . S ince the auxiliary sensor would not be deployed in the combination of gyroscope and accelerometer in real application, it is necessary to conduct som e steps of calibration offli ne . T he external information source applied in this research was an encoder providing exact angles of rotations. T he output of encoder was considered to be ground truth in experiment. The parameters requiring offline calibration steps in this document are offset and sensitivity of gyroscope and zero G bias of accelerometers. O nly the sensitivity of gyroscope was calibrated w ith the help of encoder, gyroscope bias and accelerometer s offsets were obtained without e ncoder information . The coarse estimation of offset of gyroscope was simply determined by averaging gyroscope data sequence obtained in static state for a sufficiently long period. However, the gyroscope offset can be checked and corrected in real measurement if only some particular mo tion style could be conducted, making it possible to cancel the calibration for gyroscope offset before deployment , but use the value provided by manufactory or the averaging as the initial estimation instead. T he correction procedure is discussed in Offset and Long term Drift Compensation . T he sensitivity of gyroscope was obtained by comparing encoder output with integral of gyroscope signal given a rotational motion. T he zero G bias of each accelerometer was d etermined using an eight orientation calibration method in which the accelerometer was placed approximately horizontal and in 8 orientations, then the 8 outputs of accelerometer were recorded and applied for error reduction. T he parameter s that can be estimat ed online are sensitivit ies of accelerometers, because the external information provided for this calibration could be the integral of output of calibrated gyroscope, that is, the rotation angle. A dual position calibration was PAGE 28 28 applied to find out b oth X Y sensitivities as well as one cross sensitivity of each accelerometer. T he calibration is discussed in detail in Dual position Calibration Method . O bviously this calibration can be conducted with enc oder also, however, the use of gyroscope reduces the reliance of calibration for accelerometer sensitivities on auxiliary equipment, and make it possible to correct the sensitivities while the joint to which the sensors are attached is rotating, improving the accuracy of measurement and robustness of system . I n conclusion, we have the process of calibration shown in Figure 2 2 . The implementations of calibration are to be described in detail in Chapter 3 and Chapter 4 . Data P reprocessing Coordinate System Establishment O nce the gyroscope and offsets of accelerometers were calibrated, the combination of sensors could be attached to th e link whose orientation is to be tested. H owever, before positioning, some information is required from output data to establish the necessary global coordinate system. Figure 2 2 . Calibration pro cedure for gyroscope and accelerometer . PAGE 29 29 T o determine the global coordinate system based on A1 and the sensitive axes system of A2, the vector of gravitational force was used as the reference. Outputs of both accelerometers were obtained in static condition, in which the angles between sensitive axes and G vector can be found. T he angles of A1 established the global coordinate system, while the angles of A2 give the rotation with those of A1. N ote that the accelerometer A 2 was fixed relative to A1 all through measur ement , thus the angle wa s constant during measurement. T o rotate the readings and of A2 to the equivalent accelerations in global coordinates for further positioning process , the following rotation matrix was applied, T he preprocessing of accelerometer signal was applied to the one which is not the base of global coordinate system only. U sing the rotation matrix R obtained above , all signals from accelerometers were rotated to the global system. Obviously in static situation, since the only applied acceleration is the G force, the rotated outputs of accelerometers should be all the same. H owever, in dynamical case, they may not be identical due to the existence of radial and tangential acceleration components shown in Figure 2 1 . Since the radial and tangential acceleration both depend on the distance between sensor and the joint, the differences of output of accelerometers on global system are one of the key point of positioning. Reducing Error in Data While positioning the accelerometers, the data obtained by sensors should be first processed before the positioning algorithm . PAGE 30 30 T he purpose of pre processing of gyroscope signal was mainly reducing the random noise and long term drift in data sequence , giving a better performance of sensed angular velocity for further calculation and for the integral of the s ignal obtaining angle. T he output of gyroscope was pre processed by a complementary filter which was a combination of Savitzky Golay filter and Kalman filter in this document. I t is discussed in detail in Chapter 5 . Positioning Algorithm Once accelerometers readings are transformed to the global coordinate system, we have the following sum and difference of readings , w here a A1 and a A2 are the sensed acceleration vector of A1 and A2, x A1 and y A1 forms the global coordinate system based on accelerometer A1, and are angular acceleration and angular velocity respectively, r x,y are the projection of distance vector r shown in Figure 2 1 on x A1 and y A1 , similarly, l x,y are projection of vector l in Figure 2 1 . a 0 is the acceleration of joint. is angle between vector a 0 and axis x A1 , and is the one between gravitational vector and x A1 . Then w e have four equations from the two equations, PAGE 31 31 where in which a 0 is the acceleration of joint in inertial space, g is the gravity acceleration, 1 is the angle between X axis of global system and the a 0 vector at this moment (Mom ent 1), and 1 is the angle between X axis of global system and the g vector at the same time. Theoretically, there are four unknowns , the four position components r x,y and l x,y . However, due to the significant effect of noise in differentiation , it is not available to obtain the a ngular acceleration by simply differentiating angular rate signal from gyroscope. T hus the angular acceleration should be the fifth unknown, and at least five equations are needed. To this end , difference of acceleration at a second time (Moment 2) was introduced , leading to the following two equations, T o find analyti cal solutions of five unknowns, we conducted mathematical manipulation as follows. M ultiplying l x and l y to both sides of Equation 2 6 and 2 7 respectively gives, T hus we have . S imilarly, we have S ince, PAGE 32 32 W e have the ratio of two distance components, S ubstituting the ratio above to Equation 2 12 gives, Now firstly we have the two distances, and , To find the angular acceleration, w e have from Equation 2 10 and 2 11 . T hus, T o obtain coordinates r x and r y , we have, By Equation 2 1 6 we know that, PAGE 33 33 T hus, W e have, and, T he position of joint, that is, the rotation center, is located at ( 0.5 l x + r x , 0.5 l y + r y ) in the global coordinate system giv en origin at accelerometer A1 . E quivalently, the positions of A1 and A2 can be determined relative to the joint, since location of A2 relative to A1 has also been obtained. Note that the parameters C 1 and S 1 are both function o f a 0 , , g and 1 , where g is known, 1 can be determined by integral of gyroscope signal together with the last value of angle . T he value of a 0 , could be determined in cascade style in practice . S ince the acceleration of base joint which connects the base link to ground can be assumed to be zero, values of C 1 and S 1 of the base link depend on g and 1 only. T hen the five unknowns of the base link can be obtained from accelerometers and gyroscope attached to it . Usually it is also assumed that the distance between adjacent links are PAGE 34 34 known, since they are fixed and provided by designer or manufactory of the objects. B y the given distance from base joint to the next joint and obtained positions of sensors, angular a cceleration and angular velocity of the base link, the acceleration vector of the next joint can be easily determined using the DCMR principle [ 15 ], including the information of a 0 and orientation of vector a 0 in inertial space . S ince orientation of a global coordinate s system of one link in inertial space can be determined from all integral s of gyroscopes from the base to that link, the relative orientation between vector a 0 and the corresponding global X axis can be derived directly. T he cascade for m of positioning and determination of C 1 and S 1 are shown in Figure 2 3 . I n this research , only one link was applied in experiment for simplicity , in other words, the acceleration of joint and its corresponding angle are set zero in experimental validation . I n conclusion, we have the five solutions of unknowns for each link , Equation 2 1 3 , 2 1 4 , 2 1 5 , 2 1 8 and 2 1 9 for l x,y , , r x,y respectively. These are the basic analytical solution s for accelerometer s positioning. Fig ure 2 3. Positioning program in cascade form for links chain starts from the base link . PAGE 35 35 CHAPTER 3 CALIBRATION OF GYROSCOPE AND DRIFT REDUCTION Two parameters of a gyroscope, offset voltage and sensitivity, are to be determined in calibration process. I n our project, an encoder wa s required in obtaining sensitivity of gyroscope. Since encoder wa s not the sensor deployed in inclinometer, calibration of gyroscope wa s implemented offline . A single axis MEMS gyroscope ( ADXRS614 ) was calibrated, compensated and applied in further experiment. A m agnetic encoder ( US Digital MA3 A10 125 N ) w as used for calibration of gyroscope. The measurement of the encoder angle was assumed to be the ground truth. Usually the calibration of gyrosco pe offset and sensitivity is carried out on a turntable with sufficiently high accuracy . However, due to high cost and relatively large volume , turntable may not be available in all calibration case s . Instead, w e calibrated the gyroscope using output of an encoder only, and implemented online compensation for long term drift of offset in a way that is not difficult to be carr ied out in real application. T he calibrating processes of offset and sensitivity were implemented separately . T he initial value of co nstant offset voltage was obtained by averaging the output signal over a relatively long period under static situation. T hen the sensitivity was obtained by comparing integra l of the gyroscope signal (initially estimated offset subtracted) in a static rota ting static motion style with the corresponding output of an encoder. S ince the initial offset was obtained without considering temperature effect, the temperature compensation is necessary for online application. The way of using a steady dynamic detector to identify static period during working and then finding average value in static situation is applicable [ 55, 99 ] . B ut once two rotations in opposite direction with same PAGE 36 36 angle are implemented, temperature effect could be reduced by this pair of motions a nd consequently more accurate offset could be obtained. Offset and Long term Drift Compensation T he offset of gyroscope is initially obtained by averaging the output signal V 0 in static situation for a certain period. B y the datasheet of ADXRS614, the typical bandwidth of this gyroscope is 80Hz at zero angular rates . I t had been proved by AVAR method that the sampling rate should be at least 3~5 times the bandwidth of gyroscope to obtain reliable performance of error ana l ysis [ 100 ] . T hus the required sampling rate should be no less than 240~400Hz. At each task of averaging for offset here, output of gyroscope is recorded for 1 hour with a sample rate 5 00Hz under static condition. T he recorded signal sequence of one test i s shown in Figure 3 1 . T he initial value of gyroscope offset voltage is obtained directly as, where V 0 (i) is the i th measured value in static condition. However, when the offset obtained above was applied in determination of sensitivity and further online measurement, it is usually problematic due to existence of error which is not taken into account in averaging method . S ince noise at high frequency i s eliminated in averaging process, the major part deteriorates performance of determined offset is random drift in relatively lower frequency, which can be detected by AVAR at higher averaging time . T he primary long term components of offset error in static situation include temperature drift, floating power supply and offset instability. Since effect of temperature and the one of drift in power supply voltage on the offset i s relatively significant in some cases and is accumulated in integration, they should be considered and compensated. PAGE 37 37 Figure 3 1 . Gyroscope output sequence recorded in 1 hour . Figure 3 2. Long term drift in gyroscope output sequence. By the datasheet of ADXRS614 gyroscope, the typical temperature scale factor is 9 mV/ at around 25 . To clearly demonstrate the effect of temperature and floating power supply on offset voltage, an one hour signal sequence shown in Figure 3 1 wa s smo othed by a 4 th order Savitzky Golay filter with 2001 one side points. T he filtered signal is shown in Figure 3 2 . I n this case, i t is clear that a drift caused by temperature or floating power supply makes a voltage deviation at about 5 mV in 1 hour measurement . N ote that it is not necessary to separate effect of temperature from the one of power supply on drift. They are both compensated together using the following method. PAGE 38 38 T o compensate for the long term deviation of offset in calibration procedure, two potential ways could be applied for both online or offline cases. One is traditional averaging method with dynamic static detector. The other one is an integral correction method. T o implement the first ap proach, a dynamic static detector [ 55, 99 ] using variance and a previously set threshold was applied to signal from one accelerometer to determine the periods in which the gyroscope accelerometer system was in static situation. W ithin the determined static period, averaging was applied to gyroscope signal to obtain the current estimation of offset. T o implement the second one, at least two rotations of the gyroscope to be compensated are needed . T he requirement of the two rotations is the same angle of rota tion but opposite directions. A ngular rate distribution over time during rotations is not required here. O ne example of gyroscope and encoder outputs obtained in the required pair of rotations is shown in Figure 3 3 , A B Figure 3 3 . Outputs of gyroscope and encoder during a pair of rotations . A) Gyroscope signal sequence, B) encoder signal sequence. PAGE 39 39 S ince the angle between the two static positions is integral of angular velocity measured by gyroscope, we have, where is the angle between the two static positions , S gyro is sensitivity of gyroscope , t 0 is an arbitrarily selected time in initial static period, and t 1 in the second static period. I n calibration, time between t 0 and t 1 is restricted so that long term drift in gyroscope signal does not quite affect the integral result. Note that output of gyroscope is a sequence of voltages, given the sample interval t , Equation 3 2 should be changed to discrete form, F or the case of the pair of rotation shown in Figure 3 3 , the angles of both rotations are equal but in opposite direction, thus, after both subtracted by offset, the discrete integrals of the two gyroscope signals should give same result if the applied offset O gyro.accu is accu rate . From Equation 3 3 , we have the following equation, w here superscript s cw and ccw of V stand for the output voltages obtained in clockwise and counterclockwise rotation respectively, t 0 and t 1 denote the beginning and ending points of the clockwise rotation respectively, t 2 and t 3 denote the ones of the co unterclockwise rotation. the applied offset is not drift compensated, usually we don t have the equality. A ssume that the uncompensated offset is O gyro , the needed is O , we have O gyro.accu = O gyro + O and the following equation, PAGE 40 40 T hus the difference between two discrete integrations is O(N + M) t (N + M) t O Ã” gyro.accu = O gyro + O Figure 3 4 . Integral of gyro scope signal using original offset vs . the one using offset compensated by integral correction method . PAGE 41 41 t he original offset was obtained by averaging previous static signal, which was 2.40338V . The correction O found by the mean difference between three peaks and three valleys points was + 0.00395 V. Thus the estimated accurate offset was 2.40733V . The integrals of six rotations using the corrected offset are shown in Figure 3 4 in solid li ne and marked by . C learly results of the six offset compensated integrals are flatter and thus closer to the corresponding one same angle. However, t he integral correction method has a requirement of a fixed angle of rotation where the exact value of a ngle is not necessary . In practice, with the help of position limit switches which are widely used in industries, this requirement is not difficult to be satisfied. Therefore t his compensation method can be conducted online intermittently to reduce the effect of offset deviation on the sensing performance. Sensitivity The sensitivity of gyroscope is the only parameter obtained with the help of auxiliary sensor which is not included in combination of gyroscope and accelerometers. H ere an en coder is used to find this sensitivity, thus it is the only parameter that should be determined offline. T o obtain the sensitivity, the gyroscope is set to a static position and holds for a while initially, then rotates to a new static position. From Equation 3 3 , we have, N ote that the angle was provid ed by the encoder. S ince the offset O gyro has been determined in the previous section, once the rotation angle by encoder and output sequence by gyroscope are both obtained, we can find out the sensitivity S gyro by way of discrete integral. Thirty four tests were implemented to calculate sensit ivity of gyroscope, PAGE 42 42 the result are shown in Figure 3 5 . The variance of obtained sensitivities was 0.0806 , which was determined as the gyroscope parameter used in further experiment . B ased on the mean value of the 34 test results and the maximum, minimum values among them, t he range of deviation of sensitivities fr om their mean value was calculated to be 2.12% , showing a good accuracy performance of the result . Figure 3 5 . G yroscope sensitivity results in 34 tests . PAGE 43 43 CHAPTER 4 C ALIBRATIONS OF ACCELEROMETERS Five parameters of each accelerometer, zero g bias voltages of X and Y axis, sensitivities of X and Y axis and the cross sensitivity between X/Y axes , are to be determined in calibration process. I n our project, the zero g bias voltages were obtained off line, and all sensitivities were calibrated when the accelerometers were attached to object and doing the measurement . It means that accelerometer sensitivities could be calibrated with the help of a calibrated gyroscope onl y and thus can be done online, no extra sensor was required . T he biases of both axes of each accelerometer were measured off line independently. F or each accelerometer, it was laid on a flat surface to make both X and Y axes closely approaching zero g posit ion. T hen the same accelerometer was placed upside down obtaining the other approximate zero g position. T he average value of the two positions was considered real bias . O nce bias obtained, the accelerometer is ready to be fixed to link, manipulator, etc. to make a VDI system with gyroscope. Then, to calibrate for sensitivities, the gyroscope signal is to be processed to reduce noise and then integ rated for determining orientation of accelerometers. A dual position method was then implemented online to obtain three sensitivities for each accelerometer. Two dual axis MEMS accelerometers ( ADXL320 J ) were calibrated, compensated and applied in further experimen t. Least Square Method for Calibration T o calibrate for the biases and sensitivities of accelerometers, least square method was commonly conducted [ 55, 56, 10 1 , 10 2 ] . T he advantage of this method is that it does not require external information except for the gravity vector. PAGE 44 44 T he normally used least square error criterion for calibration of accelerometer is where k denotes the k th test of N tests, and, Note that the accelerometer outputs a X ( k, ) and a Y ( k, ) are normalized since for acceler ation sensors the magnitude of gravity is usually considered as unit 1. The proper estimation of the vector which contains all parameters to be calibrated is the estimation which minimize s the e rror . T o solve the least square problem above, two conventional methods were usually applied, the Gauss Newton method or the Newton Raphson iteration method [ 55, 56, 101, 102 ] . If the Gauss Newton method is applied, Jacobian matrix and error vector are introduced as , where is the estimation of at each time step. F or the first step in iteration, is an initial guess of . T he Gauss Newton method works iteratively to give a series of estimations . E ach corresponds to an error , and the magnitude of decreases as the subscript n increases . T he solution of Gauss Newton method converges once the magnitude of falls in a predefined threshold, giving the solved estimation . I ntroduce the error of the i th estimation of vecto r , PAGE 45 45 B y Gauss Newton method, since the Equation 4 1 can be rewritten as, the estimation of the error is obtained as, where , . Then the updated estimation of the parameter vector further reduce s ( ) to ( ) . Once converged, the final estimation of vector can be considered as the calibrated parameters. I f Newton Raphson iteration is applied, to minimize the least square error criterion, we have a new estimation of vector based on its previous estimation , w here H ( k, ) is the Hessian matrix of h ( k, ) , is a damping parameter smaller than 1 . H owever, the least square method using either Gauss Newton method or Newton Raphson iteration may have a problem of ill conditioned Jacobian matrix. I n the current accelerometer calibration case, t he acceleration components a X and a Y can be solve d from sens itivities, biases O X,Y and measured signals V X,Y as , By mathematical manipulation, five elements of Jacobian matrix h ( k, p ) / p are, PAGE 46 46 A mong the five elements in Jacobian vector shown above, though not directly clear, the magnitude of the last two elements , that is, h ( k, p ) / O X,Y , were about thousands times greater than the first two elements h ( k, p ) / S XX,YY and hundreds times greate r than the third element h ( k, p ) / S XY in typical conditions, making the Jacobian vector as well as the matrix to be ill conditioned. F or example, given a typical value of vector p in experiment p = [145 145 0.5 1.5 1.5] T and typical measured values V X = 1.27 V, V Y = 1.23 V , w e have the following numerical Jacobian vector, T he Jacobian matrix above is ill conditioned, since too large ratio s among elements exist . O bviously in the matrix consists of the Jacobian vectors, elements in PAGE 47 47 the first three rows are significantly smaller than ones in the lower two rows. B ecause of this, the matrix is also ill conditioned. Since ill conditioned Jacobian matrix or vector may cause the solution fou nd by either Gauss Newton way or Newton Raphson way being in correct, the least square method is not applied in this research. Zero acceleration Offset Unlike the way in which some papers have done, in this research , the zero g bia s es of all active axes were calibrated separately from the calibration of sensitivities. T o determine the bias at zero input, both accelerometers A1 and A2 were placed approximately in horizontal situation. T hen the accelerometer in calibration was set to 8 o rientations in horizontal plane. T he angle between two sequential orientations was approximately 45 degrees. A fter all 8 orientations completed, the accelerometer was inverted to be upside down. A gain, eight orientations were implemented. T he orientations of accelerometer axes in test s are shown in Figure 4 1 . Figure 4 1 . Eight accelerometer o rientations in z ero g bias c alibration . PAGE 48 48 A B C D Figure 4 2 . O utput of accelerometer A1 for bias calibration . A) X axis at position A, B) X axis at position B, C) Y axis at position A, D) Y axis at position B . PAGE 49 49 A B C D Figure 4 3 . O utput of accelerometer A 2 for bias calibration . A) X axis at position A, B) X axis at position B, C) Y axis at position A, D) Y axis at position B . PAGE 50 50 Tab le 4 1 . Result of accelerometers offset calibration . Accelerometer X axis zero g b ias (V) Accelerometer Y axis zero g b ias (V) A1 1.3582 A1 1.3477 A2 1.3594 A2 1.3550 T he calibration process was done at two different places A and B trying to eliminate error in local geometric shape of surface. N ow we have two accelero meters, two places, two axes for each accelerometer, then the eight combinations are shown in Figure 4 2 and Figure 4 3 describing approximate zero g bias characters of all axes of the two accelerometers. I n Figure 4 2 and Figure 4 3 , the c urves with marks are the measurement data, while the thicker curves between them are the averaged data. N ow that we have obtained eight curves of four axis at places A and B. Base on analysis on the curves, for each axis, the mean value of peak and valley along the curve of eight orientations was found as the offset voltage of corresponding axis. T he result of offset calibration is shown in Table 4 1 . Dual position Calibration Method A ll four offset values have been obtained and then can be applied in determination of sensitivities . For each accelerometer, the sensitivities to be calculated include the sensitivities S XX , S YY on X axis and Y axis respe ctively, and the cross sensitivity S XY of X direction acceleration on Y axis, as well as the one S YX of Y direction acceleration on X axis. T he outputs of X and Y axes are, PAGE 51 51 N ormally the two cr oss sensitivities are not identical. However, if the cross sensitivity was assumed due to misalignment of two sensitive axes, S XY c an be set equal to S YX artificially by using a proper orthogonal coordinate to express components of sensed acceleration. A s the case shown in Figure 4 4 , the sensed acceleration is decomposed into two components a X and a Y on axes X and Y respectively. We have approximately , where the 0 in subscripts denote the sensitivities of real sensitive axes. Figure 4 4 . Sensor Axes and its equivalent orthogonal coordinate system . Since the X Y coordinate system is used for decomposition of acceleration only, it can be arbitrarily selected. I f we have S XX 0 sin 2 = S YY 0 sin 1 , approximately 2 / 1 = S YY 0 / S XX 0 since the two angles 1,2 are small, then the sensitivities are, A nd the output of accelerometer in matrix form is, PAGE 52 52 Figure 4 5 . Outputs of four accelerometer axes in a pair of static situations . T o determine the three values of S XX , S YY , S XY , at least three equations are required. F or each static position, the accelerometer has two outputs corresponding to the two gravitational components exerted on X Y coordinate system defined by 1 and 2 determined previously. T hus two static positions are required for the three parameters. One example of test of a pair of static situations is shown in Figure 4 5 . For each static case, only gravitational acceleration 1 g applies to the sensor. U sually unit of the sensitivity of accelerometer is mV /g , thus the acceleration applied on X or Y coordinate axis can be expressed by the angle between the axis and gravity vector. S ince two static positions were used for calibration, one of them can be denoted by angle and the other one by ( + ) . is the angle change between the two positions. T he Equation 4 9 can be rewritten as, where subscript s and ( + ) denote the corresponding positions. U se the inverse matrix of the second term on left, PAGE 53 53 where . Since, We have, F inally, Thus we have, Sinc e has a range from 180 to 180 , arctangent function has ambiguity. The signs of ( V X O X ) and ( V Y O Y ) are needed for determining correct value o f . If, and , satisfies If, and , satisfies If, and , satisfies If, and , satisfies Now we have , then the Sensitivity matrix can be obtained as, PAGE 54 54 where , B y this way, the sensitivities of the two accelerometers were calibrated in 40 tests. The sensitivity of each axis is the mean value of the 40 results. The calibrated results together with deviation errors of sensitivities in the 40 tests are shown in Table 4 2 . The error in the table is the ratio of range of results calibrated of one axis to their mean value. Tab le 4 2 . Results of calibrated sensitivities of accelerometers . Accelerometer S XX (mV/g) S XX error S YY (mV/g) S YY error S XY (mV/g) A1 1 43.378 0.113% 145.368 0.420% 0.160 A2 145.674 0.150% 144.904 0.473% + 0.316 Validation of the Calibrated Parameters T o validate the calibrated sensitivities and offsets, the two accelerometers were applied in measurement simultaneously then cross validation was conducted. F or each accelerometer in validation test, the residual acceleration was foun d by, where a X and a Y were obtained by, T he residual based on the calibrated parameters were tested for 10 times. T he typical results of A1 and A2 are shown in Figure 4 6 . N ote that the wave parts of the residual correspond to the rotating moments of system, while flat sections mean the static moments. C learly during static moments, due to non existence of extra acceleration, the magnitude of sensed acceleration should be PAGE 55 55 pur ely 1 g , causing residual equaling to zero. I n fact, the residuals of static periods shown in Figure 4 6 did lie on 0 g . These results show that the calibrated sensitivities and offsets are acceptable for further application. H owever, one thing noticeable in the following figure which is worthy of being mentioned is that the vibrations shown during rotating moments provide information of extra acceleration exerted on the accelerometers which contains information of both position s of sensors and dynamic characteristics. T his vibration form is actually the base of positioning algorithm discussed in this research. A B Figure 4 6 . Residual of accelerometer based on the calibrated parameters . A) Case 1, B) C ase 2. PAGE 56 56 CHAPTER 5 ONLINE ERROR REDUCTION OF GYROSCOPE T he angular velocity measured by gyroscope i s to be integ rated in determining positions of accelerometers. Since error in output of gyroscope would be accumulated by integral, it i s of great importance to be compensated for. A complementary filter consist ing of a low pass filter and a Kalman filter is applied here to reduce noise and drift in gyroscope output, and Allan Standard Deviation is used to validate compensation. S ince Alla n Standard Deviation as well as Allan Variance plays important rule in filter parameter determination and performance validation, the basic principle of the two concepts is introduced first in this chapter. Allan Variance and Allan Deviation E rror in gyros cope output signal plays significant rule in calibration and further measurement application. For better performance of the gyroscope, noise and drift should be eliminated or compensated. T o estimate error components in original signal, and to validate the effect of online error elimination by some technical ways, a variable characterizing statistical deviation of sensor output and the filtered signal is needed. A commonly used variable describing statistical dispersion of sequence is variance or standard deviation [ 103 ]. H owever, angular velocity obtained by gyroscope is usually integrated to give angle orientation. Drift in gyroscope signal has multiple components in terms of different characteristic time lengths. F or example, changing trend of random ci rcuit noise in MEMS gyroscope chip may varies within several milliseconds , while trend of temperature drift of output signal may varies in hours. In integration over a certain time range , these components have different performance in error accumulation an d thus have different effect on the result angle integrated. F or PAGE 57 57 example, white noise component and the noise with short characteristic time length usually have relatively smaller effect on integration due to the averaging of integral , while random walk co mponent which has all most same or even longer characteristic time lengths compared with time range of integration may have significant effect by accumulation. T hus, variable characterizing statistical deviations with different characteristic time lengths is needed for gyroscope output signal. Obviously variance is problematic here since it provides global dispersion only. Allan Variance (AVAR) method recommended by IEEE standards was initially put forward by David W. Allan in 1966. I t was initially used to estimate frequency stability of cesium clock due to noise processes . The principle of AVAR bases on the fact that frequency or phase of output signal would be deviated or modulated by noise, making instability of frequency. T he nois e components from different sources may cause different effect on this instability which can be distinguished in time domain, that is, instability at various time sizes. U sually in application, instability on short time range shows fluctuation in signal, w hile the one on longer time range stands for drift. Because vibration or resonant system are widely used in gyroscopes today, making gyroscope a special type of vibrator , AVAR can be used in test and analysis of gyroscopes . Today, AVAR is commonly used in testing of inertial measuring equipment because that AVAR has the capability of identifying various error sources and their statistic characteristics, making it easy to see noise or drift components in signal. AVAR would be applicable in case noise or drif t components need to be distinguished. M oreover, Allan Standard D eviation (ADEV) is defined as the square root of A VAR. I n this research, ADEV was applied in parameter determination and denoising validation. PAGE 58 58 T o derive the ADEV or AVAR of a given gyroscope signal sequence with N points and certain sampling time interval , t he entire sequence should be separated into K groups with M points in each group first . H ere we have M (N 1) 2 . The time length of each group is averaging time = M . T he average of each group is, where ( k 1) M + i is the measured angular velocity at time ( k 1) M + i . T hen, AVAR of averaging time is defined as, ADEV is defined as as we mentioned previously. AVAR describes drift intensity in different averaging time , thus can be treated as a spectrum of variance. H ere the averaging time corresponds to the characteristic time length of corresponding AVAR component of noise. Actually we have a relationship between AVAR and Power Spectrum Density (PSD) S (f) of the output signal [ 100 ] , It was assumed that the process of original signal was stationary while deriv ing the equation above. Eq uation 5 3 provi des another way to obtain AVAR . I t shows that AVAR at averaging time is always the same as passing power of original signal x through a filter whose transfer function is given as (sin 4 x ) x 2 . O bviously here the transfer function of filter or even pass band of the filter is determined by . T his means that adjusting changes pass band of the filter, then making different components passing through the filter, showing different random processes in the filtered signal. PAGE 59 59 H owever, accuracy of AVAR at every averaging time depends on every indivi dual case. T he degree of confidence of AVAR determined by a sequence at a given averaging time depends on both lengths of sequence and of groups. A longer original sequence gives higher accuracy of AVAR. O nce the length of entire original sequence was fixe d, groups with shorter averaging time have greater reliability. I n practice , lengths of groups have a supremum which is determined by the drift component need to be detected that changes most slowly. F or instance, at least one group with 10 seconds averagi ng time is required if AVAR is used to identify a drift whose character istic time period is around 10 seconds. T hus, the most slowly changing drift to be detected and its corresponding degree of confidence required give a lower bound of length of entire si gnal. However, length of sequence means computation burden, memory and time cost, a proper length of sequence depends on both required accuracy and time range of drift in application is needed. T he error percentage estimat or of AVAR with length N at group length M is [ 104 ] , Figure 5 1 . Error percentage of Allan Variance vs. ratio of N and M . PAGE 60 60 O bviously the percentage above is a function of ratio N M . To clearly demonstrate the relationship between N M and error percentage, we ha ve Figure 5 1 . Obviously this error percentage grows as M grows from 1 to its maximum ( N 1) 2 given a certain sequence length N . By Equation 5 4 , the supremum number of groups in sequence which guarantees certain error percentage % error is, W e require an error percentage lower than 10% in experiment, that is, % error 0.1 . By the equation above, we have N M .To determine proper value of M whose degree of confidence reaches certain upper limit, two aspects should be considered. On one hand, an error component whose correlation time is longer than a certain time length M 0 0 may not be detected by A VAR at group length M < M 0 , and on the other hand, the supremum of erro r component averaging time to be detected should at least be the same length of integration time range. Thus we set value of whose degree of confidence reaches certain upper limit to be the same length of integration time range. I n experiment, gyroscope output was usually integrated in 1 minute, thus a sequence of at least 51 minutes are required to obtain an error less than 10% for 1 minute process. I n the calibration for offset and sensitivity, the output of gyroscope was tested for 1 hour for each task. Procedure Outline of Online Error Reduction T he signal from gyroscope is used both directly and integrated in calibration of accelerometers and further measurement to obtain angular velocity and angle. T he short term or high frequ ency components in noise or error of gyroscope output may have negative effect on direct use of the signal, and the relatively longer term or lower PAGE 61 61 frequency components in error as drift may cause deviation in result of angle by integrating of signal. T hus , it is necessary to reduce or eliminate error with less loss of information carried by original signal. The signal V from gyroscope can be expressed as sum of two parts in different information domains, where V m is signal part carrying specific information of each meas uring task, n is the noise and drift part whose statistic character depends on the gyroscope itself instead of each individual measurement task . D ue to the existence of bandwidth of sensor as well as nature of motion measured by gyroscope, the frequency of V m is usually relatively low, and higher frequency components in output signal are dominated by noise which ha s less dependence o n macroscopic motion. Fig ure 5 2. S chematic power density frequency distribution of information signal and noise/drift signal in output of gyroscope. T he purpose of online error compensation is to reduce the noise components in signal as much as possible, while to keep the components carrying useful information for further analysis. A n intuitional thought is using frequency characteristic of both PAGE 62 62 information and noise signals to tear them apart by applying a proper designed frequency domain filter. U nfortunately , no obvious frequency boundary distinguishing information and noise signal exists in output of gyroscope, the V m part in Equation 5 6 usually overlaps with the n part in frequency domain as shown in Figure 5 2 . Thus it is not applicable of using a frequency domain filter directly in this case. T o solve the current problem where components in signal have overlap in frequency domain, some characteristics of the noise part which are distinctively different from those of the information part in original signal should be firstly noticed. O n one hand, the ADEV curve of gyroscope measured in static condition is absolutely governed by noise and drift. On the other hand, the noise signal in low domain can be considered as white noise, while the drift in lower frequency domain influence primarily the offset of gyroscope which can be s imply compensated by the i ntegral c orrection m ethod mentioned previously. To apply the characteristics above, the Equation 5 6 can be rewritten in higher and lower frequency domains as, where subscript L denotes lower frequency (LF) part, H higher frequency (HF) part. Based on the Equation 5 7 , a complementary filter was designed and applied to gyroscope signal. The filter f irstly separate signal into LF and HF parts , both contain corresponding V m and n signals. S ince the LF part is contaminated by some long term drift, while the HF part also contains some components carrying fa st changing information, filter s for both part s individually are needed to give V m L and V m H to recover both components bearing information. F inally find summation of filtered sequences of both parts as the recreated, complementary filtered signal . PAGE 63 63 U nfortunately, for the LF part, statistic characters and frequency characters of both information signal and drift in this part were unknown prior to data processing. T hus it is not that easy to remove drift signal from original sequence using either a real time statistical filter or a frequency domain filte r . I nstead, as mentioned previously, t he long term drift including temperature effect or power supply floating was detected and compensated intermittently by the integral correction way mentioned previously in I ntegral correction Method in this research . On the contrary , for the HF part in original signal , usually it is assumed that noise component can be treated as white noise whose statistic character is different from that of fast changing information. H oweve r, the frequency band of information may have overlap with the band of noise , making the filters based on frequency domain manipulation in applicable here. B ut it is possible to separate signal carrying information in this part from noise using filters base d on statistical principle, for example, Kalman filter, Wiener filter, etc. Figure 5 3 . Complementary filter structure for gyroscope signal processing. PAGE 64 64 N ow we have the outline of an online complementary filter for error of output of gyroscope. F irstly, a low pass filter is used to do the separation , obtaining LF and HF parts of original signal dominated by motion information and noise respectively. Then for the two parts, a statistical filt er and integral correction method were applied to the one changing faster and slower respectively. The obtained information from the two parts was finally combined again to recover the useful signal. T o clearly demonstrate the complementary filter process, we have Figure 5 3 . D etermination of cutoff frequency for complementary filter So far, determining of the cutoff frequency of the low pass filter to separate the original signal , that is, the frequency boundary between LF and HF parts, is needed . T o find a proper frequency threshold to set gyroscope signal ap art, ADEV method was applied here. Since ADEV shows a spectrum of standard deviation of signal, by using ADEV of signals of various type of motions, the standard deviation components that do not depend on specified motion, that is, the noise independent on individual measured signal, can be determined, corresponding to the HF part to be separated . T o this end, ADEV of 20 tests were calculated and plotted to one fi gure, shown in Figure 5 4 . W ithin the 20 tests, one wa s obtained under static condition, showing the lowest bound of ADEV, and another one is obtained under dynamic condition in which maximum permitted angular rate ( 90 degrees second for ADXRS614) was applied to the gyroscope. T he ADEV curves of both bound ary cases are marked by two arrows here. Besides the two ADEV curves above, the other 18 curves correspond to 18 arbitrarily implemented rotational motion sequences. C learly i n Figure 5 4 , all curves converge at the left part of figure to the curve section which has slope 1 2 at averaging time 0.01 seconds and less. I t shows that all PAGE 65 65 components with averaging time < 0.01 seconds, that is, frequency gre ater than 100 Hz, are independent on the motion in test, but determined by gyroscope itself. Actually this section of curve corresponds to angle random walk which is created by integration of white noise process of angular rate in time [ 10 0] . B y datasheet of ADXRS614, the gyroscope has a bandwidth of about 80 Hz, thus the part of ADEV whose averaging time < 0.01 seconds could be treated as pure noise carrying no information of motion, In other words, components above 100 Hz could be safely consi dered noise only. Figure 5 4 . ADEV curve s of gyroscope signals in 20 tests . T he curve sections at averaging time 0.01 seconds and above give the correlated noise and sinusoidal parts of ADEV. H ere the both parts have relationship with certain motion measured and output by gyroscope [ 10 4 ] . T he error parts of bias instability, angular rate random walk and rate ramp are not shown in Figure 5 4 since the test time length is not sufficient for demonstration of the two components. PAGE 66 66 It can be seen clearly from the curve of static case that sinusoidal component shows up between 0.2~4 seconds of averaging time. T his sinusoidal part may due to the inaccurate gain matching in close loop circuit [ 10 5 ]. F or the cases in which the inputs to gyroscope were not zero, ADEV of the lowest five curve s (the five most slowly changing cases) were still influenced by the mismatch of circuit, as well as standard deviation by faster motion of the other cases covered that sinusoidal component. T he sinusoi dal component in ADEV curve always shows up in corresponding time domain sequence as periodic waves. T o determine how to deal with this sinusoidal part, the corresponding frequency in time domain sequence should be obtained. F or this error part, we hav e the sinusoidal ADEV component magnitude described as [ 104 ] , where 0 is sinusoidal error coefficient, is averaging time, and f the corresponding frequency in original signal. F rom the static case curve, we have the first valley at 0.4 s econd , corresponding to the first time sin( f ) = 0 . Thus we have f 2.5 Hz . T his periodic component can be clearly observed in original sign al shown in Figure 5 5 , Figure 5 5 . Periodic component ( 2.5 Hz) in gyroscope signal . PAGE 67 67 In principle, since this is an error produced by the gyroscope itself, it could not be manipulated by integral correction method for LF components carrying information, it should be taken into account to the HF part to be dealt with by statistical filter . However, on the one hand, this component of error is not just noise following statistic principle, but the one that has certain time domain shape and can be reduced directly by a frequency domain band stop filter, on the other hand, usually a statistical filter requires its input to be stationary, having no obvious periodic components . T herefore the sinusoidal should b e filtered into the part with information, and then be eliminated by a frequency domain band stop filter. H ere, the right bound of cutoff frequency of low pass filter is determined. H owever, the capability of eliminating drift in LF part was limited because of constraint composed by the integral correction method that it could be deployed intermittently only instead of continuous effect of drift reducing . T hus t he cutoff frequency should be relatively low. A lso there wa s no significant periodic component detected at a frequency > 2.5 Hz obtained previously. A ll curve sections with +1 2 slope in Figure 5 4 were correlated components due to the input motion to gyroscope [ 105 ], and no significant tr end term exist in static curve. T hus we can safely assume that apart from the correlated components due to the rotation applied to gyroscope, the curve sections with averaging time < 0.1 second have no periodic and trend term, and may satisfy the requireme nt for statistical filter, making them applicable to be filtered. E ven further, the cutoff frequency can be selected around the averaging time range of 0.1~0.4 second. T o fully utilize the advantage of statistical filter s in high frequency domain and avoid the disadvantage of integral correction method in low frequency domain, the cutoff frequency was finally set at around 0.2 second, that is, about 5 Hz. PAGE 68 68 Low Pass Filter for Obtaining LF and HF Components Since the output of gyroscope is heavily impaired by measurement and process noise, a complementary method was put forward previously to separate the original gyro scope signal into two parts , and then apply different filters to both parts respectively to suppress the erro r to the largest extent possib le in order to use the data for comfort evaluation. To do the separation , two commonly used digital filters were to be selected : the Butterworth filter and the Savitzky Golay filter . Both filters have the advantage of computational simplicity and are well suited to represent time varying features in acceleration measurements. Butterworth Filter The Butterworth filter is a type of signal processing filter widely applied to experimental data in a broad range of discipli nes by low pass filtering. I t is a typical frequency domain filter operating based on the difference between band of information and the one of noise. The Butterworth filter is designed to have as flat a frequency response as possible in the pass band , and is usually referred to as a maximally flat magnitude filter. The Butterworth filter was first proposed in 1930 by the British engineer and physicist Stephen Butterworth [ 106 ] . T he frequency response (gain) G( ) of amplitude of an n th low pass Butterworth filter is, where c is cutoff frequency. T he advantage of Butterworth filter includes the extreme flat response in pass band and good integrative performance. I t has bett er impulse response than that of Chebyshev filter, and better attenuation of pass band than the one of Bessel filter. PAGE 69 69 However, the Butterworth filter is a frequency filter which works primarily in frequency domain. I f the frequency band of noise or the one of information is not quite clear or not fixed, or even if only statistic character istic of noise or information is provided , the frequency filter such as Butterworth filter may not be applicable here, since the pur e cutoff frequency may cut some useful frequency components off without considering any statistic characteristic of the i nput to filter , that is, the original signal . I n these cases, some statistical concepts should be introduced in filter application. T o implement low pass filtering or smoothing for this field considering the statistic characteristic of signal , least square method is introduced and a commonly used smoother based on this method is the Savitzky Golay filter. Savitzky Golay filter The Savit zky Golay filter is a special type of low pass FIR filter which is usually applied in data smoothing. T oday it is adopted in many engineering fields [ 107 110 ]. Unlike some commonly used filters, Savitzky Golay filter does not have previously defined exact properties in frequency domain. Instead, it processes the signal sequence in time domain directly, based on theory of least squares method. Savitzky Golay filter was popularized by Abraham Savitzky and Marcel J. E. Golay in 1964. They published tables of convolution coefficients which significantly reduce computation cost for polynomial coefficients in the filter and thus left it with better speed performanc e [ 111 ] . T o smooth data sequenc e, value at each time point is replaced by a new one which is obtained by a predefined function of some measured values within certain time range around the point to be replaced. S ince true signal value does not change rapidly in typical measurement case, replacement base on a correctly defined function would not introduce much deviation from true value but increase S ignal to N oise R atio (SNR). PAGE 70 70 F or Savitzky Golay filter, a polynomial with predefined order is applied as the function. Coefficients of polynomial are determined by look up tables of convolut ion coefficients or matrix method using all measured values in each time range, that is, the time window. T he length of time window is also predefined as a parameter of the filter. Now assume a data sequence x ( i ) where i = m, m +1, , 1, 0, 1, , m 1, m , a nd an n order polynomial f i = a nk i k is used to fit the 2m+1 points. T o this end, we should determine all a nk to minimize the quadratic sum of error, T o calculate values of a nk , the table of convolution coefficients initially created by Savitzky and Golay is used. S ince we have the following equations given that the quadratic sum of error above is minimized [ 109 ], F or each k , the value of the k th derivative of f i a t i = 0 can be obtained using the 2m+1 points and the table with corresponding order of derivative and order of polynomial provided by Savitzky and Golay . Then a nk can be determined by dividing d k f i d i k a t i = 0 by factorial of k . O nce all coefficients of polynomial are obtained, the value at time point 0, that is, x (0) which lies at the middle of current time window, is then replaced by the value of polynomial at that time point. T hen the time window step one point forward, making the current x (1) to be the new x (0) , and repeat the pr ocess of fitting and replacing. N ote that for each time window, the values from x ( m ) to x ( 1) used for polynomial fitting are the original signal values, not the previously replaced values. PAGE 71 71 H owever, the size of table is limited. T he table provided by Savi tzky and Golay gives coefficients of polynomial up to 6 th order, by up to 25 points each time window [ 111 ]. Obviously it is far from sufficiency since the length of sequence and time resolution of signal is relatively high today, and the length needed of time window is usually much larger than 25. T hus the way to calculate each d k f i d i k at i = 0 online is required here. A ssume again a sequence g with k points, and an other n order polynomial f ( x ) = a n p x p to fit the sequence g . I ntroduce the following matrix B , D efine, Obviously we have Ba = g , thus a = ( B T B ) 1 B T g . S ince the c th derivative of f ( x ) is , W e have the following equation for d c f ( x ) d x c , where, O nce the derivatives d k f ( x ) d x k are obtained as shown above, the Equation 5 1 1 can be applied here to find all co efficients of the fitting polynomial. PAGE 72 72 T he filtering or smoothing performances of the Butterworth filter and the Savitzky Golay filter are not quite distinctive to each other. H owever, t he latter one was selected as the filter to separate original signal in this research . T he reasons are described as follows . F irstly, the advantages of Savitzky Golay filter include smoothing the signal without phase shift, as well as better mean square de viation performance compared with Butterworth filter or Chebyshev filter, making the signal after filtering closer to true value sequence, thus guarantee the accuracy and performance of following fitting, calibration and measurement process [ 112 ] . Secondly, the Savitzky Golay filter is essentially a smoother using generalized moving averag ing process . F or the signal to be filtered by statistical filter in noisy part, it is required to be zero mean, stationary and normal. T o satisfy the zero mean req uirement , usually average values of signal were calculated and subtracted from data to obtain satisfactory data. T his is automatically fulfilled by the Savitzky Golay filter while the Butterworth filter does not have the averaging nature. F urthermore, the Savitzky Golay filter helps to preserve the peaks while have a good smoothing performance compared with Butterworth filter for data that have abrupt changes [ 107 ] . I t was found that the Savitzky Golay filter showed better performance to deal with shocks or peaks in signal as low pass filter than the Butterworth filter [ 113 ] . B y observing the data from gyroscope, the vast amount of random peaks or sudden changes are remarkable, making it a large signal floating relative to the average (true) value. F or this case, it is appropriate to use Savitzky Golay filter to gain the approximate average signal. Savitzky Golay filter deals with signal in time domain directly. H owever, since it is applied as low pass smoother, it has implicit frequency character istic , which can be PAGE 73 73 recovered by comparing PSD of original signal and that of filtered signal. Given order of polynomial N and time window length M , t he normalized cutoff frequency f c at 3dB of a Savitzky Golay filter is [ 112 ] , T he order of polynomial was set to N = 4 , and sampling rate of signal was 1kS sec. To obtain a cutoff frequency of 3dB at around 5 Hz, that is, normalized frequency at around 0.005 Hz, the length of time window should be M = (5 0.005 + 4.6) 3.2 313 . In experiment, this length was set to 301, making a 5.2 Hz cutoff frequency which is acceptable here. T o validate the cutoff frequency of the Savitzky Golay filter by its determined parameters, the PSD of original and filtered signal were calculated and plotted in Figure 5 6 . Obviously the curves share an almost same spectrum part at frequency < 6 Hz, while the components beyond 6 Hz were significantly reduced, verifying the applicability of Savitzky Golay filter for this gyroscope. Figure 5 6 . PSD of original gyroscope signal vs . PSD of Savitzky Golay f iltered gyroscope signal . PAGE 74 74 Original HF Signal Preconditioning T o eliminate noise in HF part of gyroscope output obtained by the Savitzky Golay filter selected previously, three steps are implemented. The first step is to precondition the original signal from gyroscope so that the proce ssed HF data satisfies requirement of further procedure. Then an error model for signal needs to be established for err or estimation. Finally the established error model is combined with preconditioned signal to obtain compensated output by a selected type of filter or algorithm . Once collected from sensor and given by the Savitzky Golay filter , the HF part of original sig nal should be tested then and manipulated if necessary to provide a normal, zero average, stationar y input sequence for error model identification and filter. F or a static output sequence of gyroscope, requirement of zero average can be obtained by simply subtract ing average constant sequence from original signal. T o confirm the stationar ity , the periodic and long term trend components should be first detected and then eliminated from original signal. The periodic components were discovered by checking PSD . T he PSD of static signal from gyroscope is shown in Figure 5 6 in dotted line . Obviously there is an outstanding peak marked by an arrow at frequency 2.3~2.8 Hz. This peak was detected in every test of gyroscope and should be fil ter ed by a filter before identifying model parameters. A band stop Butterworth filter was applied here to obtain signal sequence which has almost no significant periodic components. S ince the HF part of gyroscope output whose stationar ity is being tested here is provided as the result of the output from the 4 th order Savitzky Golay filter with 301 side sequence length mentioned previously being subtracted from the original signal, long term trend components should be contained in the output of Savit zky Golay filter and erased from the HF part signal. T herefore, there should be no long term trend PAGE 75 75 in the current case. T o confirm this, a 4 th order Savitzky Golay filter with large side sequence length (3001 points) was applied to the HF part of gyroscope signal . As expected, n o trend was discovered . A fter identifying and eliminating the periodic and long term trend components in original signal, stationar ity of processed sequence was tested as follows . T o test stationar ity of discrete series, there are several methods including data sequence observation method, autocorrelation observation method, Dickey Fuller test (DF), A ugmented Dickey Fuller test (ADF), Phillips Perron test (PP), reverse test and runs test. H ere DF test, ADF test and PP test can also be categorized as unit root test. The d ata sequence observation method is plotting the data sequence and see if there is periodic or trend component in the signal. Obviously it is direct but inaccurate . The a utocorrelation o bservation method is to plot the autocorrelation sequence of the original signal, and see if the autocorrelation values damps slowly or periodically. O nly quickly damping autocorrelation sequences means stationary original signal. B ut since it is not easy to remove all periodic components in signal practically, making it difficult to absolutely reduce the periodic wave in autocorrelation sequence, a slight vibration part in signal may cause difficulty in determining if it is really stationary in practice . U nit root test is widely used in confirming stationarity . H owever, to apply the unit root test, constant and linear trend terms as well as lag number of the signal to be tested should be estimated prior to the test. Thus it is relatively complex and difficu lt to apply the unit root test in practice. R everse test is applicable in testing for monotone trend only. T he runs test judges the stationarity of a data sequence using values in the sequence only . No previous estimation of signal parameter is required. I t also provides PAGE 76 76 much more reliable result than that obtained by data sequence observation method, autocorrelation observation method and reverse test. I n conclusion , runs test was applied here and implemented as follows [114] . 1. Separate the original sequen ce to be tested into N groups, 2. RMS of each group was calculated and then we have an RMS sequence with N points, 3. F ind the median of RMS sequence and set the threshold RMS threshold = RMS median + where is a positive number, 4. C lassify each value in RMS sequence. I f the value is greater than RMS threshold , it is marked + , or else it is marked . Finally we have a sequence of plus and minus signs , 5. S can the plus and minus signs sequence and count number of runs. Consecutive same signs between the other signals make one run. For instance, a sequence , here three + signs make one run, 6. O nce total number of runs r is obtained , check the look up table using r and N or use p robability distribution calculation to see if r lies in certain acceptance domain of N . I f r lies in the domain, the original sequence is stationary. I n application, usually we have look up runs test table including N up to 200 runs. However, in our cases, number of points in the sequence to be tested was large (3.6 10 5 ), and the length of each group was set to 100, that is, we have N = 3600 , which has no corresponding entry in normal look up runs test table. T hus we have to obtain the probabilities for N = 3600 at different confidence level . I f N is suff iciently large, the runs number r follows normal probability distribution N( r , r ) , where, PAGE 77 77 H ere , N + and N are number of plus signs and minus signs respectively . Obviously we have N + + N = N . S ince the threshold to classify sign of each group is the median, and N is sufficiently large, usually we could assume that approximately N + N N , and then we have simplified expectation and variance, T hus we have statistic r which follows N(0,1) , Now we have a look up Table 5 1 for runs test given N = 3600 . G iven the look up table, 3 gyroscope output sequences obtained under static condition were tested using runs test the result were compared with the elements in Table 5 1 . N umber of runs of each sequence is shown i n Table 5 2 . W e reach the conclusion that sequence after periodic terms remov ed is stationary. T o confirm the normality of sequence, there are also several methods applicable here, in cluding Pearson's 2 tes t, Jarque Bera test, or Kolmogorov Smirnov test. We applied a simple and intuitional way to test if the gyroscope output is normal using a MATLAB Â® function normplot which operates in a traditional way as normal probability paper . A typical result is shown in Figure 5 7 . Tab le 5 1 . L ook up table for runs test given N = 3600. Statistical significance 0.99 0.975 0.95 0.05 0.025 0.01 r 1698 1712 1723 1879 1890 1904 Tab le 5 2 . N umber of runs of gyroscope signals . Test (3.6 10 5 samples each) 1 2 3 R uns 1796 1818 1792 PAGE 78 78 Figure 5 7 . Result of normplot of output of gyroscope . T he output of normplot shows that data sequence fits the line very well . T he gyroscope output follows normal probability distribution . S o far, the requirement of being zero average has been satisfied by removing average consequence. T he normality and stationary have been confirmed by normplot function and runs test respectively. The preconditioned HF part signal is now ready for model identification for further statistical filtering process. ARMA Model for Preconditioned HF Signal T o eliminate the noise in gyroscope signal, the error sequence model should be firstly determined to describe character of error changing in discrete time domain. In practice , a normal, stationary random series can be described by Auto Regressive Moving Average (ARMA) model. A n ARMA( p,q ) model can be expressed as, PAGE 79 79 where a i < 1 ( i = 1,2, , p ) ar e auto regressive coefficients , b j < 1 ( j = 1,2, , q ) are auto regressive moving coefficients, k ~ W (0, 2 ) ( k = n q , n q+ 1 , , n 1 , n ) are white noise sequence which has zero expectation and a variance of 2 . I f q = 0 , ARMA model reduces to p th Auto Regressive (AR) model, that is, I f p = 0 , ARMA mo del reduces to q th Moving Average (MA) model, that is, A B Figure 5 8 . A utocorrelation function and P artial autocorrelation function of gyroscope noise . A) Autocorrelation, B) partial autocorrelation. PAGE 80 80 ARMA model identification is to find proper value of p and q so that ARMA( p,q ) best fits the signal series to be identified. I f the proper model is eithe r AR( p ) or MA( q ) , determination of the appropriate values p or q can be facilitated by plotting the autocorrelation functions for an estimate of q , and the partial autocorrelation functions for an estimate of p . I f the autocorrelation sequence shows a cutoff lag q while the partial autocorrelation sequence do es n t have apparent end but just decays to small values, we have MA( q ) . I f the partial autocorrelation sequence shows a cutoff lag p while the autocorrelation sequence do es n t have apparent end but just decays to small values, we have AR( p ) , otherwise, we ha ve ARMA( p,q ) and values of p and q should be determined by some other ways. W e calculated the autocorrelation function and the partial autocorrelation function shown in Figure 5 8 , O bviously we do not have any cutoff style, implying an ARMA( p,q ) model to fit the signal sequence. To determine the two parameters, there are three commonly used criterions for model identification, Akaike Information Criterion (AIC), Bayesian Information Criterion (BIC) and H annan Q uinn C riterion (HQC). The AIC was applied here using the residual of each model while applied to original real signal sequence to find m aximum likelihood solution of the two parameters. T he AIC value is defined as, where (p,q ) is loss function, d is number of model parameters, and N is the number of values in the estimation data set . T he loss function is defined as, w here (n,p,q) is the residual vector of model ARMA( p,q ) at time n . PAGE 81 81 Figure 5 9 . AIC values of 14 candidates of ARMA, AR, MA models . The order of the model describing the error behavior of gyroscope is usually no greater than 4 b ased on experience . M oreover, since the ARMA model is equivalent to a linear system, the transfer function of the system is usually a rational fraction . I n other words, paramete r p is no less than q in most cases. T herefore, the candidate models include AR(1) ~ AR(4), ARMA(1,1), ARMA(2,1~2), ARMA(3,1~3) and ARMA(4,1~4). T o determine the parameters of the model, the AIC value of each model was calculated and shown in Figure 5 9 . C learly in Figure 5 9 , the model that has the lowest AIC value is the ARMA(3,2) . T hus the parameters p and q are identified to be 3 and 2 respectively. The coefficients in the model were then estimated , giving the specific ARMA(3,2) model as shown in Equation 5 2 5 describing changing of error of gyroscope HF signal , Kalman Filter for Preconditioned HF Signal The Kalman filter is a widely used method of formulating the linear minimum mean square error filtering problem which utilizes state space methods. Th e principle of PAGE 82 82 Kalman filter and Extended Kalman filter is discussed in Kalman Filter (KF) and Extended Kalman Filt er (EKF) . I n this gyroscope case, since the model for state process , that is, the ARMA(3,2) model, is linear, a normal Kalman filter was applied to reduce error in gyroscope signal. T he process state vector was set as, A nd output o f Kalman filter was Y n = x n . W e have the state equation and measurement equation , where we have the expectations of noise terms n and w n , A ccording to the ARMA(3,2) identified previously, we have the matrices, By commonly used form, the following five governing equations for Kalman filter were ap plied to establish recursive calculation process , Estimation of state: Estimation of error covariance: Find current Kalman Gain: C orrection of state based on measurement: C orrection of error covariance: PAGE 83 83 wh ere is the covariance matrix of priori estimation error of state at time step n , is the covariance matrix of posteriori estimation error of state, is the priori estimation of state, while is the poste riori estimation. Result of Kalman Filter and Complementary Filter The Kalman filter proposed here was applied in experiments. T o demonstrate the performance of Kalman filter in reducing noise in gyroscope HF s ignal, the output signal of the Savit zky Golay filter , that is, the LF signal, wa s subtracted from the original gyroscope data and the consequent signal was used directly as input to the filter. O ne typical comparison of the input HF signal and Kalman filtered signal is shown in Figure 5 10 A ) . Both the input and output of the Kalman filter were filtered by a Savitzky Golay filter for clearer demonstration in Fig ure 5 10 B ) . It is obvious that the Kalman filter provides a good reduction of the random noise component in input signal. F or further demonstration of the performance of Kalman filter, t ypical curve of ADEV of original offset signal is presented as the upper curve in Figu re 5 1 1 . C learly the static signal almost represents the HF part of signal obtained in dynamic case, except that it also contains extra but weak LF drift components. Since parameters of one gyroscope are time variance, A DEVs of measurements would seem a little different. A fter the original static signal of the ADEV above was filtered by the Kalman filter established previously, the ADEV of filtered signal is shown as the lower curve in Figure 5 1 1 . Obviously the ADEV of filtered s ignal is apparently lower than the original signal, showing again the performance of the Kalman filter. It is clear that the variance reduction effect provided by the Kalman filter is almost same on the average time domain; the Kalman filter did not change the relative distribution of ADEV on frequency. PAGE 84 84 A B Figure 5 10 . High frequency part of o riginal gyroscope si gnal (in black) and its Kalman filtered s equence (in gray) . A ) T he original signal s , B) the Savitzky Golay filter ed signal s . Figure 5 1 1 . ADEV of original static gyroscope signal vs . ADEV of Kalman filtered static signal . PAGE 85 85 A B C Figure 5 1 2 . Result of Kalman filter to process the gyroscope signal, Case 1 . A ) Original gyroscope signal , B ) f iltered gyroscope signal , C ) ADEV of both signals . PAGE 86 86 A B C Figure 5 1 3 . Result of Kalman filter to process the gyroscope signal, Case 2 . A ) Original gyroscope signal , B ) f iltered gyroscope signal , C ) ADEV of both signals . PAGE 87 87 T he complementary filter which is combination of a Kalman filter, a Savitzky Golay filter and the integral correction method was also applied to measurement of gyroscope. O nce the HF part signal was preconditioned and filtered by the Kalman filter, and the LF part signal was compensated by integral correction method, the two processed parts were then combined to recreate the entire gyroscope signal in which noise parts in both HF and LF part were eliminated to some extent. T wo typical results are shown as fo llows in Figure 5 1 2 and Figure 5 1 3 . O bviously in both cases shown above, the noise components are significantly reduced, giving smoother and more accurate curves. Note that in ADEV sub figures, o nly the variance components located at high frequency domain are reduced, showing the performance of Kalman filter in eliminating the HF noise. However, the normal motion located at lower frequency is not significantly influenced by the complementary filte r , showing the maintenance of useful information determined by gyroscope while the signal was passing through the filter . Since the two cases shown in Figure 5 12 and Figure 5 13 were both obtained in minutes in which the drift effect is not sign ificant, the integral correction method applied in dealing with LF part shows no apparent effect in long term drift reduction. PAGE 88 88 CHAPTER 6 EXPERIMENT OF ACCELEROMETERS P OSITIONING The experiment was set with two dual axis linear MEMES accelerometers and one single axis MEMS gyroscope. T he accelerometers whose positions were to be determined were attached to a virtual body which was determined by an inverted pendulum, and the gyroscope was attached to the inverted pendulum. T he position of accelerometers could be adjusted relative to the pen d ulum or the base joint, making different cases for positioning process. H owever, the position of gyroscope does not have influence on out experiment result. I nt erfacing of the analog voltage signal was done using NI DAQmx and LabVIEW . A magnetic encoder was fixed at the base joint of the pendulum. S ignal s from two accelerometers (4 channels), one gyroscope (1 channel) and one encoder (1 channel) were sent to a laptop for data processing via two NI USB6009 ports using difference connection . T he sensed data were obtained and recorded by LabVIEW . The data were then processed by Matlab . All algorithm s, digital filters applied in this research were implemented in .m file. T he architecture of the experiment equipment combination is shown in Figure 6 1 . Fig ure 6 1 . Architecture and data flow in experiment of accelerometers positioning. PAGE 89 89 Analytical Result W e have obtained the analytical solutions Equation 2 1 3 , 2 1 4 , 2 1 5 , 2 18 and 2 19 of four distance components r x,y and l x,y in global coordinate system as well as the angular acceleration given the readings from four axes of calibrated accelerometers A1 and A2 and the filtered gyroscope signal. For cleara nce, the five solutions are rewritten here in a sequential form which can be easily implemented in computer programming , where l x,y are firstly obtained, then 1 using l x,y , finally r x,y using 1 . U sing the equations above, several tests were conducted. T he example motion is shown here i n Figure 6 2 . PAGE 90 90 A B C Figure 6 2 . Result of EKF of kinematic and distance parameters . A ) Output of accelerometers , B ) o utput of gyroscopes , C ) c alculated angular acceleration . PAGE 91 91 A B C D Figure 6 3 . Calculated value sequences of four distance parameters r x,y and l x,y . A) r x , B) r y , C) l x , D) l y . PAGE 92 92 A B C D Figure 6 4 . Probability distribution of four distance parameters r x,y and l x,y . A) r x , B) r y , C) l x , D) l y . PAGE 93 93 N ote that the result curves shown in Figure 6 3 were noisy, especially during the static moments. S ince there were no radial and tangential acceleration components in static condition, the difference between readi ngs of the two accelerometers approaches to zero, making the calculated distances nonsense in the moment. I n other words, the static moments provide no information about the accelerometers positions. T o deal with this, a conventional way is use of a static dynamic detector applied to output of sensors to deactivate the algorithm while the system is in static [ 55, 99 ] . H owever, all elements in the calculated sequences of r x,y and l x,y are obtained independently from their corresponding readings of sens ors at that time. It is nature assum ing that the calculated result would be the real value at the most probability. T hus it is possible to find the real value by plotting the probability density figure of the sequences of r x,y and l x,y , the peak of each pl ot lies at the real value of corresponding parameter in maximum likelihood. F or the case shown in Figure 6 3 , the probability densities of r x,y and l x,y were calculated and shown in Figure 6 4 . F rom the result above, the most possible estimated values of the four distance components are, B y this, we have the estimation of distance between the two accelerometers as well as the distance between each accelerometer and the joint. PAGE 94 94 I n fact, the distance between accelerometers A1 and A2 was 0.440 m, the estimation error was 0.6%. The distance between accelerometer A1 and joint was 0.652 m, the estimation error was +4.9%. The distance between accelerometer A2 and joint was 0.736 m, the estimation error was 5.6%. S everal tests were co nducted in another sensor configuration which is different from the one tested above. T he calculated results using the analytical solutions and probability density method are shown in the Table 6 1 , including the estimated length of each radial vector as well as the relative deviation of the estimator from the true value . F or the third configuration with changed orientation of both accelerometers, several tests were also conducted. The results are shown in Table 6 2 . H owever, the application of probability density is not that convenient in practice, since it is less accurate due to the fuzziness in identification for peak in probability density plot, and time consuming due to the fact that for each probabil ity density plot, a sufficiently long signal sequence is required, making it quite slow in updating rate if the sensor is moving relative to the object. Tab le 6 1 . Results of the analytical solutions, probability density method, Case 1 . Real distances Estimated in Test1 0.504 m +2.2% 0.640 m 3.6% 0.438 m 7.7% Estimated in Test2 0.536 m +8.6% 0.669 m +0.8% 0.423 m 10.9% Estimated in Test3 0.521 m +5.7% 0.675 m +1.6% 0.442 m 6.9% Estimated in Test4 0.525 m +6.5% 0.668 m 0.6% 0.436 m 8.2% Tab le 6 2 . Results of the analytical solutions, probability density method, Case 2 . Real distances Estimated in Test1 0.508 m +3.1% 0.659 m +1.0% 0.470 m 3.4% Estimated in Test2 0.516 m +4.7% 0.666 m +2.2% 0.462 m 5.1% Estimated in Test3 0.528 m +7.2% 0.647 m 0.7% 0.473 m 2.9% Estimated in Test4 0.502 m +1.9% 0.670 m +2.8% 0.452 m 7.2% PAGE 95 95 N ote that we have a probability distribution character for each distance parameter solved in the analytical way mentioned previously. I t is possible to filter the calculated parameters if statistical filter is applied here. Statistic al Filters Discrete Wiener Filter S ince it is discrete case here , the continuous Wiener filter is not discussed here. The Wiener filter for discrete sequence is a weighting function approach. T he filter weights all past data trying to find the best estimation of curr ent value. G iven an n element noisy input sequence x i , i = 1,2, ,n , the optimal estimation of y n is, T hen the mean square error is, T o determine c oefficients a i to minimize the error , the partial derivative s with respect to each a i are calculated and should be all equal to zero. G iven that all of the autocorrelation and crosscorrelation functions of y n and past data x i are known , the c oefficients vector a = [ a 1 a 2 a n ] T can be obtained using matrix algebra [ 6 ], then the c oefficients vector gives the best estimation proposed in Equation 6 6 . Kalman Filter (KF) and Extended Kalman Filter (EKF) I n 1960, Kalman published his celebrated paper in which discrete data sequence was filtered in a recursive way to estimate the state vector in process. F rom then on, benefited from development of digital computation technologies, Kalman filte r is wide spreading in a lot of fields as a subject of research and application, especially in the field of inclinometer or navigation system based on inertial sensors. Kalman filter is another type of statistical filter considered in this research. The op timality of this filter is PAGE 96 96 based on Bayes theorem and the use of conditional probability density functions [ 6 ] . It operates on the base that the spectral characteristics of the processes involved and the statistic characteristics of error are previously known. Normally all noise are considered or should be preprocessed to be white and normal before processing by Kalman filter. T wo governing equations used in a Kalman filter are the state equation and the measurement equation. The state equation describes the changing of state vector in process involved. I t is the base of prediction of state. The measurement equation shows the relationship between the state vector and the measured signal obtained by observer outside the process sys tem. It is the background of correction to the prediction made by the state equation earlier in current time step. T he equation of state is, w here X n and X n 1 are the state ve ctor of the system at time step n and n 1 respectively. u n 1 is the input to system at n 1 , n 1 is the process noise vector, A and B are constant state transition matrix and input matrix describing projection of input onto the state. T he measurement equation is, w here Y n is the measurement vector obtained at time step n . H is the measurement matrix rel ates state and observed signal. w n 1 is the measurement noise vector. T wo covariance matrices R and Q of measurement error and process error respectively are introduced in Kalman filter to control the change of Kalman gain and estimation error covariance in recursive process, which are, PAGE 97 97 Note that the two covariance matrices are actually time various . However, for a stationary process, they are usually considered constant. T he covariance of posteriori estimation error of the state at time step n is, A Kalman gain K n is applie d to obtain discrete linear posteriori estimation of state, where is the priori estimation of state obtained by prediction of state equation. I n Equation 6 12 , the factor gives the difference between real measurement and the expected measurement. A nd the term plays as a rule of updater to correct the priori estimation of the state. I n the updating equation above, Kalman gain is obtained in the following way , T he equation for Kalman gain at each time step was derived in order to minimize the covariance of posteriori estimation error of the state mentioned in Equation 6 11 . T he recursive process of a Kalman filter is shown in Figure 6 5 . Figure 6 5 . Outline of Kalman filter . PAGE 98 98 H owever, the state updating process and the relationship between state and measurement are both linear in the normal Kalman filter described above. I n practice, we may have a nonlinear process involved, where either the dynamic of state or measurement relat ions may be nonlinear . T hus, the state equation and the measurement equation may have the nonlinear form as, where either or both f and h are nonlinear functions. T o deal with this situation, a variant of Kalman filter which is called Extended Kalman filter (EKF) was developed. Both nonlinear functi ons f and h may be linearized by taking partial derivatives of the nonlinear functions to obtain the linear form utilized by normal Kalman filter . Obviously once partial derivatives are conducted, the obtained transition or measurement matrices will no lon ger be constant and must be updated at each time step. I f the input to process is not taken in consideration, we have the time variant transition and measurement matrices at time step n as, where is the posteriori estimation of state X at time step n 1 . Tilde s above vectors denote approximation of the vector without considering noise or w . S ubscripts [ i ] and [ j ] denotes the i th and j th element in state X or the i th row and j th column in matrix A or H . O nce the matrices A and H are obtained by partial derivatives , the other part of processing in EKF is the same as that of normal Kalman filter. D ue to the nonlinear equations in process, optimality of Kalman filter lost in EKF, since the normality of random signal sequence does not hold once past through a nonlinear system. This PAGE 99 99 disadvantage makes it risky of applying EKF in some cases. H owever, for most cases, EKF would lead to the suboptimum solution which is not far away from the best one. One advantage of Kalman filter compared to the previous probability method is that the filter could be easily implemented in computer programming, ho wever, there may be difficult and ambiguous by identifying peaks using probability density method. A nother advantage of Kalman filter is that it is recursive filter, while the Wiener filter as well as the probability density method is non recursive method. Unlike the discrete time Wiener filter discussed previously, recursive processing of Kalman filter allows an updated estimate to be made using only the results from the previous estimate. B y this, Kalman filter can be applied online with a lower memory re quirement. Once the measurement process starts, only a short initial period in the measurement is needed for convergence of the filter. T he filter would have an updating output of parameters as soon as the filter is converged. Unlike the probability densit y method and the Wiener filter, the updating rate of Kalman filter could be usually as fast as sampling rate. Positioning using Extended Kalman Filter T he first step is to establish state system. T he state vector was selected as, There was no input signal to the EKF , t hus the state update equation was, where, PAGE 100 100 T hen to establish the measurement equation, a ssume that the joint is static, we have already obtained first four of the following measurement equations by the analysis mentioned previously, F or the measurement matrix equation used in EKF, Jacobian matrix in EKF was applied here since the relationship between the state and the measurement vector was nonlinear. T he measurement equation was, M oreover, the covariance matrix R of measurement noise vector W n was obtained by calculating the covariance of residual signals corresponding to elements in state vector while average values were removed from measurement, PAGE 101 101 And the covariance matri x Q of process noise vector was adjusted to be, U sing the matrices above, the EKF was implemented by computer programming, and measured signal sequences were applied to the filter for obtaining converged result of the state vector. T ake the test sequence shown in analytical algorithm for example again for comparing purpose; we have the filtered elemen ts in state vector as follows, Figure 6 6 . EKF filtered angular velocity signal . Figure 6 7 . EKF filtered angular acceleration signal . PAGE 102 102 Figure 6 8 . EKF filtered values of four distance parameters r x,y and l x,y . C ompare the filtered curves in Figure 6 6 and Figure 6 7 with the curves shown in Figure 6 2 B ) and Figure 6 2 C ) respectively, obviously the EKF filtered curve are much more s mooth and less noisy. N ote that all four curves representing the four distance components to be found in Fig ure 6 8 converge after a period of vibration. T he obtained values of the four parameters in this test are, Again, we have the estimation of distance between the two accelerometers as well as the distance between each acce lerometer and the joint. PAGE 103 103 The distance between accelerometers A1 and A2 was 0.440 m, the estimation error was 6.6%. The distance between acceleromete r A1 and joint was 0.652 m, the estimation error was 7.8%. The distance between accelerometer A2 and joint was 0.736 m, the estimation error was 7.1%. F rom the EKF converged result obtained above, it is clear that the output of EKF converged successfully to the real positions of the accelerometers. T he validity of positioning using EKF was thus confirmed. All estimated position parameters fell in 10% deviation range from the true values. H owever, the result around 7% is not a very accurate solution compar ed to the current way using ruler or other measuring equipment. T he governing problem so far is that the relatively less accuracy of inertial sensors compar ed to the relatively small dimension of the inverted pendulum deployed in this experiment does not t hat match. T he small dimension, that is, the relatively near of the accelerometers to the joint, gave relatively small extra acceleration signal in output of accelerometers. T he extra acceleration is not that large so that may cause relatively large error in estimation of locations. T hus, the positioning algorithm should perform better in the case of larger dimension of link, for example, cranes, industry robots. Since the maximum angular rate in our research was limited by the range of gyroscope, better ac curacy may be obtained if a more capable gyroscope is deployed in application, because larger angular rate gives larger extra acceleration exerted on accelerometers , making a better SNR for the algorithm. PAGE 104 104 CHAPTER 7 C ONCLUSION AND FU TURE WORK The position ing algorithm presented here is a method aiming at reducing reliance of in ertial inclinometers on the pre determined information of positions of accelerometers. The combination of inertial sensors obtains the distance parameters of the accelerometer from the joint center by themselves . O nce the positions of all accelerometers are determined automatically by this way, the obtained position information as well as the filtered angular acceleration, angular rate can be applied in further application of inclino meter. The method gives both analytical solution and statistic solution of the positions to be determined . F or the analytical solutions, a probability distribution method was applied to find parameters. For the statistic solution, an E KF converges estimat ions of parameters directly to their values during online measurement. I n experiment, the converging process costs 1~2 minutes , which is usually acceptable in real application. T here is no special configuration needed for the positioning process. T he simpl est combination which can provide the positions is two dual axis MEMS linear accelerometers and one single axis MEMS gyroscope. A number of experiments were conducted using a pair of MEMS dual axis accelerometers and one MEMS single axis gyroscope attached to a link with a rotating joint to validate the performance of positioning method proposed here. The results show that the algorithm is really applicable in practice to obtain the distances from accelerometers to the joint, and reach the degree of positio ning error less than 10%. The positioning algorithm for inertial inclin ometers is applicable in types of circumstances without special requirement such as lower magnitude of acceleration or PAGE 105 105 specified position of sensors. O nce both types of sensors are cal ibrated, they can be attached directly to anywhere on the surface of object, and provides necessary information by the combination of sensors itself. D ue to lower cost and better performance of MEMS sensors, the applications of inertial based inclinometer extend to various circumstances, making a wider potential application field for positioning algorithm. Potential cases in which the self positioning method may be required including human body tracking, gait analysis, robot and manipulator working in aggre ssive condition, remote control, etc. The formulation and experiment done in this research was limited for one rotation joint case in which the link moves in vertical plane only. T he future work requires to attempt to extend the positioning algorithm to 3 D motion cases where spatial joints are used or to multi link manipulator cases where the method can be applied in cascade way. H owever, the accuracy of the solution so far is not quite satisfactory, making room for improving the precision. Once the positioning process gains sufficiently high accuracy, i t is also desired to apply the method to detect slow motion of sensors relative to its object using appropriate updating strategy to refresh information of po sitions . PAGE 106 106 LIST OF REFERENCES [1] S. Billat, H. Glosch, M. Kunze, F. Hedrich et al ., based micromachined in The 14th IEEE International Conference on Micro Electro Mechanical Systems , 2001, pp. 159 161. [2] D . Lapadatu , S . Habibi, B . Reppen, G . Salomonsen et al . , axes capacitive inclinometer/low in IEEE International Conference on Micro Electro Mechanical Systems , 2001, pp. 34 37. [3] R.A. Yotter, R.R. Baxter, S. Ohno , S.D . Hawley et al . , in 12th International Conference on Solid State Sensors, Actuators and Microsystems , vol. 2 , 2003, pp. 1279 1282. [4] R. Poppe , , Computer Vision and Image Understanding , vol. 108 , no. 1 2, pp. 4 18, 2007. [5] G. Baselli, G. Legnani, P. Francoc, F. Brognoli et al . , , Journal of Biomechanics , vol. 34 , no. 6, pp.821 826, 2001. [6] E.R. Bach mann , , Ph.D. Dissertation, N aval P ostgraduate S chool , 2000. [7] F . Alonge, E . Cucco, and F . , and , in Preceedings of 2013 IEEE International Conference on Mechatronics and Automation , 2013, pp. 939 944. [8] S . K urata and M . M akzkawa , Joint motion monitoring by accelerometers set at both near sides around the joint , in The 20th Ann ual International Conference of the IEEE Engineering in Medicine and Biology Society , 1998, pp. 1936 1939. [9] B . Najafi , K . Aminian, F . Loew, Y . Blanc et al ., Measurement of stand sit and sit stand transitions using a miniature gyroscope and its application in fall risk evaluation in the elderly , IEEE Transactions on Biomedical Engineering, vol. 49, no. 8, pp. 843 851, 2002. [10] H. Zeng and Y. Zhao, Sensing Movement: Microsensors for Body Motion Measurement , Sensors , vol. 11, no. 1, pp. 638 660, 2011. [11] J. Baerveldt and R . Klang , cost and Low weight Attitude Estimation System , in IEEE International Conference on Intelligent Engineering Systems , 1997 , pp. 391 395 . [12] L . Hertig, D . Schindler , M . Bloesch, C. David Remy et al ., , in I EEE International Conference on Robotics and Automation (ICRA ) , 2013, pp. 2471 2476 . PAGE 107 107 [13] A . J. van den Bogert, L . Read and B . M. Nigg , A method for inverse dynamic analysis using accelerometry , Journal of Biomechanics , vol. 29, no. 7, pp. 949 954, No. 7, 1996. [14] M . El Gohary and J . McNames, Shoulder and Elbow Joint Angle Tracking With Inertial Sensors , IEEE Transactions on Biomedical Engineering , vol. 59, no. 9, pp. 2635 2641, 2012. [15] P. Cheng and B. Oelmann , Joint angle measurement using accelerometers and gyroscopes A survey , IEEE Transactions on Instrumentation and Measurement, vol. 59, no. 2, pp. 404 414, 2010. [16] W . Quan , H . Wan g and D . T. Liu , A Multifunctional Joint Angle Sensor with Measurement Adaptability , Sensors , vol. 13, pp. 15274 15289, Nov. 2013. [17] J . M. Lin and C . H. Lin , Thermal Convection Non Floating Type Inclinometer with Hemispherical Chamber , in Proceedin gs of the 2nd International Conference on Intelligent Technologies and Engineering Systems (ICITES2013) , Lecture Notes in Electrical Engineering 293, 2013, pp. 223 230. [18] P. T. Gibbs and H. Harry , axis human joint , Journal of Neuroengineering and Rehabilitation , vol. 2, no. 7, pp. 1 18 , 200 5 . [19] L. Dipietro, A. M. Sabatini and P. Dario , Evaluation of an instrumented glove for hand movement acquisition. Journal of Rehabilitation , Journal of Rehabilitation Research & Development, vol. 40, no. 2, pp. 179 190, 2003. [20] R . Jafari, F . Dabiri, P . Brisk and M . Sarrafzadeh , Adaptive and fault tolerant medical vest for life critical medical monitoring , in Proceedings of the 2005 AC M symposium on Applied computin g , 2005, pp. 272 279. [21] System for in Proceedings of International Conference on Body Sensor Networks (BSN) , 2011 , pp. 35 40. [22] Saba Bakhshi et al ., Dev elopment of a Body Joint Angle Measurement System Using IMU Sensors , in 33rd Annual International Conference of the IEEE EMBS Boston , 2011, pp. 6923 6926. [23] S. Corazza , Mohammad H. Mahoor, and B . S. Davidson , A markerless motion capture system to study musculoskeletal biomechanics: Visual hull and simulated annealing approach , Annals of Biomedical Engineering , vol. 34, no. 6, pp. 1019 1029, 2006. PAGE 108 108 [24] V. Macellari, CoSTEL: a computer peripheral remote sensing device for 3 dim ensional monitoring of human motion , Medical and Biological Engineering and Computing , vol. 21, no. 3, pp. 311 318, 1983. [25] T.B. Moeslund , A . Hilton, and V . Kru o ger , A survey of advances in vision based human motion capture and analysis. , Computer vision and image understanding , vol. 104, no. 2 3, pp. 90 126, 2006. [26] P . Roan, N . Deshpande, Y . Z. Wang, and B . Pitzer , with Low Cost Accelerometers and , in IEEE/RSJ International Conference on Intelligent Robots and Systems , 2012 , pp. 4822 4827. [27] Carmen M. N. Brigante, N . Abbate, A . Basile, A . C . Faulisi et al ., Miniaturization of a MEMS Based , IEEE Transactions on Industrial E lectronics , vol . 58, no . 8, pp. 3234 3241, 2011 . [28] S. Miyazaki, Long term unrestrained measurement of stride length and walking velocity utilizing a piezoelectric gyroscope , IEEE Transactions on Biomedical Engineering , vol. 44, no. 8, pp. 753 759, 1997. [29] R. Moe Nilssen, A new method for evaluating motor co ntrol in gait under real life environmental conditions. Part 1: The instrument , Clinical Biomechanics , vol. 13, no. 4 5, pp. 320 327, 1998. [30] I.P.I. Pappas , Milos R. Popovic, T . Keller, V . Dietz et al ., A reliable gait phase detection system , IEEE Transac tions on Neural Systems and Rehabilitation Engineering , vol. 9, no. 2, pp. 113 125, 2001. [31] R. Moe Nilssen and J.L. Helbostad, Trunk accelerometry as a measure of balance control during quiet standing , Gait and posture , vol. 16, no. 1, pp. 60 68, 2002. [32] T . Kyriacou , the parameters of a biologically inspired model , J . Comput Neurosci , vol. 32 , pp. 281 295 , Jul. 2012. [33] Emma Beanland , Luana C. , Maina, Brad Aisbetta, Paul Gastin et al ., Val idation of GPS and accelerometer technology in swimming , J Sci Med Sport , 2013, to be published. [34] Michelle Norris , Ross Anderson and Ian C Kenny et al ., Method analysis of accelerometers and gyroscopes in running gait: A systematic review , J Sports Engi neering and Technology , 2013, to be published. [35] Yo Shih , Chih Shan Hob and Tzyy Yuang Shiang , Measuring kinematic changes of the foot using a gyro sensor during intense running , Journal of Sports Sciences , 2013, to be published. PAGE 109 109 [36] Y . W. Guo , G .R. Zhao , Q .Q. Liu , Z .Y. Mei et al ., evaluation of hemiplegic gait using an inertial body sensor , BioMedical Engineering OnLine , vol. 12, p . 83, 2013 . [37] V . Marotto , A . Serra, D . Carboni , M . Sole et al ., Orientation Analysis through a Gyroscope Sensor for Indoor Navigation Systems , in SENSORDEVICES 2013 : The Fourth International Conference on Sensor Device Technologies and Applications , 2013, pp. 85 90. [38] S . Moel and I . Schjplbergl, Real Time Hand Guiding of Industrial Manip ulator in 5 DOF using Microsoft Kinect and Accelerometer , in The 22nd IEEE International Symposium on Robot and Human Interactive Communication , 2013, pp. 644 649. [39] A.D. King, Inertial navigation forty years of evolution , General Elecric Company Review , vol. 13, no. 3, pp. 140 149, 1998. [40] P . G. Savage , Navigation , J ournal of G uidance , C ontrol , and D ynamics , vol. 36, no. 3, pp. 637 655, 2013. [41] N. Barbour and G. Schmidt, Inertial sensor technology trends , IEEE Sensors Journal , vol. 1, no. 4, pp. 332 339, 2001. [42] H . S. Zeng and Y . Zhao , Sensing Movement: Microsensors for Body Motion Measurement , Sensors , vol. 11, pp. 638 660 , 2011. [43] D. Yamane , T. Konishi , T. Matsushima , G. Motohashi et al ., n A rrayed MEMS A ccelerometer with a W ide R ange of D etection , in Transducers , 2013 , pp. 22 25. [44] Walther , B. Desloges , C . Lejuste, B . Coster et al ., gyroscope , J. Micromech. Microeng. , vol. 23, pp. 1 8, 2013. [45] K . Y. Tong and Malcolm H. Granat , A practical gait analysis system using gyroscopes , Medical Engineering & Physics , vol. 21 , pp, 87 94, 1999. [46] H . Huang , V . Agafonov , and H . Y. Yu , Motion , Sensors , vol. 13, pp. 4581 4597 , 2013. [47] C.V.C. Bouten , Karel T. M. Koekkoek, Maarten Verduin, Rens Kodde et al ., A triaxial accelerometer and portable data processing unit for the assessment of daily physical activity. IEEE Transactions on Biomedical Engineering , vol. 44, no. 3, pp. 136 147, 1997. [48] Y . F. Chen , C . J. Liu, and Q.J. Chen , A Vestibular System Model for Robots and Its Application in Environment Perception , in International Conference on Computing, Control and Ind ustrial Engineering , 2010, pp. 230 235. PAGE 110 110 [49] E . Bernmarka and C . Wiktorin , , Applied Ergonomics , vol. 33 , pp. 541 547 , 2002. [50] S . Thiemjarus , P . Poomchoompol, I. Methasate, and T . Theeramunkong , Constraints of Accelerometer based Range of Motion Estimation , in Proceedings of the IEEE EMBS International Conference on Biomedical and Health Informatics , 2012, pp. 551 554. [51] P.H. Veltink , Hans B. J. Bussmann, Wiebe de Vries, Wim L. J. Martens et al ., Detection of static and dynamic activities using uniaxial accelerometers , IEEE Transactions on Rehabilitation Engineering , vol. 4, no. 4, pp. 375 385, 1996. [52] A . Caroselli , Fabio BagalÃ , and Angelo Cappello , Quasi Real Time Estimation of Angular Kinematics Using Sin gle Axis Accelerometers , Sensors , vol. 13, pp. 918 937, Jan. 2013. [53] B. Kemp , Ad J.M.W. Janssen, and Bob van der Kamp , Body position can be monitored in 3D using miniature accelerometers and earth magnetic field sensors , Electroencephalography and Clinical Neurophysiology , vol. 109, no. 6, pp. 484 488, 1998. [54] G . X . Lee and Kay Soon Low, A Factorized Quaternion Approach to Determine the Arm Motions Using Triaxial Accelerometers With Anatomical and Sensor Constraints , IEEE Transactions on Instrumenta tion and Measurement , vol. 61, no. 6, pp. 1793 1802, Jun. 2012. [55] F . Ghassemi , S . Tafazoli, Peter D. Lawrence , and K . Hashtrudi Zaad , Design and calibration of an integration free accelerometer based joint angle sensor , IEEE Transactions on Instrumentation and Measurement , vol. 57, no. 1, pp. 150 159, 2007. [56] F. Ghassemi , S . Tafazoli, Peter D. Lawrence , and K . Hashtrudi Zaad , An accelerometer based joint angle sensor for heavy duty manipulators , i n IEEE In ternational Conference on Robotics and Automation . vol. 2, 2002, pp. 1771 1776. [57] R . E. Mayagoitia , Anand V. Nene , and Peter H. Veltink ., gyroscope measurement of kinematics: an inexpensive alternative to optical motion , Journal of Biomechanics , vol. 35 , pp. 537 542 , 2002. [58] H . Dejnabadi , Brigitte M. Jolles, and Kamiar Aminian , Measurement of Uniaxial Joint Angles Based on a Combination of Accelerometers , IEEE T ransactions on Biomedical Engineering , vol . 52, no . 8, pp. 1478 1 484, 2005. [59] A.T.M. Willemsen , Carlo Frigo, and Herman B. K. Boom , Lower extremity angle measurement with accelerometers error and sensitivity analysis , IEEE Transactions on Biomedical Engineering , vol. 38, no. 12, pp. 1186 1193, 2002. PAGE 111 111 [60] A.T.M. Willemsen , J . A. van A lste, a nd H. B. K. B oom , Real time gait assessment utilizing a new way of accelerometry , Journal of Biomechanics , vol. 23, no. 8, pp. 859 863, 1990. [61] Toshiaki Otani , Takateru Urakubo , Satoshi Maekawa , Hisashi Tamaki et al ., Position and Attitud e Control of a Spherical Rolling Robot Equipped with a Gyro , in AMC 06 , 2006, pp. 416 421. [62] Y.F. Jiang and Y.P. Lin, Improved strapdown coning algorithms , IEEE Transactions on Aerospace and Electronic Systems , vol. 28, no. 2, pp. 484 490, 1992. [63] Y . X. Wu and X . F. Pan, Velocity/Position Integration Formula Part I: Application to In Flight Coarse Alignment , IEEE Transactions on Aerospace and Electronic Systems , vol. 49, no. 2, pp. 1006 1023, 2013. [64] E. R. DAMIANOT and R. D. RABBITT , ation model of fluid , J. Fluid Mech. , vol. 307, pp. 333 372 , 1996. [65] R. D. RABBITT and E. R. DAMIANO , , J. Fluid Mech. , v ol. 238, pp. 337 369 , 1992. [66] T. Mergner , G . Schweigart, L . Fennell , Vestibular humanoid postural control , Journal of Physiology , vol. 103, no. 3 5, pp. 178 194, 2009. [67] T. Mergner , G. Schweigart , C. Maurer , A. Blumle , Human postural responses to motion of real and virtual visual environments under different support base conditions , Experimental brain research , vol. 167, no. 4, pp. 535 556, 2005. [68] J. Laurens and J. Droulez, Bayesian processing of vestibular information , Biological Cybernetics , vol. 96, no . 4, pp. 389 404, 2007. [69] F. Patane , C. Lasch i, H. Mi w a , E. Guglielmelli et al ., Design and development of a biologically inspired artificial vestibular system for robot heads , in IEEE/RSJ International Conference on Intelligent Robots and Systems . vol. 2, 2004, pp. 1317 1322. [70] Peter H. Veltink , Henk J. Luinge , Bart J. Kooi , Chris T.M. Baten et al ., vestibular system design of a tri axial inertial sensor system and its application in , Control of Posture and Gait, Proceedings of the International Society for Postural and Gait Research ISP G , 2001. [71] Giovanni Passetti , Federico Corradiyz, Marco Ragliantiz, Davide Zambrano et al ., , In Biomedical Circuits and Systems Conference (BioCAS), 2013 , pp. 174 177. PAGE 112 112 [72] R. Williamson , , Medical & Biological Engineering & Computing , v ol. 39 , pp. 294 302, 2001. [73] G . To and Mohamed R. Mahfouz , Quaternionic Attitude Estimation for Robotic and Human Motion Tracking Using Sequential Monte Carlo Methods With von Mises Fisher and , I EEE Transactions on Biological Engineering , vol . 60, no . 11, pp. 3046 3059, 2013. [74] H. J. Luinge and P.H. Veltink, Inclination measurement of human movement using a 3 D accelerometer with autocalibration , IEEE Transactions on Neural Systems and Rehabilitation Engineering , vol. 12, no. 1, pp. 112 121, 2 004. [75] H. J. Luinge and P.H. Veltink, Measuring orientation of human body segments using miniature gyroscopes and accelerometers , Medical and Biological Engineering and Computing , vol. 43, no. 2, pp. 273 282, 2005. [76] J . Honkakorpi , Juho Vihonen , and Jouni M attila , MEMS based State Feedback Control of Multi Body Hydraulic Manipulator , in IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS) , 2013, pp. 4419 4425. [77] J . Vihonen , R. Porter, and B. Shirinzadeh , Aided MEMS Motion State Estimation for Multi , in IEEE/ASME International Conference on Advanced Intelligent Mechatronics (AIM ) , 2013, pp. 341 347. [78] S . Tadano et al ., and Gyro Sensors Based on , Sensors , vol. 13, pp. 9321 9343 , 2013. (1314) [79] Z . Q. Zhang , J . Z. Shang, C . Seneci and Guang Zhong Yang , Sensing Using Micro , in IEEE/RSJ International Conference on Intelligent Robots and Syste ms (IROS) , 2013, pp. 831 836 . [80] A.R. Schuler , A . G rammatikos , and K . A. Fregly , Measuring rotational motion with linear accelerometers , IEEE Transactions on Aerospace and Electronic Systems , vol. 3, no. 3, pp. 465 472, 1967. [81] Sebastian O.H. Madgwick , Andrew J.L. Harrison , Paul M. Sharkey , Ravi Vaidyanathan et al ., Measuring motion with kinematically redundant accelerometer arrays: Theory, simulation and implementation , Mechatronics , vol. 23, pp. 518 529, May 2013. [82] J . Vihonen , R. Porter, and B. Shirinz adeh , Geometry Aided MEMS Motion State Estimation for Multi Body Manipulators , in IEEE/ASME International Conference on Advanced Intelligent Mechatronics (AIM) , 2013, pp. 1349 1354. PAGE 113 113 [83] G . Ligorio and A . M . Sabatini, Extended Kalman Filter Based Methods fo r Pose Estimation Using Visual, Inertial and Magnetic Sensors: Comparative Analysis and Performance Evaluation , Sensors, vol. 13, pp. 1919 1941, Feb. 2013. [84] K.A. Tahboub, Optimal estimation of body angular velocity based on otolith canal interaction , in Mediterranean Conference on Control and Automation , 2008, pp. 848 853. [85] E. Foxlin, Inertial head tracker sensor fusion by a complimentary separate bias kalman filter , in Proceedings of Virtual Reality Annual International Symposium , 1996, pp. 185 194. [86] A . J. Petruska and S . G. Meek, Non Drifting Limb Angle Measurement Relative to the Gravitational Vector During Dynamic Motions Using Accelerometers and Rate Gyros , in IEEE/RSJ International Conference on Intelligent Robots and Systems , 2011, pp. 3632 3637 . [87] P . Goel , Stergios I. Roumeliotis and Gaurav S. Sukhatme , Robust Localization Using Relative and Absolute Position Estimates , in Proceedings of the 1999 IEEWJ International Conference on Intelligent Robots and Systems , 1999, pp.1134 1140. [88] Y . R. Wang , H . Zhang, and Q . F. Zhou , , TELKOMNIKA , v ol.10, n o.7, pp. 1869 1878 , Nov . 2012 . [89] Carlos Rodriguez Donate , Roque Alfredo Osornio Rios, Jesus Rooney Rivera Guillen , and Rene de Jesus Romero Troncoso , Multi , S ensors , vol. 11, pp. 4335 4357 , 2011. [90] Y . Teruyama and T . Watanabe , Gain Kalman Filter based on A ngle Error Calculated from Acceleration Signals for Lower Limb Angle , in 35th Annual International Conference of the IEEE EMBS , 2013, pp. 3423 3426. [91] L . Peppoloni , A . Filippeschi , E . Ruffaldi , and Carlo A . Avizzano , A novel 7 degrees of freedom model for upper limb kinematic reconstruction based on wearable sensors , in IEEE 11th International Symposium on Intelligent Systems and Informatics , 2013, pp. 105 110. [92] Y . Wang , Chieh Chien, James Xu, Gregory Pottie et al ., Ga it Analysis Using 3D Motion Reconstruction With an Activity Specific Tracking Protocol , in ICASSP, 2013, pp. 1041 1045. [93] J. Vaganay , M. J. Aldon, and A. Fournier , Mobile robot attitude estimation by fusion of inertial data , in IEEE International Confere nce on Robotics and Automation , 1993, pp. 277 282. PAGE 114 114 [94] B. Barshan and H.F. Durrant Whyte, Inertial navigation systems for mobile robots , IEEE Transactions on Robotics and Automation , vol. 11, no. 3, pp. 328 342, 1995. [95] X .L. Men , Z . Q. Zhang , J .K. Wu, and W .C. Wong , , I EEE T ransactions on Biomedical Engineering , vol . 60, no . 7, pp. 2052 2063, Jul. 2013 . [96] Vishesh Vikas , V estibular D ynamic I nclinometer and M easurement of I nclination P arameters , Ph.D. Dissertation, University of Florida , 20 11 . [97] K . Shimamoto , K . Watanabe, and K . Ohnishi , Force Projection for 2 DOF Tendon , i n IEEE International Conference on Mechatronics (ICM) , 2013, pp. 286 291. [98] K . Kunze and P . Lukowicz, Dealing With Sensor Displacement In Motion Based Onbody Activity Recognition Systems , in UbiComp 08 , 2008, pp. 20 29. [99] M . Glueck , D . Oshinubi , and Y. Manoli , Automatic realtime offset calibration of gyroscopes , Sensors Applications Symposium (SAS) 2013 IEEE , pp. 214 218, 2013. [100] H . Y. , M.S. Thesis, U niversity of C algar y, 2004. [101] J.C. Lu tters , J. Schipper , P.H. Veltink , W. Olthuis et al ., Procedure for i n use calibration of triaxial accelerometers in medical applications , Sensors and Actuators A , vol. 68, pp. 221 228 , 1998. [102] I . Frosio , F . Pedersini, and N. Alberto Borghese , Autocalibration of MEMS Accelerometers , IEEE T ransactions on I nstrumentation and M easurement , vol . 58, no . 6, pp. 2034 2041, J un. 2009 . [103] J. Li, Research on the Application of the Time Serial Analysis Based Kalman Filter in MEMS Gyroscope Random Drift Compensation , Chin ese Journal of Sensors and Actuators , vol. 10, no. 5, pp. 2215 2219, 2006. [104] N . El Sheim , H . Y. Hou, and Xi . J. Niu , , IEEE T ransactions on I nstrumentation and measurement , vol . 57, no.1, pp. 140 149 , J an . 2008 . [105] V. V. Chikovani , Coriolis Vibratory and Fiber Optic Gyros Based on in Proceedings of , 2013, pp. 153 156. [106] S. Butterworth, On the theory of filter , Experimental wireless and the wireless engineer , vol. 7, pp. 536 541, Oct. 1930. PAGE 115 115 [107] M . Kaur , B Singh , and Seema , of , in International Conference and Workshop on Emerging Trends in Technology , 2011, pp. 1290 1294. [108] N. P askaranandavadivel , , P . Du , and Leo K Cheng , sl , Neurogastroenterol . and Motil . , vol. 25, pp. 79 83, 2013 . [109] T. Pulecchi , A. Manes , M. Lisignoli , and M. Giglio , data acquired during the intervention , Measurement , vol. 43 , pp. 4 55 468, Dec. 20 09. [110] M . Dyck and M . Tavakoli , Human Arm , in IEEE International Conference on Rehabilitation Robotics , 2013, pp. 1 8 . [111] A . S avitzky and M arcel J. E. G olay, by , A nalytical C hemistry , vol . 36, no.8, pp. 1627 1639, J ul. 1964 . [112] R . W. Schafer , O n the Frequency Domain Properties of Savitzky Golay Filters , Hewlett Packard Laboratories , 2010. [113] H . Nozato , T. Bruns , H . Volkers , and A . Oota , A study of Savitzky Golay filters for derivatives in primary , ACTA IMEKO , vol.2, no. 2, pp. 41 47, Dec . 2013 . [114] A. Wald and J. Wolfowitz , , The Annals of Mathematical Statistics , v ol. 11, no.2, pp. 147 162, Jun. 1940 . PAGE 116 116 BIOGRAPHICAL SKETCH Linghan Zhao was born in Hangzhou, China in 1986 . He received his B achelor of Science in mechanical engineering from Zhejiang University, Hangzhou in 2008. After that he worked at Zhige Intelligence Technology Co., Ltd . He received his first Master of Science in mechatronic engineering in Zhejiang University of Technology in summer 2012. H e completed his second Master of Science in mechanical engineering in summer 2014. His research interests include robotics, artificial intelligence and nonlinear control. PAGE 1 136IEEETRANSACTIONSONBIOMEDICALENGINEERING,VOL.44,NO.3,MARCH1997ATriaxialAccelerometerandPortableData ProcessingUnitfortheAssessment ofDailyPhysicalActivityCarlijnV.C.Bouten,*KarelT.M.Koekkoek,MaartenVerduin,RensKodde,andJanD.JanssenAbstractâ€“ Thepresentstudydescribesthedevelopmentof atriaxialaccelerometer(TA)andaportabledataprocessing unitfortheassessmentofdailyphysicalactivity.TheTAis composedofthreeorthogonallymounteduniaxialpiezoresistive accelerometersandcanbeusedtoregisteraccelerationscovering theamplitudeandfrequencyrangesofhumanbodyacceleration.Interinstrumentandtestretestexperimentsshowedthat theoffsetandthesensitivityoftheTAwereequalforeach measurementdirectionandremainedconstantontwomeasurementdays.Transversesensitivitywassignicantlydifferentfor eachmeasurementdirection,butdidnotinuenceaccelerometer output(<3%ofthesensitivityalongthemainaxis).Thedata unitenablestheon-lineprocessingofaccelerometeroutputto areliableestimatorofphysicalactivityovereight-dayperiods. Preliminaryevaluationofthesystemin13malesubjectsduring standardizedactivitiesinthelaboratorydemonstratedasignicantrelationshipbetweenaccelerometeroutputandenergy expenditureduetophysicalactivity,thestandardreferencefor physicalactivity(r =0.89).Shortcomingsofthesystemareits lowsensitivitytosedentaryactivitiesandtheinabilitytoregister staticexercise.Thevalidityofthesystemfortheassessmentof normaldailyphysicalactivityandspecicactivitiesoutsidethe laboratoryshouldbestudiedinfree-livingsubjects. IndexTermsâ€“ Assessmentofdailyphysicalactivity,dataprocessingunit,energyexpenditure,physicalactivity,piezoresistive accelerometer,portableunit,triaxialaccelerometerI.INTRODUCTIONTHEquantitativeassessmentofdailyphysicalactivityin humansrequiresanobjectiveandreliabletechniquetobe usedunderfree-livingconditions.Fromaphysiologicalpoint ofview,physicalactivitycanberegardedasanymovement orposturethatisproducedbyskeletalmusclesandresults inenergyexpenditure[1].Currently,theenergyexpenditure duetophysicalactivityiswidelyacceptedasthestandard referenceforphysicalactivity[2],butmeasurementofthis variableunderconditionsofdailylivingisimpracticaland notfeasibleforpopulationstudies.Therefore,theinterest forestimatesofenergyexpenditurebasedonobservations, questionnaires,heartraterecordings,ormovementregistration isgrowing.Atpresent,movementregistrationwithbody-ManuscriptreceivedDecember15,1994;revisedSeptember30,1996. Asteriskindicatescorrespondingauthor. *C.V.C.BouteniswiththeUniversityofTechnology,Divisionof ComputationalandExperimentalMechanics,P.O.Box513,Eindhoven5600 MBTheNetherlands(e-mail:carlijn@wfw.tue.nl). K.T.M.Koekkoek,M.Verduin,R.Kodde,andJ.D.Jansenarewith theUniversityofTechnology,DivisionofComputationalandExperimental Mechanics,Eindhoven5600MBTheNetherlands. PublisherItemIdentierS0018-9294(97)01466-3.xedmotionsensorsoffersthebestalternativefordaily physicalactivityassessment.Variousmotionsensorshave beendesignedforthispurpose,rangingfrommechanical pedometers[3]andactometers[4]toelectronicaccelerometers [5][7].Unfortunately,motionsensorscannotbeusedto determinethestaticcharacteristicsofphysicalactivity.Itis assumed,however,thatthecontributionofstaticexerciseto thetotallevelofphysicalactivityisnegligibleundernormal dailylivingconditions[8],[9]. Electronicaccelerometersarethemostpromisingmotion sensorsforphysicalactivityassessmentinfree-livingsubjects. Thesesensorsrespondtobothfrequencyandintensityof movement,andinthiswayaresuperiortopedometersand actometers,whichareattenuatedbyimpactortiltandonly countbodymovementifacertainthresholdispassed.Dueto thecurrentstateofartinintegratedcircuittechnologythereis alsogoodopportunitytobuildverysmallandlightweightaccelerometersystemsthatcanbewornfordaysorevenweeks. Accelerometershavebeenusedforseveraldecadestostudy gaitandothermovements[10],[11]orforthemeasurement oftremorandmotoractivityinneurologicalpatients[12], [13].Theuseofaccelerometersfortheassessmentofphysicalactivityisbasedondemonstratedrelationshipsbetween accelerometeroutputandenergyexpenditureinstudieson gaitanalysisandergonomics.Instudyingtheverticalforces resultingfrombodymovementinindustrialworkerswitha forceplatform,Brouha[14]foundasignicantcorrelation betweentheintegraloftherectiedforce-timecurveandenergyexpenditure( 0.83.96).Ismail etal .[15]measured forcesduringwalkinginthreeorthogonalplanesbyusing aforceplatformplacedunderthebeltofatreadmill.By integrationoftherectiedforce-timecurvesitwaspossibleto predictenergyexpenditurefromeachmeasurementdirection ( 0.73.92).Thebestpredictionofenergyexpenditure, however,wasachievedbysummationoftheabsolutevalues ofthethreeorthogonalforces.Reswick etal .[16]used ahead-mountedaccelerometerduringwalkingonalarge walkway.Theyconcludedthattheintegralofthemodulus ofaccelerometeroutputwaslinearlyrelatedtotheenergy expenditureduringwalking.Theseresultshaveledseveral researcherstohypothesizethattheintegralofthemodulus ofaccelerationmeasuredonthehumanbodyâ€“especiallyin verticaldirectionâ€“canbeusedtopredictenergyexpenditure duetophysicalactivity.In1981,Wong etal .[7]developed aportablemotionsensorwithapiezo-electricaccelerometer0018/97.00 1997IEEE PAGE 2 BOUTEN etal. :ACCELEROMETERANDDATAPROCESSINGUNITFORASSESSMENTOFDAILYPHYSICALACTIVITY137toestimateenergyexpenditureduringvariousactivities.This sensor,namedCaltrac,couldbeattachedtothewaistand registeredaccelerationsparalleltotheverticalaxisofthebody. Theabsolutevalueofaccelerometeroutputwasintegrated andsummedforthetotaltimeitwasworn.Thismeasure wasreferredtoasaccelerationcount.Reproducibilityof accelerationcountsmeasuredatthelowbackduringvarious exerciseswasgood( 0.94)andapooledcorrelation coefcientof0.74wasfoundfortherelationshipbetween accelerationcountandenergyexpenditure[6].Frommore recentstudiesitcanbeconcludedthatCaltracoutputshows goodcorrelationwithmeasuredenergyexpenditureduring separatewell-denedactivitiesinthelaboratory,liketreadmill walking( 0.68.94),butgenerallyoverestimatesenergy expenditureunderthesecircumstances[17],[18].UnderfreelivingconditionstheuniaxialCaltraccanbeusedtodistinguish amonginterindividuallevelsofdailyphysicalactivity,but tendstosystematicallyunderestimatetheenergyrequirements whencomparedtowholebodyindirectcalorimetry[19]or thedoublylabeledwatermethod[20],[21].Itseemsobvious thattheuseoftheuniaxialCaltraclimitstheregistrationof multidirectionalbodymovementandhencethepredictionof energyexpenditure.ThiswasdemonstratedbyAyenandMontoye[22]whoshowedthatenergyexpenditureduringexercises likewalking,running,andsquat-thrustsisbetterpredicted byusingthreeseparateCaltracs,mountedatrightangleson thewaist( 0.75),thanbyusingasingleverticalCaltrac ( 0.65).Recently,Meijer etal .[5]developedaportable triaxialaccelerometer(TA)fortheassessmentofdailyphysical activity.Thisaccelerometerwasbaseduponabendedpiezoelectricplate,sensitivetoaccelerationsinthreedirections.A data-acquisitionunitwithsolid-statememorywasusedfor lowpassltering,amplication,rectication,andintegrationof theresultingaccelerometeroutputover1-mintimeintervals.In thelaboratoryalinearrelationshipbetweenaccelerometeroutputandenergyexpenditurewasfoundforthepooleddataof16 subjectsperformingspeciedexercisesresemblingactivitiesof dailyliving[23].However,athighlevelsofenergyexpenditure,especiallyinrunning,theaccelerometeroutputsystematicallyunderestimatedenergyexpenditure.Duringevaluations underfree-livingconditionsover50%ofthemeasureddata werelostduetobadperformanceofthedevices[24].Other shortcomingsofthedeviceusedbyMeijer etal .arethehigher sensitivityinverticalmeasurementdirectioncomparedtothe horizontalmeasurementdirections,thelackofinformation fromseparatemeasurementdirections,andthelackofadc responseforadequatecalibrationoftheaccelerometer.This lastitemismanifestinallpiezo-electricaccelerometers. Theprimaryaimofthepresentstudywastodevelopan improvedelectronicTAdevicefordailyphysicalactivityassessment.Initialinvestigationwasaimedatspecifyingthetype ofaccelerometertobeused.Forthispurpose,thefrequency andamplitudecharacteristicsofhumanbodyacceleration, whichwillultimatelydeterminethetechnicalspecicationsof theaccelerometer,weredescribed.Next,aTAwasdesigned. Inordertoinvestigatetherelativecontributionofdifferent measurementdirectionstotheassessmentofphysicalactivity, thisTAwascomposedofthreeseparate,orthogonallymounted uniaxialaccelerometers.TheuniaxialsensorsintheTAwere evaluatedconcerningtheirinterinstrumentandtestretestvariabilityinabench-test.Finally,aportabledataunitforthe on-lineacquisition,processing,andstorageofaccelerometer outputoverprolongedperiodsoftimewasdevelopedtostudy patternsofdailyphysicalactivity.Preliminaryevaluationof theTAandthedataunitagainstmeasurementsofenergyexpenditureduetophysicalactivity( )wasperformedduringastandardizedactivityprotocolinarespirationchamber. II.DEVELOPMENTOFATRIAXIALACCELEROMETERWhenchoosinganappropriateaccelerometerfortheassessmentofdailyphysicalactivityoneshouldconsiderthe specicationsofthevariousavailableelectronicaccelerometersaswellasthecharacteristicsofhumanmovementwhich determinetheoutputofabody-xedaccelerometer.These aspectsaredescribedbelow. A.ElectronicAccelerometers Electronicaccelerometersarebasedonpiezo-electricor piezoresistiveproperties.Piezo-electricaccelerometerscanbe consideredasdampedmass-springsystems,inwhichapiezoelectricelementactsasspringanddamper.Thiselement generatesanelectricalchargeinresponsetothemechanical force,andhence,theaccelerationappliedtoitbyasmall mass.Inpiezoresistiveaccelerometersspringanddamperare replacedbysiliconresistorswhichchangeelectricalresistance inresponsetotheappliedmechanicalload.Theresistorsare electricallyconnectedinaWheatstonebridgetoproducea voltageproportionaltotheamplitudeandfrequencyofthe accelerationofthesmallmassinthesensor.Piezoresistive accelerometersaresmallerthanpiezo-electricaccelerometers, butrequireanexternalpowersource.Furthermore,piezoresistiveaccelerometershavedcresponse,whilepiezo-electric accelerometersdonotrespondtoconstantacceleration. B.SourcesofAccelerometerOutput Theoutputofanidealaccelerometerwornonthehuman bodyoriginatesfromseveralsources:1)accelerationdue tobodymovement;2)gravitationalacceleration;3)external vibrations,notproducedbythebodyitself(e.g.,resulting fromvehicles);4)accelerationsduetobouncingofthesensor againstotherobjectsorjoltingofthesensoronthebodydueto looseattachment,eventuallyresultinginmechanicalresonance [25].Ofthese,onlythersttwoaredirectlyrelatedto intentionalmovementofthebody.Gravitationalacceleration mayvarybetween 1and1g,dependingontheorientation ofthemeasurementdirectionofthesensorinthegravitational eld.Thissourceofaccelerometeroutputisoftendescribed asthegravitationalcomponent.Theoutputduetobody movement,usuallyreferredtoasthekinematiccomponent, isdependentonthetypeofactivityperformed,thepartofthe bodywhereaccelerationsaremeasured,andthemeasurement direction(antero-posterior,medio-lateral,orvertical),asis discussedbelow.Sources3and4,whichmayaddnoise totheaccelerometeroutput,canbeattenuatedbyadequatelteringtechniques(ifthefrequencyrangeofthenoisedoesnot PAGE 3 138IEEETRANSACTIONSONBIOMEDICALENGINEERING,VOL.44,NO.3,MARCH1997interferewiththefrequencyrangeofhumanbodyacceleration) andproperattachmentofthesensortothebody. C.HumanBodyAcceleration Accelerometersusedfortheassessmentofphysicalactivityinhumansmustprovideaccurateregistrationsofthe frequenciesandamplitudesofaccelerationsinvolvedinhuman movement.Oncethebroadestpossiblespectraoffrequencies andamplitudesisknown,therightsensorcanbechosen.In general,frequenciesandamplitudesofaccelerationsinvolved inhumanmovementarerelativelylow.Thelargestaccelerationswiththehighestfrequenciesmightbeexpectedduring runningandjumping. 1)Frequency: Duringlocomotion,frequenciesare generallyhigherintheverticalthaninthemedio-lateralor theantero-posteriordirectionsandthefrequencyspectrum shiftstowardhigherfrequenciesfromcranialtomorecaudal partsofthebody.Inwalkingatnaturalvelocitythebulkof accelerationpowerintheupperbodyrangesfrom0.8Hz, whereasthemostabruptaccelerationsoccuratthefootin verticaldirectionduringheelstrikeandsometimesamountup to60Hz[26].Bymeasuringtheseworstcaseaccelerations atheelstrikewithaforceplatformAntonnsonandMann [27]demonstratedthat99%oftheaccelerationpowerduring walkingwithbarefeetisconcentratedbelow15Hz.Higher frequenciesarecausedbytheimpactbetweenfootand walkingsurfaceanddonotdirectlyresultfromvoluntary muscularwork.Byusingpiezoresistiveaccelerometers (range: 20g,frequencyresponse:0Hz),tapedtothe skin,Bhattacharya etal .[28]foundthemajorityoffrequency componentsduringrunningtovarybetween1Hzin verticaldirectionattheankle.Atthelowbackandthehead, thefrequencycontentofaccelerationproleswassmaller. Duringtrampolinejumping,wheretheimpactbetweenfoot andjumpingsurfaceislesspronouncedthaninrunning, thefrequencycontentofverticalaccelerationsattheankle, thelowback,andtheheadwasapproximatelyequal,ranging from0.7Hz.Thefrequencycontentofdailyactivities hasrecentlybeenstudiedbySunandHill[29].FastFourier analysisofdailyactivitiesperformedonaforceplatform revealedthemajorenergybandtobebetween0.3.5Hz. 2)Amplitude: Likethefrequencycharacteristics,theamplitudesofaccelerationsinvolvedinlocomotionareusually higherintheverticaldirectionthaninbothhorizontaldirectionsandincreaseinmagnitudefromcranialtowardcaudal bodyparts.Duringwalking,forinstance,theaccelerationsof theupperbodydeterminedfromstereophotogrammetryrange from 0.3to0.8gintheverticaldirection,whereasinthe horizontaldirectionstheyrangefrom 0.2to0.2gatthe headandfrom 0.3to0.4gatthelowback[26].Atthetibia, theamplitudeofaccelerationsduringwalkingvariesbetween 1.7and3.3gintheverticaldirectionandbetween 2.1 and2.3ginthehorizontaldirections,asmeasuredwithbonemountedpiezoresistiveaccelerometers(range: 25g, :1000 Hz)[30].Duringrunning,Bhattacharya etal .[28]observed absoluteverticalpeakaccelerationsrangingfrom0.8.0gat thehead,from0.9.0gatthelowback,andfrom3.0.0g attheanklebyusingtheirskin-mountedaccelerometers. Duringtrampolinejumping,differencesinabsolutevertical peakaccelerationsmeasuredatthehead,thelowback,and theanklewerelesspronouncedthaninrunning:theyvaried between3.0.6gatthehead,between3.9.0gatthelow back,andbetween3.0.0gattheankle. 3)DesirableFrequencyandAmplitudeRanges: Despite thedifferenceinmeasurementtechniquesusedtodetermine thefrequencyandamplitudespectraofhumanbodyaccelerations,theabove-mentionedstudiesdoprovideinsightinto thefrequenciesandamplitudesthatmightbeexpectedduring normaldailyactivities.Consideringthendingsofthese studies,body-xedaccelerometersmustbeabletoregister accelerationswithintheamplituderangeof 12to 12gand withfrequenciesupto20Hzinordertoassessdailyphysical activity.Ingeneral,body-xedaccelerometersforphysical activityassessmentareplacedatwaistlevel[5],[23],[31].At thissiteanamplituderangeofabout 6to 6gwillsufce. Althoughthemajorkinematicaccelerationcomponentduring humanmovementisusuallyfoundintheverticaldirection, accelerationsinotherdirectionsdocontributetothetotal, complexmovementpattern.Thus,foracompleteregistration ofmultidirectionalhumanmovement,accelerationsshould bemeasuredinthreedirections.However,toinvestigatethe relativecontributionofseparatemeasurementdirectionsto theestimationof ,itshouldbepossibletoanalyzethe outputofthethreemeasurementdirectionsindependently. D.TheTriaxialAccelerometer Anaccurateregistrationofthefrequenciesandamplitudes ofaccelerationsinvolvedinhumanmovementisdependent onthetypeofaccelerometerused.Asismentionedabove, piezoresistiveaccelerometersyieldadcresponse,whereas piezo-electricaccelerometersdonot.Amajoradvantageof dcresponseisthatitenablesthecalibrationofpiezoresistiveaccelerometers,forinstancebyrotationwithinthe gravitationaleld.Humanmovement,however,willnever correspondtoadcresponse.Therefore,itwasdecidedto useapiezoresistivesensorincombinationwithahighpass ltertoeliminateanyresultingdccomponentduringthe registrationofhumanmovement,whileretainingaverylow frequencycutoff(about0.1Hz).Anumberofdifferentuniaxial piezoresistiveaccelerometerswasinvestigatedandasmall, lightweightaccelerometer(ICSensors,type3031-010,size: 4 4 3mm;weight0.3g;range: 10g,frequencyresponse: 0Hz; :1200Hz)wasselectedforfurtherresearch. Threeuniaxialaccelerometers(A ,A ,andA )weremounted orthogonallyontoa12 12 12-mmcubemadeofCeleron toaccomplishaTAwithindependentmeasurementdirections. III.EVALUATIONOFTHETRIAXIALACCELEROMETERInordertoassesstheusefulnessoftheTAforthemeasurementofaccelerationsinvolvedinhumanmovement,thedevice wastestedinalaboratoryexperiment.Theoutputofthethree uniaxialaccelerometersA ,A ,andA asaconsequenceof amechanicallyappliedaccelerationwasstudiedundervarious conditions.Interinstrumentandtestretestcomparisonswere PAGE 4 BOUTEN etal. :ACCELEROMETERANDDATAPROCESSINGUNITFORASSESSMENTOFDAILYPHYSICALACTIVITY139 (a)(b) Fig.1.Experimentalsystemtoevaluateinterinstrumentandtestretestvariabilityofaccelerometercharacteristics.(a)Lateralviewofequipm ent( = 0). (b)FrontalviewoftherotationalmotionappliedtoaTA( 6 =0).System-xedmeasurementdirections(~ e1,~ e2,~ e3)areindicated.RadiusRrepresents thedistancebetweentheTAandthecenterofrotation,~ gisthegravitationalaccelerationvector,andtheanglebetweenRand~ g.R25andR100refertoradiiof25and100mm,respectively.madetotestwhethertheoutputsofA ,A ,andA were similarunderthesameexperimentalconditionsandremained constantondifferentoccasions. A.Methods Aschematicrepresentationoftheexperimentalsetupis giveninFig.1.TheTAwasmountedontoacounterweighted armatdistancesof25mm( )or100mm( )from thecentralshaft.Aftermounting,thepositivemeasurement directionsoftheTAwereparalleltotheaxesofthelocal systemofreference( , , ).ByturningtheTAaround 90 anglestheuniaxialaccelerometerscouldbetestedin differentdirections.Bridgeampliers(PMI,type:OP9OGS) fortheaccelerometerswereattachedtotheoppositesideof thecounterweightedarm.Connectionsbetweentheampliers andtheTAwereestablishedviaa12-conductorshieldedcable. Thecentralshaftofthearmwassecuredtotheelectricmotor ofalathe(LaaglandRotterdam,type:Celtic14)withan adjustablenumberofrevolutionspersecond(rps)tocause arotationalmotionofthearmwithconstantangularvelocity (accuracy 0.1rps).Aslipringassembly(HottingerBaldwin Messtechnik,type:SK12)wasusedforthetransmission oftheoutputvoltagesofthethreeseparateaccelerometers toacomputer.Thesupplyvoltagetotheaccelerometers andamplierswasalsoprovidedbymeansoftheslipring assembly. Therotationalmovementofthearmresultedinaconstant radialaccelerationinthe directionandgravitationalcomponentsalong and ,varyingsinusoidallybetween 1and 1g(g 9.812m s attheexperimentalsite)withabasic frequencyequaltotheappliednumberofrps.Noaccelerations wereappliedin direction.Apartfromtheradial(kinematic) andgravitationalcomponents,theaccelerometeroutputduringtheexperimentsresultedfromtheoveralloffsetofthe accelerometerandtheexperimentalsystem,andâ€“sincethe accelerometerswereassumedtobenonidealâ€“fromtransverse sensitivity.Thislastcomponentisdenedasthesensitivityof theaccelerometerstoaccelerationsfromdirectionsotherthan themeasurementdirection.Thus,accelerometeroutputs and in and direction,includingoffsetandtransverse sensitivity,canbedescribedby (1) (2) withtheindex ( 1,2,3)denotingthenumberof theuniaxialaccelerometertestedinthedenedmeasurement direction. representstheoffset, thesensitivityalong themainaxis,and thetransversesensitivityofthe th accelerometer. representstheanglebetweenradius and thegravitationalvectoroftheearth( ), thetimederivative of ,and theradialacceleration.Accelerometeroutput ismeasuredinmV, isexpressedasmV g ,and is expressedasmV g orasapercentageof (% ).Note thatalong and thetransversesensitivityduetooutput fromthe directionwasneglectedsincenoaccelerations wereappliedin direction.Accelerometeroutputalong isonlyproducedbyoffsetandtransversesensitivitydueto accelerationsin and direction (3) PAGE 5 140IEEETRANSACTIONSONBIOMEDICALENGINEERING,VOL.44,NO.3,MARCH1997TABLEI OFFSET(V0 ;i),SENSITIVITY(Gi),ANDTRANSVERSESENSITIVITY(ki,GIVENASABSOLUTEVALUEANDASAPERCENTAGEOFGi)FOREACHOFTHETHREEUNIAXIALACCELEROMETERS(A1,A2,ANDA3)INTHETRIAXIALACCELEROMETER.VALUESAREGIVENASMEANSDFROM33OBSERVATIONS Fig.2.TypicalexampleofTAoutputintimein~ e1,(|),~ e2({),and~ e3()directionduringrotationalmotionat4rpswitharadiusof100mm.where ,thetransversesensitivityduetoaccelerationalong ,and ,thetransversesensitivityduetoaccelerationalong ,areassumedtobeequal. Byusing11combinationsof1rpswiththetworadii, and ,theconstantaccelerationalong wasvaried from0.1.6g(1.2.7m s ).Accelerometeroutputs duringcombinationsof withrps 10and with rps 6couldnotbemeasured,sincetheoutputvoltage exceededthemaximalallowedoutputvoltageoftheampliers( 15V).Thetotalamplituderangeofhumanbody acceleration,however,wasincludedintheexperiments.For eachtrialtheoutputsfromA ,A ,andA weresampled andanalyzedsynchronouslyusingadata-acquisitionboard (DIFAMeasuringSystems).Signalsweresampledat800Hz toproducearepresentationofaccelerometeroutputintime. Atypicalexampleofaccelerometeroutputintimealong , ,and isshowninFig.2forrotationat4rpsat aradiusof100mm.Next,theamplitudespectrumofthe signalswasdetermined.Theaveragespectrumoffoursubsequentmeasurementsof2048samples,sampledat25Hz,was usedforfurthercalculations.Thefrequencyofthemeasured accelerationscouldbedeterminedwithanaccuracyof0.01Hz. Thedccomponentsoftheamplitudespectrawereusedto determine plustheoutputduetoconstantradialaccelerationin direction.Thedifferencebetweenaccelerometer outputsfromtrialsperformedatthesameangularvelocitybut atdifferentradiiwastakentoeliminate andtodetermine thesensitivityduetotheradialacceleration.Inthesameway thetransversesensitivityduetotheradialaccelerationsin and directionweredetermined.Frequencycomponentsof theamplitudespectrawerewindowedwithaat-topwindow beforetheywereusedtocomputethesensitivityandtransverse sensitivityduetogravitationalaccelerationin and directionandthetransversesensitivityduetogravitational accelerationin direction.Besidesthedeterminationof , ,and ,itwasinvestigatedwhetheraccelerometeroutput along wasproportionaltotheappliedradialacceleration withintheamplituderangeofhumanbodyacceleration. Theexperimentalprotocolof11trialswasrepeatedthree timesondayone,withtheTApositionedinanotherorientation oneachrepetition.Thus,33combinationsof ,rps,and orientationoftheTAwereperformedtocalculate , , and underdifferentconditions.Dataarepresentedasmean andstandarddeviation(sd).Interinstrumentvariabilityin , ,and wasstatisticallyanalyzedwitharepeated-measures analysisofvariance(ANOVA),followedbyapost-hoctest (ScheffeF-test)toindicatepossiblesignicantdifferences.A levelof5%( 0.05)wastakenaslevelofsignicance. Partoftheexperimentalprotocolwasperformedontwo separateoccasions(daysoneandtwo)todeterminetestretest variabilityin , ,and .Onlytrialsat2,4,and6rps wereincludedisthisanalysis.Statisticalcomparisonsbetween themeasuredparametersondaysoneandtwoweremadeby usingapairedT-test. Accelerationswithfrequenciesabove10Hzwerenotincludedintheexperiments.Itwasassumed,however,that frequencieswithintherangeofhumanbodymovement(up to20Hz)couldbemeasuredaccuratelywiththeuniaxialaccelerometerswhichhaveafrequencyresponseof0Hz.To verifythisassumption,oneoftheuniaxialaccelerometers(A ) wastestedusingavibrationexcitator(LingDynamicSystems, type:201)atfrequenciesbetween0Hzandamplitudes of0.50.25mm.Theamplicationofaccelerometeroutput wassimilartothatduringtheinterinstrumentandtestretest experiments. B.Results Valuesfor , ,and ondayonearegiveninTable I. and weresimilarforA ,A ,andA ,whereasthe absolutevaluesof ,aswellas expressedasapercentage PAGE 6 BOUTEN etal. :ACCELEROMETERANDDATAPROCESSINGUNITFORASSESSMENTOFDAILYPHYSICALACTIVITY141of ,weredifferentforA ,A ,andA .Inallthree accelerometers,however, didnotsignicantlyinuence accelerometeroutput( 5%ofthesensitivityalongthemain axis)anddidnotexceedthemaximalvaluefor given bythemanufacturer(3% ).Nodifferencesinparameter valueswerefoundfortheseparatemeasurementdirections. Forinstance,thetransversesensitivityofaccelerometerA , wassimilarwhenA ,wastestedalong , ,or . Fig.3showsthedccomponentsofaccelerometeroutputresultingfromconstantangularvelocityalong for accelerometersA ,A ,andA plottedagainsttheappliedconstantradialacceleration.Accelerometeroutputwascorrected foroffsetvalues.Theguresshowthatwithintheamplitude rangeofhumanbodyacceleration( 12g),thedccomponent ofaccelerometeroutputislinearlyrelatedtotheapplied radialacceleration.Thesensitivityoftheaccelerometers,as determinedfromtheslopeoftheregressionlinesinFig.3, was583mV g inA ,573mV g inA ,and578mV g inA . Valuesfor , ,and (mean sd)ondaysoneandtwo areindicatedinTableII.Notethataveragevaluesondayone maydifferfromthoseinTableI,sinceonly18observations wereusedfortestretestanalysis,while33observationswere includedintheinterinstrumentanalysis. Nodifferencesinparametervalueswerefoundbetweendays oneandtwo.Valuesfor ,eitherexpressesasmV g or as% ,weresignicantlydifferentfromeachotheronboth days.Again,inallthreeaccelerometers waswellbelow 5%of .Accelerometeroutputasafunctionoftheapplied radialaccelerationalong ondaytwowascomparabletothe datagiveninFig.3. Duringvibrationalmotionwithfrequenciesupto35Hz theamplitudeoftheoutputfromA (averageofthreetrials, correctedforoffset)waslinearlyrelatedtotheamplitudeofthe appliedacceleration(Fig.4).Thesensitivity, ofA ,during thevibrationexperiments,asdeterminedfromtheslopeofthe regressionline,was580mV g . FromtheseresultsitwasconcludedthattheuniaxialaccelerometersoftheTAcouldbeusedfortheregistrationof accelerationswithintheamplitudeandfrequencyrangeof humanbodyacceleration.Theaccelerometerswerebuiltin aathousingofCeleron(50 30 8mm,16gr)withtwo slitsforanelasticbelt.ByusingthisbelttheTAcouldeasily beattachedtoseverallocationsonthehumanbody. IV.DEVELOPMENTOFADATAPROCESSINGUNITInordertocorrelateanalogaccelerometeroutputtodiscrete dataonenergyexpenditure,theoutputmustbeprocessedto ausefulquantity.Recently,Bouten etal .[32]reportedonthe mostoptimalwayofdataprocessingfortheestimationof fromaccelerometeroutput.Fromsynchronousrecordingsof andaccelerometeroutput,measuredwiththe TAattachedtothelowback,itwasconcludedthat wasbetterpredictedfromthree-directionalthanfromunidirectionalaccelerometeroutputwhendifferenttypesofactivity (sedentaryactivities,walking)wereperformed.Summation ofthetimeintegralsofthemoduliofaccelerometeroutput fromtheseparatemeasurementdirections( ),resulted (a) (b) (c) Fig.3.Accelerometeroutput(dc-component),correctedforoffset,asa functionoftheappliedradialaccelerationinthreeuniaxialaccelerometers (A1,A2,andA3).inthemostaccuratepredictionof ( 0.95).Thedata processingreferredtoisgivenby (4) PAGE 7 142IEEETRANSACTIONSONBIOMEDICALENGINEERING,VOL.44,NO.3,MARCH1997TABLEII OFFSET(V0 ;i),SENSITIVITY(Gi),ANDTRANSVERSESENSITIVITY(ki,GIVEN ASABSOLUTEVALUEANDASAPERCENTAGEOFGi)OFTHETHREEUNIAXIALACCELEROMETERS(A1,A2,A3)INTHETRIAXIALACCELEROMETERONDAYSONEANDTWO.VALUESAREGIVENASMEANSDFROM18OBSERVATIONS *: ,signicantdifferencesbetween , ,and ondayoneas wellasondaytwo. #: ,signicantdifferencesbetween (% ), (% ),and (% )ondayoneaswellasondaytwo. Nosignicantdifferencesinparametervaluesbetweendaysoneandtwowere found.with thetimeperiodforintegrationand , ,and theaccelerometeroutputsalongthe , ,and directions ofabody-xedsystemofreference.Thisdataprocessing wassuperiorto,forinstance,theintegrationofthetotal accelerationvector. Whenusing fortheassessmentofphysicalactivityinfree-livingsubjects,theambulatoryrecordingsof thisvariableshouldnotinterferewiththesubjects'physical activities.Therefore,aportabledataunitwithminimalweight anddimensionsfortheon-lineacquisition,processing,and storageofaccelerometerdataisrequired.Inaddition,thisdata unitmustenabletheamplicationandlteringofacceleration signalsfromtheTAaswellasthestorageof over periodsofdaysorweekstostudypatternsofdailyphysical activity.Consideringthesedemands,afreeprogrammable dataunitwasdevelopedfortheassessmentofdailyphysical activitywithminimaldiscomforttosubjects. A.TheDataProcessingUnit TheblockdiagraminFig.5showstheprocessingofthe outputfromtheTA,whichisimplementedinthedataunit. TheconnectionbetweenthedataunitandtheTAisestablished viaa0.5mexible12conductorshieldedcable.Individual outputsfromthethreemeasurementdirectionsoftheTAare ampliedandhighpass(0.11Hz,5.6dB/octave)andlowpass (20Hz,9dB/octave)lteredtoattenuatedc-responseand frequenciesthatcannotbeexpectedtoarisefromvoluntary humanmovement.Next,accelerationsignalsaredigitized(100 Fig.4.Amplitudeofmeasuredaccelerations,correctedforoffset,asafunctionoftheamplitudeoftheappliedaccelerationinauniaxialaccelerometer (A1)duringvibrationalmotion.Hz)andfurtherprocessedusingaminiaturizeddatalogger (OnsetComputerCorporation,type:Tattletale5F).Thisdataloggerenablesafreeprogrammabledataprocessingandis programmedandstartedfromacomputerviaaserialinterface (OnsetComputerCorporation,type:TC-1).Fortheassessment ofdailyphysicalactivitythedataloggerwasprogrammedto calculate ,(4),usingitsTX-Basicsoftwarepackage. Thetimeperiodforintegrationisvariableandcanbeadjusted atthestartofameasurementperiod.Afterprocessing,the obtaineddataarestoredina512kB,16bitdatamemorychip thatcanbereadoutwiththeserialinterfacetoacomputer. Thereisalsothepossibilitytostartthedataloggerandtoread outthememorychipbymodem.Thememoryisresetby disconnectingthesupplyvoltagetotheTAandthedataunit, whichisprovidedbybatteries.Two9V,1200mAhbatteries arerequiredtoregisterandprocessaccelerationsignalsovera periodofeightdays.Forameasurementperiodofthreedays two9V,500mAhbatteriescanbeused.Batteries,datalogger, andotherelectroniccomponentsfordataprocessingarebuilt inahousingofPVC.Thishousingcanbeopenedbythe investigatorforreplacementofbatteriesandcalibrationof theaccelerometers.Theac/dcswitcheswithinthehousing areusedtoomitthehighpasslters.Thedcresponsesof theuniaxialaccelerometerscanthanbeappliedtodetermine theirsensitivitybyalteringtheorientationoftheTAwith respecttothegravitationalvectoroftheearth.Thegainofthe accelerometersaswellasthebalanceoftheWheatstonebridge intheaccelerometers(zeroing)canbeadjustedmanually.In thiswaytheseparatemeasurementdirectionsoftheTAcan becalibratedequally.Fortheassessmentofdailyphysical activitythegainisusuallysettoproduceanaccelerometer outputof1.5V g ,correspondingtoanoutputcountof 1000overanintegrationperiodof1mininthedataunit. Thedataunitmeasures110 70 35mmandweighs170gr withoutbatteries(250grincludingbatteries).Itcanbeworn aroundthewaistinasmallbag(fannypack),butitcanalso beattacheddirectlytoawaistbeltbyusingtwoslitsonboth sidesofthePVChousing. PAGE 8 BOUTEN etal. :ACCELEROMETERANDDATAPROCESSINGUNITFORASSESSMENTOFDAILYPHYSICALACTIVITY143 Fig.5.Block-diagramoftheTAandthedataunitforacquisition,processingandstorageofaccelerometeroutput.Thecompletedataprocessingis indicatedforoneoftheuniaxialaccelerometers(A2).B.EvaluationoftheTriaxialAccelerometer andDataProcessingUnit PreliminaryevaluationoftheTAanddataprocessingunit wasperformedinagroupof13youngmalesubjects(age: 27 4yr,bodymass:77 12kg,height:1.83 0.07m) duringastandardizedlong-termactivityprotocolinarespirationchamber[33].Thischamber(14m )isfurnishedwitha bed,table,chair,toilet,washing-bowl,radio,andtelevision andisprovidedwithequipmentforthedeterminationof metabolicenergyexpenditurefromrespiratorygasexchange. Theventilationratethroughthechamberwasmeasuredwith adry-gasmeter(Schlumberger,G4)andanalyzedbyaparamagneticO analyzer(Hartmann&Braun,Magnos6G)and aninfraredCO analyzer(Hartmann&Braun,Uras3G). Ingoingairwasanalyzedonceevery15minandoutgoing airtwiceevery5min.FromtheventilationrateandO and CO concentrationsinin-andout-goingair,O consumption andCO productionwerecalculatedon-lineonacomputer. Totalenergyexpenditurewascalculatedat5-minintervals fromO consumptionandCO production,accordingtoWeir [34].Thesubjectsstayedinthechamberfor36h:twonights andtheinterveningday.Duringday-time(0830)they performedstandardizedactivities,resemblingnormaldaily activities(sedentaryactivities,householdactivities,walking). Eachactivitywasperformedfor30min.Exceptforsteppingandcarryingloads,theactivitieswereperformedatthe subjects'preferredrate.Steppingwasperformedat5-min intervals,alternatelywith5-minrestperiods,onabench33 cmhighandatarateof30steps min or60steps min . Duringcarryingloadsthesubjectscarried1-kgirondisksfrom onesideoftheroomtotheothersideat5-minintervals. Foreach30-minactivityperiodtheaveragetotalenergy expenditurewasdetermined.Next, foreachactivity wascalculatedastheaveragetotalenergyexpenditureminus sleepingmetabolicrate,whichwasmeasuredovera3-h intervalbetween0230withthelowestlevelofactivity asindicatedbyaDopplerradarsysteminthechamber. wasexpressedinwattsandcorrectedforbodymass(W kg ).DuringtheentireactivityprotocoltheTAwasattached tothelowbackofthesubjectsatthelevelofthesecond lumbarvertebrabyusinganelasticbeltaroundthewaist. Thetimeintervalforintegrationofaccelerometeroutputwas setto1minand (counts min )wasaveragedfor each30-minactivityperiod.Associationsbetweenaverage and weredeterminedwithlinearregression analysis,accordingtotheleastsquaresprinciple.Correlation coefcients(Pearson's )andstandarderrorsofestimatewere calculated.Fig.6showsthemeansandstandarddeviations for and forthe13subjectsduringtheactivity protocol.Individualcorrelationsbetween and variedbetween0.87.97(mean:0.92),whereasstandard errorsofestimaterangedfrom0.2.5W kg (mean: 0.3W kg ).Whenusingthedataofallsubjectsandall activities,apooledcorrelationcoefcientof0.89wasfound (pooledstandarderrorofestimate:0.4W kg ).Comparison betweenestimatedandmeasuredvaluesof byusingthe pooledregressionequationfor versus resulted inanaverageoverestimationof by7.5%whenthe PAGE 9 144IEEETRANSACTIONSONBIOMEDICALENGINEERING,VOL.44,NO.3,MARCH1997 Fig.6.Accelerometeroutput(IMAtot)andenergyexpenditureforphysicalactivity(EEact)duringastandardizedactivityprotocolinarespiration chamber(meansdfrom13malesubjects).entireactivityprotocolwasconsidered.Regardingtheseparate activities,differencesbetweenestimatedandmeasuredvalues of rangedfrom 40%(underestimation)duringdish washingto 33%(overestimation)duringlying.Intensive activities,likestepping,carryingloads,walking,cleaning,and makingthebed,wereunderestimatedonaverageby6.2%, whereaslessdemandingactivities,likesitting,lying,and deskworkwereoverestimatedonaverageby6.6%usingthe accelerometerdevice. V.DISCUSSIONThepresentstudydescribesthedevelopmentandevaluation ofaTAanddataprocessingunitfortheassessmentofdaily physicalactivityintermsofenergyexpenditure.TheTAis composedofthreeseparateorthogonallymounteduniaxial accelerometersinordertoinvestigatewhether isbetterpredictedfromthree-directionalthanfromunidirectional accelerometeroutputandtostudytherelativecontribution oftheseparateuniaxialmeasurementdirectionstotheestimationof .Theseaspectscouldnotbestudiedbythe three-directionaldevicedescribedbyMeijer etal .[5],[23] whichconsistsofasinglesensor,sensitivetomultidirectional acceleration.Besidesthelackofinformationfromseparate measurementdirections,thedeviceofMeijer etal .isalso moresensitivetoaccelerationsintheverticaldirectionthan toaccelerationsinbothhorizontaldirections.Thismightbe adisadvantagewhenassessingactivitieslikewalkingand running,sinceithasbeenshownthatalthoughthemajor accelerationcomponentduringtheseactivitiesisinthevertical direction, isbetterpredictedfromaccelerometeroutput intheantero-posteriordirection[32].Asdeterminedfrom interinstrumentexperiments,thesensitivity( )ofourTA wassimilarforeachmeasurementdirectionâ€“i.e.,foreach uniaxialaccelerometer.Nodifferencesinoffsetvalues( ) werefound,andthetransversesensitivity( )didnotinuence accelerometeroutput.Testretestexperimentsshowedthat valuesfor , ,and weresimilarontwoconsecutive measurementdays.Furthermore,itwasconcludedthat,within PAGE 10 BOUTEN etal. :ACCELEROMETERANDDATAPROCESSINGUNITFORASSESSMENTOFDAILYPHYSICALACTIVITY145theamplituderangeofhumanbodyacceleration,theoutput ofallthreeuniaxialaccelerometerswasproportionaltothe appliedacceleration.Foroneoftheaccelerometers(A )it wasshownthattheoutputwasproportionaltotheapplied accelerationduringvibrationalmotionatfrequencieswithin thefrequencyrangeofhumanbodyacceleration.Itisassumed thatthisisalsotrueforA andA .Ascanbeobservedfrom TablesIandII,relativelylargedeviationsin werefound. Thisiscausedbytheincreaseof withthemagnitudeof theappliedacceleration.Duringaccelerationsabove10g, reachedvalues 5% .However,asaccelerationsabove 10garerelativelyscarceinhumanmovement,itisassumed that willnotaffecttheassessmentofdailyphysicalactivity. Amajoradvantageofthepiezoresistiveaccelerometersin theTAistheirdcresponse,whichfacilitatescalibration.Yet, theoutputvoltageoftheseaccelerometersisinuencedbyoffset.Inordertoavoidover-orunderestimationofthemeasured acceleration,accelerometeroutputmustbecorrectedforoffset values.Theoffsetmaydrift,however,asaconsequenceof temperaturechanges,resultinginamorecomplicatedcorrectionforoffsetvalues.Inourlaboratoryexperimentstheoverall offsetconsistedofaccelerometeroffsetandtheoffsetofthe ampliers.Signicantoffsetdrift,typically1%ofthefullscale overatemperaturerangeof0 Cinreferenceto25 Cin theaccelerometersand1.2 V C intheampliers,was notobservedduringtheseshort-termexperiments.During thelong-termassessmentofdailyphysicalactivitythe offsetâ€“andpossibleoffsetdriftâ€“isattenuatedbythehighpass lterinthedataunit.Nevertheless,offsetandoffsetdrift shouldbeconsideredduringcalibration,whenthehighpass lterisomitted.Thesensitivityofaccelerometersmayalso driftduetotemperaturechanges.Liketheoffsetdrift,the sensitivitydriftofthepiezoresistiveaccelerometersusedin thisstudyistypically1%ofthefullscaleandmayaffect theassessmentofphysicalactivity.However,nodriftin sensitivitywasobservedduringthelaboratoryexperiments. Fromcalibrationdatabeforeandaftertheactivityprotocolin therespirationchamber,itwasconcludedthatsensitivityand offsetdidalsonotdriftduringtheselong-termmeasurements. Apartfrombodyacceleration,theoutputofbody-xedaccelerometersresultsfromgravitationalaccelerationandnoise duetoexternalvibrationsorinadequateattachmentofthe accelerometer.Thegravitationalcomponentisdependenton theorientationoftheaccelerometerwithrespecttothegravitationalvectoroftheearth.Itmayinuencetotalaccelerometer outputconsiderably,especiallywhentheangle between themeasurementdirectionandthegravitationalvectorof theearthisrelativelylargeandthekinematiccomponentof accelerometeroutputissmall.Consequently,itmayaffectthe assessmentof andhencethepredictionof . Correctionforthegravitationalcomponentunderdailyliving conditionsispracticallyimpossible.Inordertominimizethe effectofthegravitationalcomponentonaccelerometeroutput, Servais etal .[31]arguedthattheattachmentofaccelerometers atlocationswhere(thevariationin) issmallâ€“e.g.,the waistorthelowbackâ€“issuperiortolocationswhere(the variationin) islargeâ€“e.g.,thelimbs.Thepreciseeffectof thegravitationalcomponenton andtherelationship between and ,however,isunknownandshould bestudied. Externalvibrationsmayalsoconsiderablyinuence underdailylivingconditions.Alowpasslter(20 Hz)isbuiltintothedataunittoattenuatefrequenciesthat cannotbeexpectedtoarisefromvoluntarymovement.Yet, contactoftheaccelerometerâ€“orthesubjectwearingitâ€“with vibratingexternalsources,likevehiclesormachinery,may poseamajorproblemwhenfrequenciesoftheexternalsources areinterferingwiththefrequencyrangeofhumanmovement. Forinstance,thevibrationofapowerlawnmower(about5.5 Hz)orthevibrationofgroundvehicles[25],[35]willaffect accelerometeroutputandhence . Theaccelerometershouldbexedproperlytothehuman bodytoavoidthesensorfrommovingorjoltingonthe skin.Attachmentdirectlytotheskinisessential,sincemovementsofclotheswillcauseartifactsinaccelerometeroutput. Preferably,theaccelerometerisxedwithadhesivetapeor elasticstraps.Inthiswayarmattachmentwithminimal discomforttosubjectsisachieved.Still,thesofttissuelayer underthemountedaccelerometermayaffectaccelerometer output,possiblyleadingtoresonanceoftheaccelerometer ontheskin.Recently,KitazakiandGrifn[36]studiedthe resonantbehaviorofskin-mountedaccelerometersatthelow back.Accelerometerswithdifferentmassesanddimensions wereattachedtotheskinoverthespinousprocessofthethird lumbarvertebrawithdouble-sidedadhesivetape.Duringa freevibrationtestineightmalesubjectstheyfoundresonant frequenciesaslowas15Hz(range:15Hz)whenusing anaccelerometerwithsimilarmeasuresastheTA(contact surface:40 35mm;totalmass:16gr).Thus,likeexternal vibrationswithfrequencies 20Hz,theseresonantfrequencies mayinuence andshouldbeconsideredduringdaily physicalactivityassessment. Theplaceofattachmentofaccelerometersonthehuman bodyisanimportantissue[37],[38].First,theaccelerometer onthehumanbodymaynotinterferewiththesubjects' activities.Second,thekinematicandgravitationalcomponents ofaccelerometeroutputaredependentonthemeasurement location.WechoosetoattachtheTAatthetrunk(secondlumbarvertebra)asthissegmentrepresentsthemajor partoftotalbodymassandismovingduringmostdaily activities.AttachmentoftheTAtothelowbackcaused minimaldiscomforttothesubjectsanddidnotinuence theperformanceoftheiractivities.Furthermore,asdiscussed above,theeffectofthegravitationalcomponentontotal accelerometeroutputatthislocationissmall.However,further researchshouldprovideevidenceregardingthebestplaceof attachmentofaccelerometersforphysicalactivityassessment. Inaddition,theinuenceofsmallvariationsinplaceof attachmentshouldbestudiedtotestwhetherintersubjectand intrasubjectvariabilityinplacementaffectstheassessmentof physicalactivity. Aportabledataunitwasdevelopedfortheon-lineacquisition,processing,andstorageofTAoutput.Thisdataunitwas programmedtocalculate ,avariablethatishighly relatedtothestandardreferenceofphysicalactivity: [32].EvaluationoftheTAandthedataunitagainst , PAGE 11 146IEEETRANSACTIONSONBIOMEDICALENGINEERING,VOL.44,NO.3,MARCH1997duringastandardizedactivityprotocolresemblingnormal dailyactivity,showedstrongindividual( 0.87.97)and pooled( 0.89)correlationsbetween and . Thesecorrelationsarehigherthanthosefoundforstandardized activitiesinthelaboratoryusingtheuniaxialCaltrac[6], [39][41]ortheTAofMeijer etal .[5],[23].Discrepancies betweenmeasuredandestimatedvaluesof during high-intensityactivitiesandlowintensityactivitieswerealso smallerthanthosereportedfortheCaltrac[17][19]andthe accelerometerofMeijer etal .[23]. Atleastpartofthediscrepanciesbetweenmeasuredand estimatedvaluesof canbeattributedtotheperformance ofactivitiesthatinvolvestaticexercise,likestepping.During staticexercisetheincreaseinaccelerometeroutputisnot proportionaltotheincreasein andtheactual willbeunderestimated.Althoughitisassumedthattheeffects ofstaticexerciseonthetotallevelofdailyphysicalactivityare negligible[8],[9],predictionsof fromaccelerometer outputduringactivitieslikewalkinguphill,carryingloads, orcyclingwithheadwindshouldbeevaluatedwithcaution. Thegravitationalcomponentofaccelerometeroutputand noisefrominadequateattachmentoftheTAmayalsohave inuencedthepredictionof from .Inaddition, itshouldberealizedthatthevaluesof inthisstudy arenotcorrectedforthethermogeniceffectoffoodconsumption.Therefore,relativelyhighvaluesof ,compared to ,maybefoundduringactivitiesfollowingfood consumption,ascanbeseenfordishwashing.Duringthis activitytheactual wasunderestimatedby40%when usingthepooledregressionequationbetween and . Itisobviousthattherelativelyhighvaluesfor and duringintensiveactivitiesmayconsiderablyinuence thecorrelationbetweenthesevariables.Aftereliminationof values 1500counts min ,correspondingtomakingthebed,stepping,cleaning,walking,andcarryingloads, individualcorrelationsrangedfrom0.57.83,withameanof 0.80(standarderrorofestimate:0.1.3W kg ,mean:0.1 W kg ).Thepooledcorrelationcoefcientdecreasedto0.77 (standarderrorofestimate:0.3W kg ).Thedecreasein correlationisnotonlycausedbythesmallerrangesin and duringlessintensiveactivities,butalsobythe interindividualdifferenceinperformanceoftheseactivities. Inaddition,arelativelylargecontributionofthethermogenic effectoffoodconsumptionaswellasarelativelylarge contributionofthegravitationalcomponentofaccelerometer outputmayhaveaffectedtherelationshipbetween and duringtheless-intensiveactivities. VI.CONCLUSIONThedevelopmentofaTAandportabledataunithas resultedinanewaccelerometerdeviceforphysicalactivity assessment.Majoradvantagesofthedevicearetheuseof piezoresistiveaccelerometers,whichfacilitatescalibration,the abilitytomeasureandanalyzeaccelerationsfromthreedifferentmeasurementdirections,andtheon-linedataprocessing toquantifyaccelerometeroutputasafunctionofphysical activity( ).Thecalculationof overadjustable short-termintervalsâ€“forinstanceoneminâ€“enablestheinvestigationofpatternsofphysicalactivityintime.From mechanicaltestingitwasconcludedthattheTAwasreliable andvalidforthemeasurementofaccelerationswithinthe frequencyandamplituderangeofhumanbodyacceleration. Preliminaryevaluationoftheaccelerometerdeviceunder laboratoryconditionsshowedsignicantrelationshipsbetween and .Theserelationships,however,maybe inuencedbyseveralfactors,liketheperformanceofstatic exercise,gravitationalacceleration,andnoisefromexternal sourcesorresonantbehavioroftheaccelerometerontheskin. Theseaspects,aswellastheplaceofattachmentoftheTA, shouldbeconsideredwhenevaluatingthedeviceinfuture applications.Inaddition,thevalidityandusefulnessofthe devicefortheassessmentof undernormaldailyliving conditionsorduringspecicactivitiesoutsidethelaboratory shouldbestudiedinfree-livingsubjects. REFERENCES [1]C.J.Caspersen,K.E.Powell,andG.M.Christenson,Physicalactivity, exerciseandphysicaltness:Denitionsanddistinctionsforhealth relatedresearch, PublicHealthRep., vol.110,pp.126,1985. [2]R.E.Laporte,H.J.Montoye,andC.J.Caspersen,Assessmentof physicalactivityinepidemiologicalresearch:Problemsandprospects, PublicHealthRep., vol.200,pp.131,1985. [3]H.G.C.KemperandR.Verschuur,Validityandreliabilityofpedometersinhabitualactivityresearch, Eur.J.Appl.Physiol., vol.37,pp. 71,1977. [4]P.Avons,P.Garthwaite,H.L.Davies,P.R.Myrgatroyd,andW.P.T. James,Approachestoestimatingphysicalactivityinthecommunity: Calorimetricvalidationofactometersandheartratemonitoring, Eur. J.Clin.Nutr., vol.42,pp.185,1988. [5]G.A.L.Meijer,K.R.Westerterp,F.M.H.Verhoeven,H.B.M. Koper,andF.tenHoor,Methodstoassessphysicalactivitywithspecial referencetomotionsensorsandaccelerometers, IEEEBiomed.Eng., vol.38,pp.221,1991. [6]H.Montoye,R.Washburn,S.Servais,A.Ertl,J.G.Webster,andF.J. Nagle,Estimationofenergyexpenditurebyaportableaccelerometer, Med.Sci.SportsExerc., vol.15,pp.403,1983. [7]T.C.Wong,H.J.Webster,H.Montoye,andR.Washburn,Portable accelerometerdeviceformeasuringhumanenergyexpenditure, IEEE Trans.Biomed.Eng., vol.28,pp.467,1981. [8]W.H.M.SarisandR.A.Binkhorst,Theuseofpedometerand actometerinstudyingdailyphysicalactivityinman.PartII:Validity measuringthedailyphysicalactivity, Eur.J.Appl.Physiol., vol.37, pp.229,1977. [9]R.VerschuurandH.C.G.Kemper,Adjustmentofpedometerstomake themmorevalidinassessingrunning, Int.J.SportsMed., vol.1,pp. 95,1980. [10]G.Currie,D.Rafferty,G.Duncan,F.Bell,andA.L.Evans,Measurementofgaitbyaccelerometerandwalkway:Acomparisonstudy, Med.Biol.Eng.Comput., vol.30,pp.669,1992. [11]J.R.W.Morris,Accelerometryâ€“Atechniqueforthemeasurementof humanbodymovements, J.Biomech., vol.6,pp.729,1973. [12]J.D.Frost,Triaxialvectoraccelerometry:Anewmethodforquantifyingtremorandataxia, IEEETrans.Biomed.Eng., vol.BME-25,pp. 17,1978. [13]M.S alzer,Three-dimensionaltremormeasurementsofthehand, J. Biomech., vol.5,pp.217,1972. [14]L.Brouha,Evaluationofrequirementsofjobs,in Physiologyin Industry. Oxford,U.K.:Pergamon,1960,pp.94. [15]A.H.Ismail,J.W.Baraby,andC.B.Smith,Relationshipsbetween mechanicalforceandphysiologicalcostduringgaitinadultmen,in Biomechanics:ProceedingoftheC.I.C.SymposiumBiomechanics, J.M. Cooper,Ed.Chicago:TheAthleteInst.,1971,pp.99. [16]J.Reswick,J.Perry,D.Antonelli,N.Su,andC.Freeborn,Preliminary evaluationoftheverticalaccelerationgaitanalyzer(VAGA),in Proc. 6thAnnu.Symp.ExternalControlExtremities, Dubrovnik,1978,pp. 305. PAGE 12 BOUTEN etal. :ACCELEROMETERANDDATAPROCESSINGUNITFORASSESSMENTOFDAILYPHYSICALACTIVITY147[17]J.A.Balogun,D.A.Martin,andM.A.Clendenin,Calorimetric validationoftheCaltracaccelerometerduringlevelwalking, Phys. Ther., vol.69,pp.501,1989. [18]E.M.HaymesandW.C.Byrnes,Comparisonofwalkingandrunning energycostusingtheCaltracandindirectcalorimetry, Med.Sci.Sports Exerc., vol.25,pp.1365,1993. [19]Y.Schutz,F.Froidevaux,andE.J equier,Estimationof24henergy expenditurebyaportableaccelerometer, Proc.Nutr.Soc., 1988,vol. 47,p.23A. [20]R.J.Gretebeck,H.Montoye,andW.Porter,Comparisonofthedoubly labeledwatermethodformeasuringenergyexpenditurewithCaltrac accelerometerrecordings, Med.Sci.SportsExerc., vol.23,suppl.356, 1991. [21]M.B.Heyman,P.Fuss,V.R.Young,W.J.Evans,andS.B. Roberts,PredictionoftotalenergyexpenditureusingtheCaltrac activitymonitor, Int.J.Obestet., vol.15,p.21,1991. [22]T.G.AyenandH.J.Montoye,Estimationofenergyexpenditurewitha simulatedthree-dimensionalaccelerometer, J.AmbulatoryMonitoring, vol.1,pp.293,1988. [23]G.A.L.Meijer,K.R.Westerterp,H.Koper,andF.tenHoor, Assessmentofenergyexpenditurebyrecordingheartrateandbody acceleration, Med.Sci.SportsExerc., vol.21,pp.343,1989. [24]G.A.L.Meijer,Physicalactivity.Implicationsforhumanenergy metabolism,Ph.D.thesis,Univ.ofLimburg,TheNetherlands,1990. [25]D.P.RedmondandF.W.Hegge,Observationsonthedesignandspecicationofawrist-wornhumanactivitymonitoringsystem, Behavior. Res.Methods,Instrum.,Comput., vol.17,pp.659,1985. [26]A.Cappozzo,Lowfrequencyself-generatedvibrationduringambulationinnormalmen, J.Biomech., vol.15,pp.599,1982. [27]E.K.AntonssonandR.W.Mann,Thefrequencycontentofgait, J. Biomech., vol.18,pp.39,1985. [28]A.Bhattacharya,E.P.McCutcheon,E.Shvartz,andJ.E.Greenleaf, BodyaccelerationdistributionandO2uptakeinhumansduringrunning andjumping, J.Appl.Physiol., vol.49,pp.881,1980. [29]M.S.SunandJ.O.Hill,Amethodformeasuringmechanicalwork andworkefciencyduringhumanactivities, J.Biomech., vol.26,pp. 229,1993. [30]M.A.Lafortune,Three-dimensionalaccelerationofthetibiaduring walkingandrunning, J.Biomech., vol.24,pp.877,1991. [31]S.B.Servais,J.G.Webster,andH.J.Montoye,Estimatinghuman energyexpenditureusinganaccelerometerdevice, J.Clin.Eng., vol. 9,pp.159,1984. [32]C.V.C.Bouten,K.R.Westerterp,M.Verduin,andJ.D.Janssen, Assessmentofenergyexpenditureforphysicalactivityusingatriaxial accelerometer, Med.Sci.SportsExerc., vol.26,pp.1516,1994. [33]P.F.M.Schoffelen,W.H.M.Saris,K.R.Westerterp,andF.tenHoor, Evaluationofanautomaticindirectcalorimeterformeasurementof energybalanceinman,in HumanEnergyMetabolism:PhysicalActivity andEnergyExpenditureMeasurementsinEpidemiologicalResearch BaseduponDirectandIndirectCalorimetry, A.J.H.vanEs,Ed.The Hague,TheNetherlands:KoninklijkeBibliotheek,1984;Eur.Nutrition Rep.5,pp.51. [34]J.B.deV.Weir,Newmethodsforcalculatingmetabolicratewith specialreferencetopredictproteinmetabolism, J.Physiol., vol.109, pp.1,1949. [35]A.O.Radke,Man'snewenvironment:Vehiclevibration, Mech.Eng., vol.80,pp.38,1958. [36]S.KitazakiandM.J.Grifn,Adatacorrectionmethodforsurface measurementofvibrationonthehumanbody, J.Biomech., vol.28, pp.885,1995. [37]J.A.Balogun,L.O.Amusa,andI.U.Onyewadume,Factorsaffecting CaltracandCalcountaccelerometeroutput, Phys.Ther., vol.68,pp. 1500,1988. [38]R.A.WashburnandR.E.Laporte,Assessmentofwalkingbehavior: Effectofspeedandmonitorpositionontwoobjectivephysicalactivity monitors, Res.Quart.Exerc.Sport, vol.59,pp.83,1988. [39]J.F.Sallis,M.J.Buono,J.J.Roby,D.Carlson,andJ.A.Nelson, TheCaltracaccelerometerasaphysicalactivitymonitorforschool-age children, Med.Sci.SportsExerc., vol.22,pp.698,1990. [40]J.F.Nichols,P.Patterson,andT.Early,Avalidationofaphysical activitymonitorforyoungandolderadults, Can.J.SportSci., vol.17, pp.299,1992. [41]M.S.Bray,W.W.Wong,J.R.Morrow,N.F.Butte,andJ.M.Pivarnik, Caltracversuscalorimeterdeterminationof24-henergyexpenditurein femalechildrenandadolescents, Med.Sci.SportsExerc., vol.26,pp. 1524,1994. CarlijnV.C.Bouten wasborninVleuten,The Netherlands,onAugust31,1967.Shereceived theM.Sc.degreeinhumanmovementsciences (functionalanatomy,exercisephysiology)fromthe FreeUniversityofAmsterdam,TheNetherlands,in 1991andthePh.D.degreefromtheEindhovenUniversityofTechnology,Eindhoven,TheNetherlands, in1995.ThepresentpaperispartofherPh.D. thesisTheassessmentofdailyphysicalactivityby registrationofbodymovement. Atpresent,sheisapostdoctoralresearchfellow attheDepartmentofMechanicalEngineeringoftheEindhovenUniversityofTechnology.Herresearchincludessofttissuemechanicsandtissue metabolismwithspecialreferencetothepressuresoreproblem. KarelT.M.Koekkoek wasborninBreda,The Netherlands,onFebruary26,1942.Hegraduated withthedegreeinelectricalandelectronicengineeringfromtheCollegeforTechnicalEngineering, Breda,TheNetherlands,in1965. In1967hebecameaTechnicianattheDepartmentofMechanicalEngineeringoftheEindhoven UniversityofTechnology,Eindhoven,TheNetherlands.Hisworkinvolvesthedevelopmentofmechanicalinstrumentsandtechniques,withemphasis ondataacquisitionanddigitalandanaloghardware. MaartenVerduin wasborninMiddelharnis,The Netherlands,onDecember24,1929. In1970hewasregisteredasaProfessionalEngineerinelectricalandelectronicengineering.From 1952to1961hewasanAssistantatthePhilipsResearchLaboratoriesandtheDepartmentofScientic InstrumentsofPhilipsinEindhoven,TheNetherlands.From1961to1972hewasaSeniorAssistant attheDepartmentofMechanicalEngineeringofthe EindhovenUniversityofTechnology,Eindhoven, TheNetherlands.In1972hebecameLecturerat thesamedepartment. RensKodde wasbornonJanuary5,1946in Biggekerke,TheNetherlands.Hegraduatedwith theM.Sc.degreeinelectricalengineeringfromthe EindhovenUniversityofTechnology,Eindhoven, TheNetherlands. HeiscurrentlyaLecturerattheDepartmentof MechanicalEngineeringoftheEindhovenUniversityofTechnology.Hiseldofresearchincludes experimentalmechanics. JanD.Janssen wasborninSpaubeek,TheNetherlands,onAugust31,1940.Hegraduatedinmechanicalengineering(1963)andreceivedthePh.D. degreefromtheEindhovenUniversityofTechnology,Eindhoven,TheNetherlands,in1967. HewasappointedProfessorofContinuumsMechanicsatthesamedepartmentin1968.Since1986, hehasalsobeenpart-timeProfessorofBiomechanicsattheDepartmentofHealthSciencesofthe MaastrichtUniversity,TheNetherlands.Hismain researchinterestslieintheareasofcomputational andexperimentalmechanics,biomechanics,andbiomedicaltechnology. |