"Applied Science, Faculty of"@en . "Engineering, School of (Okanagan)"@en . "DSpace"@en . "UBCO"@en . "Sadi, Fazle"@en . "2012-04-12T16:55:39Z"@en . "2011"@en . "Master of Applied Science - MASc"@en . "University of British Columbia"@en . "Algorithms to detect athletic jumps and to determine in real-time performance parameters such as jump air time (AT), horizontal distance, height\nand drop, are developed in this work. These algorithms are customized to be implemented onboard action sports goggles developed by Recon Instruments Ltd. These goggles are equipped with low cost micro-electro-mechanical inertial sensors and a single point GPS receiver which feed raw data to the algorithms. The micro-LCD display system in the goggles displays jump statistics to the user wearing the goggles. Two novel methods, namely WMCM (Windowed Mean Canceled Multiplication) and PFAD (Preceding and Following Acceleration Difference), are introduced for jump detection using accelerometer data. Four characteristic points in the resultant acceleration data are selected as the AT defining epochs for a jump. A novel threshold independent, probabilistic method using MADM (Multiple\nAttribute Decision Making) and the Closest Peak method are proposed to detect these characteristic points and determine the corresponding AT of a jump. A GPS/INS integration algorithm is developed to determine jump horizontal distance, height and drop. A novel sensor error compensation scheme is developed using sensor fusion and Linear Kalman Filters (LKF). The LKF parameters are varied to address the fluctuating dynamics of the athlete during a jump. The Extended Kalman Filter (EKF) used for GPS/INS integration has an observation vector augmented with sensor error measurements derived from sensor fusion.\n\nThe performance of the proposed algorithms was evaluated through experimental field tests. The proposed jump detection algorithm successfully\ndetected 92% of the jumps performed by a snowboarder wearing the goggles whereas the current Recon algorithm only detects 60%. The AT determination algorithm exhibited an average error of 0.033 s (4.8%) which is well within the accuracy requirement of Recon, \u00B10.1 s, and betters the current Recon algorithm which has an average error of 0.111 s (8.4%). For determination\nof jump horizontal distance, height and drop, the proposed algorithm has an error of 14.34 cm (5.55%), 1.56 cm (38.21%) and 6.71 cm (9.43%)\nrespectively. The accuracy achieved is deemed to fulfill expectations of both recreational and professional athletes."@en . "https://circle.library.ubc.ca/rest/handle/2429/41971?expand=metadata"@en . "Jump Parameter Estimation with Low Cost MEMS Sensors and GPS for Action Sports Goggles by Fazle Sadi B.Sc., Bangladesh University of Engineering and Technology, 2009 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF APPLIED SCIENCE in The College of Graduate Studies (Electrical Engineering) THE UNIVERSITY OF BRITISH COLUMBIA (Okanagan) August 2011 c\u00C2\u00A9 Fazle Sadi 2011 Abstract Algorithms to detect athletic jumps and to determine in real-time per- formance parameters such as jump air time (AT), horizontal distance, height and drop, are developed in this work. These algorithms are customized to be implemented onboard action sports goggles developed by Recon Instruments Ltd. These goggles are equipped with low cost micro-electro-mechanical in- ertial sensors and a single point GPS receiver which feed raw data to the algorithms. The micro-LCD display system in the goggles displays jump statistics to the user wearing the goggles. Two novel methods, namely WMCM (Windowed Mean Canceled Multiplication) and PFAD (Preced- ing and Following Acceleration Difference), are introduced for jump detec- tion using accelerometer data. Four characteristic points in the resultant acceleration data are selected as the AT defining epochs for a jump. A novel threshold independent, probabilistic method using MADM (Multiple Attribute Decision Making) and the Closest Peak method are proposed to detect these characteristic points and determine the corresponding AT of a jump. A GPS/INS integration algorithm is developed to determine jump horizontal distance, height and drop. A novel sensor error compensation scheme is developed using sensor fusion and Linear Kalman Filters (LKF). The LKF parameters are varied to address the fluctuating dynamics of the athlete during a jump. The Extended Kalman Filter (EKF) used for GP- S/INS integration has an observation vector augmented with sensor error measurements derived from sensor fusion. The performance of the proposed algorithms was evaluated through ex- perimental field tests. The proposed jump detection algorithm successfully detected 92% of the jumps performed by a snowboarder wearing the goggles whereas the current Recon algorithm only detects 60%. The AT determi- nation algorithm exhibited an average error of 0.033 s (4.8%) which is well within the accuracy requirement of Recon, \u00C2\u00B10.1 s, and betters the current Recon algorithm which has an average error of 0.111 s (8.4%). For determi- nation of jump horizontal distance, height and drop, the proposed algorithm has an error of 14.34 cm (5.55%), 1.56 cm (38.21%) and 6.71 cm (9.43%) respectively. The accuracy achieved is deemed to fulfill expectations of both recreational and professional athletes. ii Table of Contents Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vi List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii Nomenclature . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . x Notations and Abbreviations . . . . . . . . . . . . . . . . . . . . xvii 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1 Context . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Recon Action Sports Goggles . . . . . . . . . . . . . . . . . . 3 1.2.1 Key Performance Variables for Jumps . . . . . . . . . 4 1.2.2 Jump Categories . . . . . . . . . . . . . . . . . . . . . 7 1.3 Objective . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 1.4 Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 1.5 Thesis Outline . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2 Literature Review . . . . . . . . . . . . . . . . . . . . . . . . . 12 2.1 In Field Sports Performance Monitoring . . . . . . . . . . . . 12 2.1.1 General Strategies . . . . . . . . . . . . . . . . . . . . 12 2.1.2 Potential of Micro-Technology . . . . . . . . . . . . . 13 2.2 Jump Detection and Air Time Calculation . . . . . . . . . . 14 2.2.1 Inertial Sensor Based Techniques . . . . . . . . . . . 15 2.2.2 Non-Inertial Sensor Based Techniques . . . . . . . . . 17 2.2.3 Application of Wavelets . . . . . . . . . . . . . . . . . 18 2.2.4 Recon Application: Feasibility of Current Methods . 19 2.3 Jump Horizontal Distance, Height and Drop Detection Strate- gies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 iii Table of Contents 2.3.1 Imagery Based . . . . . . . . . . . . . . . . . . . . . . 21 2.3.2 Satellite Based . . . . . . . . . . . . . . . . . . . . . . 23 2.3.3 Satellite and Inertial Navigation Based . . . . . . . . 23 2.3.4 Alternative Techniques . . . . . . . . . . . . . . . . . 24 2.4 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 3 Jump Detection and Air Time Determination . . . . . . . . 25 3.1 Sensor Selection . . . . . . . . . . . . . . . . . . . . . . . . . 25 3.2 Accelerometer Data . . . . . . . . . . . . . . . . . . . . . . . 25 3.3 Fast Fourier Transform . . . . . . . . . . . . . . . . . . . . . 26 3.4 Sensor Correlation and Covariance . . . . . . . . . . . . . . . 31 3.5 Windowed Mean Canceled Multiplication . . . . . . . . . . . 36 3.6 Preceding and Following Acceleration Difference . . . . . . . 37 3.7 Jump Detection . . . . . . . . . . . . . . . . . . . . . . . . . 39 3.7.1 Detection through WMCM . . . . . . . . . . . . . . . 40 3.7.2 Validity Check through PFAD . . . . . . . . . . . . . 41 3.7.3 Jump Indication Categories . . . . . . . . . . . . . . . 43 3.8 Air Time Determination . . . . . . . . . . . . . . . . . . . . 51 3.8.1 Characteristic Acceleration Points . . . . . . . . . . . 51 3.8.2 Primary and Secondary Detection . . . . . . . . . . . 53 3.8.3 Window Selection . . . . . . . . . . . . . . . . . . . . 55 3.8.4 Positive and Negative Peaks . . . . . . . . . . . . . . 55 3.8.5 Primary Detection . . . . . . . . . . . . . . . . . . . . 57 3.8.6 Secondary Detection . . . . . . . . . . . . . . . . . . 57 3.9 Sensor and Visual Air Time . . . . . . . . . . . . . . . . . . 65 3.10 Experimental Results and Comparison . . . . . . . . . . . . 66 3.10.1 Field Data Collection . . . . . . . . . . . . . . . . . . 66 3.10.2 Visual AT from SAT . . . . . . . . . . . . . . . . . . 66 3.10.3 Performance Comparison . . . . . . . . . . . . . . . . 68 3.11 Implementation Aspects for Other Devices/Applications . . 71 4 Trajectory Determination . . . . . . . . . . . . . . . . . . . . 75 4.1 Coordinate Frames . . . . . . . . . . . . . . . . . . . . . . . 75 4.2 Rotation Matrix . . . . . . . . . . . . . . . . . . . . . . . . . 76 4.3 Body Frame Attitude Computation . . . . . . . . . . . . . . 79 4.3.1 Rotation Matrix Update . . . . . . . . . . . . . . . . 79 4.3.2 Orientation Vector Computation . . . . . . . . . . . . 80 4.3.3 Coriolis Effect . . . . . . . . . . . . . . . . . . . . . . 81 4.4 Inertial Navigation Equations . . . . . . . . . . . . . . . . . 82 4.5 Rotation Angle Computation . . . . . . . . . . . . . . . . . . 84 iv Table of Contents 4.5.1 Roll and Pitch from Accelerometer . . . . . . . . . . 84 4.5.2 Yaw from Magnetometer . . . . . . . . . . . . . . . . 87 4.5.3 Rotation Angles from Rotation Matrix . . . . . . . . 88 4.6 Sensor Errors . . . . . . . . . . . . . . . . . . . . . . . . . . . 89 4.7 Proposed Algorithm . . . . . . . . . . . . . . . . . . . . . . . 90 4.7.1 Acceleration Error Compensation . . . . . . . . . . . 92 4.7.2 Angular Rate Error Compensation . . . . . . . . . . . 94 4.7.3 INS Mechanization . . . . . . . . . . . . . . . . . . . 97 4.7.4 Integration with GPS . . . . . . . . . . . . . . . . . . 98 4.7.5 Period Segmentation and Parameter Customization . 101 4.8 Field Experiments . . . . . . . . . . . . . . . . . . . . . . . . 104 4.9 Experimental Results . . . . . . . . . . . . . . . . . . . . . . 106 4.10 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116 5 Conclusions and Recommendations . . . . . . . . . . . . . . 117 5.1 Objectives and Algorithm Performance . . . . . . . . . . . . 117 5.2 Methods and Key Findings . . . . . . . . . . . . . . . . . . . 119 5.3 Recommendations . . . . . . . . . . . . . . . . . . . . . . . . 120 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122 Appendix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128 Appendix A Multiple Attribute Decision Making . . . . . . . . 128 v List of Tables 3.1 Two sample t-Test for low impact JIs . . . . . . . . . . . . . 49 3.2 Two sample t-Test for high impact JIs . . . . . . . . . . . . . 50 3.3 Thresholds for secondary detection of JEP and JEN . . . . 60 3.4 Thresholds for secondary detection of JSP and JSN . . . . . 60 3.5 Attribute type and weight for probabilistic secondary detec- tion of JEP and JEN . . . . . . . . . . . . . . . . . . . . . . 63 3.6 Attribute type and weight for probabilistic secondary detec- tion of JSP and JSN . . . . . . . . . . . . . . . . . . . . . . 63 3.7 Comparison of AT among Recon, Ripxx and the proposed algorithm. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 3.8 Comparison of number of false processor wake up calls for Recon algorithm and the proposed algorithm. . . . . . . . . . 73 3.9 Comparison of the number of false detected jumps + true undetected jumps among Recon, Ripxx and the proposed al- gorithm. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 4.1 Quadrant information of roll from acceleration signs. . . . . . 86 4.2 Quadrant information of pitch from acceleration signs. . . . . 87 4.3 Quadrant information of yaw from magnetic field component signs. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88 4.4 Comparison of jump horizontal distance determination per- formance between the proposed algorithm and Catapult GPS solution. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112 4.5 Comparison of jump height determination performance be- tween the proposed algorithm and Catapult GPS solution. . . 113 4.6 Comparison of jump drop determination performance between the proposed algorithm and Catapult GPS solution. . . . . . 114 vi List of Figures 1.1 Snowboard enthusiast catching some \u00E2\u0080\u0098big air\u00E2\u0080\u0099 [1]. . . . . . . . 1 1.2 Recon Goggle. Click on the image to see promotional video [1]. 4 1.3 Micro-LCD display in Recon Action Sports Goggles [1]. . . . 5 1.4 Orientation of the three dimensional body/sensor frame. . . . 6 1.5 Real-time feedback in Recon goggles [1]. . . . . . . . . . . . . 6 1.6 Categories of ski and snowboard jumping. . . . . . . . . . . . 8 2.1 StroMotion highlights the essential moments in an athletic performance. . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 2.2 The movement and form of two performances is compared with SimulCam. . . . . . . . . . . . . . . . . . . . . . . . . . . 22 3.1 Typical acceleration in the body frame. . . . . . . . . . . . . 27 3.2 Resultant acceleration ab of the body frame during jump. . . 28 3.3 Disturbance and head motion in resultant acceleration ab of the body frame during jump. . . . . . . . . . . . . . . . . . . 28 3.4 Sliding windowed FFT of ab (resultant acceleration) near jump start region. Jump starts at 19.93 s. . . . . . . . . . . . 29 3.5 Sliding windowed FFT of ab (resultant acceleration) near jump end region. Jump ends at 21.62 s. . . . . . . . . . . . . 30 3.6 Sliding GT of ab (resultant acceleration) near jump start re- gion. Jump starts at 19.93 s. . . . . . . . . . . . . . . . . . . 31 3.7 Sliding GT of ab (resultant acceleration) near jump end re- gion. Jump ends at 21.62 s. . . . . . . . . . . . . . . . . . . . 32 3.8 Cross-correlation at zero lag of windowed abx and a b y. . . . . . 34 3.9 Covariance at zero lag of windowed abx and a b y. . . . . . . . . 35 3.10 WMCM of abx, a b y and a b z. . . . . . . . . . . . . . . . . . . . . 37 3.11 Windows for PFAD around a point of interest near take-off. . 38 3.12 PFAD value of the resultant acceleration ab. . . . . . . . . . . 39 3.13 Proposed jump detection method. Note: mxyz is shifted by -2 g3 for visual convenience. . . . . . . . . . . . . . . . . . . . 42 3.14 High disturbance JI. . . . . . . . . . . . . . . . . . . . . . . . 44 vii List of Figures 3.15 Low disturbance JI. . . . . . . . . . . . . . . . . . . . . . . . 45 3.16 High impact JI. . . . . . . . . . . . . . . . . . . . . . . . . . . 46 3.17 Low impact JI. . . . . . . . . . . . . . . . . . . . . . . . . . . 47 3.18 Proposed jump detection method. . . . . . . . . . . . . . . . 52 3.19 U-shaped region and characteristic points for an ideal jump. . 53 3.20 U-shaped region and characteristic points for a practical jump. 54 3.21 Window for characteristic acceleration points detection and potential positive peaks. . . . . . . . . . . . . . . . . . . . . . 56 3.22 Primary detection. . . . . . . . . . . . . . . . . . . . . . . . . 58 3.23 (a) Potential positive peaks for JEP . (b) Attribute value of the potential positive peaks. . . . . . . . . . . . . . . . . . . . 59 3.24 Secondary detection of JEP and JEN with threshold based method. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 3.25 D\u00E2\u0088\u0097D\u00E2\u0088\u0092 plane for the secondary detection of JEP . . . . . . . . 64 3.26 DOIRPj and rank of the top candidate peaks for JEP . . . . . 64 3.27 Secondary detection of JEP and JEN with probabilistic method. 65 3.28 Experimental video and captured sensor resultant acceleration. 67 3.29 Video AT vs SAT. . . . . . . . . . . . . . . . . . . . . . . . . 68 3.30 Error in AT comparison among the Recon, Ripxx and pro- posed algorithm. . . . . . . . . . . . . . . . . . . . . . . . . . 70 4.1 Coordinate frames. . . . . . . . . . . . . . . . . . . . . . . . . 76 4.2 Coordinate frame before rotation. . . . . . . . . . . . . . . . . 78 4.3 Coordinate frame after rotations in ZYX sequence where roll = pitch = yaw = 90\u00E2\u0097\u00A6. Click on the image to see sequential rotations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 4.4 Coordinate frame after rotations in ZXY sequence where roll = pitch = yaw = 90\u00E2\u0097\u00A6. Click on the image to see sequential rotations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79 4.5 Roll and pitch calculation from accelerometer data. . . . . . . 85 4.6 Error compensation block. . . . . . . . . . . . . . . . . . . . . 95 4.7 INS mechanization block. . . . . . . . . . . . . . . . . . . . . 95 4.8 Flow diagram of jump trajectory determination. . . . . . . . 99 4.9 Period segmentation for trajectory determination. . . . . . . 102 4.10 Catapult Minimax unit affixed to the biker\u00E2\u0080\u0099s helmet. . . . . . 105 4.11 Catapult Minimax unit strapped to the Novatel GPS receiver antenna pole. [Photograph by Darren Handschuh] . . . . . . 106 4.12 Test jump performed by the athlete. [Photograph by Darren Handschuh] . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107 viii List of Figures 4.13 Resultant acceleration of Catapult unit strapped to the an- tenna pole carried in backpack. . . . . . . . . . . . . . . . . . 108 4.14 First example of jump trajectory determination with the pro- posed algorithm, with the Catapult GPS solution and with the Novatel receiver GPS solution. . . . . . . . . . . . . . . . 109 4.15 Second example of jump trajectory determination with the proposed algorithm, with the Catapult GPS solution and with the Novatel receiver GPS solution. . . . . . . . . . . . . . . . 110 4.16 Error comparison in KPV determination among Catapult GPS solution, GPS/INS integration excluding the error compen- sation block and proposed algorithm [Jump no. 1-7]. . . . . . 115 4.17 Error comparison in KPV determination among Catapult GPS solution, GPS/INS integration excluding the error compen- sation block and proposed algorithm [Jump no. 8-14]. . . . . 115 A-1 TOPSIS method for MADM. . . . . . . . . . . . . . . . . . . 131 A-2 M-TOPSIS method for MADM. . . . . . . . . . . . . . . . . . 132 ix Nomenclature a\u00CC\u0084bIfw Average of the signal ab within the PFAD window following tI a\u00CC\u0084bIpw Average of the signal ab within the PFAD window preceding tI a\u00CC\u0084bJIfw Average of the signal ab within the PFAD window following tJI a\u00CC\u0084bJIhigh Maximum between a\u00CC\u0084 bJI pw and a\u00CC\u0084 bJI fw a\u00CC\u0084bJIlow Minimum between a\u00CC\u0084 bJI pw and a\u00CC\u0084 bJI fw a\u00CC\u0084bJIpw Average of the signal ab within the PFAD window preceding tJI y\u00CC\u00841 Mean of the sample of a\u00CC\u0084 bJI high y\u00CC\u00842 Mean of the sample of a\u00CC\u0084 bJI low \u00CE\u00B1 Integral of the angular rate vector \u00CE\u00B2\u00CF\u0089 Reciprocal of the process correlation time \u00CF\u0084\u00CF\u0089 \u00CE\u00B2a Reciprocal of the process correlation time \u00CF\u0084 a \u00CF\u0089 Angular rate vector \u00CF\u0089k True angular rate of body frame with respect to navigation frame at epoch tk \u00CF\u0083 Orientation vector \u00CF\u0083\u00CF\u0089 Amplitude of the angular rate bias power spectral density \u00CF\u0083a Amplitude of the acceleration bias power spectral density \u00CF\u0084\u00CF\u0089 Process correlation time of angular rate bias \u00CF\u0084 a Process correlation time of acceleration bias \u00E2\u0088\u0086t Sampling interval x Nomenclature \u00CE\u00B4xk+1 State vector of the Extended Kalman Filter (EKF) at epoch tk+1 \u00CE\u00B4zk+1 Observation vector of the EKF at epoch tk+1 \u00CE\u00B6\u00CC\u0082k estimated attitude error vector in Inertial Navigation Systems (INS) mechanization for update cycle starting at epoch tK a\u00CC\u0082k+1 Final estimate of the acceleration for update cycle starting at epoch tk+1 a\u00CC\u0082ck+1 Initial estimate of the acceleration for update cycle starting at epoch tk+1 b\u00CC\u0082 \u00CF\u0089E k Gyroscope bias error estimate from the EKF for GPS/INS integration in cycle starting at tk b\u00CC\u0082 aE k Acceleration bias error estimate from the EKF for GPS/INS integra- tion in cycle starting at tk b\u00CC\u0082 \u00CF\u0089L k+1 Gyroscope bias error estimate derived via sensor fusion in the cycle starting at tk+1 b\u00CC\u0082 aL k+1 Acceleration bias error estimate derived via sensor fusion in cycle starting at tk+1 R\u00CC\u0082k final estimation of rotation matrix for update cycle starting at tk r\u00CC\u0082nk+1 Estimated position vector in navigation frame from GPS/INS inte- gration at epoch tk+1 \u00CE\u00BB Longitude \u00C2\u00B51 Mean of the parameter a\u00CC\u0084 bJI high \u00C2\u00B52 Mean of the parameter a\u00CC\u0084 bJI low \u00E2\u0088\u0087tlce Extension parameter of the characteristic window end \u00E2\u0088\u0087tlcs Extension parameter of the characteristic window start \u00E2\u0084\u00A6E Earth\u00E2\u0080\u0099s rotation rate \u00CF\u0086 Roll \u00CF\u0088 Yaw xi Nomenclature ab Acceleration of body frame observed from body frame ab The resultant acceleration in the body frame an Acceleration of body frame observed from navigation frame ak True acceleration at epoch tk b\u00CF\u0089k Bias in angular rate at epoch tk bak Bias in acceleration at epoch tk gn Normal gravity gnl Plumb bob gravity mb Earth\u00E2\u0080\u0099s magnetic field vector observed from body frame qak Process noise of acceleration with covariance matrix Q a k QEk Covariance matrix of the EKF at epoch tk Rnb Rotation matrix to convert from body frame to navigation frame R\u00CF\u0086\u00CE\u00B8 Rotation matrix to make the x-y plane of body frame coplanar with the x-y plane of local level frame rnk+1 True position vector in navigation frame at epoch tk+1 rnINSk+1 Position vector in navigation frame derived from INS equations at epoch tk+1 ubk Incremental velocity of body frame observed from body frame for update interval starting at tk unk Incremental velocity of body frame observed from navigation frame for update interval starting at tk vn Velocity of the body frame with respect to the earth expressed in local level frame vnk+1 True velocity vector in navigation frame at epoch tk+1 vnINSk+1 Velocity vector in navigation derived from INS equations at epoch tk+1 xii Nomenclature w\u00CF\u0089k Measured noise in angular rate at epoch tk with covariance matrix R\u00CF\u0089k wak Measured noise in acceleration at epoch tk with covariance matrix Rak w\u00CF\u0089Lk+1 Measurement noise of LKF 3 with covariance matrix R\u00CF\u0089Lk+1 waLk+1 Measurement noise of LKF 2 with covariance matrix RaLk+1 zack+1 Observation vector for LKF 1 zbaLk+1 Observation vector for LKF 2 z\u00CF\u0089\u00CF\u0089Lk+1 Observation vector for LKF 3 \u00CE\u00B8 Pitch \u00CF\u0089\u00CC\u0083k Measured body frame angular rate with respect to navigation frame at epoch tk \u00CF\u0089\u00CC\u0083bibk+1 Gyroscope reading for the update cycle starting at tk+1 a\u00CC\u0083k Measured acceleration by accelerometer at epoch tk m\u00CC\u0083k Measured magnetometer readings at epoch tk R\u00CC\u0083k initial estimation of rotation matrix for update cycle starting at tk \u00CF\u0095 Latitude \u00CE\u00B6 Significance level of the t\u00E2\u0088\u0092Test A\u00E2\u0088\u0097 Positive ideal solution of TOPSIS A\u00E2\u0088\u0092 Negative ideal solution of TOPSIS Aj j th alternative of the MADM process abx The forward/backward (anteroposterior) axis acceleration in the body frame aby The right/left (mediolateral) axis acceleration in the body frame abz The vertical (down/up) axis acceleration in the body frame C\u00E2\u0088\u0097j Similarities to positive ideal solution of j th alternative xiii Nomenclature cwxy(kT ) Covariance of windowed a b x and a b y with center epoch w and at lag k D Decision matrix D\u00E2\u0088\u0097j Distance from the positive ideal solution to j th alternative D\u00E2\u0088\u0092j Distance from the negative ideal solution to j th alternative DOIRPj Distance from the OIRP to j th alternative dab(tI) PFAD of a b at tI F Sampling frequency JEN Jump end negative peak JEP Jump end positive peak JSN Jump start negative peak JSP Jump start positive peak l Extension counter index of the characteristic window LIfw Length of the PFAD window following tI LIpw Length of the PFAD window preceding tI LJIfw Length of the PFAD window following tJI LJIpw Length of the PFAD window preceding tJI mxyz(t j s + nT ) WMCM of abx, a b y and a b z for a window starting at t j s N Number of samples in the data capture window njmax Index of the maximum WMCM value for the jth window n1 Sample size of a\u00CC\u0084 bJI high n2 Sample size of a\u00CC\u0084 bJI low NLIfw Number of samples within a window of length L I fw NLIpw Number of samples within a window of length L I pw xiv Nomenclature rwxy(kT ) Cross-correlation of windowed a b x and a b y with center epoch w and at lag k R0(\u00CF\u0095) Earth\u00E2\u0080\u0099s radius at latitude \u00CF\u0095 Re Earth\u00E2\u0080\u0099s equatorial radius Rp Earth\u00E2\u0080\u0099s polar radius rji Vector normalized decision matrix element for j th alternative and ith attribute S21 Sample variance of a\u00CC\u0084 bJI high S22 Sample variance of a\u00CC\u0084 bJI low S2p Estimate of the common variance of a\u00CC\u0084 bJI high and a\u00CC\u0084 bJI low T Sampling period tje Ending epoch of the jth data capture window tjs Beginning epoch of the jth data capture window tNc Epochs of negative peak in characteristic window tPc Epochs of positive peak in characteristic window tjumpe Jump end epoch tAGPs Start epoch of Acceleration Gain Period (AGP) tjumps Jump start epoch tce End epoch of the characteristic window tcs Start epoch of the characteristic window tJI Epoch of detected JI thm Threshold for WMCM to detect JI vji Weighted normalized rating for j th alternative and ith attribute w Center epoch of the data capture window wci Weight of i th attribute xv Nomenclature Xh Earth\u00E2\u0080\u0099s magnetic field component along X-axis of the realigned body frame xbtilt The angle between horizontal plane and X b Yh Earth\u00E2\u0080\u0099s magnetic field component along Y -axis of the realigned body frame ybtilt The angle between horizontal plane and Y b I Identity matrix tNJE Epoch of JEN tPJE Epoch of JEP tNJS Epoch of JSN tPJS Epoch of JSP xvi Notations and Abbreviations Convention Matrices are represented by capital case bold letters and vectors are rep- resented by lower case bold letters. A superscript in a vector represents the particular frame in which the vector is represented. Rotation matrices between two coordinate systems are defined with a su- perscript and a subscript denoting the two coordinate systems., e.g. Rnb transforms from body frame b to navigation frame n. Angular rate between two frames is expressed by a angular rate vector, e.g. \u00CF\u0089bib represents the angular rate between the inertial frame and the body frame observed from the body frame. The skew symmetric matrix form of any vector \u00CF\u0089 = [\u00CF\u00891\u00CF\u00892\u00CF\u00893] is represented as [\u00CF\u0089\u00C3\u0097] = \u00EF\u00A3\u00AE\u00EF\u00A3\u00B0 0 \u00E2\u0088\u0092\u00CF\u00893 \u00CF\u00892\u00CF\u00893 0 \u00E2\u0088\u0092\u00CF\u00891 \u00E2\u0088\u0092\u00CF\u00892 \u00CF\u00891 0 \u00EF\u00A3\u00B9\u00EF\u00A3\u00BB . (1) Abbreviations KPV Key Performance Variable AT Air Time TDR Total Degrees of Rotation INS Inertial Navigation Systems GPS Global Positioning System xvii Notations and Abbreviations MEMS Micro-Electro-Mechanical Systems IMU Inertial Measurement Unit FFT Fast Fourier Transform PSD Power Spectral Density HMM Hidden Markov Model MADM Multiple Attribute Decision Making WMCM Windowed Mean Canceled Multiplication LKF Linear Kalman Filter EKF Extended Kalman Filter FFT Fast Fourier Transform PFAD Preceding and Following Acceleration Difference JI Jump Indication OIRP Optimized Ideal Reference Point SAT Sensor Air Time IAP Initial Acceleration Period AGP Acceleration Gain Period APP Aero Phase Period xviii Chapter 1 Introduction 1.1 Context Many sports enthusiasts crave the sensation of \u00E2\u0080\u0098catching the air\u00E2\u0080\u0099 while participating in such sports as snowboarding and skiing. \u00E2\u0080\u0098Catching air\u00E2\u0080\u0099 is a phrase often used to describe what happens when a snowboarder or skier travels fast enough over a jump to allow their board or skis to lift off the snow and into the air. Well executed and stylish routines consisting of com- plex aerial acrobatic maneuvers, as in Figure 1.1, are also highly appreciated in both sports. During competitions, judges often score competitors based upon the elevation, distance and time elapsed during the loft motion when \u00E2\u0080\u0098catching air\u00E2\u0080\u0099. The interest in witnessing and maneuvering complex jumps becomes apparent from the regular exclamation after a jump that he \u00E2\u0080\u0098caught\u00E2\u0080\u0099 some \u00E2\u0080\u0098big sky\u00E2\u0080\u0099, \u00E2\u0080\u0098big air\u00E2\u0080\u0099 or other assorted remarks when referring to the ele- vation, distance or time elapsed of the loft motion. Unfortunately, both the spectators and athletes themselves typically only have a qualitative sense as to the performance variables (e.g. jump height) without ever quantitatively knowing the statistics. However, many believe that objectivity in evaluation of these acrobatic performances could significantly enhance the experience for these athletes. Enhancement of experience and entertainment is not the only motiva- Figure 1.1: Snowboard enthusiast catching some \u00E2\u0080\u0098big air\u00E2\u0080\u0099 [1]. 1 1.1. Context tion for finding objective information for evaluating sports performance. In competition environments, performances are also evaluated based upon sub- jective measures, referred to as \u00E2\u0080\u0098overall impression\u00E2\u0080\u0099, by a panel of judges [2]. During training sessions for these competitions, in the absence of appropriate measurement tools, athletes and coaches much rely largely upon evaluating performances based on their own experiences and feelings. Recent analysis has revealed that sport specific Key Performance Vari- ables (KPVs), e.g. loft duration and elevation for snowboard and ski jump- ing, strongly correlates with the subjectively judged scores in competition. Therefore, judges, coaches, and athletes could be greatly assisted by any complementary tools that objectively measure the sport specific KPVs. Even though video based assessment is widely used for training and evalu- ation, it has often proved to be misleading. Furthermore, video recording based performance analysis is generally time consuming and requires consid- erable ground infrastructure (e.g. video setup). It is also often very difficult to acquire and analyze information on the KPVs through the labor intensive post processing of the video data [3]. For instance, sophisticated video pro- cessing software is needed to determine the KPVs from the recorded data which is both resource and time consuming. In more recent years, satellite based positioning, such as Global Posi- tioning System (GPS), has been developed as a performance assessment technique to quantitatively measure sport specific variables. There are a number of commercial products available on the market, e.g. [4, 5], which provide performance parameters, such as distance, velocity and height pro- files derived from satellite based position fixes. While these products may be effective for certain parameters, coaches and athletes are often also in- terested in other parameters such as orientation, degree of rotation, and epochs of take-off and landing which cannot be fully measured by satellite based monitoring systems. For example, in many instances the athlete\u00E2\u0080\u0099s environment may have naturally occurring obstructions that could block or attenuate satellite signals. To over come many of these variables, which affect the ability to mea- sure KPVs, Inertial Measurement Units (IMUs) are now being used along with GPS units for sports monitoring. Conventional GPS/INS equipment comprising dual-frequency GPS receivers and tactical grade INS, can pro- vide highly accurate data (cm for positioning, cm/s for velocity and 1/100 degree for orientation). In the context of sports applications, however, such equipment is often not suitable because of adverse weight (a few kg) and expense (>= e40, 000). Recently, utilizing the advancements of micro tech- nology, Micro-Electro-Mechanical Systems (MEMS) IMUs as well as low cost 2 1.2. Recon Action Sports Goggles L1 GPS receivers has entered the world of sports monitoring [6\u00E2\u0080\u00939]. The ben- efit of such a GPS/INS system includes its global availability, light weight, low cost and fast setup. This system, therefore, diminishes the necessity for external infrastructure, such as video cameras [3]. The idea of measuring forces and position fix parameters with a GP- S/INS integrated system to provide quantitative insight into athlete\u00E2\u0080\u0099s motion dynamics is academically sound. However, the assumption then declaration that this raw information will directly help the end users in functional eval- uation and training seems to be an optimistic statement [2]. In addition, the interpretation of the system\u00E2\u0080\u0099s gathered information into something suit- able for use by athletes and coaches constitutes a major component of any research focus. It is also imperative that consistency between the sport- ing communities expectations and the focus of scientific research be reached for the successful integration of any technological innovation into sports. Consequently, strong rationale has been provided for researchers to explore algorithms and systems that can bring cutting edge technology to the service of sports communities. Recon Instruments, a Vancouver based company, has recently developed and is making Recon Action Sports Goggles, which are able to impart real- time information on KPVs to the individual wearing the goggles. A built-in MEMS GPS/INS unit is used to generate this data. The research conducted for this thesis is part of the Recon goggle development project, conducted in collaboration between The University of British Columbia and Recon Instru- ments. The goal of this collaborative study was to develop an innovative algorithm to derive sport specific KPVs in real-time onboard the goggles with the best achievable accuracy using the available sensors as well as to provide instantaneous feedback to individuals wearing the goggles. 1.2 Recon Action Sports Goggles Recon Instruments has developed the world\u00E2\u0080\u0099s first GPS/INS enabled goggles with a head-mounted display system (Figure 1.2). The micro- LCD at the lower portion of the goggles, as shown in Figure 1.3, displays performance-enhancing data and statistics on KPVs. The graphics and op- tics of this direct-to-eye communication technology has made the display unobtrusive for front and peripheral vision. Furthermore, little interaction is needed during use to view the information making it applicable for use in fast-paced environments. The goggles can be used both in entertainment situations (for general interest during recreational skiing/snowboarding), or 3 1.2. Recon Action Sports Goggles Figure 1.2: Recon Goggle. Click on the image to see promotional video [1]. as a tool to aid in training for competitions. The Recon goggles have sensors including a GPS receiver, a three axis MEMS accelerometer, a three axis MEMS gyroscope, a three axis MEMS magnetometer and an altimeter. The body frame coordinate system is shown in Figure 1.4, and the sensors axes are congruent with the body frame. The update rate of the GPS receiver position fixes is 1 Hz, and the data capture rate of the MEMS sensors is 100 Hz. The entire system is controlled with an ARM9TM processor. Currently, these goggles provide real-time feedback including speed, lat- itude/longitude, altitude, vertical distance traveled, total distance traveled, chrono/stopwatch mode, a run-counter, temperature and time as shown in Figure 1.5. The goggles can be charged by USB, which can also be used to transfer data, along with the post-processing software included. 1.2.1 Key Performance Variables for Jumps Recon Instruments has specially developed these goggles with sport per- formance quantifying capabilities specially relevant to the motions involved in snowboarding and skiing jumps. Four KPVs have been selected as the 4 1.2. Recon Action Sports Goggles Figure 1.3: Micro-LCD display in Recon Action Sports Goggles [1]. 5 1.2. Recon Action Sports Goggles Vertical Mediolateral Anteroposterior X b Y b Z b Figure 1.4: Orientation of the three dimensional body/sensor frame. Figure 1.5: Real-time feedback in Recon goggles [1]. 6 1.2. Recon Action Sports Goggles standard parameters of interest for the aerial acrobatic routines performed by skiers or snowboarders, especially during competitions. The KPVs are defined as follows: 1. Air Time. Total time traveled through the air from the time of take- off to the time of landing. 2. Horizontal Distance. The length over land traveled through the air from the time of take-off to the time of landing. 3. Height. Total vertical distance achieved from the point of take-off to the apex of the jump. 4. Drop. Total vertical distance achieved from the apex to the point of landing. 1.2.2 Jump Categories Recon Instruments has also categorized the different types of ski and snowboard jumping into five major families. Each jump category possesses distinct characteristics with respect to the KPVs previously mentioned. The jump categories are illustrated in Figure 1.6 and defined as follows: 1. Ollie. Height and drop are equal, and the goggle wearer achieves distance and Air Time (AT). 2. Step-up. Height is larger than drop due to elevation change during the jump. AT and distance are both achieved, but AT is truncated. 3. Standard. This is the most common jump, where drop is larger than height due to elevation drop (common on jumps in mountain terrain). Longer AT and distances are achieved. 4. Half-pipe. AT is very large, and distance is very minimal. Drop is larger than height. 5. Cliff-drop. This is where there is very little to no height achieved, and drop is very large. AT is also very big, but distance is minimal. 7 1.2. Recon Action Sports Goggles Drop Horizontal Distance Height Air Time Horizontal Distance Drop Height Air Time Drop Horizontal Distance Height Air Time Horizontal Distance Drop Height Air Time Horizontal Distance Drop Air Time Jump Trajectory Jump Start Jump End Jump Apex Legend1. OLLIE 2. STEP-UP 3. STANDARD JUMP 5. CLIFF-DROP4. HALF-PIPE Figure 1.6: Categories of ski and snowboard jumping. 8 1.3. Objective 1.3 Objective The objective of this research is to develop an algorithm for the on- board goggle processor to analyze all available sensor data. By doing so, the processor will then be able to successfully detect each of the above men- tioned jump categories, with a minimum of false positives, and determine the standard KPVs as listed in Subsection 1.2.1. Recon has set the maxi- mum allowable error in AT to be \u00C2\u00B10.1 s. For the other KPVs, Recon has set no specific requirements. However, errors should be negligible in magnitude in comparison with the KPVs measured such that the KPVs have practical application to the user. The key factors involved in this research are as follows: \u00E2\u0080\u00A2 Sensor Fusion. The fundamental problem is how to combine the various sensor data, each of which has its own particular sensitivity, update rate, and noise, to achieve reliable conclusions. Whereas cur- rent work has concentrated on processing data from only two types of sensors, this project will combine data from five types of sensors - accelerometers, gyroscopes, magnetometers, GPS and altimeter. This task is complicated by the variability of situations a skier or snow- boarder encounters over the course of a typical day with the goggles on. Therefore, the processing algorithm must cope with a large variety of dynamics. \u00E2\u0080\u00A2 Minimal Resources. Since battery life is limited, power manage- ment is an important consideration in the goggle design. For instance, the minimum update rate from each the sensors required to reliably detect jumps must be determined. The lower the update rate, the greater the potential for battery and memory savings. \u00E2\u0080\u00A2 Online Processing. The algorithm developed to process all sensor data and compute the parameters must be suitable for real-time oper- ation and instantaneous feedback by the ARM9TM processor onboard the Recon goggles. The limited computing resources and battery life on the Recon goggles is a challenging factor in algorithm development. \u00E2\u0080\u00A2 Head Motion. Since all sensors used are mounted within the goggles, these sensors will be subject to movement of the user\u00E2\u0080\u0099s head. In other documented work, sensors are primarily mounted on the user\u00E2\u0080\u0099s back, waist, legs or underfoot. Independent head motion could complicate (or aid) the interpretation of data from head-mounted sensors. 9 1.3. Objective \u00E2\u0080\u00A2 Unconventional Jumps. Apart from the standard jumps, there are also many different ways that a skier or snowboarder can leave the ground that must be taken into account during jump detection. For instance, in some cases, the skier may never rise above the elevation of the take-off point. Flips, half-pipe moves, and other jump types each can result in radically different sensor readings. \u00E2\u0080\u00A2 Low Data Rate. One of the major hardware limitations that the goggles have is the low rate (1 Hz) of position fixes from the GPS receiver whereas the typical ski jump duration is 0.5\u00E2\u0088\u00922.0 s. Due to the scarcity of available data in such a brief operation time, conventional GPS dependent tracking algorithms are not suitable. \u00E2\u0080\u00A2 Noise. MEMS sensors produce very noisy data. Depending upon the noise characteristics, denoising algorithms are also often computation- ally heavy. Therefore, focus should be given to the tradeoff between computational cost and denoising benefits. \u00E2\u0080\u00A2 Poor Sensitivity. The on board altimeter is not sensitive enough to accurately capture typical jump height or drop of 0 \u00E2\u0088\u0092 2.0 m. The GPS receiver also does not have enough accuracy and sensitivity in the vertical axis. Consequently, none of the sensors can be used directly to compute the KPVs. \u00E2\u0080\u00A2 Accuracy. The algorithm needs to satisfy the accuracy criteria de- manded by sport professionals and recreational athletes. Recon speci- fied that a \u00C2\u00B10.1 s error accuracy for AT must be met. Although Recon has not set any accuracy specification for jump horizontal distance, height and drop, error in the measurement should be constrained to what may render the estimated parameters to be useful. \u00E2\u0080\u00A2 Fast and Easy Operation. Coaches and athletes often lack knowl- edge about emerging technologies. In addition, skiers and snowboard- ers often prefer to execute jumps without any interruptions. Therefore, it is important that the algorithm is independent of any tedious and repetitive routines such as manual sensor calibration, bias estimation or initial waiting period for filter convergence, which could interrupt one\u00E2\u0080\u0099s performance. 10 1.4. Methodology 1.4 Methodology \u00E2\u0080\u00A2 Jump Detection. For jump detection, accelerometer data has been used by the novel Windowed Mean Canceled Multiplication (WMCM) method as well as by the novel Preceding and Following Acceleration Difference (PFAD) method. \u00E2\u0080\u00A2 Air Time. To determine the AT, two methods are proposed. One is a threshold based technique and another is a probabilistic approach that uses no threshold. The novel probabilistic approach uses Multiple Attribute Decision Making (MADM). \u00E2\u0080\u00A2 Horizontal Distance, Height and Drop. To determine the jump horizontal distance, height and drop, GPS aided and GPS integrated INS algorithms have been proposed. The GPS aided INS algorithm uses Linear Kalman Filters (LKFs) for continuous recursive error cor- rection. For, the integration of GPS with INS, the EKF has been explored. Both the aided and integrated algorithm have been devel- oped with direct and indirect state estimation techniques. 1.5 Thesis Outline The thesis is structured chronologically. Chapter 2 summarizes all the relevant current research and development in sports monitoring. A detailed discussion is presented on the feasibility of the current technology satisfy- ing the requirements of the Recon goggles. In Chapter 3, the developed jump detection and AT determination algorithm is presented with practi- cal results. Chapter 4 presents the detailed architecture of the algorithm developed for determining jump horizontal distance, height and drop with performance evaluation. Concluding remarks and potential future research is presented in Chapter 5. 11 Chapter 2 Literature Review The research in determining KPVs for Recon Action Sports Goggles comprises several fields of study. For jump detection and Air Time (AT) determination, the main focus is on the dynamics of the signals during jump start, aero phase and jump end epochs available from the sensors. On the other hand, the research in determining horizontal distance, height and drop mostly involves the kinematics of the sensor data. These two fields of study are distinct in signal processing, tools of manipulation and denoising aspects. Therefore, the relevant previous and current research to determine the men- tioned KPVs will be discussed in two main categories. Firstly, the prevailing jump detection and AT calculation strategies will be presented. Secondly, the shortcomings and feasibility of these current methods for the Recon goggles application will be discussed. Finally, the current work pertinent to horizontal distance, height and drop determination will be presented in this chapter. The application of INS and GPS is widely practiced in similar applications. Therefore, a priority will be given to the possibilities and im- plications of the present INS and GPS involved strategies to meet the Recon goggle feature requirements. Preceding these discussions, general monitor- ing strategies for acrobatic maneuvers and the advantages of MEMS sensors in this particular application will also be presented. 2.1 In Field Sports Performance Monitoring 2.1.1 General Strategies Research has shown that subjectively judged scores for athletes doing acrobatic maneuvers in sports like skiing and snowboarding has strong cor- relation with the KPVs such as AT, Total Degrees of Rotation (TDR), jump height, drop, horizontal distance, etc. As laboratory performance assess- ment for these sports is not feasible, a number of attempts have been made to develop and use current technologies to quantify sports specific perfor- mance variables in the field. Recently, sports scientists from the Australian Institute of Sport (AIS) have undertaken video based analysis of KPVs as- 12 2.1. In Field Sports Performance Monitoring sociated with elite half-pipe snowboard competitions [10]. The video based analysis is currently the most readily available method for coaches, sports scientists, competition judges and recreational users to gather the most ob- jective information from acrobatic maneuvers as possible and assess style and run execution. However, video based analysis poses intrinsic time de- lay in information feedback and is a labour intensive form of analysis. Due to the requirement for synchronization of cameras, exterior orientation and determination of the cameras\u00E2\u0080\u0099 position, all video based systems require a large amount of infrastructure as well as considerable setup time. Video based analysis has, therefore, not been adapted for everyday recreational use as well as for training purposes, despite its ability to provide objective observations. Novel strategies have been evolving to allow sensors to be more capable of interacting with the physical world, and to provide sport statistics in a convenient way. An electromagnetic tracking system has been introduced in [11] to quan- tify the kinetics and kinematics associated with snowboarding. The joint moment torque and angular displacement information related to snowboard turns is monitored by post-processing the data captured on snow by the electromagnetic tracking system. When the object is placed inside vary- ing yet controlled magnetic fields, voltages are induced in the sensor coils. These induced voltages are used by the measurement system to calculate the position and orientation of the object [12]. Additional research has been conducted to measure the kinematics of the skier and snowboarder feet and legs during maneuvers using a similar electromagnetic motion tracking sys- tem [13]. A dynamometric platform has been capitalized for field based data acquisition while snowboarding [14]. Without any interaction with current snowboard binding technology this dynamometric platform measures all the load components transmitted between the boots and snowboard bindings. Therefore, this technology facilitates the capability to acquire quantitative information on the torques and forces imparted by the boot on the bindings. 2.1.2 Potential of Micro-Technology The research in the previously mentioned field based assessment methods is a worthwhile goal for collecting information on the kinematics of technical skill based sports, such as skiing and snowboarding. Unfortunately, certain intrinsic drawbacks preclude its growth and the wide application of these methods. The bulky electromagnetic tracking system [11] attached to the athlete\u00E2\u0080\u0099s leg and the 10 kg measurement system which is carried by the researcher and tethered to the athlete is quite an obtrusive measurement 13 2.2. Jump Detection and Air Time Calculation system. The weight of this monitoring system on the subject alone would be enough to constrict and alter their normal movement patterns, removing any relevance to real life snowboarding or skiing techniques. Therefore, a heavy monitoring system with tethered requirements is unsuitable for any data collection for elite level athletes training or competing in aerial acrobatic events with elements of risk involved [2]. With the advent of small scale technologies, a promising alternative, MEMS has drawn major attention from the sports science community due to its light weight and unobtrusive nature. MEMS sensors are capable of measuring human motion thousands of times per seconds in multiple axes without the problems associated with size and weight. Micro sensors are easy to attach anywhere on the athlete and do not limit the natural mo- bility. Moreover, onboard data storage technology along with the MEMS sensors diminishes the necessity for tethering bulky equipment. Therefore, MEMS inertial sensors, such as triaxial accelerometer, triaxial gyroscope, triaxial magnetometer etc., possess enormous potential in monitoring nat- ural, unaltered and unobstructed acrobatic maneuvers of athletes in high risk sports in a realistic manner. There are challenges involved in using the capabilities of MEMS IMUs to calculate reliable and accurate objective information relevant to the end user [2]. Furthermore, limited processing capability and battery power cause the challenges to increase for small size wearable integrated systems dedicated to online processing, such as the Re- con goggles. 2.2 Jump Detection and Air Time Calculation One of the primary features of the Recon goggles is the ability to detect jumps while skiing and snowboarding as well as during acrobatic maneuvers in similar sports. These features are also essential for determining the ex- pected KPVs. Once the jump start and jump end epochs are determined, the duration of the aero phase of the athlete, i.e. AT, can be derived easily. A number of different types of sensors are being used to detect the occur- rence of the jump and determine the associated AT. Among MEMS inertial sensors, the accelerometer is significantly more popular for this purpose due to its capability of capturing the changes in an athlete\u00E2\u0080\u0099s dynamics related to jumps. 14 2.2. Jump Detection and Air Time Calculation 2.2.1 Inertial Sensor Based Techniques A significant amount of work has been done by Australian sports engi- neers to provide sport specific KPV information automatically and imme- diately post run for skiers and snowboarders. These developed technologies allow coaches and competition judges to automatically and objectively as- sess athletic performance in the environment in which half-pipe snowboard- ing takes place [2, 10, 15, 16]. For example, a signal processing technique was developed to calculate the AT associated with each individual aerial ac- robatic manoeuvre performed during all half-pipe snowboard runs detected in [2]. MEMS accelerometers were used as data capturing sensors. All ac- celerometer units underwent a static calibration in three axes (up/down, forward/backward and left/right) prior to each data collection session align- ing each axes of sensitivity with and against the direction of gravity. As reported in [2], half-pipe snowboarding generates a very distinct raw triaxial accelerometer data trace. The forward/backward and up/down acceleration axes were reported to display elevated accelerations during the aero phase of aerial acrobatic maneuvers. As a result of the substantial increases and decreases in acceleration throughout the performance, the power levels were reported to rise during the aero phase. Thus, jump locations were detected and jump windows were extracted by analyzing the entire trace for elevated levels of acceleration and utilizing a sliding Fast Fourier Transform (FFT) window and subsequent power analysis of average power levels. A threshold based algorithm was used to identify the jumps within a trace by evaluating each FFT window for frequencies ranging from 0.25 \u00E2\u0088\u0092 0.85 Hz (as aerial acrobatic maneuvers occur relatively rhythmically with a period of 1.2 \u00E2\u0088\u0092 4.0 s) [2]. For detecting the AT, a two pass method was proposed in this study. In the first pass, the accelerometer data was quan- tized into three levels/states, namely maximum, minimum and transition. In the second pass, a technique was applied to discard false jumps detected by focusing on the sport specific possible durations of aerial acrobatic ma- neuvers. The duration of any component between 0.8\u00E2\u0088\u00922.2 s was considered to be a valid aerial acrobatic maneuver. It was claimed by Harding et al. that this range would cover most aerial acrobatic ATs completed by half- pipe snowboard athletes. However, the threshold was selected heuristically for the algorithm and no general threshold selection scheme was presented in [2]. In a similar study, a preliminary automated feedback system based upon MEMS sensors was designed to calculate objective information on these sport specific KPVs [10]. The authors, Harding et al., argued that the results 15 2.2. Jump Detection and Air Time Calculation would provide practical benefits for elite half-pipe snowboard training and current subjective judging protocols. According to the method for detecting flight phase (i.e. the aero phase of the jump) in ski jumping proposed in [17], acceleration data, sampled at 200 Hz was stored in the digital memory and a FFT analysis was conducted on the overlapping Hanning windows comprising 256 points. It was reported that the maximum amplitude of the acceleration at the initial flight phase was clearly large (over 2.5 G) in each local coordinate system. In this study, angular rate/velocity of the athlete was also used to detect jumps. The peak angular velocity, captured by MEMS gyroscopes, at the initial flight phase was reported to be produced 0.2 s later than the peak acceleration. The peak power spectrums for acceleration in three directions just after take-off were observed at 9.37 Hz. Peak power spectrums for angular velocity just after take-off were reported to be produced at the frequencies of 9.37 Hz in the anteroposterior axis of the trunk, and at 8.59 Hz in the transverse axis and the longitudinal axis. These peak power spectrums were granted as the evidence of jump or flight of the athletes. In [18], an unorthodox way of sensing loft was proposed. Single and multiple accelerometers were used to sense the vibration created by a snow- board, skis or mountain bike, moving along a surface. While moving, the voltage output provided by the accelerometer generates a noisy spectrum over time. It was reported that the vibrational spectrum generated during the aero phase for the athlete is generally much smoother than the spectrum generated while in motion on the surface. This is because while in air, the sporting equipment (i.e. skis or snowboard) is not subjected to the random vibrations of the road or ski slope. Accordingly, this relatively smooth spec- trum was discerned from the rest of the spectrum by the microprocessor subsystem and evaluated to generate the AT. The time difference between the epochs of the end and start of the smooth spectrum was deemed as the airtime by the authors. Because the vibrational spectrum is affected by the particular activity of a user, e.g. standing still, included with the micro- processor subsystem of the invention was a means for assessing boundary conditions of the spectrum. For excluding certain conditions from determin- ing the air/loft time, a lower boundary of 500 ms and an upper boundary of 5 s were fixed. Any events with elapsed time duration outside this boundary condition are excluded from the determination of loft time. 16 2.2. Jump Detection and Air Time Calculation 2.2.2 Non-Inertial Sensor Based Techniques A number of non-inertial sensor based techniques for jump and airtime detection were proposed in [18]. The authors, Fleantov et al., proposed loft sensors constructed with the following technologies: \u00E2\u0080\u00A2 a microphone assembly that senses a noise spectrum \u00E2\u0080\u00A2 a switch that is responsive to the weight of the user of the vehicle \u00E2\u0080\u00A2 a voltage-resistance sensor that generates a voltage indicative of a vehicle\u00E2\u0080\u0099s speed. The loft sensors proposed in this study sense a spectrum of information, e.g. a vibrational or sound spectrum and the microprocessor subsystem determines the relative change in the spectrum of information. Further- more, the microprocessor interprets this relative change in the spectrum to determine the loft time. For example, the microphone based embodiment detects sound waves and provides a voltage output accordingly. Similar to the accelerometer based system described earlier, the microphone senses the vibration of the vehicle, e.g. snowboard or skis, moving along the surface. From the smoother spectrum during aero phase, the loft sensor senses the jump and determines the AT. Another embodiment of the loft sensing system that is proposed in [18] is a switch that rests below the boot of the skier and that senses the pressure caused by the weight of the skier. When the skier is on the ground, the boot squeezes the switch, thereby closing it. The closed switch sends a distinct input to the microprocessor subsystem. When the skier jumps into the air, the switch opens up and renders the microprocessor to detect the openness of the switch. When the skier lands on the ground the switch closes again. The duration of the switch state being open is determined by the subsystem to calculate the loft time. In the same paper [18], Fleantov et al. proposed a pad that can be placed under the skier\u00E2\u0080\u0099s boot and changes capacitance as a function of the change in applied pressure. The capacitance changing element is then connected to a circuit, incorporating a monostable multivibrator through a resistor. The value of the resistor is fixed while the variable capacitance causes the pulse train generated by the monostable multivibrator to change its frequency depending on the pressure exerted on the pad. When the pulse train repe- tition rate increases, the value of the capacitance decreases and the skier\u00E2\u0080\u0099s boot applies less pressure on the pad. This event marks the start of a jump 17 2.2. Jump Detection and Air Time Calculation and beginning of the AT measurement. When the pulse repetition rate de- creases, meaning a increase in capacitance, this indicates that the boot is applying greater pressure on the pad. This event marks the end of the jump as well as the end of the AT measurement. 2.2.3 Application of Wavelets The acceleration of a snowboarder or skier at the begining of a jump (take-off for the flight) and at the end of the jump (landing of the flight), re- sults in a sudden change in acceleration. The detection of this abrupt change in signals and/or systems is a classic problem in signal processing, especially in the vast research field of detecting singularities or non-stationarity in a stationary signal. Cusp or jump location identification techniques of any signal encompasse a wide variety of signal processing tools and demand very customized and application specific solutions. Among the various sig- nal processing tools, wavelet analysis has proven to be useful in irregularity detection of a signal for numerous areas of application. Wavelet coefficients generally exhibit high peaks near abrupt changes or cusps of the underlying signal. Therefore, even in the presence of additive noise, by checking the absolute value of wavelet coefficients it is possible to detect discontinuities in an otherwise smooth curve. The problem of detecting a mean value jump of a stationary random pro- cess was addressed in [19] by means of the wavelet transform. A matched filter based approach was also proposed in this study for which an ideal knowledge of the Power Spectral Density (PSD) of the signal was needed. However, this is not generally the case. Therefore, a second approach was proposed using the classic pattern recognition principle relying on the corre- lation between the signal and the signature. Typically, wavelet based statis- tical signal processing techniques model the wavelet coefficients as indepen- dent or jointly Gaussian. To address the non-Gaussian statistics encountered in real world signals, a new framework for statistical signal processing based on the wavelet domain Hidden Markov Model (HMM) was proposed in [20]. The efficient expectation maximization algorithm was also developed in this study for fitting the HMM to observational signal data. The new framework was claimed to be suitable for a wide range of applications, including signal estimation, detection, classification, prediction, and even synthesis. A method of detecting the cusp or an incident in the traffic flow was presented in [21]. Unlike the existing wavelet incident detection algorithm, where the wavelet technique was utilized to denoise data before the data was input into an algorithm, features that were extracted from traffic mea- 18 2.2. Jump Detection and Air Time Calculation surements by using wavelet transformation were directly utilized in detect- ing changes in traffic flow. Using the Peaks Over Threshold modeling of the noisy wavelet transformations and extreme value theory, a method for detecting an arbitrary number of discontinuities in an unknown function observed with noise was proposed in [22]. 2.2.4 Recon Application: Feasibility of Current Methods To develop the jump detection and AT calculation algorithm for the Re- con goggles, extensive insight can be achieved from the prevailing relevant studies as previously discussed, even though most of the techniques were devoted to particular applications different from the application considered here. However, the majority of the inventions were developed using com- mon platforms, e.g. FFT. Unfortunately, none of the current techniques meet Recon\u00E2\u0080\u0099s requirements entirely. The significant characteristics of cur- rent methods which are not conducive with the application of Recon goggles are as follows: 1. The Recon goggles demand a online (real-time and onboard) algorithm to detect jumps which should be computationally light to save battery power and instantaneous feedback. An FFT algorithm for jump de- tection, as proposed in [2, 10, 15, 17], requires continuous heavy pro- cessing of the acceleration data. Therefore, FFT analysis is not likely to be a good option for any online algorithm required by Recon even if it is useful for offline processing . 2. The use of PSD extracted from a sliding FFT window [2, 10, 15, 17] for detecting jump and AT is heavily dependent upon the window size and its position relative to the abrupt change in acceleration. The expected high power may not even be captured if the window size is not precisely accurate. However, no criterion was proposed for choos- ing the proper window length in the related papers. Moreover, high power level, due to the sudden shift in acceleration, may continue for a number of consecutive FFT windows. In such cases, determination of the proper epoch for jump start or jump end becomes ambiguous. These phenomena will be discussed in detail in Chapter 3. 3. In most of the accelerometer data based jump detection and AT deter- mination techniques, sensor calibration and bias correction is neces- sary. However, sensor calibration before using a device like the Recon goggles is very inconvenient for recreational users, even for coaches and professional athletes. 19 2.2. Jump Detection and Air Time Calculation 4. Almost all the current methods use stringent thresholds for the ac- celeration or post processed data, which causes the algorithm to lose generality. For fixed thresholds to work, bias correction and calibra- tion is also necessary. This is not suitable for devices like the Recon goggles which are targeted for users in different environmental condi- tions with ease of operation. Some of the studies, e.g. [2], used a very low threshold so that no event is missed. However, implications such as too many false detections are there for these sort of strategies. 5. Most of MEMS accelerometer based techniques [2, 10, 15, 17] work with raw acceleration data captured with individual accelerometers aligned in specific directions with respect to the athlete\u00E2\u0080\u0099s body. At- taching the accelerometers in exact alignment with direction is itself fastidious. Moreover, the orientation of the body part to which the sensors are attached is not the same for all body parts even for any particular acrobatic maneuver. The magnitude of the readings from the sensor may drastically change with a slight change from the ex- pected alignment of the sensors. This phenomenon poses high risk of failure to detect events for the fixed threshold and/or individual axis acceleration based algorithms. 6. Strategies like quantization, as in [2], can also cause information loss when trying to find the exact epoch for jump start or end. 7. Certain boundary conditions, e.g 0.8\u00E2\u0088\u0092 2.2 s in [2], 0.5\u00E2\u0088\u0092 5.0 s in [18], were set to distinguish jumps from other events, such as abrupt head movements. However, small jumps with AT less than these lower boundaries are reported frequently, especially for nonprofessional ath- letes. Therefore, the use of these boundary conditions will limit the use of the Recon goggles. 8. Examples shown in [2, 10, 15, 17, 18] are mostly for very well behaved data. However, it can be observed that low cost IMU sensors produce very noisy readings. It is quite common for the low cost IMU sensors to produce readings with true features buried under the noise. Movement artifacts also arise frequently from the rapid movements of athletes to maintain balance or to add an extra component of style to their acrobatic performance. All these undesired signal attributes should be accounted for to develop any robust algorithm. 9. In [18], it was assumed that the spectrum of acceleration is always smooth during the flight phase of the athlete. However, this is often 20 2.3. Jump Horizontal Distance, Height and Drop Detection Strategies not the case for the jumps with complex acrobatic maneuvers such as spinning and flipping. The extra acrobatic components during the aero phase will render the IMU sensors to lose their spectral smoothness. 10. The biggest challenge in implementing any wavelet based algorithm is the computational demand of the wavelet decomposition. For online processing of the KPVs it is required that the computation tool be power efficient with very little time lag in producing feedback which is not always the case for wavelet analysis. 2.3 Jump Horizontal Distance, Height and Drop Detection Strategies The task of determining the horizontal distance, height and drop of jumps for outdoor sports disciplines actually falls under the research for trajectory determination of the athlete during the loft period. Trajectory determination in sports application is relatively more explored than the AT detection discussed earlier. Even though every sport has its own method for performance analysis depending on the unique dynamics, ergonomic re- quirements, and trajectory determination techniques of the sport, most have developed with similar techniques. The general strategies of trajectory de- termination are discussed below. 2.3.1 Imagery Based Dartfish [23] has developed a powerful tool using cameras to provide qualitative information about athletic trajectories. Their StroMotion tech- nology is able to create trajectory video footage revealing the evolution of an athlete\u00E2\u0080\u0099s movement, technique, execution and tactics. It is possible to superimpose one athletic movement over another for precise comparisons using the SimulCam technology developed by Dartfish. The developed sys- tem utilizes invariant image retrieval techniques, which detect the common features in the images captured for two different athletes, so that the images of the two different trajectories can then be superimposed and compared. Figure 2.1 and Figure 2.2 show examples of these two technologies. These methods are only useful for relative analysis as quantitative data can not be extracted. However, stereo-camera systems are able to provide highly accurate, quantitative information about the 3D trajectory and motion of athletes [3]. In [24], trajectory analysis was performed using single cam- era and stereo-camera systems. Codamotion has developed systems with 21 2.3. Jump Horizontal Distance, Height and Drop Detection Strategies Figure 2.1: StroMotion highlights the essential moments in an athletic performance. Figure 2.2: The movement and form of two performances is compared with SimulCam. infrared LEDs which work as active markers that are detectable by the cameras [25]. Passive marker systems for optical or infrared cameras have also been developed [26]. Only the markers are visible on the image in an infrared system. An advanced markerless system has been developed by Organicmotion [27]. For this system, volume models are fitted specially to the athlete. Each of these unique systems (active marker, passive marker and markerless) [25\u00E2\u0080\u009327], are able to capture the 3D motion of the subject which can be used to derive the movement and trajectory of the athlete. 22 2.3. Jump Horizontal Distance, Height and Drop Detection Strategies 2.3.2 Satellite Based KPVs such as jump horizontal distance, height and drop are directly related to the position information of the athlete. Therefore, satellite based positioning has been widely used in sports applications. Many commercial products are available that can provide the satellite navigation tracking information and/or KPVs derived from satellite positioning [4, 5, 28\u00E2\u0080\u009333]. GPS has been used for trajectory determination in multiple sports disciplines such as skiing, ski jumping and snowboarding [3, 34, 35], water sports [36, 37] and many others. Unfortunately, GPS signal reception is not ensured for athletes perform- ing maneuvers in adverse environments where satellites may get masked due to naturally occurring obstructions, such as cliffs, and hills. Therefore, in most sport disciplines, continuous observation of the athlete\u00E2\u0080\u0099s performance cannot be guaranteed. Moreover, the accuracy of low cost single point re- ceivers does not always meet the requirements. For example, measuring jump height or drop with the Recon goggles is not possible because of the poor accuracy of the integrated GPS receiver along the vertical axis. Highly accurate dual-frequency GPS receivers are also not feasible for most sports applications because of their larger size and cost. In addition, the low up- date rate of the GPS receiver precludes the measurement of KPVs during relatively short acrobatic maneuvers in such sports as ski jumping or snow- board half-pipe. For example, the 1 Hz update rate of the single point GPS receiver available in the Recon goggles is not suitable for measuring the horizontal distances for jumps with small AT such as \u00E2\u0088\u00BC 0.5 s. Even for long jumps, the low update rate of the GPS receiver is unable to guarantee to record position fixes near the jump start or the jump end epochs with acceptable accuracy. 2.3.3 Satellite and Inertial Navigation Based Although a MEMS IMU cannot provide standalone navigation, INS with GPS aiding or GPS integration can provide continuous reliable tracking for sports applications. Integration of INS with GPS potentially bridges the GPS outages along with compensating for the systematic inertial sensor er- rors. The combination of the high short-term (relative) accuracy of the INS results and longterm (absolute) accuracy of GPS results smears the varia- tion in positioning performance and provides data at a high rate including good orientation estimates [3]. Therefore, GPS/INS integration has been widely used for reliable tracking such as in [38\u00E2\u0080\u009342]. However, GPS/INS in- 23 2.4. Summary tegration incurs a heavy computational load on the microprocessor of the system which is not always affordable. 2.3.4 Alternative Techniques There are a few radio technologies that can be used for position fixing instead of satellite navigation such as Ultra-Wide Band (UWB) [43, 44] and Wireless Local Area Networks (WLAN) [45]. These systems can perform well for indoor applications where the GPS signal is unavailable. Unfortu- nately, these technologies only work for a short range and the infrastructure needed for them is considerable. Moreover, shadowing by the human body and multipath propagation limit the use of these technologies for perfor- mance evaluation of athletes [3]. 2.4 Summary Upon reviewing the above-mentioned systems it becomes apparent that none of these systems or their methods can be seen as a complete solution for the Recon goggle requirements. Indeed, every system and strategy has its advantages and drawbacks for sport application. In terms of available re- sources, observed parameters and desired accuracies, however, a GPS aided inertial sensor based power efficient algorithm is the most suitable solution for the Recon goggles. 24 Chapter 3 Jump Detection and Air Time Determination In this chapter, sensor signals are investigated to develop an algorithm to detect jumps and the corresponding Air Time (AT) involved. The com- monly applied FFT is also explored to detect jump indicative changes in the sensor signals. A novel Windowed Mean Canceled Multiplication (WMCM) method to detect jump occurrence is proposed. Another novel method to determine the AT is presented in the later part of the chapter. Experimental results and comparative analysis among the proposed algorithm, the current Recon goggle algorithm and the Ripxx [4] algorithm are presented. 3.1 Sensor Selection Among all the sensors available in the Recon goggles, the accelerometers and gyroscopes measure the dynamics involved in the motion of the athlete. Accelerometers capture the change in acceleration involved in jump take- off and jump landing. The change in acceleration is very consistent and significant from jump to jump. On the other hand, the gyroscopes are able to record the high rate of angular movement of the head that is likely to occur during the start and end of the loft motion. However, the high angular rate of head movement is not guaranteed for all jumps. Moreover, the peak angular rate tends to vary its relative position with respect to jump end or start from jump to jump. Therefore, the accelerometers are deemed to be the most suitable sensor for jump detection and AT determination. 3.2 Accelerometer Data Ideally, in static conditions such as standing still on the ground, the sensor body frame experiences antigravity (vertically upward), which has the magnitude of 1 g, where g = 9.81 ms\u00E2\u0088\u00922. While in free fall, no acceleration, including antigravity, is experienced by the sensor body frame. Therefore, 25 3.3. Fast Fourier Transform in an ideal noiseless condition the resultant body frame acceleration should be 1 g in static condition and 0 g in free fall. These ideal scenarios will help understanding the accelerometer data recorded while ski jumping. An example of typical acceleration data captured by the triaxial MEMS accelerometer during a ski jump is shown in Figure 3.1. This data was recorded with Recon goggles which were mounted on the skier\u00E2\u0080\u0099s head. The acceleration along the forward/backward (anteroposterior), right/left (medi- olateral) and down/up (vertical) axes in the body frame are termed as abx, aby and a b z respectively. It is significant in the figure that a b x and a b y exhibit a relatively low value (\u00E2\u0088\u00BC 0 g) during the aero phase, i.e. while the athlete is in the air. Except for the aero phase these sensors report significant accel- erations which cause the resultant acceleration to be near 1 g. The sudden change in accelerations are also noticeable at the start and end of the aero phase. It should be noted that the accelerations exhibited by the three axes of the body frame depend entirely on the orientation of the body frame. It is not guaranteed that any specific body frame axis will capture the sudden change in acceleration during the loft motion. However, the abrupt change in acceleration is always captured in the total resultant acceleration ab as shown in Figure 3.2, where ab = \u00E2\u0088\u009A (ab2x + a b2 y + a b2 z ) . The resultant acceleration falling near 0 g while the skier is in the air, provides an indication of the free fall state or, in other words, jump occur- rence. However, due to the presence of noise, head motion and disturbance signals from other activities, the detection of this drop in acceleration is not straight forward. A typical signal due to disturbance and head motion is shown in Figure 3.3. The similarities among the signal patterns during the time of these unwanted motions and jump are clearly noticeable. A signifi- cant amount of noise is also present in these low cost MEMS sensor signals. 3.3 Fast Fourier Transform The sliding windowed Fast Fourier Transform (FFT) is widely proposed to detect the abrupt change in acceleration as in [2, 10, 15, 17]. Figure 3.4 shows the sliding windowed FFT of the jump start region of the resultant acceleration shown in Figure 3.2. Here, the FFT window length is 0.6 s and each window overlaps with the following by 0.4 s. The sampling rate of the data is 100 Hz. The time stamp indicated above each plot represents the epoch at which the sliding window begins. The FFT for several windows preceding and following the jump start region are shown. Even though fluctuating spectrums at higher frequencies are noticeable in the windows 26 3.3. Fast Fourier Transform \u00E2\u0088\u00922 0 2 4 6 A m pl itu de [g ] Acceleration a x b \u00E2\u0088\u00921 0 1 2 3 A m pl itu de [g ] Acceleration ay b 6 8 10 12 14 16 18 20 22 24 26 \u00E2\u0088\u00922 0 2 4 Time [s] A m pl itu de [g ] Acceleration a z b Jump occurence region Figure 3.1: Typical acceleration in the body frame. 27 3.3. Fast Fourier Transform 6 8 10 12 14 16 18 20 22 24 26 \u00E2\u0088\u00921 0 1 2 3 4 5 Time [s] A m pl itu de [g ] Jump end region Jump start region Jump occurrence region Figure 3.2: Resultant acceleration ab of the body frame during jump. 87 88 89 90 91 92 93 94 95 96 97 0 1 2 3 4 Time [s] A m pl itu de [g ] Signal due to disturbance and head motion Signal due to jump Figure 3.3: Disturbance and head motion in resultant acceleration ab of the body frame during jump. 28 3.3. Fast Fourier Transform 0 20 40 0 0.5 1 1.5 2 Frequency (Hz) A m pl itu de [g ] Starts at 19.19s 0 20 40 Frequency (Hz) Starts at 19.39s 0 20 40 Frequency (Hz) Starts at 19.59s 0 20 40 Frequency (Hz) Starts at 19.79s 0 20 40 Frequency (Hz) Starts at 19.99s 0 20 40 0 0.5 1 1.5 2 Frequency (Hz) A m pl itu de [g ] Starts at 20.19s 0 20 40 Frequency (Hz) Starts at 20.39s 0 20 40 Frequency (Hz) Starts at 20.59s 0 20 40 Frequency (Hz) Starts at 20.79s 0 20 40 Frequency (Hz) Starts at 20.99s Jump start window Figure 3.4: Sliding windowed FFT of ab (resultant acceleration) near jump start region. Jump starts at 19.93 s. neighboring that which contains the jump start, no distinct pattern in the spectrum can be observed for the jump start window. Similar observations can also be made for the sliding windowed FFT of the jump end region as shown in Figure 3.5. The Gabor Transform (GT) [46] is explored as an alternative for the regular windowed FFT. In the GT, the data is multiplied with a Gaussian window before applying the FFT, whereas in the regular windowed FFT the data can be multiplied with a square window with the unity magnitude before applying the FFT. Figure 3.6 and 3.7 show the results of the sliding GT of the jump start and end region of the resultant acceleration shown in Figure 3.2. The spectrums are smoother than what is observed for the FFT. Still, the ripples at higher frequencies exhibit a random behavior near the jump start/end region without showing any consistent detectable pat- 29 3.3. Fast Fourier Transform 0 20 40 0 0.5 1 1.5 2 A m pl itu de [g ] Starts at 20.99s Frequency (Hz) 0 20 40 Starts at 21.19s Frequency (Hz) 0 20 40 Starts at 21.39s Frequency (Hz) 0 20 40 Starts at 21.59s Frequency (Hz) 0 20 40 Starts at 21.79s Frequency (Hz) 0 20 40 0 0.5 1 1.5 2 Frequency (Hz) A m pl itu de [g ] Starts at 21.99s 0 20 40 Frequency (Hz) Starts at 22.19s 0 20 40 Frequency (Hz) Starts at 22.39s 0 20 40 Frequency (Hz) Starts at 22.59s 0 20 40 Frequency (Hz) Starts at 22.79s Jump end window Figure 3.5: Sliding windowed FFT of ab (resultant acceleration) near jump end region. Jump ends at 21.62 s. 30 3.4. Sensor Correlation and Covariance 0 20 40 0 0.1 0.2 0.3 0.4 0.5 0.6 A m pl itu de [g ] Starts at 19.19s Frequency (Hz) 0 20 40 Starts at 19.39s Frequency (Hz) 0 20 40 Starts at 19.59s Frequency (Hz) 0 20 40 Starts at 19.79s Frequency (Hz) 0 20 40 Starts at 19.99s Frequency (Hz) 0 20 40 0 0.1 0.2 0.3 0.4 0.5 0.6 Frequency (Hz) A m pl itu de [g ] Starts at 20.19s 0 20 40 Frequency (Hz) Starts at 20.39s 0 20 40 Frequency (Hz) Starts at 20.59s 0 20 40 Frequency (Hz) Starts at 20.79s 0 20 40 Frequency (Hz) Starts at 20.99s Jump start window Figure 3.6: Sliding GT of ab (resultant acceleration) near jump start re- gion. Jump starts at 19.93 s. tern. Even if the window for jump start/end is detected through the FFT or GT, the jump start/end epoch can only be detected with the very poor accuracy of the sliding window length, which, in this case is 0.6 s. The fluc- tuation spectrum patterns exhibited in and neighboring the jump start/end windows also vary significantly with the window size chosen and jump du- ration. Therefore, besides the heavy computational load, the inconsistent spectrum and poor accuracy provides enough rationale to discard the sliding windowed FFT or GT as a tool for jump detection for the Recon goggles. 3.4 Sensor Correlation and Covariance In steady motion or static condition of the athlete the recorded acceler- ations in the three orthogonal axes of the body frame are understandably 31 3.4. Sensor Correlation and Covariance 0 20 40 0 0.1 0.2 0.3 0.4 0.5 0.6 Frequency (Hz) A m pl itu de [g ] Starts at 20.99s 0 20 40 Frequency (Hz) Starts at 21.19s 0 20 40 Frequency (Hz) Starts at 21.39s 0 20 40 Frequency (Hz) Starts at 21.59s 0 20 40 Frequency (Hz) Starts at 21.79s 0 20 40 0 0.1 0.2 0.3 0.4 0.5 0.6 Frequency (Hz) A m pl itu de [g ] Starts at 21.99s 0 20 40 Frequency (Hz) Starts at 22.19s 0 20 40 Frequency (Hz) Starts at 22.39s 0 20 40 Frequency (Hz) Starts at 22.59s 0 20 40 Frequency (Hz) Starts at 22.79s Jump end window Figure 3.7: Sliding GT of ab (resultant acceleration) near jump end region. Jump ends at 21.62 s. 32 3.4. Sensor Correlation and Covariance uncorrelated due to the random nature of noise, head motion and fluctuation in the accelerations. However, it is virtually guaranteed that all three ac- celerometers in three orthogonal directions will experience a sudden impact if there is an abrupt change in motion dynamics. During the take-off for a jump, forces exerted by the athlete and the subsequent free fall state cause a rapid change in the athlete\u00E2\u0080\u0099s dynamics. Hence, all the accelerometers are expected to experience a peak (positive or negative) value during this brief period of sudden change in acceleration. With similar reasoning, all the accelerometers are also expected to record a local maximum or minimum during the landing of the jump. Assuming this rationale to be valid, the following analysis of the accelerometer data is conducted. The cross-correlation of two discrete time signals, e.g. abx and a b y, can be defined as rxy(kT ) = +\u00E2\u0088\u009E\u00E2\u0088\u0091 n=\u00E2\u0088\u0092\u00E2\u0088\u009E abx(nT )a b y(nT \u00E2\u0088\u0092 kT ), k = 0,\u00C2\u00B11,\u00C2\u00B12, .... (3.1) Here, the sampling period T = 1F and the sampling frequency F = 100 Hz. The index k is the discrete lag parameter and kT is the time shift. For online signal processing, the data is captured as a window. Therefore, cross- correlation of the two signals captured in a window from tjs to t j e can be expressed as rwxy(kT ) = N\u00E2\u0088\u00921\u00E2\u0088\u0091 n=\u00E2\u0088\u0092(N\u00E2\u0088\u00921) abx(t j s + nT )a b y(t j s + nT \u00E2\u0088\u0092 kT ), k = 0,\u00C2\u00B11, ...,\u00C2\u00B1(N \u00E2\u0088\u0092 1) (3.2) where tjs = beginning epoch of the jth data capture window, (integer multiple of T) tje = ending epoch of the jth data capture window, (integer multiple of T) N = number of samples in the data capture window = (tje \u00E2\u0088\u0092 tjs)/T + 1 w = center epoch of the window = (tje + t j s)/2. The purpose here is to capture the epoch of the highest correlation be- tween any of the two accelerometer signals, given that one expects two or more accelerometers to be highly correlated at jump start and end. The 33 3.4. Sensor Correlation and Covariance 6 8 10 12 14 16 18 20 22 24 26 \u00E2\u0088\u009220 \u00E2\u0088\u009210 0 10 20 30 40 50 60 Center epoch of the cross\u00E2\u0088\u0092correlation window, w [s] Co rre la tio n va lu e at z er o la g, rw x y(0 ) [ g2 ] Peak near landing Jump occurrence region Peak near take\u00E2\u0088\u0092off Figure 3.8: Cross-correlation at zero lag of windowed abx and a b y. region of simultaneous high correlation in all the acceleration signals is very brief in span. Therefore, the cross-correlation at zero lag, i.e. rwxy(0), is the main focus of interest for each window. Figure 3.8 shows the crosscorre- lation of abx and a b y of the data shown in Figure 3.1. For convenience, the cross-correlation values are plotted against the center epoch of the pertinent window. The window size, i.e. tje\u00E2\u0088\u0092tjs, is 0.6 s. The window size is empirically chosen to produce best results. The peak values at the jump start and end region are clearly noticeable in this cross-correlation plot. Hence it proves the hypothesis that the acceleration signals are highly correlated during the start and end of the aero phase due to the abrupt change in the dynamics involved. In Figure 3.8, undesired fluctuations with considerable magnitude are present before and after the jump region. The nonzero average value of the sensor signal causes this high cross-correlation value. Moreover, as ac- celerometer values are very low during a jump, it is possible to miss a cross- correlation peak if one of the accelerometer values is essentially zero. To overcome these difficulties, covariance of the sensor signals may be used. 34 3.4. Sensor Correlation and Covariance 6 8 10 12 14 16 18 20 22 24 26 \u00E2\u0088\u00928 \u00E2\u0088\u00926 \u00E2\u0088\u00924 \u00E2\u0088\u00922 0 2 4 6 8 Center epoch of the covariance window, w [s] Co va ria nc e va lu e at z er o la g, cw x y(0 ) [ g2 ] Peak near take\u00E2\u0088\u0092off Relatively smoother region neighbouring jump Jump occurrence region Peak near landing Figure 3.9: Covariance at zero lag of windowed abx and a b y. Covariance of the two discrete time signal abx and a b y can be defined as cwxy(kT ) = N\u00E2\u0088\u00921\u00E2\u0088\u0091 n=\u00E2\u0088\u0092(N\u00E2\u0088\u00921) { abx(t j s + nT )\u00E2\u0088\u0092 1 N N\u00E2\u0088\u00921\u00E2\u0088\u0091 n=0 abx(t j s + nT ) } \u00C2\u00B7 { aby(t j s + nT \u00E2\u0088\u0092 kT )\u00E2\u0088\u0092 1 N N\u00E2\u0088\u00921\u00E2\u0088\u0091 n=0 aby(t j s + nT ) } , (3.3) k = 0,\u00C2\u00B11, ...,\u00C2\u00B1(N \u00E2\u0088\u0092 1). (3.4) Figure 3.9 shows the covariance of abx and a b y at zero lag, i.e. c w xy(0) versus the center epoch of the window. Similarly, in this plot, the peaks near jump start and end are prominently distinguishable. Moreover, it is apparent that the region neighboring the jump is smoother than what is seen for the cross- correlation plot. Therefore, the covariance of the sensor signals would be a better tool for jump detection. The presence of a peak with high magnitude is also noticed after the peak near take-off. These scenarios may appear due to the two extremes (one minimum and one maximum) experienced by the accelerometers involved in jump take-off and landing. The acceleration experienced by the body frame depends entirely on the orientation and direction of the forces involved in the jump. Consequently, 35 3.5. Windowed Mean Canceled Multiplication the covariance will differ significantly depending on the selection of two sensor signals for any particular jump. Moreover, using the covariance value, it is only possible to predict the jump start and end epochs with an accuracy of the covariance window length. Window length cannot be reduced, as it will diminish distinct features of covariance near the jump region. Exploiting the concept of covariance, a novel intuitive method is presented in the next section, which can subdue the above mentioned problems. 3.5 Windowed Mean Canceled Multiplication From the previous discussion, it is understood that all three accelerome- ter signals should be involved in the jump detection algorithm to overcome the problem pertaining to the orientation uncertainty of the body frame. An- other key factor is that the high peak value of the covariance is contributed to predominately by a very small number of epochs which fall within the brief time span of the event \u00E2\u0080\u009Cjump start\u00E2\u0080\u009D or \u00E2\u0080\u009Cjump end\u00E2\u0080\u009D. Based on these ratio- nale, the novel Windowed Mean Canceled Multiplication (WMCM) method is proposed. WMCM of the three acceleration signals abx, a b y and a b z captured in the jth window spanning from tjs to t j e is defined as mxyz(t j s + nT ) = { abx(t j s + nT )\u00E2\u0088\u0092 1 N N\u00E2\u0088\u00921\u00E2\u0088\u0091 p=0 abx(t j s + pT ) } \u00C2\u00B7 { aby(t j s + nT )\u00E2\u0088\u0092 1 N N\u00E2\u0088\u00921\u00E2\u0088\u0091 p=0 aby(t j s + pT ) } \u00C2\u00B7 { abz(t j s + nT )\u00E2\u0088\u0092 1 N N\u00E2\u0088\u00921\u00E2\u0088\u0091 p=0 abz(t j s + pT ) } , n = 0, 1, 2, ..., N \u00E2\u0088\u0092 1. (3.5) Since the multiplied values are not added together for a particular win- dow in WMCM, the information on epochs with the most correlations is not lost. The involvement of all three axes accelerations also ensures that the change in dynamics is captured regardless of sensor orientation. Figure 3.10 shows the WMCM values of the acceleration data set shown in Figure 3.1. It is evident from the plot that near the take-off and landing regions the WMCM value exhibits a very large magnitude response. WMCM values other than in the jump region are very low in magnitude without any con- siderable fluctuation. Therefore, with a safe threshold for WMCM, jump 36 3.6. Preceding and Following Acceleration Difference 6 8 10 12 14 16 18 20 22 24 26 \u00E2\u0088\u00920.5 0 0.5 1 1.5 2 2.5 3 3.5 t s + nT [s] W M CM v al ue , m w x yz (t s + n T) [g 3 ] High WMCM value due to other activities Take\u00E2\u0088\u0092off region Jump occurrence region Landing region Figure 3.10: WMCM of abx, a b y and a b z. occurrence can be detected with a good estimate of the jump start and end epochs. However, due to disturbance activities and abrupt head motion a high value in WMCM can be observed as also shown in Figure 3.10. To verify the validity of any high value of WMCM indicating jump occurrence, another complementary method is presented in the next section. 3.6 Preceding and Following Acceleration Difference In a static state or in motion with little change in the dynamics, the resultant acceleration ab of an athlete exhibits noisy characteristics with a steady mean. Therefore, the average values of the resultant acceleration over consecutive windows should not differ by much if there is no abrupt change in the resultant acceleration. However, at the jump occurrence region the average value of the resultant acceleration falls drastically. If a time epoch tI is selected near the take-off or landing region for a b, the neighboring accelerations preceding and following tI are likely to differ significantly. This difference in the preceding and following accelerations should serve as an effective tool for jump detection. Preceding and Following Acceleration Difference (PFAD) is computed considering a certain window around the point of interest. Figure 3.11 shows the PFAD windows of ab for a point of interest at epoch tI near the take- 37 3.6. Preceding and Following Acceleration Difference 6 8 10 12 14 16 18 20 22 24 26 \u00E2\u0088\u00921 0 1 2 3 4 5 Time [s] A m pl itu de o f a b [g ] Point of interest at epoch tI Following window of length LIfw a bI pw a bI fwPreceding window of length LIpw Figure 3.11: Windows for PFAD around a point of interest near take-off. off region. LIpw and L I fw are the lengths of the windows preceding and following tI respectively. The average resultant acceleration captured within the windows preceding and following tI is termed a\u00CC\u0084 bI pw and a\u00CC\u0084 bI fw respectively. PFAD at tI is the difference of these two averages. Therefore, the PFAD of ab for any point of interest at epoch tI , where tI is an integer multiple of T, can be defined as dab(tI) = NLIpw\u00E2\u0088\u0091 n=1 ab(tI + nT ) NLIpw \u00E2\u0088\u0092 NLIfw\u00E2\u0088\u0091 n=1 ab(tI \u00E2\u0088\u0092 nT ) NLIfw = a\u00CC\u0084bIpw \u00E2\u0088\u0092 a\u00CC\u0084bIfw (3.6) where NLIpw = L I pw/T and NL I fw = L I fw/T are the number of samples within the windows of length LIpw and L I fw respectively. It is clear from the captured acceleration magnitude that the average resultant acceleration within the preceding window a\u00CC\u0084bIpw, will be significantly greater than the average resultant acceleration within the following window a\u00CC\u0084bIfw. The opposite scenario will be observed if tI is near the landing region. Therefore, dab(tI) should be a maximum if tI is near the take-off and a minimum if tI is near the landing. At any other position except these two, the PFAD value dab(tI) is not expected to exhibit any considerable mag- nitude. Figure 3.12 shows the PFAD values of the resultant acceleration 38 3.7. Jump Detection 6 8 10 12 14 16 18 20 22 24 26 \u00E2\u0088\u00922 \u00E2\u0088\u00921.5 \u00E2\u0088\u00921 \u00E2\u0088\u00920.5 0 0.5 1 1.5 2 tI [s] PF A D v al ue o f t he re su lta nt a cc el er at io n d ab (t I ) [ g] Maxima after landing Minima near landing Minima before take\u00E2\u0088\u0092off Jump occurrence region Maxima near take\u00E2\u0088\u0092off Figure 3.12: PFAD value of the resultant acceleration ab. shown in Figure 3.2. For this example, the window lengths LIpw and L I fw are set to 0.6 s. Selection of this window length depends on the jump category, which is explained in Subsection 3.7.3. A maximum and a minimum with large magnitude are evident near the jump start and the jump end respec- tively. Another minimum just before the take-off and another maximum just after the landing are also evident. These peak PFAD values appear due to the gain in acceleration (mostly vertically upward and forward) just before take-off and the gain in acceleration (mostly vertically upward and backward) just after landing. Apart from the jump occurrence region, the PFAD value dab(tI) is generally small and varies little. Therefore, PFAD provides a powerful tool for detecting any jump occurrence. To be more specific, from the PFAD value it is possible to determine whether tI is near the jump start or the jump end. If dab(tI) is positive with a significant mag- nitude, then tI is near the jump start or the take-off. On the other hand, if tI is negative with a significant magnitude, then tI can be determined to be near the jump end or the landing. 3.7 Jump Detection Using the novel concepts of WMCM and PFAD a jump detection al- gorithm is proposed in this section. The algorithm consists of two main 39 3.7. Jump Detection steps: \u00E2\u0080\u00A2 Step 1. Primarily detect the jump occurrence by WMCM. \u00E2\u0080\u00A2 Step 2. Check the validity of the jump detection with PFAD. 3.7.1 Detection through WMCM Data from the accelerometers in three orthogonal axes of the body frame are captured in a window of (tje - t j s) = NT = 0.6 s. For each window the WMCM value mxyz(t j s+nT ) is calculated. Then for each window, the max- imum value of mxyz(t j s + nT ) is compared with a predefined threshold thm. If, for a particular window, the maximum value of mxyz(t j s + nT ) is greater than or equal to thm, then it is considered that a Jump Indication (JI) has been detected at the corresponding time epoch tJI = t j s+n j maxT . Here, n j max is the index of the maximum WMCM value for the jth window. To ensure that the maximum WMCM value in a potential locality is considered as JI, the maximum value of mxyz(t j+1 s +nT ), i.e. maximum WMCM value of the following window, is compared with the mxyz(tJI). If the maximum value of mxyz(t j+1 s + nT ) is greater than mxyz(tJI), the epoch of JI is updated to tJI = t j+1 s + n j+1 maxT . Otherwise the previous value for tJI is preserved. Thus, indication of jump occurrence, i.e. JI, is detected primarily through WMCM and tJI is passed to the next step for verification. WMCM Threshold It is important to note that detection of either jump start or jump end is adequate for successful detection of jump occurrence or JI. Depending on where the sensors are mounted, the relative magnitude of the WMCM value near the jump start and the jump end may vary. It has been empirically found that, for the Recon goggles where the sensors are mounted on the head, the WMCM value near the jump end is consistently greater than or equal to what is observed near the jump start. For example, the magnitude of the WMCM value near landing in Figure 3.10 is much higher than the magnitude observed near take-off. Therefore, to be more resistive to false JI detection and save power consumption, the threshold thm is set to a higher value that only secures the detection of JI near the jump end or landing. Hence, any JI near the jump start is not guaranteed to be detected by this threshold. However, if any JI is detected near the jump start, the corresponding JI at the jump end is overlooked. The sign of the PFAD value at tJI , dab(tJI), is used to determine whether the JI is near the jump start 40 3.7. Jump Detection or end. In this research, the threshold thm is set to a heuristically optimum value of 0.20 g3. For the data captured with the Ripxx device [4], which is strapped on to the leg of an athlete, the scenario is opposite. That is, the WMCM value near the jump start is consistently greater than or equal to the WMCM value near the jump end. This is because the different body parts experience different changes in dynamics for maneuvers. Consequently, the acceleration sensors capture different accelerations depending on where the sensors are mounted. Therefore, the threshold for the WMCM value should be set considering the change in dynamics of the sensor attached to the specific body part. 3.7.2 Validity Check through PFAD If any JI is detected by the WMCM procedure, the resultant acceleration ab is calculated from the three axes acceleration data captured. Then the PFAD value of ab at tJI , da btJI is calculated and compared with a prede- fined threshold thd. The selection procedure of thd is demonstrated later in Section 3.7.3. If the magnitude of the PFAD value, |dabtJI |, is greater than thd the detected JI is considered as a valid jump detection. Otherwise, the detected JI is discarded as a false jump indication. It is also possible to determine whether the JI is near take-off or landing from the sign of dabtJI as mentioned earlier. This knowledge is important for the AT determination algorithm proposed later in Section 3.8. Figure 3.13 summarizes the concept of the proposed jump detection al- gorithm. The plotted WMCM value mxyz is shifted by -2 g 3 for visual convenience. The high WMCM value will primarily report the JIs at the jump start and jump end. Both the JIs will be validated by the high value of PFAD at corresponding tJI . It should be noted that another JI will be reported in the first step due to the high WMCM value near 9.5 s. This detected JI will be discarded in the second step, as the PFAD value at the corresponding tJI is negligible in magnitude. Thus, the proposed jump de- tection method will only detect true jumps by filtering the false indication in a two step method. In Figure 3.13, LIpw and L I fw are selected to be 0.6 s. However, to efficiently address a wide range of jumps with different characteristics, the span of the PFAD windows can be varied. Depending on the acceleration variation characteristics at the neighboring epochs of the detected tJI , JI can be divided into two categories, namely high disturbance JI and low disturbance JI. Both these categories can also be subdivided into two more categories, high impact JI and low impact JI. Depending on the JI category, 41 3.7. Jump Detection 6 8 10 12 14 16 18 20 22 24 26 \u00E2\u0088\u00923 \u00E2\u0088\u00922 \u00E2\u0088\u00921 0 1 2 3 4 5 6 Time [s] ab [g ], da b [g ], m x yz [g 3 ] Resultant acceleration ab WMCM value m xyz PFAD value dab Take\u00E2\u0088\u0092off regionFalse Jump Indication Landing region Figure 3.13: Proposed jump detection method. Note: mxyz is shifted by -2 g3 for visual convenience. 42 3.7. Jump Detection the threshold value, thd, for the JI validation process is selected. 3.7.3 Jump Indication Categories JIs are categorized into high disturbance, low disturbance, high impact and low impact. The high/low disturbance categories set the window length for the PFAD calculation and the high/low impact categories determine the threshold for the PFAD value. High Disturbance JI. Generally for big jumps, where greater force is involved in take-off and stronger thrust is felt at landing, abrupt head and body motion cause high fluctuations in resultant acceleration near the jump start and end regions. In such cases, the detected JI is termed as high disturbance JI. Figure 3.14 shows an example of high disturbance JI. The trough in resultant acceleration following tJI is due to the abrupt head and body motions associated with hard landings. Jumps associated with this type of JI generally have lengthy aero phase. It is understandable from the figure that the PFAD value at tJI with narrow windows will not be significant in magnitude due to the trough. On the other hand, wide windows for PFAD calculations will lessen the effect of the trough and will produce a PFAD value with greater magnitude at tJI . Therefore, relatively higher values for LJIpw and L JI fw should be selected for high disturbance JIs. For this research, LJIpw and L JI fw are heuristically selected to be 1 s for high disturbance JIs. Low Disturbance JI. In many cases, athletes try to minimize the re- action force experienced by the body by conducting smooth take-offs and landings. In such cases, neighboring resultant accelerations of tJI show less fluctuations. Therefore, JI detected for this type of take-off or landing is termed as low disturbance JI. An example of low disturbance JI is shown in Figure 3.15. From the figure, it can be observed that the narrow window for averaging is adequate to extract the high PFAD value at tJI . It is impor- tant to note that the shallow trough is not deep enough to cause the PFAD value to be insignificant. In fact, wide windows may render the PFAD value to become less significant in magnitude. This is because low disturbance JIs are generally detected for short aero phase jumps. Hence, the span of wide windows around tJI may exceed the jump occurrence region, causing higher resultant acceleration (as shown in dotted circle in Figure 3.15) to be included in the averaging process. Therefore, a relatively lower value for 43 3.7. Jump Detection 18.5 19 19.5 20 20.5 0 0.5 1 1.5 2 2.5 Time [s] A m pl itu de o f r es ul ta nt a cc el er at io n ab [g ] Jump occurrence region Epoch of JI, tJI Narrow window for PFAD (0.15s) Deep trough Wide window for PFAD (0.95s) Figure 3.14: High disturbance JI. LJIpw and L JI fw is chosen for low disturbance JIs, which is 0.30 s in this study. Depending on the WMCM value at tJI the above mentioned JI categories can also be subdivided into two subcategories: High Impact JI. If the abrupt change in acceleration during take-off or landing occurs swiftly within a very brief time period, then the correspond- ing WMCM value tends to be higher in magnitude. This is because the mean cancelation process in WMCM causes the peak accelerations in the sensors to stand out from their neighbors. Hence the resultant multipli- cation at tJI , mxyz(tJI), exhibits high magnitude. Therefore, JIs with high WMCM magnitude are defined as high impact JIs. In this study, any JI with mxyz(tJI) \u00E2\u0089\u00A5 0.50g3 is considered a high impact JI. Figure 3.16 shows a high impact JI along with the resultant acceleration and corresponding WMCM value. Understandably, the JI is detected at the highest peak of WMCM at 40.5 s. However, it can also be observed, as previously discussed, that the PFAD value for narrow windows might not achieve a high magnitude due to the presence of a high disturbance near tJI . Similar scenarios are observed for general high impact JIs. Wide windows for PFAD will resolve 44 3.7. Jump Detection 73.5 74 74.5 75 75.5 76 \u00E2\u0088\u00921 0 1 2 3 4 5 6 Time [s] A m pl itu de o f r es ul ta nt a cc el er at io n ab [g ] Wide window for PFAD (1.15s)Narrow window for PFAD (0.20s) Shallow troughJump occurrence region Higher resultant acceleration near jump start Epoch of JI, tJI Figure 3.15: Low disturbance JI. this issue for lengthy aero phase jumps such as in Figure 3.16. Also, as pre- viously noted, wide windows will not work for jumps with short duration. Therefore, due to the unsteady nature of the resultant acceleration near JI, the threshold thd is set at a relatively low value for validity check through PFAD. Fortunately, this is also logically suited to the fact that any JI with higher WMCM magnitude is a more eligible candidate than one with a lower magnitude. Consequently, it is justified to have a less stringent threshold for the validation procedure of the high impact JI primarily detected by the WMCM procedure. Low Impact JI. Take-offs and landings are often observed where the change in resultant acceleration is steady and spans a prolonged period of time. Therefore, the acceleration at tJI cannot excel itself significantly from its neighbors. This phenomenon results in a relatively low value of WMCM at tJI . In this case, the JI with a relatively low WMCM value is termed as low impact JI. In this research, JIs with 0.20g3 \u00E2\u0089\u00A4 mxyz(tJI) < 0.50g3 are considered low impact JIs. As mentioned earlier, JIs with a WMCM value less than 0.20g3 are not considered as a valid JI. Figure 3.17 shows an example of a low impact JI. It can be discerned from the figure that the PFAD value is likely to have a high magnitude due to the steady nature of 45 3.7. Jump Detection 38.5 39 39.5 40 40.5 41 41.5 42 0 1 2 3 4 5 6 Time [s] ab [g ], m x yz [g 3 ] a b m xyz Jump occurrence region Epoch of JI, tJI Narrow window for PFAD (0.25s) Wide window for PFAD (1.45s) Figure 3.16: High impact JI. 46 3.7. Jump Detection 31.5 32 32.5 33 33.5 34 34.5 \u00E2\u0088\u00920.5 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 Time [s] ab [g ], m x yz [g 3 ] a b m xyz Epoch of JI, tJI Narrow window for PFAD (0.25s) Wide window for PFAD (1.20s) Jump occurrence region Figure 3.17: Low impact JI. the resultant acceleration on opposite sides of tJI for both wide and narrow windows. Therefore, for low impact JI the thd can be set to a higher value relative to the case of high impact JI. Fortunately, this also makes logical sense, as a JI with low WMCM value is more prone to be a false indication of jump than a JI with high WMCM magnitude. Hence, a rigorous false detection procedure with high threshold is rational for a low impact JI. PFAD Threshold Selection The selection of the PFAD threshold thd is done through a statistical analysis, namely a Two-Sample t-Test [47], for both high and low impact jumps. Acceleration data was collected with the Recon goggles for 27 ski jumps. Then JIs were detected through WMCM and corresponding a\u00CC\u0084bJIpw and a\u00CC\u0084bJIfw were calculated. For the purpose of this study, for jumps with aero phase duration well over 1 s, LJIpw and L JI fw are set at 1 s. For the others, the window lengths are set at 0.30 s. Depending on the WMCM value, the jumps are then divided into sets of high impact and low impact jumps. Since interest is focused on the difference in magnitude of a\u00CC\u0084bJIpw and a\u00CC\u0084 bJI fw , the two 47 3.7. Jump Detection sample t-Test is conducted with a\u00CC\u0084bJIlow as one sample set and a\u00CC\u0084 bJI high as another sample set. For any JI, a\u00CC\u0084bJIlow and a\u00CC\u0084 bJI high are the minimum and maximum value between a\u00CC\u0084bJIpw and a\u00CC\u0084 bJI fw respectively. In the two sample t-Test the test statistic is [47] t0 = y\u00CC\u00841 \u00E2\u0088\u0092 y\u00CC\u00842 Sp \u00E2\u0088\u009A 1 n1 + 1 n2 (3.7) where y\u00CC\u00841 and y\u00CC\u00842 are the means of the samples of a\u00CC\u0084 bJI high and a\u00CC\u0084 bJI low , and n1 and n2 are the sample sizes. It is assumed that a\u00CC\u0084 bJI high and a\u00CC\u0084 bJI low have equal variances. S2p is the estimate of the common variance computed from S2p = (n1 \u00E2\u0088\u0092 1)S21 + (n2 \u00E2\u0088\u0092 1)S22 n1 + n2 \u00E2\u0088\u0092 2 (3.8) and S21 and S 2 2 are two individual sample variances, i.e. variances of a\u00CC\u0084 bJI high and a\u00CC\u0084bJIlow respectively. The hypotheses being tested by the t-Test are \u00E2\u0080\u00A2 Null hypothesis H0: \u00C2\u00B51 = \u00C2\u00B52 \u00E2\u0080\u00A2 Alternative hypothesis H1: \u00C2\u00B51 > \u00C2\u00B52 where \u00C2\u00B51 and \u00C2\u00B52 are the means of a\u00CC\u0084 bJI high and a\u00CC\u0084 bJI low . The observation models are: a\u00CC\u0084bJIhigh = \u00C2\u00B51 + error and a\u00CC\u0084 bJI low = \u00C2\u00B52 + error. To determine whether to reject H0: \u00C2\u00B51 = \u00C2\u00B52, t0 is compared with the t distribution with n1 + n2 \u00E2\u0088\u0092 2 degrees of freedom. If t0 \u00E2\u0089\u00A5 t\u00CE\u00B6,n1+n2\u00E2\u0088\u00922, where t\u00CE\u00B6,n1+n2\u00E2\u0088\u00922 is the upper \u00CE\u00B6 percent point of the t distribution with n1 +n2\u00E2\u0088\u0092 2 degrees of freedom, H0 will be rejected. Here, \u00CE\u00B6 is the significance level, which represents the possibility of rejecting H0 while it is in fact true. For this test, \u00CE\u00B6 = 0.01 or 1%. Table 3.1 lists the two samples collected for the low impact JIs. For a 1% significance level the calculated test statistic for low impact jumps is t0 = 23.7975 > t.01,17+17\u00E2\u0088\u00922 = 2.4487. As expected, the null hypothesis is rejected and it is confirmed that \u00C2\u00B51 > \u00C2\u00B52. The extremely low P\u00E2\u0088\u0092value of 3.2184\u00C3\u009710\u00E2\u0088\u009222, where P\u00E2\u0088\u0092value is the smallest level of significance that would lead to the rejection of the null hypothesis [47], proves the strong rejection of the null hypothesis. However, the major goal of conducting this t-Test is to know the range within which (\u00C2\u00B51\u00E2\u0088\u0092 \u00C2\u00B52) is expected to lie for the low impact JIs. Therefore, interest is focused on the confidence interval of the test, which gives the range of the difference in means of the parameters (a\u00CC\u0084bJIhigh, 48 3.7. Jump Detection Table 3.1: Two sample t-Test for low impact JIs Jump No. a\u00CC\u0084bJIlow a\u00CC\u0084 bJI high 1 0.2939 1.8681 2 0.4915 1.8792 3 0.5255 1.5765 4 0.4384 1.8636 5 0.3136 1.6200 6 0.4369 1.7972 7 0.3926 1.4994 8 0.6478 1.2860 9 0.3205 1.6463 10 0.5037 1.8141 11 0.5072 1.6335 12 0.3574 1.6188 13 0.2552 1.5931 14 0.1568 1.6714 15 0.3652 1.6318 16 0.4095 1.6547 17 0.6945 2.0160 a\u00CC\u0084bJIlow). The interval 1.1377 \u00E2\u0089\u00A4 \u00C2\u00B51\u00E2\u0088\u0092\u00C2\u00B52 <\u00E2\u0088\u009E is the calculated 100(1\u00E2\u0088\u0092\u00CE\u00B6) = 99% confidence interval of the low impact JIs. In this regard, if a large number of (a\u00CC\u0084bJIhigh\u00E2\u0088\u0092 a\u00CC\u0084bJIlow), i.e. |dab(tJI)|, is gathered, 99% of them will be \u00E2\u0089\u00A5 1.1377. The calculated confidence interval provides a reasonable estimate of the value for the threshold thd for the low impact JI verification procedure. Since a value of |dab(tJI)| = 1.1377 is copiously large for an acceleration difference, a safer lower value can be selected for thd so that possibility of declaring any valid JI as invalid is minimized. Therefore, to validate numerous possible types of JIs, thd is selected to be 0.80 g for low impact jumps. A similar test is done for the samples of high impact JIs as shown in Table 3.2. Here, t0 = 6.2140 > t.01,10+10\u00E2\u0088\u00922 = 2.5524 and P\u00E2\u0088\u0092value is 3.6460\u00C3\u0097 10\u00E2\u0088\u00926. As previously discussed, the null hypothesis is also rejected here. The 99% confidence interval is 0.4531 \u00E2\u0089\u00A4 \u00C2\u00B51\u00E2\u0088\u0092\u00C2\u00B52 <\u00E2\u0088\u009E. Again, the lower boundary of the confidence interval gives an estimate of the PFAD threshold. Due to reasoning similar to that in the case of the low impact JIs, thd for high impact JIs is set to a relatively lower and safer value of 0.25 g for this project. 49 3.7. Jump Detection Table 3.2: Two sample t-Test for high impact JIs Jump No. a\u00CC\u0084bJIlow a\u00CC\u0084 bJI high 18 0.7547 2.2298 18 0.6225 1.2437 20 0.8521 1.3228 21 0.7997 1.2297 22 0.8587 1.3239 23 0.7796 1.6223 24 0.5607 1.6255 25 0.6478 1.6595 26 0.8489 1.0170 27 0.9048 1.9652 False JI Detection After any JI is detected by the WMCM procedure, the validity checking steps through PFAD are as follows: 1. From mxyz(tJI) it is determined whether the JI is high impact or low impact and the corresponding value for thd is set. 2. dab(tJI) is calculated considering L JI pw and L JI fw to be 1 s, i.e. as- suming high disturbance JI, because high disturbance JI requires long windows. 3. If dab(tJI) \u00E2\u0089\u00A5 thd the JI is considered as valid and the AT calculation algorithm is invoked. Otherwise the assumption that the JI is high disturbance is discarded and one proceeds to the next validation step. 4. After the assumption that JI is categorized as high disturbance is discarded, LJIpw and L JI fw are set to 0.30 s and da b(tJI) is calculated again assuming JI to be low disturbance. 5. If dab(tJI) \u00E2\u0089\u00A5 thd, the JI is considered as valid low disturbance JI and the AT calculation algorithm is used. Otherwise, the JI is completely discarded and considered to be a false detection. 6. After any false detection, the algorithm returns to the data acquisition mode and the WMCM value is calculated for the new data set. 50 3.8. Air Time Determination The assumption that the detected JI is high disturbance is made earlier than the low disturbance assumption to reduce computation. Most of the regular jumps will comply with the high disturbance criteria and skip the low disturbance JI verification computation. On the other hand, jumps with small AT and disturbance activities will only be discarded in the first assumption and the algorithm for low disturbance JI verification will be invoked. The flow-chart in Figure 3.18 summarizes the entire jump detection pro- cedure. 3.8 Air Time Determination 3.8.1 Characteristic Acceleration Points To determine the Air Time (AT) accurately, it is imperative to detect the precise epochs of the events \u00E2\u0080\u0098jump start\u00E2\u0080\u0099 or \u00E2\u0080\u0098take-off\u00E2\u0080\u0099 and \u00E2\u0080\u0098jump end\u00E2\u0080\u0099 or \u00E2\u0080\u0098landing\u00E2\u0080\u0099. \u00E2\u0080\u0098Take-off\u00E2\u0080\u0099 is defined as the event of momentary transition from ground to air. Similarly, \u00E2\u0080\u0098landing\u00E2\u0080\u0099 is the event of momentary transition from air to ground. In noiseless and undisturbed ideal conditions, the resultant acceleration pertaining to any jump should resemble a U-shaped pattern as shown in Figure 3.19. The characteristic points that define the U-shaped acceleration region are the highlighted points, namely jump start positive peak (JSP ), jump start negative peak (JSN ), jump end positive peak (JEP ), and jump end negative peak (JEN ). It should be noted that the local maxima of any cusp is termed as positive peak and the local minima of any notch is termed as negative peak. As shown in Figure 3.19, JSP is the latest point before aero phase when the reported acceleration indicates that the athlete is completely on the ground. JSN is the first point in aero phase when the reported acceleration suggests that the athlete is completely in the air. Therefore, the event \u00E2\u0080\u0098take-off\u00E2\u0080\u0099 must have occurred in the time span between these two points. Similarly, \u00E2\u0080\u0098landing\u00E2\u0080\u0099 must have occurred between JEN and JEP . Depending on the specific sport maneuvers, the weighted average of the epochs of JSP and JSN , and JEP and JEN should theoretically produce the exact epochs of take-off and landing. Consequently, the difference between these epochs should produce exact AT. Since the time spans for both take-off and landing are very brief in time (\u00E2\u0088\u00BC 0.05 s), normal averages of the start and end points should also yield AT with a reasonably high accuracy. The argument presented above describes why accurate AT determina- tion is directly dependent on the precise detection of the four characteristic 51 3.8. Air Time Determination Triaxial acceleration data Calculate WMCM max{ ( )} ?m t + nT thxyz s \u00E2\u0089\u00A5 j m No max{ ( )}> ) ?m t + nT txyz s JI j+1 mxyz( t = t + n TJI s max j j Yes t = t + n TJI s max j+1 j+1 Yes No th =m 0.20 g 3 thd = 0.25 g m txyz JI( ) 0.50 g ?\u00E2\u0089\u00A5 3 Yes No Calculate PFAD L JI pw fw= = 1 sL JI da t th b ( ) ?JI d\u00E2\u0089\u00A5 Yes No Calculate PFADL JI pw fw= = 0.30 sL JI da t th b ( ) ?JI \u00E2\u0089\u00A5 d Yes No Valid JI detected. Transmit for AT determinationtJI thd = 0.80 g Figure 3.18: Proposed jump detection method. 52 3.8. Air Time Determination Time R es u lt an t ac ce le ra ti o n Typical increase in acceleration before take-off Typical decrease in acceleration after landing 0 g U-shaped region JSP JSN JEP JEN Transition from ground to air Transition from air to ground Zero acceleration in aero phase Figure 3.19: U-shaped region and characteristic points for an ideal jump. points. Hence, the focus of the AT determination algorithm is to detect the four characteristic points as accurately as possible. Of course, for practical noisy jump acceleration signals the characteristic points are not as well de- fined as in the ideal case. The U-shaped region and the possible localities of JSP , JSN , JEP and JEN are shown in Figure 3.20 for a typical noisy data set. Due to noise and other disturbances, the transition periods as well as the characteristic points for the jump are not precisely obvious. Hence, in the proposed AT determination algorithm, positive or negative peaks pos- sessing the most eligibility for the U-shape defining points are detected as the characteristic points. 3.8.2 Primary and Secondary Detection Depending on the position of a valid JI, the detection techniques for the characteristic points differ. In this research, the detection of the two characteristic points near JI is termed Primary Detection and detection of the other two is termed Secondary Detection. The relevant peaks are consequently termed primary and secondary, positive/negative peaks. For example, if the JI is near take-off, then the detection of JSP and JSN is 53 3.8. Air Time Determination 11 11.5 12 12.5 13 13.5 14 14.5 15 \u00E2\u0088\u00921 0 1 2 3 4 5 Time [s] A m pl itu de o f r es ul ta nt a cc el er at io n ab [g ] U\u00E2\u0088\u0092shaped region Possible region of JSP Possible region of JSN Possible region of JEN Possible region of JEP Figure 3.20: U-shaped region and characteristic points for a practical jump. 54 3.8. Air Time Determination primary detection and the detection of JEP and JEN is secondary detec- tion. In this case, primary positive and negative peaks are JSP and JSN respectively, and secondary positive and negative peaks are JEP and JEN respectively. Since the challenges in primary and secondary detection are dissimilar, distinct methods are proposed for each. For the sake of convenience, methods for primary and secondary de- tection are demonstrated for a particular example. As the features of the jump start characteristic points are the reciprocal of those of the jump end characteristic points, the demonstration will not lose any generality. 3.8.3 Window Selection After the detection and validation of any JI with WMCM and PFAD, a characteristic window is selected within which the AT determination algo- rithm is applied to detect all four characteristic points. This window spans from tcs to tce where tcs = tJI \u00E2\u0088\u0092\u00E2\u0088\u0087tlcs, l = 1, 2, .... (3.9) tce = tJI +\u00E2\u0088\u0087tlce, l = 1, 2, ..... (3.10) Primarily, for any JI near the jump start, \u00E2\u0088\u0087t1cs = 0.50 s and \u00E2\u0088\u0087t1ce = 3 s. For a JI near the jump end these parameter values are reciprocated, i.e. \u00E2\u0088\u0087t1cs = 3 s and \u00E2\u0088\u0087t1ce = 0.50 s. Since the jump aero phase can vary in length, the primarily selected window size may not prove to be adequate. On the other hand, selecting a very large window each time to accommodate the largest possible jump will excessively increase computation time and power. Therefore, window length is extended only if any of the characteristic points are not detected in the current window span. Each time, the window expansion parameter is increased following \u00E2\u0088\u0087tl+1ce = \u00E2\u0088\u0087tlce + 3 s for JI near jump start and \u00E2\u0088\u0087tl+1cs = \u00E2\u0088\u0087tlcs+3 s for JI near jump end. The window is expanded until tlce\u00E2\u0088\u0092tlcs > 9 s. If any of the characteristic points are not detected within the maximum span of the window (9 s), the JI is discarded as invalid. 3.8.4 Positive and Negative Peaks After selecting the characteristic window, the epochs of the positive and negative peaks within the window are determined. The positive peak epochs 55 3.8. Air Time Determination 12 12.5 13 13.5 14 14.5 15 \u00E2\u0088\u00920.5 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 Time [s] A m pl itu de o f r es ul ta nt a cc el er at io n ab [g ] Resultant acceleration Positive peaks Epoch of JI, tJI \u00E2\u0088\u0087 t1 ce \u00E2\u0088\u0087 t1 cs Figure 3.21: Window for characteristic acceleration points detection and potential positive peaks. can be mathematically represented as tPc = arg tcs+nT { (ab(tcs + nT )\u00E2\u0088\u0092 ab(tcs + nT \u00E2\u0088\u0092 1)) > 0 } \u00E2\u0088\u00A9 (3.11) arg tcs+nT { (ab(tcs + nT )\u00E2\u0088\u0092 ab(tcs + nT + 1)) \u00E2\u0089\u00A5 0 } (3.12) where arg tcs+nT {\u00C2\u00B7} is the set of epochs, tcs + nT , for which the argument {\u00C2\u00B7} is valid. The operator for intersection of the set of epochs satisfying pertinent arguments is \u00E2\u0088\u00A9. Similarly the epochs for negative peaks are tNc = arg tcs+nT { (ab(tcs + nT )\u00E2\u0088\u0092 ab(tcs + nT \u00E2\u0088\u0092 1)) < 0 } \u00E2\u0088\u00A9 (3.13) arg tcs+nT { (ab(tcs + nT )\u00E2\u0088\u0092 ab(tcs + nT + 1)) \u00E2\u0089\u00A4 0 } . (3.14) For the jump shown in Figure 3.20 the characteristic window and corre- sponding positive peaks are shown in Figure 3.21. The JI is detected at tJI = 12.09 s and is validated afterwards by the jump detection algorithm. This example will be used for illustrating the primary and secondary detection procedures. 56 3.8. Air Time Determination 3.8.5 Primary Detection From the PFAD value of the jump detection algorithm it is already known that the JI is near the jump start. Therefore, primary detection means the detection of JSP and JSN for this JI and subsequent detection of the pertinent epochs tPJS and t N JS . For the Recon goggles, primary detection is done by the straight-forward Closest Peak Method. Closest Peak Method In this method, the JI is exploited as a reference point. Within the characteristic window, the positive peak in the resultant acceleration nearest JI with a magnitude greater than or equal to 1 g is considered as JSP . The nearest negative peak with a resultant acceleration less than 1 g is detected as JSN . The time epochs for JSP and JSN can be expressed as tPJS = min tI { |tJI \u00E2\u0088\u0092 arg tI { ab(tI) \u00E2\u0089\u00A5 1g }|}, tI \u00E2\u0088\u0088 tPc (3.15) tNJS = min tI { |tJI \u00E2\u0088\u0092 arg tI { ab(tI) < 1g }|}, tI \u00E2\u0088\u0088 tNc (3.16) where min tI {\u00C2\u00B7} yields the epoch tI which yields the minimum value of the argument {\u00C2\u00B7}. Figure 3.22 shows the detected JSP and JSN . The detected peaks are at tPJS = 12.09 s and t N JS = 12.24 s. It can be observed from the figure that both the detected peaks reasonably estimate the two characteristic points on one side of the U-shaped region. The threshold of 1 g is selected theoretically since the ideal minimum resultant acceleration of the athlete is 1 g on the ground. 3.8.6 Secondary Detection Each potential secondary peak possesses four characteristic attributes or criteria on which its eligibility can be assessed. The attributes for a peak at tI are 1. Proximity to JI, tI \u00E2\u0088\u0092 tJI 2. PFAD value, dab(tI) 3. Average magnitude of the preceding or following window, a\u00CC\u0084bIpw or a\u00CC\u0084 bI fw 4. Self amplitude (resultant acceleration), ab(tI). 57 3.8. Air Time Determination 12 12.5 13 13.5 14 14.5 15 \u00E2\u0088\u00920.5 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 Time [s] A m pl itu de o f r es ul ta nt a cc el er at io n ab [g ] Epoch of JI, tJI Theoretical threshold 1 g JSP JSN Figure 3.22: Primary detection. For each candidate peak, these four criteria values are calculated. To calcu- late the preceding and following window average a\u00CC\u0084bIpw and a\u00CC\u0084 bI fw, the window spanning from tcs to tI and from tI to tce is considered respectively. Two methods are proposed to assess the criteria/attributes and determine the most eligible positive and negative peak. The first method is a threshold based method and the second method is a probabilistic method, which is independent of any threshold. Both of these methods are demonstrated in the next section. Similar to primary detection, JI serves as a reference for the secondary detection. For the JI of the given example, secondary detection means the detection of JEP and JEN and their epochs t P JE and t N JE . Figure 3.23a and 3.23b show the potential positive peaks for JEP and the corresponding attribute values for each candidate. Note that for the detection of JEP , all the positive peaks following the JI are the only viable candidates. A similar argument is valid for the negative peak candidates for JEN . Threshold Based Method As shown in Figure 3.23b, the attribute values for the candidate peaks for JEP vary significantly. It can be understood that a potential candidate for JEP at tI should have a negative da b(tI) with magnitude as high as possible, 58 3.8. Air Time Determination 12 12.5 13 13.5 14 14.5 15 \u00E2\u0088\u00920.5 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 Time [s] A m pl itu de o f r es ul ta nt a cc el er at io n ab [g ] Resultant acceleration Candidate peaks for JEP Epoch of JI, tJI Discarded region in secondary detection (a) 12 12.5 13 13.5 14 14.5 15 \u00E2\u0088\u00921 0 1 2 3 4 5 Time [s] A ttr ib ut e va lu e (tI \u00E2\u0088\u0092 tJI) [s] dab(tI) [g] a bI pw [g] a b(tI) [g] Epoch of JI, tJI (b) Figure 3.23: (a) Potential positive peaks for JEP . (b) Attribute value of the potential positive peaks. 59 3.8. Air Time Determination a\u00CC\u0084bIpw should be as low as possible, and a b(tI) should have a moderately high value. Therefore, three thresholds are set for these criteria. Among the set of candidate positive peaks that satisfies these thresholds, the peak nearest to JI is selected as JEP . On the account of a reciprocal rationale, the farthest peak from JI is selected as JEN from the set of the candidate negative peaks satisfying all the thresholds. The thresholds for the attribute values of the candidates for JEP and JEN at tI are given in Table 3.3. Furthermore, Table 3.4 shows the thresholds if JSP and JSN are to be detected as secondary detection. Table 3.3: Thresholds for secondary detection of JEP and JEN Attribute Threshold for JEP Threshold for JEN tI \u00E2\u0088\u0092 tJI Lowest Highest dab(tI) < \u00E2\u0088\u00920.25 g < \u00E2\u0088\u00920.25 g a\u00CC\u0084bIpw < 1 g < 1 g ab(tI) > 1.50 g < 0.50 g Table 3.4: Thresholds for secondary detection of JSP and JSN Attribute Threshold for JSP Threshold for JSN tJI \u00E2\u0088\u0092 tI Lowest Highest dab(tI) > 0.25 g > 0.25 g a\u00CC\u0084bIfw < 1 g < 1 g ab(tI) > 1.50 g < 0.50 g Using the predefined thresholds, detection of the epochs of JEP and JEN , i.e. t P JE and t N JE , can mathematically be represented as tPJE = min tI {{ arg tI { dab(tI) < \u00E2\u0088\u00920.25g } \u00E2\u0088\u00A9 arg tI { a\u00CC\u0084bIpw < 1g }\u00E2\u0088\u00A9 arg tI { ab(tI) > 1.50 }}\u00E2\u0088\u0092 tJI}, tI \u00E2\u0088\u0088 tPc (3.17) tNJE = max tI {{ arg tI { dab(tI) < \u00E2\u0088\u00920.25g } \u00E2\u0088\u00A9 arg tI { a\u00CC\u0084bIpw < 1g }\u00E2\u0088\u00A9 arg tI { ab(tI) < 0.50 }}\u00E2\u0088\u0092 tJI}, tI \u00E2\u0088\u0088 tNc . (3.18) 60 3.8. Air Time Determination 12 12.5 13 13.5 14 14.5 15 \u00E2\u0088\u00920.5 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 Time [s] A m pl itu de o f r es ul ta nt a cc el er at io n ab [g ] Epoch of JI, tJI JEP JEN Figure 3.24: Secondary detection of JEP and JEN with threshold based method. Similar expressions can be derived for tPJS and t N JS for the threshold based method. Figure 3.24 shows the detected JEP at t P JE = 13.78 s and JEN at tNJE = 13.76 s. It is evident from the figure that the detected peaks make reasonable sense as two of the characteristic points of the U-shaped region. Probabilistic Method The threshold based method is an effective and simple technique for sec- ondary detection. However, the loss of generality cannot be avoided as is true for all threshold based techniques. Use of several thresholds makes the algorithm susceptible to error, especially for jumps with various aerial acro- batic maneuvers. As previously mentioned, a potential candidate for JEP at tI should have the highest da b(tI) magnitude and the lowest a\u00CC\u0084 bI pw. Un- fortunately, a situation may occur whereby the candidate peaks possessing the highest dab(tI) and the lowest a\u00CC\u0084 bI pw are not the same. In such a case, the most eligible peak should be selected through a decision making procedure. To accommodate all the above mentioned issues, a probabilistic method using MADM [48] is proposed. This method uses no thresholds and selects the best option among the alternatives, i.e. the candidate peaks. The idea of secondary detection by probabilistic means stems from the fact that it is possible to make an educated conjecture on the position of the characteristic 61 3.8. Air Time Determination points only from the shape of the acceleration data. For example, without even looking at the y-axis magnitude, an educated eye should be able to de- tect the possible regions of the characteristic points as shown in Figure 3.20. Rather than directly discarding any point, the probabilistic method assesses and scores the possibility of all candidate peaks to be a characteristic point. Depending on the score, the optimum solution is selected. The most de- sirable advantage of this method is that it renders the characteristic point detection process to be very generalized and to be able to cope with a large variety of situations. A decision making problem can be formulated for the secondary detec- tion process considering the candidate peaks as alternatives to choose from. The decision making can be conducted using the four attribute values of each alternative. M-TOPSIS [49] is a modified version of the popular TOP- SIS [48] method for MADM. This method effectively addresses the \u00E2\u0080\u0098Rank Inversion\u00E2\u0080\u0099 [49] problem of the TOPSIS method. Both the TOPSIS and M- TOPSIS methods, along with the rank inversion problem, are demonstrated in Appendix A. Due to the advantage over TOPSIS, M-TOPSIS is selected as the decision making technique for secondary detection. To detect JEP , the positive peaks shown in Figure 3.23a are fed to the MADM method as alternatives. The four attribute values for each alterna- tive, as shown in Figure 3.23b, are used to construct the decision matrix. Preference information about the attributes, i.e. weights wci , i = 1, 2, 3, 4, are predefined. The information regarding whether the attributes are cost type or benefit type is also provided to the MADM algorithm. Benefit at- tributes are those which are expected to be maximized and cost attributes are those which are expected to be minimized. Depending on the specific characteristic point, the attributes swap between benefit type and cost type. Table 3.5 and Table 3.6 lists all the attribute types and weights for secondary detection of JEP , JEN , JSP and JSN . Using the attribute type and weight information, the decision matrix is normalized and the weighted normalized ratings are calculated. The distances from the alternatives to the positive ideal solution (D\u00E2\u0088\u0097j ) and to the negative ideal solution (D \u00E2\u0088\u0092 j ) are measured af- terwards. Subsequently, the D\u00E2\u0088\u0097D\u00E2\u0088\u0092 plane is constructed and the Optimized Ideal Reference Point (OIRP) is determined. Figure 3.25 shows the D\u00E2\u0088\u0097jD \u00E2\u0088\u0092 j plane for JEP detection. From this plane the distances between the OIRP and the alternatives (DOIRPj ) are calcu- lated. The color intensity of the alternatives shown in the figure is propor- tional to its distance from the OIRP. DOIRPj is sorted in ascending order and all the alternatives, i.e. candidate peaks, are ranked against it correspond- 62 3.8. Air Time Determination ingly. The peak possessing the lowest DOIRPj is ranked 1 and considered as the optimal solution, i.e. JEP , for the MADM problem. Figure 3.26 shows the top ranked peaks for secondary JEP detection. Again, the color intensity of the peaks is proportional to DOIRPj . The epoch of the ranked 1 peak is determined as tPJE . In a similar manner, JEN is detected and tNJE is determined. The result of secondary detection with the probabilistic method is shown in Figure 3.27. The detected tPJE = 13.80 s and t N JE = 13.76 s. Table 3.5: Attribute type and weight for probabilistic secondary detection of JEP and JEN For JEP For JEN Attribute Type Weight Attribute Type Weight tI \u00E2\u0088\u0092 tJI Cost 0.05 tI \u00E2\u0088\u0092 tJI Benefit 0.08 |dab(tI)| Benefit 0.50 |dab(tI)| Benefit 0.50 a\u00CC\u0084bIpw Cost 0.33 a\u00CC\u0084 bI fw Benefit 0.38 ab(tI) Benefit 0.12 a b(tI) Cost 0.04 Table 3.6: Attribute type and weight for probabilistic secondary detection of JSP and JSN For JSP For JSN Attribute Type Weight Attribute Type Weight tJI \u00E2\u0088\u0092 tI Cost 0.05 tJI \u00E2\u0088\u0092 tI Benefit 0.08 |dab(tI)| Benefit 0.50 |dab(tI)| Benefit 0.50 a\u00CC\u0084bIfw Cost 0.33 a\u00CC\u0084 bI pw Benefit 0.38 ab(tI) Benefit 0.12 a b(tI) Cost 0.04 A very safe threshold can be useful from the computational point of view for the implementation of this MADM based probabilistic method even though it does not require a threshold. For example, a safe threshold of > 1 g for secondary JEP or JSP detection will eliminate a significant number of candidate positive peaks from being considered as alternatives for the decision making process. The reduced number of alternatives will cause reduced computation resulting in efficient use of resources. 63 3.8. Air Time Determination 0 0.05 0.1 0.15 0.2 0.25 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 0.22 0.24 Distance from the positive ideal solution Dj * D ist an ce fr om th e ne ga tiv e id ea l s ol ut io n D j \u00E2\u0088\u0092 0.0004 0.023 0.053 0.064 0.1 0.14 Dj * Dj \u00E2\u0088\u0092 plane Optimized ideal reference point Alternatives/Candidate peaks Distances from the reference point Figure 3.25: D\u00E2\u0088\u0097D\u00E2\u0088\u0092 plane for the secondary detection of JEP . 12 12.5 13 13.5 14 14.5 15 0 0.5 1 1.5 2 2.5 3 3.5 4 Time [s] A m pl itu de o f r es ul ta nt a cc el er at io n ab [g ] 0.0004 0.023 0.053 0.064 0.1 0.14 Rank 4 Rank 1 Rank 2 Rank 3 Figure 3.26: DOIRPj and rank of the top candidate peaks for JEP . 64 3.9. Sensor and Visual Air Time 12 12.5 13 13.5 14 14.5 15 \u00E2\u0088\u00920.5 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 Time(s) A m pl itu de o f r es ul ta nt a cc el er at io n ab [g ] Epoch of JI, tJI JEP JEN Figure 3.27: Secondary detection of JEP and JEN with probabilistic method. 3.9 Sensor and Visual Air Time For snowboarding or skiing, Recon defines the jump start epoch as the moment when the skis/board are no longer in contact with the ground. Recon also defines the jump end epoch as the moment when any part of the skis/board come into contact with the ground again. The difference between the end and the start epochs defined in this manner leads to an AT which can be visually realized. The procedure for calculating AT from the JSP , JSN , JEP and JEN is described later in Subsection 3.10.2. However, up to this point, the jump start and end epochs have been discussed from the point of view of the sensor signals. The visual jump start and end epochs are not necessarily same as what the sensors detect as the jump start and end epochs. This is because the changes in the athlete\u00E2\u0080\u0099s dynamics as detected by the sensors, does not always reflect what can be visually observed. For example, while landing, the athlete may have been observed to touch the ground even though no real impact was sensed by the sensors. Similarly, during take-off, a change in acceleration is typically sensed by the sensors prior to what is defined as the take-off by Recon. Logically, it is expected to have a constant bias in the visually realized AT compared to the AT determined from the sensor signal, which is termed as Sensor Air 65 3.10. Experimental Results and Comparison Time (SAT). From the experimental data, a simple relationship is developed to convert the determined SAT into visual AT. For the rest of the discussion, AT will refer to what is visually observed unless otherwise stated. 3.10 Experimental Results and Comparison 3.10.1 Field Data Collection Field data was collected to test the performance of the proposed jump detection and AT determination algorithm. The athlete had a Ripxx [4] unit strapped to his leg and the Recon goggles on his head. Both the devices were set to capture and log raw data from the available sensors. A video camera with 30 fps (frame per second) capture capability was set to record the jumps. Snowboard and ski jumps with different aero phase duration and including various aerial maneuvers, e.g. spinning, were conducted. The frame numbers of the skier leaving the ground and returning to the ground were determined from video playback. Since the frame rate is known, (dif- ference in frame numbers / frame rate in fps) provides the AT in seconds as determined from the video. This video AT serves as the benchmark for comparing the performance of current Recon, Ripxx and proposed algo- rithms. However, the video AT is only accurate within a maximum error of (0.033 s + 0.033 s) = 0.067 s. The error is twice (1/frame rate) since a single frame error may occur both at the jump start and end. 3.10.2 Visual AT from SAT Data captured with the Recon goggles is fed to the proposed algorithm to detect the jump and determine the corresponding jump characteristic points. For secondary detection, the probabilistic method is used. The difference between tNJS and t N JE is used to determine the SAT, since JSN/JEN falls closest to the take-off/landing points defined by Recon. An example of the bias between visual AT and SAT is illustrated in Figure 3.28. The athlete\u00E2\u0080\u0099s position and corresponding resultant acceleration are linked in the figure. The characteristic points detected by the proposed algorithm are also shown. Image 1 (frame 3409) is the first epoch when the accelerometers sense free fall. However, it is clear from the image that the snowboard is still in contact with the ground. Visually, the take-off is detected in image 2 (frame 3411). Hence, sensor take-off is at an earlier time than visual take-off. Similarly, images 4 and 5 exhibit the bias in sensor and visual landing. Therefore, it is evident that the visual AT would always be shorter than SAT. 66 3.10. Experimental Results and Comparison 73.6 73.8 74 74.2 74.4 74.6 74.8 75 75.2 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 A m p li tu d e o f re su lt an t ac ce le ra ti o n a [g ] b Time [s] JS P JS N JE P JE N 1 2 4 Frame 3409 Frame 3411 Frame 3430 3 5 Frame 3432 Frame 3422 Figure 3.28: Experimental video and captured sensor resultant accelera- tion. 67 3.10. Experimental Results and Comparison 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 \u00E2\u0088\u00920.5 0 0.5 1 1.5 2 Sensor AT (SAT) [s] V id eo A T [s] AT = 0.9438 \u00C3\u0097 SAT \u00E2\u0088\u0092 0.0138 Figure 3.29: Video AT vs SAT. To determine the AT from SAT, the determined SATs are plotted against the ATs from video. The plot is shown in Figure 3.29. A first order linear curve is fitted to develop a relation between the visual AT and SAT. The developed relation, AT = 0.9438\u00C3\u0097SAT\u00E2\u0088\u00920.0138 s, can be used to determine the visual AT from determined SAT. Of course, the relationship could be refined with acquisition of more jump data for any specific sport application. Furthermore, the interpretation or use of the characteristic points along with the SAT will vary with different definitions of take-off/landing. The devel- oped relation between visual AT and SAT is most suitable for snowboard and ski jumps according to Recon\u00E2\u0080\u0099s definition. 3.10.3 Performance Comparison The algorithm for jump detection and AT detection are compared in different aspects. The most important aspect is the accuracy of the AT. Similar to the proposed algorithm, the Recon jump detection and AT deter- mination algorithm, which is a threshold based method, was implemented in Matlab c\u00C2\u00A9 for off-line processing. The Recon algorithm is proprietary and, therefore, cannot be described in detail here. The recorded accelerometer data sets collected with the Recon goggles are fed to both the Recon algo- rithm and the proposed algorithm to generate AT. ATs for the proposed 68 3.10. Experimental Results and Comparison algorithm are determined by converting the SATs. AT for the test jumps are also available from the Ripxx log sheet. These three ATs are compared with the ATs deduced from the video. The error percentages with respect to the true AT (video) are also calculated. The results are shown in Table 3.7. Ripxx data was not available for the first two data sets. Note that an error of 100% indicates that a jump that in fact occurred, was not detected. From Table 3.7 it can be seen that the proposed algorithm is more effi- cient in jump detection and has better accuracy in AT determination than the other two. Among the 25 jumps, the proposed algorithm successfully detected 23 of them (92% accuracy). On the other hand, the Recon algo- rithm detected only 15 jumps out of the 25 (60% accuracy). Among the available 23 jumps, the Ripxx algorithm was only able to successfully detect 10 of them (43.5% accuracy). Moreover, the Recon and Ripxx algorithms have shown an AT determination average error of 8.4% and 39.7% respec- tively where as the proposed algorithm has an average error of 4.8%. It is very important to note that in the average error calculations the undetected jumps for each algorithm were excluded. Therefore, the effect of undetected jumps does not reflect on the average error. Figure 3.30 shows, for all jumps, the measured AT error in seconds for all three algorithms. Apart from very few exceptions, error for the proposed jump detection and AT determination algorithm is well within the maximum allowable error limit of \u00C2\u00B10.1 s, whereas the other two exhibit considerable error and often exceed the specification of \u00C2\u00B10.1 s. The second important aspect when comparing algorithms is their use of resources. Depending on the jump detection method the number of false wake up calls to the processor varies from algorithm to algorithm. A false wake up means the realization of a false jump detection after most of the computation is conducted. For example, the proposed algorithm will cause a false wake up call if any of the characteristic points are not detected for the largest characteristic window. No information on Recon false wake up criteria can be provided in this study due to proprietary issues. For a given data set, if the number of false wake up calls is significantly large, more power will be unnecessarily drained, leading to a lower battery life. Therefore, the number of false wake up calls to the processor is also determined and used as an algorithm comparison criterion. The comparison is shown in Table 3.8. Since the Ripxx algorithm for AT detection is unknown, the comparison for Ripxx on this criterion could not be done. From Table 3.8 it is clear that the proposed algorithm generates far fewer false wake up calls as compared to the Recon algorithm. The Recon algorithm generates a total of 21 false wake up calls while operating on the 14 data sets containing 25 jumps in 69 3.10. Experimental Results and Comparison 0 5 10 15 20 25 \u00E2\u0088\u00921 \u00E2\u0088\u00920.5 0 0.5 1 1.5 Jump Er ro r i n ai r t im e [s] Ripxx Recon Proposed 0.1 s \u00E2\u0088\u00920.1 s Figure 3.30: Error in AT comparison among the Recon, Ripxx and pro- posed algorithm. total. Therefore, the current Recon jump detection algorithm unnecessarily wakes the processor 21 \u00E2\u0088\u0097 100/(21 + 25) = 45.6% of the time. The 11th data set contains 12 jumps which is far more than any of the other data sets. This is the reason for the large number of false wake up calls for this set for the Recon algorithm. On the other hand, the proposed algorithm generated no false wake up calls for any of the 14 data sets. The third criterion for algorithm performance comparison is the number of false jumps detected and true jumps undetected. As shown in Table 3.7, any true jump undetected also results in a 100% error in the corresponding AT calculation. The results are shown in Table 3.9 for the collected data set. In the 14 data sets (containing 25 jumps), the Recon algorithm had 11 wrong detections (44%), i.e. true jumps undetected plus false jump detected. For Ripxx, the number of wrong detections is 13 (52%). On the other hand, the proposed algorithm has shown significantly better performance by only reporting 2 wrong detections, (8%) while operating over the 14 data sets. From these results, it is evident that the proposed method has com- prehensive supremacy in AT determination performance over the other two algorithms. The false wake up calls are also reduced significantly by the proposed method. Except for the 11th data set, the proposed method also 70 3.11. Implementation Aspects for Other Devices/Applications exhibits efficiency in avoiding false detection or missing true jumps. More- over, the two missed jumps in the 11th set had extremely brief AT duration and are likely to be insignificant to the goggle user. Therefore, it can be con- cluded from the available experimental results that the proposed algorithm performs significantly better than the current algorithm Recon is using and the algorithm of the commercially available Ripxx unit. 3.11 Implementation Aspects for Other Devices/Applications The proposed methods for primary and secondary detection can be al- tered according to the sensor signal characteristics. As mentioned earlier, the jump detection thresholds are specifically customized for Recon goggles only to guarantee the detection of landing. This is done to minimize the mandatory computation as well as the possibility of false detection. There- fore, the closest peak method for primary detection and the probabilistic method for secondary detection were proposed. However, if the scenario is such that both the JI near jump start and end are guaranteed to be de- tected, the closest peak method can be used for both primary and secondary detection. On the other hand, if scenarios are observed where the closest peak method is not as accurate as expected, the probabilistic method can be used for both primary and secondary detection. In such a case, after secondary detection, any of the detected secondary characteristic points can be selected as reference points to serve as the JI for the primary detection process. 71 3.11. Implementation Aspects for Other Devices/Applications Table 3.7: Comparison of AT among Recon, Ripxx and the proposed algo- rithm. Data set Video AT(s) Recon AT(s) [%error] Ripxx AT(s) [%error] Proposed AT(s) [%error] 1 1.4333 1.28 [10.6] - [ - ] 1.4585 [1.8] 2 1.4333 1.43 [0.2] - [ - ] 1.4207 [0.9] 3 1.3333 1.10 [17.5] 1.10 [17.5] 1.3735 [3.0] 4 1.2000 1.28 [6.7] 1.35 [12.5] 1.2131 [1.1] 5 1.2333 1.36 [10.3] 1.05 [14.9] 1.1281 [8.5] 6 1.3333 1.39 [4.2] 0.60 [55.0] 1.3452 [0.9] 7 1.3333 1.34 [0.5] 0.62 [53.5] 1.3263 [0.5] 8 1.3333 0.00 [100.0] 0.00 [100.0] 1.3452 [0.9] 9 1.2667 0.72 [43.1] 0.00 [100.0] 1.1848 [6.5] 10 1.4000 1.30 [7.1] 0.94 [32.8] 1.4207 [1.5] 11-a 1.3333 1.39 [4.2] 0.65 [51.2] 1.3169 [1.2] 11-b 1.5000 1.35 [10.0] 0.00 [100.0] 1.4962 [0.2] 11-c 0.6667 0.00 [100.0] 0.00 [100.0] 0.7034 [5.5] 11-d 0.6667 0.67 [0.5] 0.00 [100.0] 0.7223 [8.3] 11-e 0.4333 0.00 [100.0] 0.00 [100.0] 0.4203 [3.0] 11-f 0.1333 0.00 [100.0] 0.00 [100.0] 0.00 [100.0] 11-g 0.4333 0.00 [100.0] 0.00 [100.0] 0.3165 [26.9] 11-h 0.3000 0.00 [100.0] 0.00 [100.0] 0.3448 [14.9] 11-i 0.1667 0.00 [100.0] 0.00 [100.0] 0.0000 [100.0] 11-j 0.3333 0.00 [100.0] 0.00 [100.0] 0.3637 [9.1] 11-k 0.5333 0.00 [100.0] 0.00 [100.0] 0.5619 [5.4] 11-l 0.6000 0.00 [100.0] 0.00 [100.0] 0.6374 [6.2] 12 1.3333 1.30 [2.5] 0.55 [58.7] 1.3263 [0.5] 13 1.7000 1.63 [4.1] 0.96 [43.5] 1.7227 [1.3] 14 1.2300 1.28 [4.1] 1.93 [56.9] 1.2131 [1.4] Avg. error(s) [%error] 0.111 [8.4] 0.538 [39.7] 0.033 [4.8] 72 3.11. Implementation Aspects for Other Devices/Applications Table 3.8: Comparison of number of false processor wake up calls for Re- con algorithm and the proposed algorithm. Data set Recon Proposed 1 1 0 2 0 0 3 1 0 4 0 0 5 1 0 6 0 0 7 0 0 8 2 0 9 0 0 10 1 0 11 9 0 12 3 0 13 0 0 14 3 0 Total 21 0 73 3.11. Implementation Aspects for Other Devices/Applications Table 3.9: Comparison of the number of false detected jumps + true unde- tected jumps among Recon, Ripxx and the proposed algorithm. Data set Recon Ripxx Proposed 1 0 - 0 2 0 - 0 3 0 1 0 4 0 0 0 5 0 0 0 6 0 0 0 7 0 0 0 8 1 1 0 9 0 1 0 10 0 0 0 11 9 10 2 12 0 0 0 13 0 0 0 14 1 0 0 Total 11 13 2 74 Chapter 4 Trajectory Determination The determination of jump horizontal distance, height and drop corre- sponds to the task of determining the trajectory of the athlete during a jump. Once the trajectory is determined, the calculation of the desired parameters is straight forward. In this chapter a GPS/INS integrated algorithm is pro- posed for trajectory determination with the Recon goggles. First, the coor- dinate frames involved in the algorithm are presented. Afterwards, the basic INS algorithm and pertinent important factors are demonstrated. Finally, the proposed algorithm with detailed performance analysis is presented. 4.1 Coordinate Frames The local level frame has been chosen as the navigation frame (index n). The local level frame is defined as the local North-East-Down (NED) axes. The north axis is labeled X, the east is labeled Y , and the down is labeled Z. As the reference frame for positioning, the Earth Centered Earth Fixed (ECEF) frame is used, which is indexed by e. The ECEF frame is defined as [50]: \u00E2\u0080\u00A2 origin - at the mass center of the earth \u00E2\u0080\u00A2 X axis - pointing towards the Greenwich meridian, in the equatorial plane \u00E2\u0080\u00A2 Y axis - 90\u00E2\u0097\u00A6 east of Greenwich meridian, in the equatorial plane \u00E2\u0080\u00A2 Z axis - axis of rotation of the reference ellipsoid. Measurement from the gyroscope is done with respect to an inertial frame which is indexed by i. The inertial frame is considered to be rotationally fixed. Coordinates in the ECEF frame can be transformed to the inertial frame by a negative rotation about the Z axis of the frame by the amount of the Greenwich Mean Sidereal Time (GMST) [50]. The navigation frame, ECEF frame and inertial frame are depicted in Figure 4.1. The body frame is indexed by b and defined as shown in Figure 1.4. The triaxial sensors are 75 4.2. Rotation Matrix Greenwich meridian Z e , Z i Y e X e Northing Easting Down Z n X n Y n Earth Rotation Axis \u00CF\u0089ie X i Y i 360 - GMST Equator Equinox Figure 4.1: Coordinate frames. aligned entirely with the body frame and the body frame is considered to be rigidly attached to the athlete. It is important to note that all of these frames follow right handedness. 4.2 Rotation Matrix The most convenient way to represent the attitude of the body frame with respect to the navigation frame is through a set of three Euler an- gles [51]. These angles are commonly known as roll (\u00CF\u0086), pitch (\u00CE\u00B8) and yaw (\u00CF\u0088) which represent rotation along the X, Y and Z axes accordingly. Any arbitrarily rotated body frame attitude with respect to the navigation frame can be represented as three consecutive coordinate rotations. The three co- ordinate rotations along the X, Y and Z axes can be represented in matrix 76 4.2. Rotation Matrix form respectively as follows: RX(\u00CF\u0086) = \u00EF\u00A3\u00AE\u00EF\u00A3\u00B01 0 00 cos(\u00CF\u0086) sin(\u00CF\u0086) 0 \u00E2\u0088\u0092 sin(\u00CF\u0086) cos(\u00CF\u0086) \u00EF\u00A3\u00B9\u00EF\u00A3\u00BB (4.1) RY (\u00CE\u00B8) = \u00EF\u00A3\u00AE\u00EF\u00A3\u00B0cos(\u00CE\u00B8) 0 \u00E2\u0088\u0092 sin(\u00CE\u00B8)0 1 0 sin(\u00CE\u00B8) 0 cos(\u00CE\u00B8) \u00EF\u00A3\u00B9\u00EF\u00A3\u00BB (4.2) RZ(\u00CF\u0088) = \u00EF\u00A3\u00AE\u00EF\u00A3\u00B0 cos(\u00CF\u0088) sin(\u00CF\u0088) 0\u00E2\u0088\u0092 sin(\u00CF\u0088) cos(\u00CF\u0088) 0 0 0 1 \u00EF\u00A3\u00B9\u00EF\u00A3\u00BB . (4.3) Multiplication of the matrices in Equation 4.1, 4.2 and 4.3 results in a ro- tation matrix. A rotation matrix is a matrix whereby the multiplication with a vector rotates the vector while preserving its length [51]. In this research a ZY X intrinsic rotation convention is followed. This convention means that rotation is applied along the Z, Y and X axes sequentially. In- trinsic rotation means that rotations are applied with respect to the rotated frame itself, not with respect to the reference frame. Therefore, the rotation matrix to convert any vector represented in the body frame into a vector represented in the navigation frame is derived as Rnb = RZ(\u00CF\u0088)RY (\u00CE\u00B8)RX(\u00CF\u0086) = \u00EF\u00A3\u00AE\u00EF\u00A3\u00B0c\u00CE\u00B8c\u00CF\u0088 s\u00CF\u0086s\u00CE\u00B8c\u00CF\u0088 \u00E2\u0088\u0092 c\u00CF\u0086s\u00CF\u0088 c\u00CF\u0086s\u00CE\u00B8c\u00CF\u0088 + s\u00CF\u0086s\u00CF\u0088c\u00CE\u00B8s\u00CF\u0088 s\u00CF\u0086s\u00CE\u00B8s\u00CF\u0088 + c\u00CF\u0086c\u00CF\u0088 c\u00CF\u0086s\u00CE\u00B8s\u00CF\u0088 \u00E2\u0088\u0092 s\u00CF\u0086c\u00CF\u0088 \u00E2\u0088\u0092s\u00CE\u00B8 c\u00CE\u00B8s\u00CF\u0086 c\u00CE\u00B8c\u00CF\u0086 \u00EF\u00A3\u00B9\u00EF\u00A3\u00BB (4.4) where cos(arg) and sin(arg) are represented as carg and sarg. It is very important to note that the matrix multiplication sequence is not commutative. Any fixed amount of roll, pitch and yaw can result in completely different results. For example, the coordinate frame shown in Figure 4.2 is rotated in ZY X and ZXY sequences both with 90\u00E2\u0097\u00A6 roll, pitch and yaw. The resultant attitudes are entirely different due to the interchange in the sequential rotations along X and Y axis as shown in Figure 4.3 and 4.4. 77 4.2. Rotation Matrix \u00E2\u0088\u00921 \u00E2\u0088\u00920.5 0 0.5 1 \u00E2\u0088\u00921 \u00E2\u0088\u00920.5 0 0.5 1 \u00E2\u0088\u00921 \u00E2\u0088\u00920.5 0 0.5 1 X Z Y Figure 4.2: Coordinate frame before rotation. \u00E2\u0088\u00921 \u00E2\u0088\u00920.5 0 0.5 1 \u00E2\u0088\u00921 \u00E2\u0088\u00920.5 0 0.5 1 \u00E2\u0088\u00921 \u00E2\u0088\u00920.5 0 0.5 1 Z X Y Figure 4.3: Coordinate frame after rotations in ZYX sequence where roll = pitch = yaw = 90\u00E2\u0097\u00A6. Click on the image to see sequential rotations. 78 4.3. Body Frame Attitude Computation \u00E2\u0088\u00921 \u00E2\u0088\u00920.5 0 0.5 1 \u00E2\u0088\u00921 \u00E2\u0088\u00920.5 0 0.5 1 \u00E2\u0088\u00921 \u00E2\u0088\u00920.5 0 0.5 1 Y Z X Figure 4.4: Coordinate frame after rotations in ZXY sequence where roll = pitch = yaw = 90\u00E2\u0097\u00A6. Click on the image to see sequential rotations. 4.3 Body Frame Attitude Computation 4.3.1 Rotation Matrix Update To compute the rotating body frame attitude with respect to the navi- gation frame from the measured angular rate vector \u00CF\u0089, the rotation matrix must be updated. In order to update the rotation matrix a differential equation of the following form [52] needs to be solved with R\u00CC\u0087 = R[\u00CF\u0089\u00C3\u0097] (4.5) where \u00CF\u0089\u00C3\u0097 is the skew-symmetric form of the angular rate vector. Note that the subscript and superscript of the rotation matrix Rnb are left out for convenience. Over an update period from time tk to tk+1, the solution of the equation can be written as Rk+1 = Rkexp( \u00E2\u0088\u00AB tk+1 tk [\u00CF\u0089\u00C3\u0097]dt). (4.6) If it is assumed that the direction of the angular rate vector \u00CF\u0089 remains fixed over the update interval, the integral of the skew symmetric matrix of the angular rate is written as [\u00CF\u0083\u00C3\u0097] = \u00E2\u0088\u00AB tk+1 tk [\u00CF\u0089\u00C3\u0097]dt. (4.7) 79 4.3. Body Frame Attitude Computation where \u00CF\u0083 is the orientation vector. Hence Rk+1 = Rkexp[\u00CF\u0083\u00C3\u0097] (4.8) = RkAk. (4.9) Therefore, Ak is the rotation matrix that transforms a vector in the body frame at the epoch tk+1 into the body frame at epoch tk. Here, tk and tk+1 are the times associated with the start and end of the computation cycle respectively. The exponential term in Equation 4.8 is expanded afterwards, which gives Ak = I + [\u00CF\u0083\u00C3\u0097] + 1 2! [\u00CF\u0083\u00C3\u0097]2 + 1 3! [\u00CF\u0083\u00C3\u0097]3 + 1 4! [\u00CF\u0083\u00C3\u0097]4... = I + [\u00CF\u0083\u00C3\u0097] + 1 2! [\u00CF\u0083\u00C3\u0097]2 \u00E2\u0088\u0092 1 3! \u00CF\u00832[\u00CF\u0083\u00C3\u0097]\u00E2\u0088\u0092 1 4! \u00CF\u00832[\u00CF\u0083\u00C3\u0097]2... = I + (1\u00E2\u0088\u0092 1 3! \u00CF\u00832 + 1 5! \u00CF\u00834 \u00E2\u0088\u0092 ...)[\u00CF\u0083\u00C3\u0097] + ( 1 2! \u00E2\u0088\u0092 1 4! \u00CF\u00832 + 1 6! \u00CF\u00834 \u00E2\u0088\u0092 ...)[\u00CF\u0083\u00C3\u0097]2 (4.10) = I + 1 \u00CF\u0083 sin(\u00CF\u0083)[\u00CF\u0083\u00C3\u0097] + 1 \u00CF\u00832 {1\u00E2\u0088\u0092 cos(\u00CF\u0083)}[\u00CF\u0083\u00C3\u0097]2 (4.11) where I is the identity matrix. The exact representation of the rotation matrix, which relates body frame attitude at times tk and tk+1, is provided in Equation 4.11. However, this exact form can not be directly implemented onboard the goggles\u00E2\u0080\u0099 processor due to the sine and cosine functions. There- fore, Equation 4.11 is approximated by Equation 4.10 where only the first three terms of the infinite series are used. 4.3.2 Orientation Vector Computation In the previous section, it was assumed that the direction of the angular vector \u00CF\u0089 remains constant during the update cycle. However, the direction of \u00CF\u0089 does not remain fixed in space as the body frame is rotating. According to Bortz [53] the rate of change of the orientation vector \u00CF\u0083\u00CC\u0087 is \u00CF\u0083\u00CC\u0087 = \u00CF\u0089 + \u00CE\u00B5\u00CC\u0087. (4.12) Here \u00CE\u00B5\u00CC\u0087 (non-commutativity rate vector) is a component of \u00CF\u0083\u00CC\u0087 due to non- inertially measurable angular motion. The angular rate vector \u00CF\u0089 is provided by the triaxial gyroscope which has its input axes coincident with the axes of the rotating body frame. 80 4.3. Body Frame Attitude Computation Differentiating Equation 4.11, where A\u00CC\u0087k = Ak[\u00CF\u0089\u00C3\u0097], and manipulating vectors gives [52] \u00CF\u0083\u00CC\u0087 = \u00CF\u0089 + 1 2 \u00CF\u0083 \u00C3\u0097 \u00CF\u0089 + 1 \u00CF\u00832 [1\u00E2\u0088\u0092 \u00CF\u0083 sin(\u00CF\u0083) 2(1\u00E2\u0088\u0092 cos(\u00CF\u0083)) ]\u00CF\u0083 \u00C3\u0097 \u00CF\u0083 \u00C3\u0097 \u00CF\u0089. (4.13) Together, the second and third terms on the right hand side of Equa- tion 4.13 account for the non-commutativity vector. After a series of ro- tations, as a result of the non-commutativity terms, the final orientation is dependent on both the individual rotations and the order in which they have occurred [52]. Ignagni in [54] has demonstrated that for practical applica- tions, where the rotation matrix is computed in an iterative manner and integration of the orientation vector rate is conducted over a brief period of time, Equation 4.13 can be simplified to \u00CF\u0083\u00CC\u0087 = \u00CF\u0089 + 1 2 \u00CE\u00B1\u00C3\u0097 \u00CF\u0089 (4.14) where \u00CE\u00B1 = \u00E2\u0088\u00AB tk+1 tk \u00CF\u0089dt. (4.15) Ignagni has also shown that the orientation vector can be found from direct integration of Equation 4.14 rather than solving the differential equation in Equation 4.13. The second term in Equation 4.14 is the simplified form of the non-commutativity rate vector. 4.3.3 Coriolis Effect In inertial navigation calculations, the effects of the earth\u00E2\u0080\u0099s rotation with respect to the inertial frame and the rotation of the navigation frame with respect to the earth, should be accounted for. These effects are known as the Coriolis effects. The earth\u00E2\u0080\u0099s angular velocity with respect to the inertial frame and as observed from the navigation frame can be represented as \u00CF\u0089nie = \u00EF\u00A3\u00AE\u00EF\u00A3\u00B0 \u00E2\u0084\u00A6E cos(\u00CF\u0095)0 \u00E2\u0088\u0092\u00E2\u0084\u00A6E sin(\u00CF\u0095) \u00EF\u00A3\u00B9\u00EF\u00A3\u00BB (4.16) and the angular velocity of the navigation frame with respect to the earth and as observed from the navigation frame is \u00CF\u0089nen = \u00EF\u00A3\u00AE\u00EF\u00A3\u00B0 \u00CE\u00BB\u00CC\u0087 cos(\u00CF\u0095)\u00CF\u0095\u00CC\u0087 \u00E2\u0088\u0092\u00CE\u00BB\u00CC\u0087 sin(\u00CF\u0095) \u00EF\u00A3\u00B9\u00EF\u00A3\u00BB (4.17) 81 4.4. Inertial Navigation Equations where \u00CF\u0095 = latitude \u00CE\u00BB = longitude \u00E2\u0084\u00A6E = earth\u00E2\u0080\u0099s rotation rate 7.292115\u00C3\u0097 10\u00E2\u0088\u00925 rad/s. Therefore, the angular velocity of the navigation frame with respect to the inertial frame is \u00CF\u0089nin = \u00CF\u0089 n ie + \u00CF\u0089 n en. (4.18) This angular velocity \u00CF\u0089nin is important to account for because the gyroscope sensors measure the angular rate of the body frame with respect to the inertial frame rather than the navigation frame or body frame. 4.4 Inertial Navigation Equations The rate of change in velocity of the body frame with respect to the earth expressed in the navigation frame is represented as [52] v\u00CC\u0087n = an \u00E2\u0088\u0092 (2\u00CF\u0089nie + \u00CF\u0089nen)\u00C3\u0097 vn + gnl (4.19) where vn = velocity of the body frame with respect to the earth expressed in local level frame = [ vN vE vD ]T an = acceleration of body frame observed from navigation frame = [ fN fE fD ]T gnl = plumb bob gravity. Plumb bob gravity is defined as gnl = g n \u00E2\u0088\u0092 \u00E2\u0084\u00A6 2 ER0(\u00CF\u0095) 2 [ sin(2\u00CF\u0095) 0 cos(2\u00CF\u0095) ]T . (4.20) Here gn = normal gravity = [ 0 0 9.81 ]T m/s2 R0(\u00CF\u0095) = earth\u00E2\u0080\u0099s radius at latitude \u00CF\u0095 = ( \u00E2\u0088\u009A (R2e cos(\u00CF\u0095)) 2 + (R2p sin(\u00CF\u0095)) 2)/( \u00E2\u0088\u009A (Re cos(\u00CF\u0095))2 + (Rp sin(\u00CF\u0095))2) Re = earth\u00E2\u0080\u0099s equatorial radius (6378137 m) Rp = earth\u00E2\u0080\u0099s polar radius (6356752.3142 m). 82 4.4. Inertial Navigation Equations The change in velocity vn during the update interval between tk and tk+1 can be calculated from integrating Equation 4.19 which is represented as \u00E2\u0088\u0086vnk+1 = \u00E2\u0088\u00AB tk+1 tk andt\u00E2\u0088\u0092 \u00E2\u0088\u00AB tk+1 tk (2\u00CF\u0089nie + \u00CF\u0089 n en)\u00C3\u0097 vnkdt+ \u00E2\u0088\u00AB tk+1 tk gnl dt. (4.21) To solve Equation 4.21, the three integral terms need to be calculated. The first integral term calculation needs further derivation, whereas the integral of the second and the third term can be calculated directly. The acceleration vector of the body frame observed from the body frame can be transformed to an acceleration vector observed from the navigation frame by multiplication with the rotation matrix, i.e. an = Rnb a b (4.22) where ab = [ abx a b y a b z ]T is the acceleration vector of the body frame observed from the body frame. Now, the incremental velocity change of the body frame as observed from the navigation frame is defined as unk = \u00E2\u0088\u00AB k+1 k Rnb a bdt = \u00E2\u0088\u00AB k+1 k RkAka bdt. (4.23) As the rotation matrix Rnb varies continuously with time over the update interval, it can be written in terms of the matrices Rk and Ak. The matrix Rk is the value of R n b at time tk and Ak is a matrix representing the trans- formation of the body frame at time tk+1 to the body frame at the start of the update interval tk. Equation 4.23 can be revised as unk = Rk \u00E2\u0088\u00AB k+1 k Aka bdt. (4.24) Utilizing Equation 4.10 and ignoring second and higher order terms unk = Rk ( \u00E2\u0088\u00AB k+1 k abdt+ \u00E2\u0088\u00AB k+1 k [\u00CF\u0083\u00C3\u0097]abdt). (4.25) The incremental velocity change of the body frame as observed from the body frame is ubk = \u00E2\u0088\u00AB k+1 k abdt. (4.26) 83 4.5. Rotation Angle Computation Using integration by parts, Equation 4.25 and Equation 4.26 can be used to write unk = Rk ( ubk + 1 2 [\u00CF\u0083\u00C3\u0097]ubk + 1 2 \u00E2\u0088\u00AB k+1 k ([\u00CF\u0083\u00C3\u0097]ab \u00E2\u0088\u0092 [\u00CF\u0083\u00CC\u0087\u00C3\u0097]ubk)dt ) . (4.27) In this study, it is assumed that ab and \u00CF\u0083 remain constant for the update interval. Therefore, as shown in [52], the integral term in Equation 4.27 becomes identically zero. After solving the integral terms of Equation 4.27, the velocity vector of the body frame as observed in the navigation frame at time tk+1 can be written as vnk+1 = v n k + u n k + g n l (tk+1 \u00E2\u0088\u0092 tk). (4.28) In general the contribution of the Coriolis terms in Equation 4.19 is small compared to the other terms in the equation [52]. Therefore, it is sufficient to include the Coriolis corrections in Equation 4.19 at a relatively lower rate than the rotation matrix computation rate. Moreover, in this study, \u00CF\u0089nen is considered negligible due to the low velocity of the athlete and very short duration (aero phase of the jump only) of the trajectory determination. Finally, the position can be found by integrating the velocity. By apply- ing the fourth order Runge-Kutta integration as in [52], the position at time tk+1 can be derived as rnk+1 = r n k + vnk\u00E2\u0088\u00921 + 4v n k + v n k+1 6 (tk+1 \u00E2\u0088\u0092 tk). (4.29) 4.5 Rotation Angle Computation In the proposed trajectory determination algorithm, as will be described in Section 4.7, rotation angles, i.e. roll, pitch and yaw, will be computed from accelerometer and magnetometer data. Rotation angles can also be derived from the rotation matrix, which is determined from the gyroscope data. In this section, the procedures for computing the rotation angles are demonstrated. 4.5.1 Roll and Pitch from Accelerometer If the body frame is in steady motion, i.e. no specific acceleration on the body frame, then the triaxial accelerometer experiences antigravity only. In such a condition, the triaxial readings represent the vertically upward antigravity axis with respect to the body frame axes. Therefore, from the 84 4.5. Rotation Angle Computation Z (vertically downwards) n - g n a b z a b x x b tilt a (outward of the page) b y Z (vertically downwards) n - g n a b z a b y y b tilt a (outward of the page) b x Figure 4.5: Roll and pitch calculation from accelerometer data. normalized accelerometer readings, it is possible to find the opposite unit vector of the Z-axis of the navigation frame. On the other hand, the third row of the rotation matrix in Equation 4.4 represents the unit vector along the Z-axis of the navigation frame. Exploiting these facts, it is possible to find the roll and pitch angle from the accelerometer data as follows. Figure 4.5 shows the body frame acceleration in the presence of anti- gravity only. From the first element of the third row of the rotation matrix in Equation 4.4 and the normalized X-axis body frame acceleration, pitch can be computed as follows: sin(xbtilt) = abx |gnl | = \u00E2\u0088\u0092(\u00E2\u0088\u0092 sin(\u00CE\u00B8)) (4.30) Therefore, xbtilt = sin \u00E2\u0088\u00921( abx |gnl | ) (4.31) and \u00CE\u00B8 = xbtilt. (4.32) Here xbtilt is the angle between the horizontal plane and X b. In [55], it has been proven that the sensitivity of Xb of the body frame decreases as it comes closer to the vertical axis. Therefore, to increase the sensitivity of the xbtilt calculation, the tangent function is used instead of the sine function in Equation 4.31 as below: xbtilt = tan \u00E2\u0088\u00921( abx\u00E2\u0088\u009A ab2y + a b2 z ). (4.33) 85 4.5. Rotation Angle Computation From the second element of the third row of the rotation matrix in Equation 4.4 and the normalized Y -axis body frame acceleration, roll can be computed as follows: sin(ybtilt) = aby \u00E2\u0088\u0092|gnl | = cos(\u00CE\u00B8) sin(\u00CF\u0086). (4.34) Therefore, ybtilt = sin \u00E2\u0088\u00921( aby \u00E2\u0088\u0092|gnl | ) (4.35) and \u00CF\u0086 = sin\u00E2\u0088\u00921 (sin(ybtilt) cos(\u00CE\u00B8) ) . (4.36) Here ybtilt is the angle between the horizontal plane and Y b. For reason- ing similar to that in [55], the tangent function is used instead of the sine function in Equation 4.35 as follows: ybtilt = \u00E2\u0088\u0092 tan\u00E2\u0088\u00921( aby\u00E2\u0088\u009A ab2x + a b2 z ). (4.37) Sign Ambiguity in Roll and Pitch From Equation 4.36 and Equation 4.32 it is not possible to find the quadrant in which the roll and pitch fall. Therefore, lookup tables are con- structed where acceleration signs along orthogonal axes are used to extract the quadrant information of the rotation angles. Table 4.1 and Table 4.2 are lookup tables required to resolve the rotation angle ambiguity that occurs while computing roll and pitch from acceleration data. Table 4.1: Quadrant information of roll from acceleration signs. Roll quadrant Sign of aby Sign of a b z 1 - - 2 - + 3 + + 4 + - 86 4.5. Rotation Angle Computation Table 4.2: Quadrant information of pitch from acceleration signs. Pitch quadrant Sign of abx Sign of a b z 1 + - 2 + + 3 - + 4 - - 4.5.2 Yaw from Magnetometer The triaxial magnetometer gives the magnetic component information vector of the earth\u00E2\u0080\u0099s magnetic field in the body frame. This vector is repre- sented as mb = [ mbx m b y m b z ]T . From this vector, it is possible to deduce the yaw of the body frame with respect to the navigation frame. Even though the earth\u00E2\u0080\u0099s magnetic north pole and geographic north pole do not entirely coincide, the error in yaw due to this phenomenon is negligible for the purpose of this study. The magnetometer sensitivity decreases with an increase in xbtilt and ybtilt angles [56]. Therefore, it is important to realign the body frame Z-axis with the navigation frame Z-axis before computing yaw. For this purpose, rotation must be applied to first remove the roll angle followed by a second rotation that removes the pitch angle. The rotation matrix required is R\u00CF\u0086\u00CE\u00B8 = \u00EF\u00A3\u00AE\u00EF\u00A3\u00B0cos(\u00CE\u00B8) 0 \u00E2\u0088\u0092 sin(\u00CE\u00B8)0 1 0 sin(\u00CE\u00B8) 0 cos(\u00CE\u00B8) \u00EF\u00A3\u00B9\u00EF\u00A3\u00BB\u00EF\u00A3\u00AE\u00EF\u00A3\u00B01 0 00 cos(\u00CF\u0086) sin(\u00CF\u0086) 0 \u00E2\u0088\u0092 sin(\u00CF\u0086) cos(\u00CF\u0086) \u00EF\u00A3\u00B9\u00EF\u00A3\u00BB = \u00EF\u00A3\u00AE\u00EF\u00A3\u00B0cos(\u00CE\u00B8) sin(\u00CF\u0086) sin(\u00CE\u00B8) \u00E2\u0088\u0092 cos(\u00CF\u0086) cos(\u00CE\u00B8)0 cos(\u00CF\u0086) sin(\u00CF\u0086) sin(\u00CE\u00B8) \u00E2\u0088\u0092 sin(\u00CF\u0086) cos(\u00CE\u00B8) cos(\u00CF\u0086) cos(\u00CE\u00B8) \u00EF\u00A3\u00B9\u00EF\u00A3\u00BB . (4.38) The order of the matrix multiplication is important, as matrix multiplica- tion is not commutative. From the magnetic field readings and R\u00CF\u0086\u00CE\u00B8 the magnetic field components along the X and Y axes of the realigned body frame can be computed as Xh = m b x cos(\u00CE\u00B8) +m b y sin(\u00CF\u0086) sin(\u00CE\u00B8)\u00E2\u0088\u0092mbz cos(\u00CF\u0086) sin(\u00CE\u00B8) (4.39) Yh = m b y cos(\u00CF\u0086) +m b z sin(\u00CF\u0086). (4.40) Finally, the yaw angle can be found as \u00CF\u0088 = tan\u00E2\u0088\u00921( Yh Xh ). (4.41) 87 4.5. Rotation Angle Computation Similar to the roll and pitch angle, the quadrant information should be resolved from the signs of Xh and Yh. The necessary lookup table is given in Table 4.3. Table 4.3: Quadrant information of yaw from magnetic field component signs. For \u00CF\u0086 Xh < 0 pi \u00E2\u0088\u0092 tan\u00E2\u0088\u00921( YhXh ) Xh > 0, Yh < 0 \u00E2\u0088\u0092 tan\u00E2\u0088\u00921( YhXh ) Xh > 0, Yh > 0 2pi \u00E2\u0088\u0092 tan\u00E2\u0088\u00921( YhXh ) Xh = 0, Yh < 0 pi 2 Xh = 0, Yh > 0 3 2pi 4.5.3 Rotation Angles from Rotation Matrix As previously discussed, in an INS algorithm the rotation matrix is up- dated in each cycle from the angular rate readings of the triaxial gyroscope. For the proposed trajectory determination algorithm, which is described later in Section 4.7, rotation angles need to be computed directly from the rotation matrix in each update cycle. For a rotation matrix Rnb = \u00EF\u00A3\u00AE\u00EF\u00A3\u00B0R11 R12 R13R21 R22 R23 R31 R32 R33 \u00EF\u00A3\u00B9\u00EF\u00A3\u00BB (4.42) the roll, pitch and yaw angles can be computed as \u00CF\u0086 = tan\u00E2\u0088\u00921( R32 R33 ) (4.43) \u00CE\u00B8 = \u00E2\u0088\u0092 sin\u00E2\u0088\u00921(R31) (4.44) \u00CF\u0088 = tan\u00E2\u0088\u00921( R21 R11 ). (4.45) (4.46) As it will later be shown, the difference between the computed rotation angles from the rotation matrix and the computed rotation angles from the accelerometer and magnetometer serve as the observation vector for a particular Linear Kalman Filter (LKF). Therefore, to make both sets of 88 4.6. Sensor Errors rotation angles coherent with each other, quadrant ambiguities are resolved for the set computed from the rotation matrix using Tables 4.1, 4.2 and 4.3 as well. 4.6 Sensor Errors For the implementation of the INS algorithm, it is important to mitigate the error generated by the noise present in the accelerometer and gyroscope readings. Among the number of error sources, such as misalignment, scaling, bias, drift, temperature effects, etc., the accelerometer and gyroscope biases are the most significant and should be accounted for [57]. Accelerometer bias has a quadratic effect on the position derived from the INS algorithm. The gyroscope bias, which is even more common, has an equally large impact on the position estimate. Due to the bias in the gyroscope, the inertial measure- ment unit is under the impression that it is rotating and the INS equations will not account for gravity correctly. This can lead to a false movement calculation due to an acceleration of maximum 9.81 ms\u00E2\u0088\u00922 depending on how much more the body frame is rotated due to the bias. Generally, for MEMS sensors the gyroscope bias is more prominent than that of the accelerometers and it also tends to drift significantly with time. Therefore, the gyroscope bias contributes to the majority of the errors in the INS algorithm position estimate. According to the research conducted by Waegli et al., it was revealed that for given MEMS sensor characteristics the misalignment, scaling, drift and bias errors cannot be decorrelated efficiently for limited integration pe- riods [38]. Therefore, in this research, a simplified error model is considered where the bias terms are only taken into account and modeled as 1st order Gauss-Markov processes as in [38]. The inertial accelerometer raw readings in discrete form can be modeled as a\u00CC\u0083k = ak + b a k + w a k (4.47) where a\u00CC\u0083k is the measured acceleration, ak is the true body frame accelera- tion, bak is the bias in acceleration and w a k is the measurement noise at epoch tk. Note that the superscript notion of the body frame is excluded for the ease of representation. According to the 1st order Gauss-Markov model, the bias can be represented with a discrete differential equation [58] as bak+1 = (1\u00E2\u0088\u0092 \u00CE\u00B2a\u00E2\u0088\u0086t)bak + \u00E2\u0088\u009A 2\u00CE\u00B2a\u00CF\u0083a2\u00E2\u0088\u0086twak (4.48) 89 4.7. Proposed Algorithm where \u00CF\u0083a2 is the square of the amplitude of the acceleration bias power spectral density, \u00CE\u00B2a is the reciprocal of the acceleration process correlation time \u00CF\u0084 a and \u00E2\u0088\u0086t is the sampling interval. These parameters can be calculated from the widely practiced autocorrelation analysis, which is described in [58]. Similarly, the measured angular rate and the bias in angular rate are modeled as \u00CF\u0089\u00CC\u0083k = \u00CF\u0089k + b \u00CF\u0089 k + w \u00CF\u0089 k (4.49) b\u00CF\u0089k+1 = (1\u00E2\u0088\u0092 \u00CE\u00B2\u00CF\u0089\u00E2\u0088\u0086t)b\u00CF\u0089k + \u00E2\u0088\u009A 2\u00CE\u00B2\u00CF\u0089\u00CF\u0083\u00CF\u00892\u00E2\u0088\u0086tw\u00CF\u0089k (4.50) where \u00CF\u0089\u00CC\u0083k is the measured body frame angular rate with respect to the navigation frame, \u00CF\u0089k is the true body frame angular rate, b \u00CF\u0089 k is the bias in angular rate and w\u00CF\u0089k is the measurement noise at epoch tk. The bias model parameter \u00CF\u0083\u00CF\u00892 is the square of the amplitude of the angular rate bias power spectral density and \u00CE\u00B2\u00CF\u0089 is the reciprocal of the angular rate process correlation time \u00CF\u0084\u00CF\u0089. The errors in magnetometer readings only affect the computation of the yaw angle of the body frame with respect to the navigation frame. These errors are not related to the critical gravity cancelation process. Hence, magnetometer errors do not impose any large error in the position estimation for a short duration. Even if there is a considerable bias in the magnetometer reading, it can be safely ignored as the goal of this research is to find the relative jump trajectory rather than the absolute trajectory. Relative jump trajectory means the evolving position estimates of the athlete with respect to a reference point, e.g. jump start position, whereas the absolute trajectory represents the evolving position estimates in global positioning terms, e.g. latitude and longitude. Therefore, magnetometer error estimation is ignored in this work and measured magnetometer readings, i.e., m\u00CC\u0083k are directly used for yaw computation to reduce the computational load. 4.7 Proposed Algorithm In this section, the proposed algorithm for jump trajectory determina- tion is demonstrated. A loosely coupled GPS/INS integrated algorithm is proposed. A flow diagram of the overall proposed algorithm is given in Fig- ure 4.8. Three Linear Kalman Filters (LKF) are used to correct the sensor data prior to computation of the position estimates with the INS algorithm. Later, an Extended Kalman Filter (EKF) is implemented to correct the trajectory from the INS with the position updates from the GPS. Even 90 4.7. Proposed Algorithm though, only the trajectory during the jump is required, the algorithm is applied starting at an epoch prior to the jump start up to the epoch of the jump end. This is done in order to capture the dynamics undertaken as a preparation for the jump. The key and novel aspects of the proposed trajectory determination algorithm are as follows: 1. In general, the INS or GPS/INS integrated algorithm corrects the accelerometer and gyroscope readings with a bias estimate from the previous update interval, as in [59], prior to resolving the INS equa- tions. This is because the sensor error estimates for any update interval are only available once the INS or GPS/INS navigation equations are solved. Computing in this manner might not address the bias drift correctly. Moreover, misalignment error and temperature effects are also not decorrelated from the bias error, which may cause the sensor error to drift rapidly. In the proposed algorithm, sensor data fusion is exploited to estimate the bias in the readings with the sensor data itself. LKF is used for the sensor fusion and to estimate the bias in one sensor with the help of another prior to resolving the INS equations. These estimated bias terms for the current update cycle are used to correct the sensor readings, which results in a better position estimate from the INS algorithm. 2. The sensor bias error estimate from the sensor signal itself also pro- vides leverage in the GPS/INS integrated algorithm. In GPS/INS integration with augmented states, i.e. where the sensor errors are also incorporated in the state vector, typically there are no sensor error terms in the observation vector as in [60]. However, the bias es- timates made prior to the INS algorithm may serve as the observations for sensor biases in the EKF for the GPS/INS integration. Therefore, two separate sensor error estimations are conducted in a particular update cycle. The estimated error from the EKF for an update cycle is also used in the next update cycle to make primary corrections in the accelerometer and gyroscope signal. In this way, the error esti- mates from the EKF contribute a better initial estimate of the sensor bias in the next update cycle. Therefore, the augmented observations for GPS/INS integration leads to more accurate overall sensor error estimation. 3. In previous work where acceleration is used to correct the inclination estimate, e.g. [59], the body frame is assumed to be entirely stationary and the relevant Kalman filter parameters also remain static. This is 91 4.7. Proposed Algorithm never the case for an athlete performing a jump. High specific accelera- tion, i.e. non gravitational acceleration, is particularly involved before the take-off and after the landing. On the other hand, during the aero phase the body frame accelerometers feel no acceleration. Therefore, the captured sensor data during and neighboring jumps are segmented into different portions depending on the athletes dynamics. The LKF parameters are customized for each segment and varied over time to properly address the different acceleration scenarios. 4. Even though the error estimation processes for accelerometer and gy- roscope through sensor fusion are not entirely independent from each other, certain assumptions are made to decouple the states of the accelerometer and gyroscope biases. This is done to reduce computa- tional load as the state vectors for the accelerometer and gyroscope errors are updated separately using transition matrices with smaller dimensions (3\u00C3\u00973) rather than using a single matrix with large dimen- sion (6\u00C3\u0097 6). The proposed algorithm is explained chronologically below. The basic equations for the implementation of LKF and EKF are not explained in this work. A detailed demonstration of the LKF and EKF can be found in [61] and [52]. 4.7.1 Acceleration Error Compensation For an update cycle, an initial estimate of the body frame acceleration is done with a LKF, which is termed as LKF1. The bias error estimated from the GPS/INS integration in the previous cycle is used to correct the accelerometer readings. This initial estimate is the observation vector for the LKF1. As the acceleration of an athlete changes gradually with time, acceleration of the current cycle can be estimated from the previous cycle with a rotation matrix. The rotation matrix rotates the acceleration in the previous cycle to reflect into the body frame attitude of the current update cycle. To adjust to the changing dynamics, enough uncertainty is injected into the process through considerable process noise. The state and measurement equations of the LKF1 in discrete form are ack+1 = A T k a\u00CC\u0082k + q a k (4.51) zack+1 = a c k+1 + w a k+1 (4.52) where 92 4.7. Proposed Algorithm ack+1 = initial corrected acceleration for the current update cycle ATk = state transition matrix, transpose of Ak a\u00CC\u0082k = final estimate of the acceleration for the previous update cycle starting at epoch tk qak = process noise of acceleration with covariance matrix Q a k zack+1 = observation vector for LKF 1 wak+1 = measurement noise of acceleration with covariance matrix R a k+1. Note that, the transition matrix ATk is the transpose of the rotation matrix derived in Equation 4.10. This is because the acceleration is reflected for- ward in time in the state equation, whereas Ak reflects the vector backwards in time. The measurement vector zack+1 for the current cycle is computed as zack+1 = a\u00CC\u0083k+1 \u00E2\u0088\u0092 b\u00CC\u0082 aE k . (4.53) Here, a\u00CC\u0083k+1 is the accelerometer reading for the current update cycle starting at tk+1 and b\u00CC\u0082 aE k is the acceleration error estimate from the EKF for GPS/INS integration in the previous cycle starting at tk. With the initial acceleration estimate a\u00CC\u0082ck+1 from the LKF 1, the algorithm proceeds to estimate the acceleration bias for the current cycle derived via sensor fusion. Another LKF, namely LKF2, is applied for this purpose. According to the 1st order Gauss-Markov model in Equation 4.48, the state and measurement equations of LKF2 are represented as baLk+1 = (1\u00E2\u0088\u0092 \u00CE\u00B2a\u00E2\u0088\u0086t)baLk + \u00E2\u0088\u009A 2\u00CE\u00B2a\u00CF\u0083a2\u00E2\u0088\u0086twaLk (4.54) zbaLk+1 = b aL k+1 + w aL k+1 (4.55) where baLk+1 = bias error in the acceleration derived via sensor fusion for the current update cycle zbaLk+1 = observation vector for LKF 2. waLk+1 = measurement noise of acceleration bias with covariance matrix RaLk+1. The observation vector is found from the initial acceleration estimate and the sensor readings as follows: zbaLk+1 = a\u00CC\u0083k+1 \u00E2\u0088\u0092 a\u00CC\u0082ck+1. (4.56) The acceleration bias estimated from LKF2, b\u00CC\u0082 aL k+1, is claimed to be derived via sensor fusion as the gyroscope readings are used in LKF1 to calculate 93 4.7. Proposed Algorithm a\u00CC\u0082ck+1. Moreover, the bias estimate from the EKF in the previous cycle indi- rectly plays a role in b\u00CC\u0082 aL k+1 calculation as b\u00CC\u0082 aE k is used to derive a\u00CC\u0082 c k+1. Here it is assumed that the gyroscope errors are decoupled from the acceleration errors. The rationale behind this assumption is that the gyroscope readings are recursively corrected with a separate LKF, as will be later discussed, for each update cycle. Therefore, the estimated error in the accelerome- ter should become independent from the gyroscope error with time in an iterative manner. Finally, with the estimation of the acceleration error from LKF2, b\u00CC\u0082 aL k+1, the accelerometer readings are corrected for the current update cycle as a\u00CC\u0082k+1 = a\u00CC\u0083k+1 \u00E2\u0088\u0092 b\u00CC\u0082aLk+1. (4.57) This final acceleration estimate a\u00CC\u0082k+1 is used for the rest of the algorithm including in the gyroscope error estimation and in the INS equations for the current update cycle. 4.7.2 Angular Rate Error Compensation The gyroscope readings provide the angular rate of the body frame with respect to the inertial frame. The Coriolis effects are accounted for to find the angular rate of the body frame relative to the navigation frame following the equation \u00CF\u0089\u00CC\u0083k+1 = \u00CF\u0089\u00CC\u0083 b ibk+1 \u00E2\u0088\u0092 (I + [\u00CE\u00B6\u00CC\u0082k\u00C3\u0097])T R\u00CC\u0082 T k\u00CF\u0089 n ie (4.58) where \u00CF\u0089\u00CC\u0083bibk+1 = gyroscope reading for the update cycle starting at tk+1 [\u00CE\u00B6\u00CC\u0082k\u00C3\u0097] = skew symmetric matrix of the estimated attitude error vector \u00CE\u00B6\u00CC\u0082k in INS mechanization for previous update cycle starting at tk R\u00CC\u0082k = estimated rotation matrix for previous update cycle starting at tk \u00CF\u0089nie = earth\u00E2\u0080\u0099s angular velocity with respect to the inertial frame and as observed from the navigation frame. It should be noted that the angular rate of the navigation frame relative to the earth is ignored. This is because the athlete\u00E2\u0080\u0099s velocity relative to the earth is insignificant enough to have no impact on the rotation matrix. The attitude error vector is computed from the state estimation with EKF, which is discussed in Subsection 4.7.4. 94 4.7. Proposed Algorithm IMU LKF1 LKF2 Update Rotation Matrix LKF3 Rotate EKF EKF + - - Error Compensation Block *Dotted line represents feedback Figure 4.6: Error compensation block. Error Compensation Block Rotate Gravity Computation EKF INS Mechanization Block + + - + + + + *Dotted line represents feedback Figure 4.7: INS mechanization block. 95 4.7. Proposed Algorithm The rotation matrix for the current cycle is computed initially using the gyroscope error estimation from the previous cycle R\u00CC\u0083k+1 = R\u00CC\u0082k(I + [\u00CE\u00B6\u00CC\u0082k\u00C3\u0097])A\u00CC\u0083k (4.59) where A\u00CC\u0083k is computed following Equation 4.10 using the angular rate vector (\u00CF\u0089\u00CC\u0083k+1 \u00E2\u0088\u0092 b\u00CC\u0082\u00CF\u0089Ek ). The estimated error in the gyroscope computed via EKF in the previous cycle is b\u00CC\u0082 \u00CF\u0089E k . It is important to note that the transpose of the matrix A\u00CC\u0083k is actually used as the transition matrix in Equation 4.51, i.e. ATk . From R\u00CC\u0083k+1, roll, pitch and yaw angles are computed using the method described in Subsection 4.5.3. The computed set of angles are represented as \u00CE\u00B3\u00CC\u0083k+1 = [ \u00CF\u0086\u00CC\u0083k+1 \u00CE\u00B8\u00CC\u0083k+1 \u00CF\u0088\u00CC\u0083k+1 ]T (4.60) Another separate set of rotation angles can be derived from the ac- celerometer and magnetometer data. The method described in Subsec- tion 4.5.1 is used to calculate the roll and pitch angles from the acceleration data. The normalized estimated acceleration for the current update cycle, a\u00CC\u0082k+1/|a\u00CC\u0082k+1|, is used as the acceleration signal for this purpose. The acceler- ation is normalized to ensure that the amplitude of the acceleration vector is 1 g. Subsection 4.5.2 is used with the available magnetometer data to find the yaw angle. The set of rotation angles derived from the accelerometer and magnetometer is represented as z\u00CE\u00B3k+1 = [ z\u00CF\u0086k+1 z\u00CE\u00B8k+1 z\u00CF\u0088k+1 ]T . (4.61) The difference between the two separate estimates of the rotation angle sets is used to derive the observation vector of the LKF, namely LKF3, for the gyroscope error estimation. The state and measurement equations for LKF3 are b\u00CF\u0089Lk+1 = (1\u00E2\u0088\u0092 \u00CE\u00B2\u00CF\u0089\u00E2\u0088\u0086t)b\u00CF\u0089Lk + \u00E2\u0088\u009A 2\u00CE\u00B2\u00CF\u0089\u00CF\u0083\u00CF\u00892\u00E2\u0088\u0086tw\u00CF\u0089Lk (4.62) zb\u00CF\u0089Lk+1 = b \u00CF\u0089L k+1 + w \u00CF\u0089L k+1 (4.63) where b\u00CF\u0089Lk+1 = bias error in the gyroscope derived via sensor fusion for the current update cycle zb\u00CF\u0089Lk+1 = observation vector for LKF 3. w\u00CF\u0089Lk+1 = measurement noise in angular rate with covariance matrix R\u00CF\u0089Lk+1. 96 4.7. Proposed Algorithm The observation vector of LKF3 is computed as zb\u00CF\u0089Lk+1 = \u00CE\u00B3\u00CC\u0083k+1 \u00E2\u0088\u0092 z\u00CE\u00B3k+1 \u00E2\u0088\u0086t . (4.64) The estimated error in angular rate using sensor fusion, b\u00CC\u0082 \u00CF\u0089L k+1, is found as the result of LKF3. This error estimate is used to correct the rotation angle derived from the rotation matrix. Therefore, the estimated rotation angle set is \u00CE\u00B3\u00CC\u0082k+1 = \u00CE\u00B3\u00CC\u0083k+1 \u00E2\u0088\u0092 b\u00CC\u0082 \u00CF\u0089L k+1\u00E2\u0088\u0086t. (4.65) Henceforth, the estimated rotation matrix for the current update cycle R\u00CC\u0082k+1 is computed from the set of corrected rotation angles \u00CE\u00B3\u00CC\u0082k+1 using Equa- tion 4.4. Figure 4.6 shows the flow diagram of the proposed error compensation technique for a particular update cycle starting at epoch tk+1. The Inertial Measurement Unit (IMU) represents the set of body frame sensors, which provides raw data to the algorithm. Among the four, two of the final outputs after the error compensation, the estimated acceleration a\u00CC\u0082k+1 and the esti- mated rotation matrix R\u00CC\u0082k+1, are fed to the INS mechanization equations. The other two final outputs, estimated acceleration error and estimated an- gular rate error from sensor fusion, namely b\u00CC\u0082 aL k+1 and b\u00CC\u0082 \u00CF\u0089L k+1, are fed to the EKF for the GPS/INS integration as a part of the observation vector. 4.7.3 INS Mechanization The estimated body frame acceleration a\u00CC\u0082k+1 and rotation matrix R\u00CC\u0082k+1 are used to solve the INS equations derived in Section 4.4. The initial ve- locity of the athlete required for the INS equations can be directly achieved from the GPS receiver data. Finally, the position update at epoch tk+1 in relation to the navigation frame is derived from the INS equations, which are termed as the vector rnINSk+1 = [ rnINSxk+1 r nINS yk+1 r nINS zk+1 ]T . Here, the sub- scripts x, y and z denote the directions toward north, east and vertically downwards respectively. The pertinent velocity vector is represented as vnINSk+1 = [ vnINSxk+1 v nINS yk+1 v nINS zk+1 ]T . Figure 4.7 depicts the flow diagram of INS mechanization steps. The final output, i.e. the position vector rnINSk+1 , is fed forward to the GPS/INS integration. 97 4.7. Proposed Algorithm 4.7.4 Integration with GPS To estimate the optimum position updates, the output of the INS needs to be integrated with the position updates from the GPS with a suitable estimator. The Extended Kalman Filter (EKF) is one of the most popular estimators for such integrations [62]. The EKF is an adaptation of the LKF to nonlinear functions, such as the INS perturbation model provided later in Equation 4.66-4.68, that is approximated by linearization. A detailed description of the EKF can be found in [52, 61]. The type of integration, where GPS position and/or velocity updates are directly used to correct the INS output, is called a loosely coupled system. In this study, a loosely coupled system is implemented for integrating GPS position updates. The integration strategy can also be either open looped or close looped. In an open looped system, the INS mechanization operates without being concerned about the presence of the estimator or external data. The estimated state vectors are used to correct the INS output only. However, in an open loop system, the error in INS mechanization with low cost MEMS senors can grow unbounded and propagate significantly large error within a short period of time. As a result, unbounded error observa- tions are delivered to the estimator, such as the Kalman filter. This causes problems in the linear filter since only small errors are allowed due to the linearization process [63]. Therefore, a closed loop system is required to continuously compensate for the error in INS mechanization. In a closed loop system, error estimates from the estimator are fed back to correct the sensor outputs and other mechanization parameters. In a closed loop loosely coupled system, the GPS/INS integration is typ- ically done by estimating the error states in the INS mechanization. The basic estimator consists of nine navigation error states consisting of three position, three velocity and three attitude error states [41]. However, due to the error in sensors, the state vector is augmented by incorporating ac- celeration and angular rate error states resulting in a fifteen element state vector. The INS error state models are obtained by the perturbation of the mechanization equations [41]. The perturbation analysis is not within the scope of this research and, therefore, is not demonstrated here. Several well documentated perturbation analyses are available as in [64] and [65]. As given in [66], the obtained model is expressed as a series of differential 98 4.7. Proposed Algorithm IMU Error Compensation Block INS Mechanization Block GPS Receiver EKF + - + + Output *Dotted line represents feedback Figure 4.8: Flow diagram of jump trajectory determination. equations as follows: \u00CE\u00B4r\u00CC\u0087n = \u00CE\u00B4vn (4.66) \u00CE\u00B4v\u00CC\u0087n = \u00E2\u0088\u00922[\u00CF\u0089nie\u00C3\u0097]\u00CE\u00B4vn \u00E2\u0088\u0092 an\u00CE\u00B6n + RnbbaE (4.67) \u00CE\u00B6\u00CC\u0087 n = \u00E2\u0088\u0092[\u00CF\u0089nie\u00C3\u0097]\u00CE\u00B6n + Rnbb\u00CF\u0089E (4.68) where the \u00E2\u0080\u009Cdot\u00E2\u0080\u009Ds denote time derivatives and rn = position error state vector vn = velocity error state vector \u00CE\u00B6n = attitude error state vector baE = accelerometer error state b\u00CF\u0089E = gyroscope error state. Following the above continuous perturbation model and the simplified discrete time approximated model demonstrated in [63], the augmented dis- 99 4.7. Proposed Algorithm crete time fifteen state model can be represented as \u00CE\u00B4xk+1 = \u00EF\u00A3\u00AE\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00B0 \u00CE\u00B4rnk+1 \u00CE\u00B4vnk+1 \u00CE\u00B6nk+1 baEk+1 b\u00CF\u0089Ek+1 \u00EF\u00A3\u00B9\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BB = \u00EF\u00A3\u00AE\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00B0 0 I 0 0 0 0 \u00E2\u0088\u00922[\u00CF\u0089nie\u00C3\u0097] \u00E2\u0088\u0092an Rnbk+1 0 0 0 \u00E2\u0088\u0092[\u00CF\u0089nie\u00C3\u0097] 0 Rnbk+1 0 0 0 (1/\u00E2\u0088\u0086t\u00E2\u0088\u0092 \u00CE\u00B2a) 0 0 0 0 0 (1/\u00E2\u0088\u0086t\u00E2\u0088\u0092 \u00CE\u00B2\u00CF\u0089) \u00EF\u00A3\u00B9\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BB \u00EF\u00A3\u00AE\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00B0 \u00CE\u00B4r\u00CC\u0087nk \u00CE\u00B4v\u00CC\u0087nk \u00CE\u00B6\u00CC\u0087 n k baEk b\u00CF\u0089Ek \u00EF\u00A3\u00B9\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BB\u00E2\u0088\u0086t+ qEk (4.69) where qEk is the driven response at tk+1 due to the presence of the pro- cess white noise during the time interval tk - tk+1 [63]. The corresponding covariance matrix is set as [62] QEk = \u00EF\u00A3\u00AE\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00B0 Qrk 0 0 0 0 0 Qvk 0 0 0 0 0 Q\u00CE\u00B6k 0 0 0 0 0 QaEk 0 0 0 0 0 Q\u00CF\u0089Ek \u00EF\u00A3\u00B9\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BB (4.70) where Qrk = covariance matrix for positioning error Qvk = covariance matrix for velocity error Q\u00CE\u00B6k = covariance matrix for attitude error QaEk = covariance matrix for acceleration error Q\u00CF\u0089Ek = covariance matrix for angular rate error. Here, the bias errors are modeled as a 1st order Gauss-Markov process as given in Equation 4.48 and 4.50. The observation vector for the EKF is derived from the difference be- tween GPS and INS position estimates, and augmented by the two estimated errors from the error compensation block. Therefore, the nine element ob- servation vector is \u00CE\u00B4zk+1 = \u00EF\u00A3\u00AE\u00EF\u00A3\u00AF\u00EF\u00A3\u00B0r nGPS k+1 \u00E2\u0088\u0092 rnINSk+1 b\u00CC\u0082 aL k+1 b\u00CC\u0082 \u00CF\u0089L k+1 \u00EF\u00A3\u00B9\u00EF\u00A3\u00BA\u00EF\u00A3\u00BB (4.71) 100 4.7. Proposed Algorithm where rnGPSk+1 is the position update from GPS at epoch tk+1. As GPS position updates occur at lower frequency than the frequency of the INS position updates, the GPS positions are interpolated to the INS update rate within the time period selected for the jump trajectory determination. Finally, the estimated states are fed back to the error compensation block and the INS mechanization block for error correction. The estimated position error is used to correct the position output from the INS. Therefore, the final estimated position update is given as r\u00CC\u0082nk+1 = r nINS k+1 + \u00CE\u00B4r\u00CC\u0082 n k+1. (4.72) The flow diagram of the entire proposed algorithm is shown in Figure 4.8. In the figure, dotted lines represent feedback. The proposed technique is depicted for a particular update cycle and the index of the update cycle is explicitly mentioned as a subscript in the variable names. 4.7.5 Period Segmentation and Parameter Customization For jump trajectory determination, it is important to capture all the abrupt changes in dynamics undertaken by the athlete in preparation for the jump. Failure to accommodate any small portion of significant change in acceleration prior to the jump may cause a large error in KPV determination, i.e. jump horizontal distance, height and drop. Therefore, the trajectory determination algorithm should be conducted on the data from a time epoch preceding jump start. However, the computation capability of the onboard processor is limited as well as the battery power. Moreover, computation should be completed shortly after the jump to provide real-time feedback on the KPVs to the goggle user. Therefore, due to resource constraints, it is also not feasible to operate the trajectory determination algorithm starting from far behind the jump start epoch to safely accommodate all of the jump preparation dynamics. Another significant difficulty in jump trajectory determination using the Kalman filter is that the rate of change in dynamics varies abruptly within a very short period of time immediately before and after a jump. The jump\u00E2\u0080\u0099s dynamic characteristics also vary a great deal for a wide range of different types of jumps. Therefore, the Kalman filter variables involved in predicting an athlete\u00E2\u0080\u0099s dynamics should be customized and fine-tuned to address highly fluctuating scenarios. To address the above difficulties, an automated algorithm is implemented to select the total trajectory determination period and segment it depend- ing on the acceleration characteristics. As shown in Figure 4.9, the jump 101 4.7. Proposed Algorithm 32 34 36 38 40 42 44 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 Time [s] A m pl itu de o f r es ul ta nt a cc el er at io n ab [g ] Initial Acceleration Period Total trajectory determination period Aero Phase Period Acceleration Gain Period Jump endJump start Figure 4.9: Period segmentation for trajectory determination. preparation period, including the jump, has been segmented into three por- tions depending on the resultant acceleration profile, namely Initial Ac- celeration Period (IAP), Acceleration Gain Period (AGP) and Aero Phase Period (APP). In IAP, the resultant acceleration is steady and there is no abrupt change in the acceleration. In this period, the athlete has only started to prepare for the jump and the average acceleration slowly increases. This period also provides enough time for the Kalman filters to converge. In AGP, the majority of the gain in acceleration (mostly in upward and forward direction) occurs. Rapid fluctuations in the acceleration also take place during the preparation for take-off. This segment is of the most im- portance for the KPV calculation. The last segment, i.e. APP, is the time period when the athlete is in the air. During APP, the acceleration falls to near zero. Any algorithm parameter that uses gravity sensed by the accelerometer should be carefully handled during AGP and APP. This is because during AGP the accelerometer senses specific acceleration in mag- nitude much higher than gravity, whereas during APP the accelerometer senses no gravity at all. The jump trajectory determination algorithm is operated from the start of IAP to the end of APP. The IAP and AGP are automatically selected by the algorithm from the resultant acceleration. APP is predetermined from the jump start epoch 102 4.7. Proposed Algorithm to jump end epoch, which are determined by the algorithm described in Chapter 3. A threshold thAGP is selected to extract the information about the gain in acceleration. If the resultant acceleration exceeds thAGP at any epoch within a considerable period, which is set to 10 s in this study, prior to the jump start, then an indication of the start of AGP is achieved. Another 2 s of safety period is annexed before this epoch to make sure that the rise of acceleration is incorporated into AGP. Therefore, the start of AGP is calculated as, tAGPs = min t { |ab(t)| > thAGP } \u00E2\u0088\u0092 2 s, tjumps \u00E2\u0088\u0092 10 s < t < tjumps (4.73) where tAGPs = start epoch of AGP ab(t) = resultant acceleration tjumps = jump start epoch. Hence the AGP spans from tAGPs to t jump s . Furthermore, the IAP is empir- ically selected to span from tAGPs - 3 s to t AGP s and, as mentioned earlier, APP spans from tjumps to t jump e . Here, t jump e denotes the jump end epoch. Two parameters of the LKFs implemented in the error compensation block are directly dependent on the dynamics of the athlete and, therefore, need to be varied over different segments. One of the parameters is the process noise covariance matrix, Qak, of LKF 1. Due to the very abrupt nature of the acceleration during AGP, Qak is set to a higher value during AGP relative to what it is set to for IAP and APP. The second parameter, and also the more important one, which needs to be tuned during different segments is the observation noise covariance matrix, R\u00CF\u0089Lk+1, in LKF 3. As previously noted, maintaining alignment, i.e. determination of the roll and pitch of the body frame, is important because of its direct relation to gravity cancelation. Therefore, acceleration mea- surements are used to correct the estimate of the body frame attitude found from the gyroscope. For accurate estimation of body frame attitude, the accelerometer should only sense vertically directed acceleration, e.g. anti gravity. However, during AGP, high specific acceleration is involved, which renders the acceleration sensed by the accelerometer in any random direc- tion. On the other hand, during APP there is practically no acceleration to calculate the attitude from. Therefore, the reliability of the rotation an- gles calculated from the accelerometer data varies amongst the segmented periods. As the measurement vector of LKF3 directly depends on the mea- sured rotation angles from the accelerometer, the measurement noise varies 103 4.8. Field Experiments signifcantly during the trajectory determination period. Therefore, the mea- surement noise covariance vector R\u00CF\u0089Lk+1 is considered a tuning parameter for proper estimation of the angular rate error with LKF3. The covariance ma- trix is tuned to a high value during AGP relative to what it is tuned to for IAP. This is done to decrease the reliability of the measurement vector dur- ing AGP. The reliability of the observation for LKF3 becomes even worse during APP due to an absence of any acceleration and, therefore, R\u00CF\u0089Lk+1 is set to a higher value to completely rely on the data achieved from the gyroscopes. 4.8 Field Experiments The version of the Recon goggles, which are equipped with all the sensors needed for trajectory determination, was not yet ready for field data collec- tion when this research was conducted. Therefore, jump data was collected with another data capturing device with similar sensors, named \u00E2\u0080\u0098Minimax\u00E2\u0080\u0099, developed by Catapult [31]. As ski resorts were closed during the time of this research, data was collected for jumps performed on a mountain bike as this closely resembles ski jumps. For benchmark data, a Novatel DL-V3 GPS receiver was carried by the mountain biker while executing the jumps. The Novatel receiver position updates were further refined with \u00E2\u0080\u0098CSRS - Precise Point Positioning\u00E2\u0080\u0099 provided by Natural Resources Canada [67]. The Novatel GPS receiver has a 5 Hz update rate while the Minimax GPS receiver gives position updates at 10 Hz. All the inertial sensors in the Catapult unit have an update rate of 100 Hz. The Catapult unit affixed to the side of the biker\u00E2\u0080\u0099s helmet, as shown in Figure 4.10, while executing jumps to mimic the similar dynamic environ- ment felt by the Recon goggle sensors. The Novatel receiver was carried in a backpack with its antenna projected outside. Synchronization of the data achieved from these two devices is a critical issue. Even though the GPS time stamps of the data captured with the Novatel receiver are reliable, the Catapult GPS time stamps suffer from a bias, which drifts with time. This is due to a varying time delay that occurs when the Catapult microprocessor writes the time stamps. This is a common issue for microsystems such as in Recon goggles, Minimax etc. Therefore, data sets from the different devices are synchronized by cal- culating the minimum Root Mean Square (RMS) value of the difference in positional data readings (latitude and longitude) from different devices. However, this synchronization strategy renders an unavoidable implication. 104 4.8. Field Experiments Figure 4.10: Catapult Minimax unit affixed to the biker\u00E2\u0080\u0099s helmet. Due to the spatial distance between the head mounted Catapult unit and the Novatel GPS antenna in the backpack, the derived jump start and end epochs from the acceleration captured by the Catapult unit\u00E2\u0080\u0099s accelerome- ter will not represent the same jump start and end epochs for the Nova- tel receiver. Therefore, to eliminate the spatial difference, the Catapult is strapped to the antenna pole of the Novatel receiver. The entire setup and a test jump performed by the athlete is shown in Figure 4.11 and Figure 4.12. Carrying the Catapult unit in the backpack also poses some negative im- plications. First of all, the collected acceleration sensor data is more noisy than the data collected with the Catapult affixed to the helmet. This is because the antenna pole is not rigidly attached to the backpack and the backpack itself generates noise from movement due to the flexible coupling with the athlete. Secondly, as the catapult unit is not firmly coupled with the athlete, the detected jump from the Catapult acceleration represents the jump of the Catapult unit itself, not necessarily the jump of the ath- lete. Experimental data has shown that the athlete\u00E2\u0080\u0099s jump occurrence is delayed in time relative to the jump detected by the Catapult unit. The probable reason for this phenomenon is that the backpack settles on the athlete\u00E2\u0080\u0099s shoulder while he is still in the air. However, this should not cause any problem in the performance analysis of the trajectory determination algorithm as all the interest is focused on the trajectory determination, not the detection of jump start or end epochs. The jump start and end are detected with the accelerometer data as described in Chapter 3 and then the detected epochs are used to find the jump region, i.e. the aero phase 105 4.9. Experimental Results Figure 4.11: Catapult Minimax unit strapped to the Novatel GPS receiver antenna pole. [Photograph by Darren Handschuh] region for the Catapult unit. The detected jump region is moved forward suitably to resemble the athlete\u00E2\u0080\u0099s jump and relevant KPVs without any loss of generality. An example of noisy resultant acceleration of the Catapult unit strapped to the antenna pole is shown in Figure 4.13. 4.9 Experimental Results Fourteen jumps were executed for which the Catapult MEMS inertial sensor and GPS receiver data, along with the Novatel receiver data, were collected. The Catapult data was fed to the proposed trajectory determi- nation algorithm and the jump trajectory was determined. The trajectory determination begins from the start of the IAP and stops at the jump end. The trajectory was also determined from the GPS position updates of the Catapult unit. Afterwards, the relevant KPVs, namely jump horizontal dis- tance, height and drop, were calculated easily. Similarly, KPVs are found from the Novatel receiver GPS position updates which are used as the bench- mark KPVs. The Novatel receiver GPS receiver updates with precise point positioning have a horizontal positioning accuracy of \u00E2\u0088\u00BC 2 cm and verti- cal positioning accuracy of \u00E2\u0088\u00BC 7 cm. Furthermore, the Novatel GPS position updates are interpolated to the inertial sensor data rate for KPV determina- tion, which is another potential error source for the benchmark data. Hence, 106 4.9. Experimental Results Figure 4.12: Test jump performed by the athlete. [Photograph by Darren Handschuh] 107 4.9. Experimental Results 16 17 18 19 20 21 22 23 24 \u00E2\u0088\u00921 0 1 2 3 4 5 Time [s] R es ul ta nt a cc el er at io n ab [g ] Increased noisy fluctuation Jump region for the Catapult unit Jump region for the athlete Figure 4.13: Resultant acceleration of Catapult unit strapped to the an- tenna pole carried in backpack. the possibility of error in the benchmark KPV value should be considered when it is used for performance evaluation. Figure 4.14 and Figure 4.15 show two examples of jump trajectory deter- mination with the proposed algorithm, with the Catapult GPS solution and with the Novatel receiver GPS solution. The jump start and end epochs for the athlete were found as described in the previous section. It is noticeable in both figures that the trajectory from the Catapult GPS solution has very poor performance in vertical positioning updates. Not being entirely certain about the Catapult algorithm, it is speculated by looking at the data from the Catapult unit that the height profile is derived from both the altimeter and GPS solution. Therefore, the poor accuracy of the altimeter and the single point GPS receiver causes the unsatisfactory trajectory determina- tion with the Catapult GPS solution. On the other hand, the trajectory determined by the proposed algorithm is very similar to the benchmark tra- jectory from the Novatel receiver GPS solution. This is due to the successful integration of the INS with the GPS updates. The fluctuation at the begin- ning of the trajectory for the proposed algorithm is due to the time required for filter convergence. However, there is a distance (mostly in vertical di- rection) between the benchmark trajectory and the trajectory derived from the proposed algorithm. This is because the proposed algorithm depends heavily on the INS for vertical position estimates and, as seen earlier, the up- 108 4.9. Experimental Results dates for vertical positioning from the Catapult height profile are also very poor in accuracy. Therefore, the noise in low cost MEMS sensors causes a drift from the benchmark in the trajectory determined by the proposed algorithm, mostly in the vertical direction. However, this should not cause considerable error as the interest is focused on relative measurements for a brief period of time. 0 5 10 15 20 \u00E2\u0088\u009220 \u00E2\u0088\u009215 \u00E2\u0088\u009210 \u00E2\u0088\u00925 0 5 \u00E2\u0088\u00922 \u00E2\u0088\u00921.5 \u00E2\u0088\u00921 \u00E2\u0088\u00920.5 0 0.5 North [m]East [m] U pw ar d [m ] Jump start Jump end Novatel GPS trajectory Catapult GPS trajectory Proposed algorithm trajectory Figure 4.14: First example of jump trajectory determination with the pro- posed algorithm, with the Catapult GPS solution and with the Novatel receiver GPS solution. Tables 4.4, 4.5 and 4.6 lists KPV values (horizontal distance, height, and drop respectively) for all the jumps along with the benchmark values deduced from Novatel trajectory. For the calculation of average errors, error in KPV for each jump was first determined and then the errors for all the jumps were averaged. Similarly, the percentage errors were also calculated. To show the advantage of the augmented error observation, provided by the error compensation block in the GPS/INS integration, two results have been shown for the proposed algorithm. One set of KPVs is calculated by the application of the proposed algorithm excluding the augmenting sensor error observations in the EKF, whereas, the second set of KPVs is calculated 109 4.9. Experimental Results 0 10 20 30 \u00E2\u0088\u009215 \u00E2\u0088\u009210 \u00E2\u0088\u00925 0 5 \u00E2\u0088\u00922.5 \u00E2\u0088\u00922 \u00E2\u0088\u00921.5 \u00E2\u0088\u00921 \u00E2\u0088\u00920.5 0 U pw ar d [m ] North [m] East [m] Jump start Jump end Catapult GPS trajectory Novatel GPS trajectory Proposed algorithm trajectory Figure 4.15: Second example of jump trajectory determination with the proposed algorithm, with the Catapult GPS solution and with the Novatel receiver GPS solution. by the complete application of the proposed algorithm, i.e. including the augmenting sensor error observations. As shown Table 4.4, the proposed algorithm provides minor improve- ments in the horizontal distance calculation. The proposed algorithm pro- vides 1.99 % and 1.71 % improvement over the Catapult GPS solution. The percentage error improvements are calculated with respect to the Catapult unit errors in KPV determination. The reason behind this slight improve- ment in horizontal distance measurement is because horizontal distance mea- surement depends greatly on the initial velocity of the body frame. With the inertial sensors, it is not possible to derive the initial velocity at the start of the trajectory determination period. Hence, the velocity from the Catapult GPS receiver is used for the proposed algorithm. Therefore, the error in the Catapult GPS solution also propagates through the proposed algorithm for horizontal distance measurements. It is postulated that the Catapult GPS solution provides velocity readings lower than the true ve- locity. Therefore, the horizontal distance measurements from the proposed 110 4.9. Experimental Results algorithm are generally lower in magnitude than the benchmark. The error in horizontal distance by the Catapult GPS solution is due to an internal mechanism error of the unit, which is unknown. On the other hand, the vertical velocity of the athlete is generally very small at the start of the trajectory determination period and, therefore, er- ror in the Catapult velocity has little effect on the proposed algorithm. As a result, great improvements are reflected in the height and drop calcula- tions for the proposed algorithm as shown in Table 4.5 and Table 4.6. In height, the complete application of the proposed algorithm provides 60.51 % improvement over the Catapult GPS solution, whereas the improvement is even greater (88.73 %) for drop determination. In Table 4.5, it is noticeable that there is no jump height for jumps no. 7, 10 and 12. This is due to the characteristic of the jump itself where the athlete achieved no elevation from the jump start point. For the application of the proposed algorithm with- out the sensor error observations in the EKF, the improvement for height and drop relative to the Catapult GPS solution are 55.70 % and 85.40 % accordingly. The smaller improvement in height and drop determination for the proposed algorithm with no sensor error observation provides the ratio- nale behind augmenting the observation vector in the EKF for GPS/INS integration with the estimated sensor errors from the error compensation block. Of course, the addition of the augmented measurements will increase the computation load. However, this is the price to pay for the increased accuracy in KPV determination. To show the effectiveness of the proposed novel error compensation scheme, GPS/INS integration was also done but excluding the entire er- ror compensation block depicted in Figure 4.6 as opposed to just excluding the use of the augmented observation vector in the EKF, the results of which are shown in Tables 4.4, 4.5 and 4.6. However, the error estimates from the EKF are still fed back to correct the sensor readings and INS mechanization parameters. The bar plots shown in Figure 4.16 and Figure 4.17 summa- rize the error comparison in KPV determination produced by the catapult GPS solution, GPS/INS integration without the error compensation block, proposed algorithm without the augmented observation in EKF, and the complete proposed algorithm. Each block is colored differently to mark the error contribution of each KPV calculation to the total error of all the KPV determinations. The height of each block represents the total error in the jump horizontal distance, height and drop calculation for the applied method to the particular jump. It is evident from the bar plot that the error compensation block has a major impact on the accuracy of the INS mechanization as well as the GPS/INS integration. 111 4.9. Experimental Results Table 4.4: Comparison of jump horizontal distance determination perfor- mance between the proposed algorithm and Catapult GPS solution. Horizontal Distance [cm] Jump no. Novatel GPS Catapult GPS Proposed algo.(without augmented observation) Proposed algo.(with augmented observation) 1 206.25 189.46 189.40 189.37 2 215.81 198.68 198.26 198.24 3 339.44 332.82 333.17 333.17 4 290.47 273.39 274.07 274.13 5 251.52 200.45 201.39 201.27 6 313.86 296.44 296.92 296.74 7 197.96 180.57 180.62 180.67 8 231.89 225.01 225.19 225.21 9 237.45 224.27 225.07 225.07 10 206.10 200.08 200.09 199.95 11 289.50 278.61 279.01 278.79 12 238.07 221.39 221.79 221.62 13 333.62 332.79 333.14 333.16 14 262.59 256.24 256.32 256.32 Average error in cm [error%] 14.59 [5.65%] 14.30 [5.53%] 14.34 [5.55%] 112 4.9. Experimental Results Table 4.5: Comparison of jump height determination performance between the proposed algorithm and Catapult GPS solution. Height [cm] Jump no. Novatel GPS Catapult GPS Proposed algo.(without augmented observation) Proposed algo.(with augmented observation) 1 0.66 0.00 1.98 1.47 2 7.13 0.00 3.37 3.37 3 13.85 0.30 10.98 11.95 4 0.21 0.00 0.00 0.32 5 2.01 1.50 4.86 4.21 6 8.02 0.00 7.84 7.42 7 0.00 0.00 0.00 0.00 8 8.80 0.00 4.14 4.18 9 6.54 0.00 7.95 7.68 10 0.00 0.00 0.00 0.00 11 2.80 0.00 3.92 3.88 12 0.00 0.00 0.09 0.10 13 3.33 0.00 8.56 8.06 14 3.81 0.00 4.57 4.53 Average error in cm [error%] 3.95 [96.85%] 1.75 [42.85%] 1.56 [38.21%] 113 4.9. Experimental Results Table 4.6: Comparison of jump drop determination performance between the proposed algorithm and Catapult GPS solution. Drop [cm] Jump no. Novatel GPS Catapult GPS Proposed algo.(without augmented observation) Proposed algo.(with augmented observation) 1 55.72 0.70 43.14 54.61 2 30.87 1.90 60.95 59.58 3 75.49 0.00 70.68 76.20 4 106.23 60.50 110.14 107.60 5 91.53 0.00 74.67 73.98 6 83.94 7.20 82.38 82.06 7 49.03 3.20 49.71 55.09 8 62.51 4.90 78.64 74.90 9 46.19 5.00 48.87 51.12 10 64.54 6.20 63.73 64.20 11 75.68 16.70 68.68 72.54 12 90.27 9.50 78.28 78.23 13 94.10 38.00 98.13 94.18 14 70.02 8.40 78.68 73.73 Average error in cm [error%] 59.57 [83.72%] 8.70 [12.22%] 6.71 [9.43%] 114 4.9. Experimental Results 1 2 3 4 5 6 7 0 50 100 150 200 250 300 Jump number To ta l e rro r i n K PV s [ cm ] Horizontal distance error Height error Drop error 1st bar: Catapalt GPS solution 2nd bar: GPS/INS integration excluding proposed Error Compensation Block 3rd bar: Proposed algorithm without augmented observation 4th bar: Proposed algorithm with augmented observation Figure 4.16: Error comparison in KPV determination among Catapult GPS solution, GPS/INS integration excluding the error compensation block and proposed algorithm [Jump no. 1-7]. 8 9 10 11 12 13 14 0 20 40 60 80 100 120 140 160 180 200 Jump number To ta l e rro r i n K PV s [ cm ] Horizontal distance error Height error Drop error 1st bar: Catapalt GPS solution 2nd bar: GPS/INS integration excluding proposed Error Compensation Block 3rd bar: Proposed algorithm without augmented observation 4th bar: Proposed algorithm with augmented observation Figure 4.17: Error comparison in KPV determination among Catapult GPS solution, GPS/INS integration excluding the error compensation block and proposed algorithm [Jump no. 8-14]. 115 4.10. Summary 4.10 Summary In this chapter, a jump trajectory determination scheme was presented using inertial sensors and GPS updates. A novel sensor fusion scheme was proposed to mitigate the error in the inertial sensor readings. The sensor fusion scheme also provides leverage to augment the observation vector for GPS/INS integration. Experimental results have shown that the proposed algorithm calculates the jump horizontal distance, height and drop with an average accuracy of 14.34 cm, 1.56 cm and 6.71 cm respectively. Even though Recon has provided no accuracy requirement for these KPVs, the proposed algorithm accuracy should be acceptable for recreational users for whom the magnitude of jump horizontal distance, height and drop are typically small but higher accuracy is not required. For training purposes of serious athletes, the jump will typically be larger and, therefore, errors in the range of 2 to 14 cm, as demonstrated here, will be a small percentage of the true jump horizontal distance, height and drop. 116 Chapter 5 Conclusions and Recommendations The objective of this research was to develop algorithms for Recon action sports goggles to determine the Key Performance Variables (KPVs) relevant to sport maneuvers. Dedicated algorithms were developed to determine the KPVs taking into account the onboard resources available in the goggles. In this chapter, the research objectives and the performance of the proposed algorithms will be highlighted. Then, the key findings of the research will be summarized. Finally, future work will be discussed and recommendations will be made on the expansion of the current work. 5.1 Objectives and Algorithm Performance Objective The objective of this research was to develop online algorithms to de- tect ski or snowboard jumps and to determine four KPVs associated with the jump, namely Air Time (AT), horizontal distance, height and drop. The algorithm should provide real-time feed back and should be developed based on the available low cost MEMS inertial sensors and single point GPS receiver contained in the Recon goggles. Recon provided an AT accuracy requirement of \u00C2\u00B10.1 s. Although Recon set no specific accuracy requirement for jump horizontal distance, height and drop, the algorithm is expected to produce results with sufficient accuracy to meet user\u00E2\u0080\u0099s expectations. For instance, for a recreational user, an accuracy of 10 \u00E2\u0088\u0092 20 cm should be suf- ficient. A professional user would likely desire better accuracy. However, a professional user would typically perform larger jumps than a recreational user and, therefore, errors no more than 10 \u00E2\u0088\u0092 20 cm will be insignificant compared to the size of the jump itself. 117 5.1. Objectives and Algorithm Performance Jump Detection For jump detection, a novel algorithm was developed using Windowed Mean Canceled Multiplication (WMCM) and Preceding and Following Ac- celeration Difference (PFAD), two novel methods developed by the author. Experimental results showed that, the proposed jump detection algorithm successfully detects 92% of the jumps performed by a snowboarder whereas the current Recon jump detection algorithm and the commercial Ripxx unit accurately detect only 60% and 43.50% of the jumps respectively. The false wake up call rate of the proposed jump detection algorithm was found to be 0% for the experimental jumps. On the other hand, the current Recon jump detection algorithm exhibits a false wake up call rate of 45.60%. Considering the number of wrong detections, i.e. the combined number of true jumps un- detected and false jumps detected, the Recon and Ripxx algorithms possess an error rate of 44% and 52% for all jump detections respectively. Unlike these two, the proposed algorithm has a wrong detection rate of only 8%. Air Time Determination To determine AT, a threshold independent probabilistic method using Multiple Attribute Decision Making was developed. The proposed algorithm exhibits an average error of 0.033 s (4.8%), which is well within the accuracy requirement of Recon. On the other hand, the current Recon and Ripxx AT determination algorithms have respective errors of 0.111 s (8.4%) and 0.538 s (39.7%) in AT determination. Jump Horizontal Distance, Height and Drop A GPS/INS integrated algorithm using a novel sensor error correction scheme and augmented measurement vector was developed for the deter- mination of jump horizontal distance, height and drop. From experimental results it was found that the proposed algorithm has an error of 14.34 cm (5.55%), 1.56 cm (38.21%) and 6.71 cm (9.43%) in the determination of jump horizontal distance, height and drop respectively, whereas, the Catapault unit GPS solution has an error of 14.59 cm (5.65%), 3.95 cm (96.85%) and 59.57 cm (83.72%) accordingly. The proposed algorithm accuracy should be sufficient for the recreational user, where magnitude of the KPVs is small, and higher accuracy is not strictly required. On the other hand, for the train- ing purposes of serious athletes, the jumps are typically large and, therefore, errors in the range of 2 to 14 cm will be insignificant with respect to the true 118 5.2. Methods and Key Findings jump horizontal distance, height and drop. Hence, the proposed algorithm should also meet the requirements of professional athletes. 5.2 Methods and Key Findings Among the available inertial sensors, namely accelerometer, gyroscope and magnetometer, the accelerometer is the most useful for jump detection as the change in acceleration during a jump is most consistent. Even though gyroscope readings change significantly as a result of the head movements associated with a jump, these changes are not consistent. The magnetometer is not usable for jump detection as this sensor is not related to jump dy- namics. Therefore, the proposed tools for jump detection, i.e. the WMCM and PFAD methods, only exploit the accelerometer signals. Use of redun- dant accelerometers will increase the accuracy and reliability of the WMCM method. This is because the WMCM method extracts jump information from each individual accelerometer oriented in different directions. With an increased number of accelerometer sensors oriented in other directions, the WMCM method will have more information about true jump occurrence. Therefore, jump detection performance will increase with the addition of more accelerometers. However, the addition of more sensors will increase the power consumption and cost of the device. Four characteristic points in resultant acceleration relevant to jump oc- currence are detected for AT determination. For detecting two of the points, the closest peak method is proposed. For detecting the other two charac- teristic points, two methods are developed, one being threshold based and the other being a threshold independent probabilistic method using MADM. The threshold dependent method incurs a loss of generality of the algorithm. Therefore, the probabilistic method is promoted even though it causes a mi- nor increase in the computational burden. M-TOPSIS, a modified version of the TOPSIS method, is used for the implementation of MADM. From the detection of these characteristic points, Sensor Air Time (SAT) is derived and related to the visual AT later on. It is of utmost importance to detect precisely the four characteristic points in order to determine the AT with high accuracy. Depending on the accuracy requirement and the definition of jump take-off and landing, detection of two characteristic points may also suffice. A GPS/INS integrated algorithm was developed for jump horizontal dis- tance, height and drop measurement and it incorporates a novel sensor error compensation scheme. Three Linear Kalman Filters (LKFs) were used to 119 5.3. Recommendations exploit sensor fusion to compensate for the errors in the accelerometer and gyroscope. It was found to be critically important to vary the Kalman filter parameters over the relevant jump period as the athlete\u00E2\u0080\u0099s dynamics signif- icantly change before, during, and after a jump. As results have shown, the error compensation scheme is very important for the detection of the KPVs. The Extended Kalman Filter (EKF) used for GPS/INS integra- tion is executed with an observation vector augmented with the sensor error observations calculated by the error compensation block, which results in better estimation of the sensor errors. Of course, the improvement in the accuracy of the state estimation comes at the cost of an increased compu- tational load due to the expansion in vector size. As mentioned earlier, the proposed algorithm with the augmented observation vector has an error of 14.34 cm (5.55%), 1.56 cm (38.21%) and 6.71 cm (9.43%) in the determina- tion of jump horizontal distance, height and drop respectively, whereas, the proposed algorithm without the augmented observation vector has an error of 14.30 cm (5.53%), 1.75 cm (42.85%) and 8.70 cm (12.22%) accordingly. 5.3 Recommendations During the course of this research, new challenges important and worthy of investigation were encountered. Quaternion. The body frame to navigation frame rotation matrix im- plemented with Euler angles may become singular at times, which can be addressed by quaternion algebra. With the application of quaternion algebra the rotation matrix calculation will also be decreased as only four elements are calculated with quaternion algebra as opposed to nine elements with Euler angles. Kalman Filter Delay. The delay in the Kalman filter should be inves- tigated for the proposed algorithms. If the error due to this phenomenon is significant, recovery measures should be undertaken for compensation of the filter delay. Tight Coupling. GPS/INS integration tight coupling should also be in- vestigated and the performance should be analyzed. The increase in compu- tational complexity and the gain in KPV determination accuracy should be evaluated in case of the adoption of a tightly coupled GPS/INS integrated system. 120 5.3. Recommendations Algorithm Application. Even though the algorithms are developed spe- cially for ski or snowboard jumping, the developed techniques might have application aspects in a variety of fields with modifications and different tuning. As seen in Chapter 4, the athlete\u00E2\u0080\u0099s jump dynamics in skiing/snow- boarding and mountain biking are very similar. Hence, this research can also be applied to sports where maneuvering jumps plays a key role. For instance, these developed algorithms can be directly implemented in BMX biking where similar KPVs are used to judge the aerial stunts executed by professional athletes. Therefore, the scope of the developed algorithms should be widely explored, not only in outdoor sports, but also in other potential sectors such as video gaming, bio-medical engineering, etc. 121 Bibliography [1] R. Instruments, \u00E2\u0080\u009Chttp://www.reconinstruments.com/,\u00E2\u0080\u009D 2011. [2] J. W. Harding, J. W. Small, and D. A. James, \u00E2\u0080\u009CFeature extraction of performance variables in elite half-pipe snowboarding using body mounted inertial sensors,\u00E2\u0080\u009D in BioMEMS and Nanotechnology III, Proc. of SPIE, vol. 6799, 2007, pp. 679 917\u00E2\u0080\u00931\u00E2\u0080\u0093679 917\u00E2\u0080\u009312. [3] A. Waegli, \u00E2\u0080\u009CTrajectory determination and analysis in sports by satellite and inertial navigation,\u00E2\u0080\u009D Ph.D. dissertation, Department of Geomatics Engineering, University of Calgary, 2009. [4] Ripxx, \u00E2\u0080\u009Chttp://ripxx.com/,\u00E2\u0080\u009D 2011. [5] SHADOWBOX, \u00E2\u0080\u009Chttp://www.shadowboxlive.com/,\u00E2\u0080\u009D 2011. [6] N. El-Sheimy and X. Niu, \u00E2\u0080\u009CThe promise of MEMS to the navigation community,\u00E2\u0080\u009D InsideGNSS, March/April 2007. [7] S. Godha and M. Cannon, \u00E2\u0080\u009CGPS/MEMS INS integrated system for navigation in urban areas,\u00E2\u0080\u009D GPS Solutions, vol. 11, no. 3, pp. 193\u00E2\u0080\u0093203, 2007. [8] Y. Li, A. Dempster, B. Li, J. Wang, and C. Rizos, \u00E2\u0080\u009CA low-cost attitude heading reference system by combination of GPS and magnetometers and MEMS inertial sensors for mobile applications,\u00E2\u0080\u009D Journal of Global Positioning Systems, vol. 5, no. 1-2, pp. 90\u00E2\u0080\u009397, 2006. [9] Y. Li, J. Wang, C. Rizos, P. Mumford, and W. Ding, \u00E2\u0080\u009CLow-cost tightly coupled GPS/INS integration based on a nonlinear kalman filtering design,\u00E2\u0080\u009D in Proceedings of the ION NTM, 2006. [10] J. W. Harding, K. Toohey, D. T. Martin, C. G. Mackintosh, A. Lindh, and D. A. James, \u00E2\u0080\u009CAutomated inertial feedback for half-pipe snow- board competition and the community perception,\u00E2\u0080\u009D in The Impact of Technology on Sport II, 2007, pp. 845\u00E2\u0080\u0093850. 122 Bibliography [11] H. Doki, T. Yamada, C. Nagai, and M. Horaki, \u00E2\u0080\u009CDevelopment of a measurement system for snowboarding turn analysis,\u00E2\u0080\u009D In: Subic A, Ujihashi S, eds. The Impact of Technology in Sport, pp. 324\u00E2\u0080\u0093325, 2005. [12] NDI, \u00E2\u0080\u009Chttp://www.ndigital.com/medical/technology-em.php,\u00E2\u0080\u009D 2011. [13] L. Ren, Y. Zhang, Y. Wang, and Z. Sun, \u00E2\u0080\u009CKinematics of the ankle joint complex in snowboarding,\u00E2\u0080\u009D Applied Mathematics Research eX- press, 2007. [14] L. Bianchi, N. Petrone, and M. Marchiori, \u00E2\u0080\u009CA dynamometric platform for load data acquisition in snowboarding: Design and field analysis,\u00E2\u0080\u009D The Engineering of Sport 5, M. Hubbard, R. D. Mehta, J. M. Pallis, J.M., vol. 2, pp. 187\u00E2\u0080\u0093193, 2004. [15] J. W. Harding, C. G. Mackintosh, A. G. Hahn, and D. A. James, \u00E2\u0080\u009CClas- sification of aerial acrobatics in elite half-pipe snowboarding using body mounted inertial sensors,\u00E2\u0080\u009D in Proceedings of 7th ISEA CONFERENCE, June 2008. [16] J. W. Harding, K. Toohey, D. T. Martin, A. G. Hahn, and D. A. James, \u00E2\u0080\u009CTechnology and half-pipe snowboard competition-insight from elite- level judges,\u00E2\u0080\u009D in Proceedings of 7th ISEA Conference. [17] S. Tsutomu, I. Kazutoyo, T. Kazuhiko, H. Hiroshi, M. Shosuke, K. Takayuki, and O. Manabu, \u00E2\u0080\u009CDevelopment of acceleration measur- ing system with a 3-D gyroscope sensor during the flight phase in ski jumping,\u00E2\u0080\u009D in 12th Annual Congress of the ECSS, July 2007. [18] B. P. Fleantov, D. D. M. Darcy, and S. C. A. Vock, \u00E2\u0080\u009CApparatus and methods for determining loft time and speed,\u00E2\u0080\u009D Patent 5 636 146, June, 1997. [19] F. Castanie and A. Denjean, \u00E2\u0080\u009CMean value jump detection using wavelet decomposition,\u00E2\u0080\u009D in Time-Frequency and Time-Scale Analysis, 1992., Proceedings of the IEEE-SP International Symposium, oct 1992, pp. 181\u00E2\u0080\u0093184. [20] M. Crouse, R. Nowak, and R. Baraniuk, \u00E2\u0080\u009CWavelet-based statistical sig- nal processing using hidden markov models,\u00E2\u0080\u009D Signal Processing, IEEE Transactionson, vol. 46, no. 4, pp. 886\u00E2\u0080\u0093902, apr 1998. 123 Bibliography [21] H. Teng and Y. Qi, \u00E2\u0080\u009CApplication of wavelet technique to freeway incident detection,\u00E2\u0080\u009D Transportation Research Part C: Emerging Tech- nologies, vol. 11, no. 3-4, pp. 289\u00E2\u0080\u0093308, 2003, traffic Detection and Es- timation. [Online]. Available: http://www.sciencedirect.com/science/ article/B6VGJ-492707R-2/2/ced7d3a53aa7787eaf87c82c01bc9184 [22] M. Raimondo and N. Tajvidi, \u00E2\u0080\u009CA peaks over threshold model for change-point detection by wavelets,\u00E2\u0080\u009D Statistica Sinica, vol. 14, pp. 395\u00E2\u0080\u0093 412, 2004. [23] DARTFISH, \u00E2\u0080\u009Chttp://www.dartfish.com/en/index.htm,\u00E2\u0080\u009D 2011. [24] C. Nicolas and M. Tavernier, \u00E2\u0080\u009CContribution to the analysis of kinemat- ics and energy and skating cross country skiing. application to the com- parative study of three nordic disciplines: cross country skiing, biathlon and nordic combined. comparative analysis of two 3d motion analysis system using a stick and a human model volumic model and application to human case studies in cross-country skiing,\u00E2\u0080\u009D Ph.D. dissertation. [25] CODAMOTION, \u00E2\u0080\u009Chttp://www.codamotion.com/,\u00E2\u0080\u009D 2011. [26] VICON, \u00E2\u0080\u009Chttp://www.vicon.com/,\u00E2\u0080\u009D 2011. [27] organicmotion, \u00E2\u0080\u009Chttp://www.organicmotion.com/,\u00E2\u0080\u009D 2011. [28] POLAR, \u00E2\u0080\u009Chttp://www.polarusa.com/us-en/,\u00E2\u0080\u009D 2011. [29] GARMIN, \u00E2\u0080\u009Chttp://www.garmin.com/,\u00E2\u0080\u009D 2011. [30] FRWD, \u00E2\u0080\u009Chttp://www.frwd.fi/,\u00E2\u0080\u009D 2011. [31] CATAPULT, \u00E2\u0080\u009Chttp://www.catapultinnovations.com/index.html,\u00E2\u0080\u009D 2011. [32] Sportvision, \u00E2\u0080\u009Chttp://www.sportvision.com/,\u00E2\u0080\u009D 2011. [33] TracTrac, \u00E2\u0080\u009Chttp://www.tractrac.com/,\u00E2\u0080\u009D 2011. [34] J. Skaloud and P. Limpach, \u00E2\u0080\u009CSynergy of CP-DGPS, accelerometry and magnetic sensors for precise trajectography in ski racing,\u00E2\u0080\u009D in In Pro- ceedings of the ION GPS/GNSS 2003. [35] S. Ducret, P. Ribot, R. Vargiolu, J. Lawrence, and A. Midol, \u00E2\u0080\u009CAnalysis of downhill ski performance using GPS and grounding force recording,\u00E2\u0080\u009D Science and Skiing III. 124 Bibliography [36] K. Zhang, R. Deakin, R. Grenfell, Y. Li, J. Zhang, W. N. Cameron, and D. M. Silcock, \u00E2\u0080\u009CGNSS for sports-sailing and rowing perspectives,\u00E2\u0080\u009D Journal of Global Positioning Systems, vol. 3, no. 1-2. [37] K. Zhang, R. Grenfell, R. Deakin, Y. Li, Z. Jason, A. Hahn, C. Gore, and T. Rice, \u00E2\u0080\u009CTowards a low-cost, high output rate, real-time GPS rowing coaching and training system,\u00E2\u0080\u009D in In Proceedings of the ION GPS/GNSS 2003. [38] A. Waegli and J. Skaloud, \u00E2\u0080\u009COptimization of two GPS/MEMS-IMU integration strategies with application to sports,\u00E2\u0080\u009D GPS Solutions, vol. 13, pp. 315\u00E2\u0080\u0093326, 2009, 10.1007/s10291-009-0124-5. [Online]. Available: http://dx.doi.org/10.1007/s10291-009-0124-5 [39] J. Skaloud, \u00E2\u0080\u009COptimizing georeferencing of airborne survey systems by INS/DGPS,\u00E2\u0080\u009D no. UCGE Report 20126. [40] I. Skog and P. Hndel, \u00E2\u0080\u009CA low-cost GPS aided inertial navigation system for vehicle applications,\u00E2\u0080\u009D in Proc. EUSIPCO. [41] S. Godha, \u00E2\u0080\u009CPerformance evaluation of low cost MEMS-based IMU in- tegrated with GPS for land vehicle navigation application,\u00E2\u0080\u009D Master\u00E2\u0080\u0099s thesis, Department of Geomatics Engineering, University of Calgary, February 2006. [42] H. Qi and J. Moore, \u00E2\u0080\u009CDirect kalman filtering approach for GPS/INS integration,\u00E2\u0080\u009D Aerospace and Electronic Systems, IEEE Transactions on, vol. 38, no. 2, pp. 687 \u00E2\u0080\u0093693, apr 2002. [43] I. Oppermann, L. Stoica, A. Rabbachin, Z. Shelby, and J. Haapola, \u00E2\u0080\u009CUwb wireless sensor networks: UWEN-a practical example,\u00E2\u0080\u009D Commu- nications Magazine, IEEE, vol. 42, no. 12, pp. S27\u00E2\u0080\u0093S32, dec. 2004. [44] G. Opshaug and P. Enge, \u00E2\u0080\u009CGPS and UWB for indoor navigation,\u00E2\u0080\u009D in Proceedings of the ION GPS/GNSS, 2001. [45] T. Kitasuka, T. Nakanishi, and A. Fukuda, \u00E2\u0080\u009CDesign of WiPS: WLAN- based indoor positioning system,\u00E2\u0080\u009D Korea Multimedia Society, vol. 7, no. 4, pp. 15\u00E2\u0080\u009329, 2003. [46] R. T. Ogden, Essential Wavelets for Statistical Applications and Data Analysis. Birkhauser Boston, 1997. 125 Bibliography [47] D. C. Montgomery, Design and Analysis of Experiments, 5th ed. John Wiley & Sons (Asia) Pte. Ltd., 2006. [48] K. P. Yoon and C. L. Hwang, Multiple Attribute decision making: An Introduction. Sage University Paper series on Quantitative Applica- tions in the Social Sciences, 07-104. Thousand Oaks, CA: Sage, 1995. [49] L. Ren, Y. Zhang, Y. Wang, and Z. Sun, \u00E2\u0080\u009CComparative analysis of a novel M-TOPSIS method and TOPSIS,\u00E2\u0080\u009D Applied Mathematics Research eXpress, 2007. [50] R. V. Wong, \u00E2\u0080\u009CDevelopment of a RLG strapdown inertial survey sys- tem,\u00E2\u0080\u009D Ph.D. dissertation. [51] J. Diebel, \u00E2\u0080\u009CRepresenting attitude: Euler angles, unit quaternions, and rotation vectors,\u00E2\u0080\u009D Stanford University, Stanford, California 94301-9010, Tech. Rep., October 2006. [52] D. H. Titterton and J. L. Weston, Strapdown Inertial Navigation Tech- nology, 2nd ed. The Institution of Electrical Engineers and The Amer- ican Institute of Aeronautics and Astronautics, 2004. [53] J. Bortz, \u00E2\u0080\u009CA new mathematical formulation for strapdown inertial navi- gation,\u00E2\u0080\u009D Aerospace and Electronic Systems, IEEE Transactions on, vol. AES-7, no. 1, pp. 61 \u00E2\u0080\u009366, jan. 1971. [54] M. Ignagni, \u00E2\u0080\u009COn the orientation vector differential equation in strap- down inertial systems,\u00E2\u0080\u009D Aerospace and Electronic Systems, IEEE Transactions on, vol. 30, no. 4, pp. 1076 \u00E2\u0080\u00931081, oct 1994. [55] Kionix, \u00E2\u0080\u009CTilt sensing with kionix MEMS accelerometers,\u00E2\u0080\u009D Tech. Rep. [56] C. Konvalin, \u00E2\u0080\u009CCalculating heading, elevation and bank angle,\u00E2\u0080\u009D MEM- SENSE, Tech. Rep. MTD-0801. [57] K. J. Walchko and D. P. A. C. Mason, \u00E2\u0080\u009CInertial navigation,\u00E2\u0080\u009D in Florida Conference on Recent Advances in Robotics, 2002. [58] S. Nassar, \u00E2\u0080\u009CImproving the inertial navigation system (INS) error model for INS and INS/DGPS applications,\u00E2\u0080\u009D Ph.D. dissertation, Department of Geomatics Engineering, University of Calgary, 2003. [59] H. Luinge and P. Veltink, \u00E2\u0080\u009CMeasuring orientation of human body segments using miniature gyroscopes and accelerometers,\u00E2\u0080\u009D 126 Bibliography Medical and Biological Engineering and Computing, vol. 43, pp. 273\u00E2\u0080\u0093282, 2005, 10.1007/BF02345966. [Online]. Available: http: //dx.doi.org/10.1007/BF02345966 [60] E. Nebot and H. Durrant-Whyte, \u00E2\u0080\u009CInitial calibration and alignment of low-cost inertial navigation units for land vehicle applications,\u00E2\u0080\u009D Journal of Robotic Systems, vol. 2, no. 16, pp. 82\u00E2\u0080\u009392, 1999. [61] S. M. Kay, Fundamentals of Statistical Signal Processing: Estimation Theory. Prentice Hall PTR, Upper Saddle River, 1993. [62] N. El-Sheimy, \u00E2\u0080\u009CThe potential of partial IMUs for land vehicle naviga- tion,\u00E2\u0080\u009D InsideGNSS, Spring 2008. [63] E.-H. Shin, \u00E2\u0080\u009CAccuracy improvement of low cost INS/GPS for land appli- cations,\u00E2\u0080\u009D Master\u00E2\u0080\u0099s thesis, Department of Geomatics Engineering, Uni- versity of Calgary, Canada, December 2001. [64] P. G. Savage, \u00E2\u0080\u009CStrapdown analytics,\u00E2\u0080\u009D Strapdown Associates, Inc, vol. 1, 2000. [65] R. M. Rogers, Applied Mathematics in Integrated Navigation Systems. AIAA Education Series, 2000. [66] N. El-Sheimy, \u00E2\u0080\u009CInertial techniques and INS/DGPS integration,\u00E2\u0080\u009D ENGO 623 Course Notes, 2004. [67] N. R. Canada, \u00E2\u0080\u009Chttp://www.geod.nrcan.gc.ca,\u00E2\u0080\u009D 2011. [68] M. Ehrgott, Multicriteria Optimization, 2nd ed. Springer-Berlin Hei- delberg New York. [69] O. Bandte, \u00E2\u0080\u009CA probabilistic multi-criteria decision making technique for conceptual and preliminary aerospace systems design,\u00E2\u0080\u009D Ph.D. dis- sertation, Georgia Institute of Technology, September 2000. [70] A. Davey and D. Olson, \u00E2\u0080\u009CMultiple criteria decision making models in group decision support,\u00E2\u0080\u009D Group Decision and Negotiation, vol. 7, pp. 55\u00E2\u0080\u009375, 1998. 127 Appendix Appendix A Multiple Attribute Decision Making In general, a Multiple Attribute Decision Making (MADM) problem can be described as [68, 69] max x1,x2,..,xp {F1(x1, x2, .., xp), F2(x1, x2, .., xp) .., FP (x1, x2, .., xp)} (A-1) subject to a set of constraints gr(x1, x2, .., xp) \u00E2\u0089\u00A4 0, r = 1, 2, .., L (A-2) where Fi, i = 1, 2, .., P are the criteria and xi, i = 1, 2, .., n are the de- sign variables which the criteria and constraints depend on. Any set X = {x1, x2, .., xp|g(x1, x2, .., xp) \u00E2\u0089\u00A4 0} is the set of feasible design variables which complies with the design constraints. For each point in X there is an asso- ciated set {F1, F2, .., FP } such that X is mapped into a set S = {F1, F2, .., FP |x1, x2, .., xp \u00E2\u0088\u0088 X} (A-3) in the criteria space. Further, {x\u00E2\u0088\u00971, x\u00E2\u0088\u00972, .., x\u00E2\u0088\u0097p} is called the optimal solution if and only if {x\u00E2\u0088\u00971, x\u00E2\u0088\u00972, .., x\u00E2\u0088\u0097p} \u00E2\u0088\u0088 X and Fi{x\u00E2\u0088\u00971, x\u00E2\u0088\u00972, .., x\u00E2\u0088\u0097p} \u00E2\u0089\u00A5 Fi{x1, x2, .., xp}, i = 1, 2, .., P . TOPSIS The Technique for Order Preference by Similarity to Ideal Solution (TOP- SIS) is one of the most popular compromising methods among the compen- satory techniques, utilizing preference information provided in the form of weights wci , i = 1, 2, .., P for each criterion. TOPSIS defines an index called similarity (or relative closeness) to the positive-ideal solution by combin- ing the proximity to the positive-ideal solution and remoteness from the negative-ideal solution [48]. The positive-ideal solution simultaneously opti- mizes each objective. The positive-ideal solution to a multi-criteria problem 128 A. Multiple Attribute Decision Making is generally not feasible, but can serve as a standard to the alternatives [70]. Assuming that the greatest utility can be achieved from the positive-ideal solution, closeness to the positive-ideal solution can be considered analogous to maximizing utility. The concept of the decision matrix D is applicable to all other MADM techniques and can be expressed as D = \u00EF\u00A3\u00AE\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00AF\u00EF\u00A3\u00B0 f11 f12 . . . f1i . . . f1P f21 f22 . . . f2i . . . f2P ... ... . . . ... . . . ... fj1 fj2 . . . fji . . . fjP ... ... . . . ... . . . ... fM1 fM2 . . . fMi . . . fMP \u00EF\u00A3\u00B9\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BA\u00EF\u00A3\u00BB . (A-4) D is an M-by-P-matrix with elements fji that indicate the value of the criterion Fi with respect to the alternative Aj . The initial step of most MADM problems is to calculate the normalized rating [48], rji, for each element of D. The vector normalization procedure can be given as rji = fji\u00E2\u0088\u009A M\u00E2\u0088\u0091 j=1 x2ji , j = 1, ...,M ; i = 1, ..., P. (A-5) Hence the weighted normalized ratings are calculated as vji = w c i rji, j = 1, ...,M ; i = 1, ..., P. (A-6) In such a context, the positive ideal solution can be denoted as A\u00E2\u0088\u0097 = {v\u00E2\u0088\u00971, .., v\u00E2\u0088\u0097i , .., v\u00E2\u0088\u0097P } = {(max j vji|i \u00E2\u0088\u0088 I1), (min j vji|i \u00E2\u0088\u0088 I2)|j = 1, ...,M} (A-7) where v\u00E2\u0088\u0097i is the best value for the ith criterion among all available alterna- tives. I1 is the set of benefit attributes and I2 is the set of cost attributes. Benefit attributes are those which are expected to be maximized and cost attributes are those which are expected to be minimized. Consequently, the 129 A. Multiple Attribute Decision Making negative-ideal solution can be denoted as A\u00E2\u0088\u0092 = {v\u00E2\u0088\u00921 , .., v\u00E2\u0088\u0092i , .., v\u00E2\u0088\u0092P } = {(min j vji|i \u00E2\u0088\u0088 I1), (max j vji|i \u00E2\u0088\u0088 I2)|j = 1, ...,M} (A-8) where v\u00E2\u0088\u0092i is the worst value for the ith criterion among all alternatives. The separation of each alternative from the positive ideal solution, A\u00E2\u0088\u0097, can then be given by D\u00E2\u0088\u0097j = \u00E2\u0088\u009A\u00E2\u0088\u009A\u00E2\u0088\u009A\u00E2\u0088\u009A P\u00E2\u0088\u0091 i=1 (vji \u00E2\u0088\u0092 v\u00E2\u0088\u0097i )2, j = 1, ...,M. (A-9) Similarly, the separation from the negative ideal solution, A\u00E2\u0088\u0092, can be given by D\u00E2\u0088\u0092j = \u00E2\u0088\u009A\u00E2\u0088\u009A\u00E2\u0088\u009A\u00E2\u0088\u009A P\u00E2\u0088\u0091 i=1 (vji \u00E2\u0088\u0092 v\u00E2\u0088\u0092i )2, j = 1, ...,M. (A-10) Afterwards, similarities to the positive ideal solution are calculated as C\u00E2\u0088\u0097j = D\u00E2\u0088\u0092j D\u00E2\u0088\u0097j +D \u00E2\u0088\u0092 j , j = 1, ...,M. (A-11) Note that 0 \u00E2\u0089\u00A4 C\u00E2\u0088\u0097j \u00E2\u0089\u00A4 1, where C\u00E2\u0088\u0097j = 0 when Aj = A\u00E2\u0088\u0092, and C\u00E2\u0088\u0097j = 1 when Aj = A \u00E2\u0088\u0097. Finally, all the alternatives are ranked against C\u00E2\u0088\u0097i in descending order and the alternative with the maximum C\u00E2\u0088\u0097i (rank 1) can be chosen as the optimal solution for the system. For a two dimensional case, i.e. for MADM for two criteria, and six alternatives, the TOPSIS method is illustrated in Figure A-1. M-TOPSIS The TOPSIS method prone to \u00E2\u0080\u0098Rank Inversion\u00E2\u0080\u0099 problem which can be resolved by the Modified TOPSIS (M-TOPSIS) method [49]. Rank inversion means the change in relative ranking among any of the two alternatives in the event of exclusion or inclusion of any other alternative. In M-TOPSIS, the D\u00E2\u0088\u0097jD \u00E2\u0088\u0092 j plane is constructed and a new Optimized Ideal Reference Point (OIRP) is constructed. In the D\u00E2\u0088\u0097jD \u00E2\u0088\u0092 j plane, D \u00E2\u0088\u0097 is the x-axis and D\u00E2\u0088\u0092 is the y-axis. (D\u00E2\u0088\u0097j , D \u00E2\u0088\u0092 j ) represents each alternative and 130 A. Multiple Attribute Decision Making A1 A3 A2 A6 A5 A4 D5 + D4 + D4 - D5 - A + A - (Negative Ideal Solution) (Positive Ideal Solution) Attribute 1 A tt ri b u te 2 Similarities to Positive Ideal Solution, C =j * D +j + Dj - Dj - Figure A-1: TOPSIS method for MADM. 131 A. Multiple Attribute Decision Making Optimized Ideal Reference Point A4 A1 A3 A2 D D planej j * - max( D )j * max( D )j - min( D )j - min( D )j * Dj - Dj * Distance from Positive Ideal Solution D is ta n ce f ro m N eg at iv e Id ea l S o lu ti o n D2 OIRP D1 OIRP D3 OIRP A5 A6D6 OIRP D5 OIRP D4 OIRP Figure A-2: M-TOPSIS method for MADM. the OIRP is the point (min(D\u00E2\u0088\u0097j ), (max)(D \u00E2\u0088\u0092 j )). Then the distances from the OIRP to the alternatives are calculated as DOIRPj = \u00E2\u0088\u009A {min(D\u00E2\u0088\u0097j )\u00E2\u0088\u0092D\u00E2\u0088\u0097j}2 + {min(D\u00E2\u0088\u0092j )\u00E2\u0088\u0092D\u00E2\u0088\u0092j }2, j = 1, ...,M. (A-12) Afterwards, DOIRPj is sorted in ascending order and all the alternatives are ranked against its\u00E2\u0080\u0099 closeness to the OIRP. The alternative with minimum DOIRPj (rank 1) is considered as the optimal choice for the MADM problem. The D\u00E2\u0088\u0097jD \u00E2\u0088\u0092 j plane and the distances between OIRP and the alternatives are depicted in Figure A-2. 132"@en . "Thesis/Dissertation"@en . "2011-11"@en . "10.14288/1.0072682"@en . "eng"@en . "Electrical Engineering"@en . "Vancouver : University of British Columbia Library"@en . "University of British Columbia"@en . "Attribution-NonCommercial-NoDerivatives 4.0 International"@en . "http://creativecommons.org/licenses/by-nc-nd/4.0/"@en . "Graduate"@en . "Jump parameter estimation with low cost micro-electro-mechanical system sensors and global positioning system for action sports goggles"@en . "Text"@en . "http://hdl.handle.net/2429/41971"@en .