UNIVERZITET U BEOGRADU ELEKTROTEHNIČKI FAKULTET Milica D. Đurić-Jovičić METODE ANALIZE SIGNALA SA INERCIJALNIH SENZORA ZA ANALIZU HODA PACIJENATA SA OŠTEĆENIM OBRASCEM HODA doktorska disertacija Beograd, 2012. ii UNIVERSITY OF BELGRADE SCHOOL OF ELECTRICAL ENGINEERING Milica D. Djurić-Jovičić INERTIAL SENSORS SIGNAL PROCESSING METHODS FOR GAIT ANALYSIS OF PATIENTS WITH IMPAIRED GAIT PATTERNS Doctoral Dissertation Belgrade, 2012 iii Podaci o mentoru i članovima komisije Mentor: Prof dr Dejan Popović, Univerzitet u Beogradu, Elektrotehnički fakultet Članovi komisije: Prof dr Dejan Popović, Univerzitet u Beogradu, Elektrotehnički fakultet Prof dr Mirjana Popović, Univerzitet u Beogradu, Elektrotehnički fakultet Prof dr Vladimir Kostić, Univerzitet u Beogradu, Medicinski fakultet Prof dr Antonije Đorđević, Univerzitet u Beogradu, Elektrotehnički fakultet Prof dr Branko Kovačević, Univerzitet u Beogradu, Elektrotehnički fakultet Datum: ________________. iv Acknowledgment I owe many people for generously giving me their time and expertise. Foremost, I would like to express my special thanks and appreciation to my supervisor, Prof. Dejan Popović, for giving me the opportunity to work in the research area of biomedical engineering, in particular in gait analysis, and for his encouragement, support, and supervision at all levels. Most of all, I am thankful for giving me the opportunity to do what I love and enjoy, for all the travels, places I’ve visited, colleagues I’ve met, and warm hospitality in all those places because they new I was coming from his group. I warmly thank Prof. Mirjana Popović, my co-advisor, for her collaboration in preparing and supervising the clinical protocols, which enabled us to gather and process all the valuable clinical data, and guidance through the writing processes for our publications. I deeply appreciate her enthusiasm, insightful comments, and helpful advices all along the way. I wish to thank Mr Nenad Jovičić for the design and producing of all the versions of the ambulatory system and the sensor units used in this research. I kindly thank Prof. Vladimir Kostić and his team, Dr Saša Radovanović, Dr Nikola Kresojević, Dr Iva Stanković and Dr Igor Petrović from the Neurology Clinic, Clinical Center of Serbia, as well as Prof. Laszlo Schwirtlich, Dr Aleksandra Dragin, therapist Jelena Milovanović from the Rehabilitation clinic “Dr Miroslav Zotović” for helping with the design of clinical protocols, suggestions for the studies, selection of the patients, assistance during experiments and interpretation of the obtained results. I would also like to thank all the subjects who participated in my studies, for their confidence, patience and collaboration in all our ideas. I also thank Serbian Ministry of Education and Science and Tecnalia Serbia who financed my work on this thesis. Special thanks to Goran Bijelić, director of Tecnalia v Serbia, for supporting me and my research. I would also like to mention our collaborators from Tecnalia Spain, Dr Ulrich Hoffman who was my project manager and who inspected and discussed with me what it inside our gait analysis system, Pierre Baralon who tested our system independently and provided valuable results and feedback for improvement to its current state, Haritz Zabaleta who helped me with the experiments and validation of the developed methods, and Dr Thierry Keller, who supported the project as the head of Biorobotics Department in Tecnalia Spain. I express my gratitude to Dr Strahinja Došen, who taught me to use motion capture system, assisted in the experiments with this system and data processing, and from whom I learned my first serious Matlab “tricks”, but most of all I thank him and Dr Jovana Kojović for welcoming me in Aalborg and making my time spent there such a pleasant memory. I also warmly thank Dr Rita Stagni, Assistant Professor at the Department of Electronics, Computer Science and Systems, University in Bologna, Italy, for organizing the experiments for validation of our developed system and comparison with state-of-the-art commercial systems, for leading the experiments and signal processing from their systems. Most of all, I am grateful for her warm hospitality, and introducing us to many Italian colleagues from the same or similar research area. Methods developed in this thesis wouldn’t be so applicable to non-engineers without user-friendly graphic interface, designed by Đorđe Klisić whom I thank for the patience and devotion to make it as perfect as possible. To my chronologically first “boss” and mentor, Prof. Antonije Đorđević, I am grateful for so many things: for the support and encouragement during my bachelor studies, for forcing me to learn to use many engineering tools, for making me love the “engineering way”, for teaching me not to be satisfied with half-results. I am also thankful for all the projects and papers we did together and for forgiving me when I chose to leave his antennas and dedicate to high heel study, which resulted in such an enjoying diploma work which overgrew to my PhD studies and this thesis as it is. I am enormously happy that, so many years later, we wrote another paper together, this time from my field of interest. It is a great pleasure to see that wonderful paper, signed by my former and present mentors, my husband and me, accepted for one of the best journals in this scientific area. vi The most beautiful and sometimes the most difficult part of my research were not equations but working with patients in clinical environment. As engineers, we were, naturally, not used to seeing all those medical conditions and gait deformities which can happen to people. These experiments and encounters with some patients were sometimes very hard and emotional, especially in the beginning of my work experience. I owe my deepest gratitude to Prof. Dejan Popović and Prof. Mirjana Popović for teaching me how to approach and work with patients, and how to organize the experiments in the most professional manner. And special thanks to my dear collaborator Ivana Milovanović, with whom I performed most of the experiments, and who always made work more enjoyable. Good experiments or observing patient’s recovery made us smile and were rewarding by default. However, there were also some bad days at the clinics. We fought against those days by taking ice-cream upon return from the Clinic, just in front of the Faculty. Besides being wonderful friend and colleague, I also thank her for reading the first few hundred versions of my first journal paper, and preventing me to start the Introduction with “Since Aristotle…”. Speaking of wonderful colleagues, to all my BMIT co-members, Nadica, Miloš, Nebojša, Lana, Milica, Andrej, Vera, Matija and Ivana with whom I worked and traveled and shared many wonderful moments, thank you for all of that. Finally, to my husband Nenad (previously mentioned as the man who made the hardware for this research), I thank for the countless conversations and debates about my research. In many ways this has been a joint venture and a joint achievement. Most of all, thank you for being there for me all these years. Lastly, to the most loyal and biased fan club one can have- my parents Branka and Dragiša, and my brother Pavle: I love you and thank for all the support. MCA vii Metode analize signala sa inercijalnih senzora za analizu hoda pacijenata sa oštećenim obrascem hoda Rezime - Analiza hoda je postala široko rasprostranjen klinički alat koji se koristi za objektivnu evaluaciju obrasca hoda, efekata hirurških intervencija, oporavka ili efekata terapije. Sve veći broj kliničara bira pogodne tretmane za lečenje pacijenata na osnovu informacija o kinematici i kinetici hoda. Procena i kvantifikacija parametara hoda je važan zahtev u oblasti ortopedije i rehabilitacije, ali takođe i u sportu, rekreaciji i posebno u razvoju tehnologija za ljude u procesu starenja. U cilju objektivne procene obrasca hoda, razvijen je bežični senzorski sistem čije su senzorske jedinice bežične, malih dimenzija i jednostavno se montiraju na segmente nogu subjekta čiji se hoda analizira. Senzorske jedinice podržavaju 3D inercijalne senzore (senzore ubrzanja i ugaonih brzina, tj. akcelerometre i žiroskope), kao i senzore sile. Osnovni cilj istraživanja je doprinos metodologiji za obradu podataka sa inercijalnih senzora i razvoj novih metoda obrade signala sa inercijalnih senzora u procesu određivanja kinematičkih veličina koje su uobičajene u analizi hoda (uglovi u zglobovima, brzina kretanja, dužina koraka). Ova metodologija je od posebne važnosti za objektivnu procenu nivoa motornog deficita, progresa bolesti i efikasnosti terapija, kao i efikasnosti primenjene motorne kontrole (prilikom funkcionalne električne stimulacije). U toku istraživanja razvijeno je nekoliko metoda za računanje uglova segmenata nogu ili zglobova, u zavisnosti od senzorske konfiguracije i složenosti algoritma. U disertaciji su odvojeno prikazani slučajevi u kojima je neophodno posmatrati kretanje u prostoru (3D analiza) i mnogo češći slučaj kad se kinematika može redukovati na sagitalnu ravan (2D analiza). Algoritmi uključuju i kalibraciju senzora, eliminaciju viii drifta, rekonstrukciju trajektorije i izračunavanje niza drugih relevantnih podataka koji karakterišu obrazac hoda. Dobijeni rezultati su poređeni sa postojećim sistemima za analizu hoda koji su validirani za kliničke primene. (sistemi sa kamerama, goniometri, enkoderi). U toku istraživanja je, prema posebno definisanim modelima koje su odobrili etički komiteti medicinskih ustanova, bio snimljen hod većeg broja ispitanika sa raznim tipovima motornih deficita. Jedan deo ovih istraživanja je rađen na Klinici za neurologiju Kliničkog centra Srbije, a drugi deo u Institutu za rehabilitaciju „Dr. Miroslav Zotović“ u Beogradu. Cilj ovih kliničkih studija je bio da se na osnovu analize odredi koje specifičnosti moraju da se analiziraju posebnim metodama obrade signala, a koje daju značajne elemente za dijagnostiku motornog deficita. Poseban fokus teze je klinička primena razvijenih algoritama na analizu hoda pacijenata sa Parkinsonovom bolešću koji, osim uobičajene analize parametara hoda, omogućavaju prepoznavanje promena parametara hoda u okviru snimljene sekvence hoda, prepoznavanje poremećaja nastalih u toku snimljene sekvence hoda (tzv. „freezing“ epizode), određivanje njihovog trajanja i klasifikaciju podtipova svake od tih epizoda. Rezultat teze takodje obuhvata i opise kliničkog protokola za praćenje efekata terapije pacijenata nakon moždanog udara pomoću razvijenog senzorskog sistema. Rezultat ovog istraživanja je sistem koji omogućava jednostavno postavljanje i korišćenje sistema sa inercijalnim senzorima. Za ovaj sistem razvijen je korisnički interfejs koji omogućava obradu snimljenih podataka primenom razvijenih metoda sa ciljem generisanja objektivne slike motornog deficita pacijenta koji je od interesa za kliničare (prikaz vremensko-prostornih parametara hoda, uglova u zglobovima, trajektorije nogu, profila sila), kao i smeštanje dobijenih rezultata u bazu podataka ili izvoz u Excel formatu. Ključne reči: inercijalni senzori, računanje uglova, objektivna evaluacija hoda, pacijenti sa oštećenim obrascem hoda, „freezing“ eipizode pacijenata sa Parkinsonovom bolešću Naučna oblast: tehničke nauke, elektrotehnika Uža naučna oblast: biomedicinsko inženjerstvo UDK broj: 621.3 ix Inertial Sensors Signal Processing Methods for Gait Analysis of Patients with Impaired Gait Patterns Abstract - Gait analysis has become a widely used clinical tool which provides objective evaluation of the gait pattern, the effects of surgical interventions, recovery or therapy progress, and more and more clinicians are choosing therapy treatments based on gait kinematics and kinetics. Measuring gait parameters is an important requirement in the orthopedic and rehabilitation fields, but also in sports and fitness, and development of technologies for elderly population. In order to provide objective evaluation of the gait pattern, we have developed sensor system with light and small wireless sensor units, which can be easily mounted on body. These sensor units comprise 3-D inertial sensors (accelerometers and gyroscopes) and force sensing resistors, and our recommended setup includes one sensor unit per each segment of both legs. The main goal of this research is contribution to the methodology for processing of signals from inertial sensors (accelerometer pairs, or accelerometer and gyroscope sensor units). By using signal processing algorithms developed for this research, inertial sensors allow objective assessment of the quality of the gait pattern. This methodology is especially important for assessment of the motor deficit, progress of the disease and therapy effectiveness, and effectiveness of performed motor control (functional electrical stimulation). We have developed several methods for estimation of leg segment angles and joint angles, which differ in sensor configuration and algorithm complexity. Methods based only on accelerometers offer reliable angle estimations, which are limited to sagittal plane analysis, while the method using accelerometers and gyroscopes allows 3- D analysis. All this algorithms include sensor calibration, drift minimization, trajectory x reconstruction and calculation of numerous other parameters relevant to gait pattern analysis. The obtained results were compared with other commercial systems which are validated for clinical applications (camera systems, goniometers, encoders). This research also includes two clinical studies. The study with patients with Parkinson’s disease was performed at the Neurology Clinic, Clinical Centre of Serbia, in Belgrade. The other study, recording hemiplegic patients during their recovery after stroke, was performed at the Rehabilitation clinic “Dr. Miroslav Zotovic” in Belgrade. The database recorded in these two clinics comprises more than one hundred patients with different types and levels of gait impairments. The goal of these studies was to determine which specificities of the gait patterns or gait deformities require special signal processing methods, and which characteristics are essential for diagnostic of motor deficit. A special focus of the thesis is clinical application of the developed gait analysis algorithms for patients with PD, which required the development of specific algorithms for signal processing that would provide detailed assessment of their gait patterns. For this application, besides providing the usual gait parameters we developed the algorithms for recognition of the changes and deformities in gait pattern (freezing of gait), duration and classification of these episodes. This thesis also includes description of the clinical protocol for monitoring the rehabilitation effects for patients after stroke. The result of this research is a system that enables simple donning and doffing of the hardware, and simple use (both signal recording and processing). For this system, graphical user interface has been developed which enables processing of recorded files (gait sequences) by implementing the developed methods and plotting the kinematics curves, force profiles, spectrum of the movement or spatio-temporal gait parameters. As well as saving the obtained results in database or exporting to Excel format. Key words: inertial sensors, angle estimation, objective evaluation of gait pattern, patients with gait impairment, freezing of gait episodes in patients with Parkinsonean disease Scientific area: technical sciences, electrical engineering Specific scientific area: biomedical engineering UDK number: 621.3 xi Table of Contents Preface .......................................................................................................................... 1 Objectives of the thesis............................................................................................. 2 Starting hypotheses of the thesis .............................................................................. 3 Research contributions ............................................................................................. 3 Thesis outline............................................................................................................ 5 1. Introduction .............................................................................................................. 7 The significance of objective gait analysis............................................................... 7 Gait analysis systems.............................................................................................. 10 Inertial sensors.................................................................................................... 12 Other gait analysis technologies ......................................................................... 14 Our hardware: custom made wireless sensor system ............................................. 16 2. Processing data from inertial sensors ..................................................................... 19 Inertial sensors and inertial measurement units.................................................. 19 IMUs for body motion sensing........................................................................... 20 Angle estimation methods for lower limb kinematics........................................ 21 Trajectory estimation for lower limb kinematics ............................................... 25 2-D Gait Analysis: Processing Accelerometer Pairs .................................................. 26 3. “Smart” filtering algorithm for 2-D angle estimation based on accelerometers .... 27 Introduction ............................................................................................................ 27 Experimental Section.............................................................................................. 27 Sensor system ..................................................................................................... 27 Algorithm ........................................................................................................... 29 Experiments ........................................................................................................ 35 xii Processing of measured data .............................................................................. 35 Results .................................................................................................................... 36 Discussion and Conclusions ................................................................................... 39 4. Nonlinear optimization for drift removal in estimation of gait kinematics based on accelerometers ............................................................................................................ 41 Drift removal by polynomials ................................................................................ 41 Drift removal from x and y components of acceleration.................................... 42 Drift removal from z component of acceleration ............................................... 43 Application of nonlinear optimization.................................................................... 44 Method and Materials......................................................................................... 44 Results .................................................................................................................... 52 Discussion and Conclusion..................................................................................... 53 Extension to 3-D analysis ................................................................................... 55 5. 3-D Gait Analysis from Accelerometers and Gyroscopes...................................... 57 Introduction ........................................................................................................ 57 Transformation matrices..................................................................................... 57 3-D gait analysis for SENSY system...................................................................... 61 Initialization and initial angles ........................................................................... 61 Creating initial rotation matrix ........................................................................... 61 Foot segment ...................................................................................................... 62 Incremental matrix rotation ( iR ) ....................................................................... 65 Stride length estimation .......................................................................................... 66 Trajectory reconstruction........................................................................................ 66 Angle estimation..................................................................................................... 68 Angles in sagittal plane ...................................................................................... 68 Defining 3-D angles ........................................................................................... 68 Stride segmentation ................................................................................................ 72 Eliminating drifts from angles................................................................................ 73 Discussion and Conclusion..................................................................................... 75 6. Processing of signals from rigid double pendulum: comparison of angle estimation methods....................................................................................................................... 78 Experiment setup .................................................................................................... 79 xiii Results for the upper segment ................................................................................ 82 Basic methods for angle estimation.................................................................... 82 Engle estimation based on accelerometer pairs.................................................. 87 Willemsen method.............................................................................................. 89 “Smart” filtering method (integration in frequency domain) ............................. 91 Nonlinear optimization for drift removal ........................................................... 93 Results for the lower segment ................................................................................ 98 Basic angle estimation methods ......................................................................... 98 Angle estimation from accelerometer pairs........................................................ 99 Willemsen method............................................................................................ 100 Smart filtering (integration in frequency domain)............................................ 101 Nonlinear optimization for drift removal ......................................................... 102 Conclusion ............................................................................................................ 104 7. Clinical Applications ............................................................................................ 106 Clinical gait assessment: no-tech, low-tech and high-tech methods .................... 106 Gait analysis for patients with Parkinson’s disease.............................................. 110 Freezing of gait phenomenon ........................................................................... 110 Gait assessment: Recognizing gait disturbances and freezing episodes............... 113 1. Gait disturbances and FOG recognition based on FSR sensors ................... 113 2. FOG recognition and classification based on inertial sensors...................... 117 3. Gait assessment by spectrogram................................................................... 126 Discussion and conclusions .............................................................................. 129 Gait analysis for patients with hemiplegia ........................................................... 131 Recording protocol ........................................................................................... 131 Quantifying rehabilitation progress/success..................................................... 132 Overall investigated clinical applicability of the system.................................. 136 8. Graphic user interface for gait analysis software ................................................. 138 9. Conclusion and future research ............................................................................ 144 Contribution of dissertation .................................................................................. 144 Suitability of the proposed system.................................................................... 146 Perspective and future research ............................................................................ 147 Extending 2-D to 3-D ....................................................................................... 147 xiv Movement analysis ........................................................................................... 148 Gait analysis ..................................................................................................... 148 Movement control ............................................................................................ 148 Holter monitor .................................................................................................. 149 Hardware adaptations: miniaturization and minimization................................ 149 References ................................................................................................................ 151 Appendix .................................................................................................................. 159 Gait terminology....................................................................................................... 160 Gait phases........................................................................................................ 160 Spatio-temporal gait parameters........................................................................... 164 Kinematic gait parameters .................................................................................... 165 Gait event detections ............................................................................................ 166 Gait Phases Nomenclature.................................................................................... 167 Kinematics and Kinetics curves: Normal gait pattern .......................................... 169 SENSY – Data acquisition system ....................................................................... 170 Curriculum Vitae ...................................................................................................... 172 xv Table of Figures Figure 1.1. – Commercial technologies for gait analysis. Upper panel: gait analysis systems based on inertial sensors, lower panel: other gait analysis systems sorted according to requirements for recording space. ...........................12 Figure 1.2. – Inertial measurement units from Xsens (MTw).....................................13 Figure 1.3. – FAB system: mounted on a subject (left), sensor units (middle), sensor insoles (right) ............................................................................................13 Figure 1.4. – KinetiSense sensors from CleveMed. Sensor units record kinematics and electromyography. .........................................................................................14 Figure 1.5 – Architecture of the SENSY system.........................................................17 Figure 1.6. – SENSY system; (a) sensor unit, (b) sensors mounted on a subject, (c) complete equipment for gait analysis in light, compact and portable package.................................................................................................................18 Figure 2.1 – Body planes: sagittal, coronal, and transverse. .......................................22 Figure 2.2 – Schematic of the system configuration with the coordinate systems......23 Figure 3.1. Setup of the sensor system. a) photo of the sensors mounted on the body during the gait analysis, b) schematic of the system configuration with the coordinate systems..........................................................................................29 Figure 3.2. Rod with two accelerometers and the coordinate system for analysis of movement in the sagittal plane. ............................................................................30 Figure 3.3. Joint angle (bottom panel) and angular velocity (middle panel) obtained by numerical integration of the measured acceleration of the segment (top panel). Dashed line envelope on the bottom panel is fitted through the points where the knee should be fully extended with zero degree xvi joint angle. However, due to the integration drift, instead of remaining approximately constant, this line has a parabolic shape.......................................31 Figure 3.4. Normalized spectrum of the angular acceleration for knee angle (shown in Figure 3.3) and magnitudes of the function 220 / sω , transfer function of a second-order Butterworth low-pass filter (lpf), and function of this low-pass filter combined with a high-pass filter (lpf+hpf)............................33 Figure 3.5. Influence of 0f on angle estimation: filtering for several cutoff frequencies compared with angle acquired from goniometers.............................34 Figure 3.6. Knee and ankle joint angles measured with goniometers and estimated from accelerometers. Thick lines represent angles from goniometers (dashed line) and accelerometers (solid line) for optimal filter frequency. Blue areas show the range of angle values estimated from accelerometers when filtered with various frequencies in the range [ gcf /3, gcf /2]. ...........................................36 Figure 3.7. Optimal filtering frequencies for knee and ankle angles versus gait cycle frequency.....................................................................................................37 Figure 3.8. Pearson’s correlation coefficient (PCC) and RMSE between goniometers and estimated angles, for knee and ankle angles, vs. filtering frequency in the range [ gcf /3, gcf /2]. The normalization is with respect to the optimal frequency. Each curve corresponds to different gait velocity (in m/s). ..38 Figure 3.9. Boxplots presenting comparison between joint angles (for knee and ankle angles) calculated from accelerometers and goniometers. Left: Pearson’s correlation coefficient (PCC). Right: root-mean-square error (RMSE). ...............................................................................................................38 Figure 4.1 – (a) Experimental setup showing sensors, markers, and coordinate systems. (b) Sensor unit with two accelerometers. ..............................................45 Figure 4.2 – The square of the angular velocity obtained directly from Eq. 2, and by integrating raw accelerometer data (Eq. 3) and using accelerometer data corrected by the optimized polynomial. ...............................................................48 Figure 4.3 – Comparison of the angle without drift correction, angle optimized up to the current time (real-time application), angle optimized until the next xvii anchoring point (optimized full stride), and angle from reference optical system. The anchoring intervals are between square and round markers. ...........49 Figure 4.4 – Segment angles estimated from the accelerometer data and measured by the camera based system (reference)...............................................................51 Figure 4.5 – Joint angles estimated from the accelerometer data and measured by the camera based system (reference)....................................................................51 Figure 4.6 – Stick diagram of the leg and trajectories of leg joints estimated from the accelerometer data (dashed lines) and measured by the camera based system (solid lines). ..............................................................................................51 Figure 4.7 –The RMS error and Pearson correlation coefficients between angles estimated from camera system and the proposed method for the ten subjects in the study. ..........................................................................................................53 Figure 5.1 – Angles between the axes of the primed and unprimed system. .............58 Figure 5.2 – Schematic of the system configuration with the coordinate systems......61 Figure 5.3 – trajectory reconstruction, SENSY (red lines) compared with motion capture camera system (black lines).....................................................................68 Figure 5.4 – Comparison of the segment angle estimation for all angles in all three planes: sagittal, coronal, and transversal. Red line shows SENSY outputs, black line shows outputs of the camera system....................................................71 Figure 5.6 – Drift elimination for segment angles. Red line represents the original signal, and the green line shows the corrected signal. Black square markers show moments in the stance phase when the foot has minimal movements........74 Figure 6.1– Double rigid pendulum with encoders placed as joints and inertial sensors placed along segments. Accelerometer positions in sensor units are denoted by round markers and gyroscope positions by square markers. .............80 Figure 6.2 – Joint angles based on encoders, complete recorded sequence. ...............81 Figure 6.3 – Joint and segment angles based on encoders, complete recorded sequence. ..............................................................................................................82 Figure 6.4 – Angle of the upper segment obtained from encoder, inclinometer, and integrated gyroscope......................................................................................83 Figure 6.5 - Angle of the upper segment obtained from encoder, inclinometer, and integrated gyroscope (exhibiting drift), complete recorded sequence. ................83 xviii Figure 6.6 – Differences in the angles estimated by inclinometer and integrated gyroscope, compared to the encoder angle. .........................................................84 Figure 6.7 – Angular velocity obtained as derivatives of encoder and inclinometer angle, and directly from gyroscope. .....................................................................85 Figure 6.8 – Angular acceleration obtained as double derivatives of encoder and inclinometer angle, and as derivative from gyroscope.........................................85 Figure 6.9 – Segment angles obtained by integrating signals from gyroscopes positioned along the segment, sequence at the beginning of recorded file (transition state) ....................................................................................................86 Figure 6.10 – Segment angles obtained by integrating signals from gyroscopes positioned along the segment, sequence at the end of recorded file. Different sensors exhibit different amount of drift. .............................................................87 Figure 6.11 – Angular acceleration obtained from different accelerometer pairs from upper segment, and by differentiation of gyroscope signal.........................87 Figure 6.12 – Angular velocity obtained by integrating from different accelerometer pairs from upper segment, and directly from gyroscope signal. Different accelerometer combinations exhibit different amount of drift. ............88 Figure 6.13 – Angles of the upper segment obtained by double integration of different accelerometer pairs and by integration of gyroscope signal. Different accelerometer combinations exhibit different amount of drift. ............89 Figure 6.14 – Angle obtained by Willemsen method from accelerometers a12 and a13, compared to reference encoder angle. ...........................................................89 Figure 6.15 – Angular velocity obtained by differentiating the angles obtained by Willemsen method, compared to gyroscope signal (reference). ..........................90 Figure 6.16 – Angular acceleration obtained by double differentiating the angles obtained by Willemsen method, compared to angular acceleration directly from a12 and a13 (reference). .................................................................................91 Figure 6.17 – Angle obtained by smart filtering, applied on accelerometer pair, compared to encoder angle (black line), during the static sequence from the beginning of the recorded sequence. ....................................................................92 xix Figure 6.18 – Angle obtained by smart filtering from accelerometer pair, compared to encoder angle (black line), ending sequence from the recorded file.........................................................................................................................92 Figure 6.19 – Quadratic angular velocity obtained from y axis of accelerometer pair (red line), preconditioned (all negative values are clipped with 5pt moving average filter), and directly from gyroscope ...........................................93 Figure 6.20 – Input data for optimization, linear accelerations, angular acceleration, and quadratic angular velocity obtained from accelerometers. Upper panel: sequence at the beginning of file (static phase), lower panel: sequence from the oscillating phase.....................................................................94 Figure 6.21 – Quadratic angular velocity: directly from the subtraction of x axes of the accelerometer pair (blue line), after optimization (red line) ......................95 Figure 6.22 – Angular velocity after optimization (red line), compared to the gyroscope signal (blue line). ................................................................................96 Figure 6.23 – Angle after optimization, compared to the angle obtained by inclinometer and from encoder.............................................................................97 Figure 6.24 – Angle after optimization (red line), compared to the angle from encoder (blue line)................................................................................................98 Figure 6.25 – Angle of the lower segment, obtained directly from encoders (red line) and by integrating g21 and g22 (green and blue line). ...................................99 Figure 6.26 – Angular acceleration of the lower segment, measured by accelerometer pair a21 and a22, and by differentiation of gyroscope signal g21. ...99 Figure 6.27 – Angular velocity obtained from accelerometer pair, from differentiation of encoder signal, and gyroscope (reference).............................100 Figure 6.28 – Angle of the lower segment obtained by double integration of angular acceleration obtained from accelerometer pair, and by integrating gyroscope signal. ................................................................................................100 Figure 6.29 – Joint angle (angle between upper and lower pendulum segment), obtained from accelerometer pairs from upper and lower segment by Willemsen method, compared to encoder angle (reference). .............................101 Figure 6.30 – Angle of the lower segment obtained by smart filtering (red line) and from encoders (black line). Upper panel: sequence from the beginning of xx file – transition state. Lower panel: zoomed sequence from the middle of the signal (stationary state).......................................................................................102 Figure 6.31 – Quadratic angular velocity from directly from accelerometer pairs, from preconditioned accelerometer pairs, and from gyroscope signal. .............103 Figure 6.32 – Angle of the lower segment obtained from accelerometer pair by optimization method, and directly from encoder (reference).............................103 Figure 7.1. - Walking protocol. Patients were asked to stand up from the chair placed in the corridor, walk toward the room, pass a doorway, turn 180° to the left (U-turn), and walk the same route back, ending with a turn to sit back in the chair. .........................................................................................................114 Figure 7.2 –An example of FOG detection using PCC analysis. Middle panel shows the ground reaction force of the entire sequence and one “normal” stride (red line). Bottom panel: result of computed PCC between the highlighted stride and the entire recorded sequence. .........................................116 Figure 7.3. Illustration of the PCC method in different walking situations. Upper panel: ground reaction force with one “normal/regular” stride selected (thick red line). Lower panel: results of computed PCC between selected “normal” stride and the entire gait sequence. Bottom bar: FOG episodes identified from video recordings. ................................................................................................117 Figure 7.4. Choosing characteristic frequency of the spectrum. Upper panel: shank accelerometer signal (vertical axis) in time domain, stride #37. Lower panel: its power spectrum (normalized) and characteristic frequencies – median frequency, peak frequency, mean power frequency, and mean frequency ............................................................................................................118 Figure 7.5. Floor plan and the pathway along which patients are walking. ..............120 Figure 7.6 – Stride segmentation shown on FSR signals, upper panel – left leg, lower panel – right leg. Each red square marker represents one recognized stride. Horizontal bars on the bottom show FOG episodes recognized by clinician from video recording. ..........................................................................121 Figure 7.7. Variability of stride length, stride time, stride speed, and median frequencies along the recorded sequence, left leg and right leg. Horizontal xxi axes show the stride number within the recorded sequence. Upper and lower threshold for classification rules are shown with black dashed lines.................122 Figure 7.8. Ground reaction forces for both legs and color-coded stride segmentation which shows the following classification: normal strides (black markers), small strides (blue markers), shuffling strides (green markers), FOG with trembling (red markers) and FOG with motor blocks (pink markers)..............................................................................................................123 Figure 7.9. Illustration of patient’s 3-D trajectory with recognized and classified FOG episodes, shown by color coding. Comparison of (a) Left leg and (b) right leg show asymmetry of the present gait disturbances................................125 Figure 7.10. Example of FOG interpretation from spectrogram: FOG episode began as a motor block (stride#2) and was followed by festinations episode with frequencies between 5 and 7 Hz (stride#3) ................................................127 Figure 7.11. Spectrograms of power spectrums from vertical acceleration axes for all three leg segments. (a) left leg, (b) right leg. For each spectrogram, the horizontal axis shows the stride number, while the vertical axis shows short- term power spectrum of the signal for each stride. ............................................128 Figure 7.12. SENSY clinical setup, assessment of the hemiplegic patients..............132 Figure 7.13. Therapy follow-up, comparison of ground reaction forces and joint angles before therapy, three weeks later and six months since the beginning of therapy............................................................................................................133 Figure 7.14. Therapy follow-up, comparison of ground reaction forces and joint angles before and after therapy (6 months later). Horizontal axes show that time interval of one stride before therapy (grey areas) corresponds to two and a half strides after therapy (red lines). ................................................................133 Figure 8.1 – Initialization of graphic user interface ..................................................138 Figure 8.2 – Initial software panel showing original sensor data. .............................139 Figure 8.3 –Kinematics panel showing joint angles of the hip, knee, and ankle for both legs. ............................................................................................................140 Figure 8.4 –Kinematics panel showing joint angles of the hip, knee, and ankle for both legs, strides plotted over each other showing stride-to-stride variability in time and range of motion. ..............................................................................140 xxii Figure 8.5 –Kinematics panel 2-D reconstruction of leg trajectories, for left and right leg...............................................................................................................141 Figure 8.6 – Kinetics panel showing force profiles, ground reaction forces and graphical visualization of force profiles. ............................................................142 Figure 8.7 – Spectrum analysis, selected signal in time and frequency domain (left side of panel) and its spectrogram (right side of panel) .....................................142 Figure 8.8 – Statistics panel showing average values and standard deviations of spatio-temporal gait parameters (left side of the panel). Right side of panel provides export of the processed data to Excel file............................................143 Figure A.1 - Gait cycle with sub phases of gait.........................................................161 Figure A.2. Gait event detection: threshold is applied to normalized FSR signals in order to detect heel contact (triangular marker pointing downwards), heel off (square marker), and toe off* moments (triangular marker pointing upwards). ............................................................................................................166 Figure A.3 Joint angles of hip, knee, and ankle joint.Typical kinematic curves for healthy gait .........................................................................................................169 Figure A.4 Ground reaction forces: load force (vertical reaction force) and friction force (horizontal reaction force). Typical force profiles for healthy gait...........169 Figure A.5 – Running SENSY Data acquisition software. .......................................170 Figure A.6 – Data delay status bar: green status showing there is no data delay , yellow status showing there is some transmission delay and the buffer is accumulating unsent data, red status indicates that the buffer is nearly full, and circled red status showing the buffer is full and all future data will be lost. .....................................................................................................................171 xxiii Table of Tables Table 1.1. Specifications of the system .......................................................................17 Table 7.1. Typical values, lower and higher thresholds for gait parameters for rule-based classification .....................................................................................119 Table 7.2. Rule-based classification of gait disturbances for PD patients. ...............119 Table 7.3. Therapy follow-up, comparison of gait patterns before, after three weeks of therapy and after the therapy (6 months later). ...................................134 Table 7.4. Therapy follow-up, comparison of outputs for three clinical scales, before therapy, after three weeks of therapy and after the therapy (6 months later)....................................................................................................................134 Table 7.5 – List of clinical experiments performed with SENSY system.................136 1 Preface This thesis is a part of the BMIT group’s more complex research for rehabilitation improvements which comprises 1) development of new sophisticated systems based on functional electrical stimulation and soft-robotic gait assistance, and 2) improvement of methodologies for objective gait assessment which is applicable for clinical practice. Technical research activities described in this thesis were performed in the Laboratory for Biomedical Instrumentation and Technologies, School of Electrical Engineering, University of Belgrade. Validations of the developed methods were performed in Health Technologies Unit, Tecnalia Research Center, San Sebastian, Spain, and in Center for Sensory-Motor Interaction, Aalborg University, Aalborg, Denmark. Clinical research activities were performed at the Neurology clinic, Clinical Center of Serbia, and Rehabilitation clinic “Dr. Miroslav Zotovic”, Belgrade, Serbia. Research presented in this thesis was supported by the Serbian Ministry of Education and Science, grant #145041, “Functional electrical therapy for forming motor patterns after cerebrovascular insult” (2008-2010), and grant #175016, “The effects of assistive systems in neurorehabilitation: recovery of sensory-motor functions” (2011-2014). Project leader of both projects is Professor Mirjana Popovic. Research and development of the sensor system for gait analysis (SENSY) and its application was supported by Tecnalia Serbia, Belgrade, Serbia (which is a node of Tecnalia Research Center, San Sebastian, Spain). 2 Objectives of the thesis In this thesis, we have developed a wearable ambulatory system based on kinematic sensors (MEMS accelerometers and gyroscopes) attached on lower limbs which can be used for objective evaluation of the gait kinematics. Gait kinematics is recorded with wireless sensor units which are placed on each leg segment of both legs. Besides kinematics, this system also supports kinetic analysis for what we use specially designed shoe insoles with integrated force sensing resistors. The main features of the measurement device are to be light, small, inexpensive and robust, with low energy demands, and not to hinder subject’s walk. The main goal of this research is development of new signal processing methods from data acquired from inertial sensors in order to obtain estimation of kinematic quantities typical for gait analysis (segment and joint angles, trajectories). The inertial sensors, organized either as combination of accelerometer pairs, or combination of accelerometer and gyroscope, can be used in ways that eliminate the consequences of numerical integration of digital signals and provide objective assessment of the quality of the gait pattern. This methodology is especially important for assessment of the motor deficit, progress of the disease, therapy effectiveness, and effectiveness of the applied motor control. In order to evaluate the quality of our new methods, we had to validate our gait analysis in clinical environment by recording gait pattern from subjects with different types and levels of deviations from the normal gait or motor deficits. Special focus of this thesis is application of the developed methods to pathological gait recorded from patients with hemiplegia and patients with Parkinson’s disease. For patients with Parkinson’s disease who can exhibit episodic gait disturbances (freezing of gait), our goal was to provide reliable algorithm that could detect these disturbances, their duration and sub-type. Results of this research provide useful gait analysis system that allows simple (neither time nor skill demanding) donning and doffing and simple calibration of the system. Result of this research includes user interface which allows user to process recorded data and get gait parameters which objectively and fully describe subject’s gait pattern and potential motor deficit. These outputs are stored in database for further post-processing or comparison with follow up results. 3 This research follows modern trends in the development of measurement systems for objective evaluation in rehabilitation, control of assistive systems, but also in sports and fitness, and development of technologies for elderly population. Starting hypotheses of the thesis • H1: It is possible to use only accelerometers for reliable and accurate gait kinematic analysis, and to obtain drift-free angle estimation of lower limbs. • H2: Inertial sensors (accelerometers and gyroscopes) can be used for 3-D gait analysis and can provide objective estimation of a subject’s gait pattern. • H3: By applying the methods developed in this research, inertial sensors can be used to detect and classify gait disturbances of patients with Parkinson’s disease • H4: By applying the methods developed in this research, inertial sensors can provide objective evaluation of the patient’s condition or change of patient’s gait pattern to clinicians, allowing them to monitor the therapy effects during rehabilitation process. Research contributions • Analysis of inertial sensors applicability for gait assessment • New methods for angle estimation, differing in their complexity and sensor configuration • Investigating sensor limitations for angle estimations • Description of the clinical protocol for gait pattern analysis of patients with Parkinson’s disease • New methods for recognition and classification of gait disturbances according to familiar sub-types • Description of the clinical protocol for usage of the proposed sensor system for monitoring therapy effects of patients after stroke • Development of a complete gait analysis system which can be used for clinical gait analysis and for assessment of various gait disorders The result of this research is a complete sensor system with custom-made hardware which is appropriate for clinical applications and new methods for signal processing integrated in user-friendly software. The software comprises new 4 algorithms for angle and trajectory estimations, as well as the algorithm for reliable stride segmentation applicable for all levels of gait impairment. Angle estimation methods based only on accelerometers offer reliable angle estimations, which are limited to sagittal plane analysis, while the method using accelerometers and gyroscopes provides 3-D analysis. All these algorithms include sensor calibration, drift minimization, trajectory reconstruction and calculation of numerous other parameters relevant to gait pattern analysis. The thesis is based on the publications listed below. A. Journal publications: 1) Kojović J., Djurić-Jovičić M., Došen S., Popović M.B., Popović D.B., “Sensor-driven four-channel stimulation of paretic leg: Functional electrical walking therapy”, Journal of Neuroscience Methods, June 2009, Volume 181(1), pp. 100-105. 2) Popović M.B, Djurić-Jovičić M., Petrović I., Radovanović S., Kostić V., “A simple method to assess freezing of gait in Parkinson's disease patients”, Braz J Med Biol Res, September 2010, Volume 43(9), pp. 883-889. 3) Djurić-Jovičić M., Jovičić N., Popović D.B., “Kinematics of Gait: New Method for Angle Estimation Based on Accelerometers”, Sensors, November 2011, Volume 11(11), pp. 10571-10585. 4) Djurić-Jovičić M., Jovičić N., Popović D.B., Djordjević A.R., “Nonlinear Optimization for Drift Removal in Estimation of Gait Kinematics Based on Accelerometers”, Journal of Biomechanics, November 2012, Volume 45(16), pp. 2849-2854. B. International conference publications: 1) Djurić-Jovičić M., Jovičić N.S., Milovanovic I., Radovanović S., Kresojevic N., Popović M.B, “Classification of Walking Patterns in Parkinson’s Disease Patients Based on Inertial Sensor Data”, Proceedings from the 10th Symposium on Neural network Applications in Electrical Engineering, Neurel 2010, September 23-25, Belgrade, Serbia, pp. 3-6,. 2) Djurić-Jovičić, M., Milovanović I.P., Jovičić, N.S., Popović D.B., “Walkaround assisted walking of stroke patients”, Proceedings from the Medical Physics and Biomedical Engineering Conference, Sept. 7-12, 2009 Munich, Germany, pp. 299-301. 5 Thesis outline Chapter 1 explains the need for objective gait analysis. It also provides information about different gait analysis technologies which are commercially available on the market. Among these gait analysis systems, we also present our own system, called SENSY, which was developed in parallel with this thesis. SENSY is a wireless wearable gait analysis system based on inertial sensors, and all the gait recordings acquired during the research for this thesis were performed on SENSY hardware. Chapter 2 introduces inertial sensors, their applicability, their advantages and possibilities for implementation, as well as computational disadvantages. It reviews different methodologies for processing of signals obtained by these sensors, with a focus on estimation of segment and joint angles. Chapters 3 and 4 describe our new methods for estimation of the lower limb angles in the sagittal plane, based on accelerometer pairs. Chapter 3 describes a new method for estimation of lower limb segment and joint angles based on adaptive band pass filtering of the differences of signals from parallel axes from two accelerometers mounted on the same sensor unit. This method eliminates the need for double integration as well as the drift typical for double integration. The other method, described in chapter 4, is computationally more complex and uses nonlinear optimization for drift removal. The key feature of the proposed method is to model the drift by a slowly varying function of time (a cubic spline polynomial) and evaluate the polynomial coefficients by nonlinear numerical simplex optimization with the goal to reduce the drift. Besides angle estimation, this method also provides trajectory estimation. The method described in chapter 4 can be used both for real-time and off-line analysis of gait. Chapter 5 describes angle estimation method based on accelerometer and gyroscope combination. In order to provide 3-D gait analysis, it uses transformation matrices and combines human locomotion and biomechanical constrains to fuse accelerometer and gyroscope data, as well as polynomial fitting to eliminate the drift. This chapter also describes methods for stride segmentation which will be a base for some clinical applications. 6 Chapter 6 shows comparison of the methods described in chapters 3, 4, and 5. These methods are tested on a double rigid pendulum, physically analog to a leg model, and validated with digital encoders. Chapter 7 describes investigated clinical applications of the developed methods, and includes two clinical studies. One part of clinical recordings were performed at the Neurology clinic, Clinical Centre of Serbia, in Belgrade (recording patients with PD) and the other part was performed at the Rehabilitation clinic “Dr. Miroslav Zotović” in Belgrade (recording hemiplegic patients during recovery after stroke). We investigated clinical application of the developed gait analysis algorithms for patients with PD. In addition to typical gait assessment comprising gait assessment through usual gait parameters, we have developed algorithms for recognition of the changes and deformities in gait pattern (freezing of gait episodes), and classification of these episodes. This thesis also includes description of the clinical protocol for monitoring the rehabilitation effects for patients after stroke. Chapter 8 proposes new gait analysis software with user-friendly interface which does not require specially educated (technically) clinical staff. This software outputs kinematic and kinetic curves including segment and joint angles, trajectory reconstruction and ground reaction forces, performs spectrum analysis and spatio- temporal gait parameters. Upon signal processing, data can be exported to Excel and saved in database. Finally, chapter 9 summarizes the contribution of this thesis and outlines some perspectives of the proposed methods. Gait terminology, definition of gait phases and spatio-temporal gait parameters, which are necessity from transferring the outputs of our developed methods to clinical practice are described in Appendix. 7 1. Introduction The significance of objective gait analysis Human gait is one of the most difficult tasks that we learn but, once learned, it becomes almost subconscious. Only when walking is disturbed by injury, disease, degeneration or fatigue, we realize our limited understanding of this complex biomechanical process (Winter, 1984). There are three main branches in studying gait analysis: gait kinetics, dynamic electromyography (EMG), and gait kinematics. A comprehensive gait analysis usually includes all three (Vaughan et al. 1992) but obtaining this complex information is time demanding, uncomfortable for subject, and requires dedicated space. Kinematics, kinetics and electromyography are fundamental to characterize gait patterns and their underlying mechanisms (Frigo et al., 1996; Romanò et al., 1996). However, simplified kinematic analysis (e.g., spatio-temporal parameters) can also be clinically valuable, and an ambulatory device may be advantageous for these types of applications (Aminian et al., 2004b). Gait kinetics is defined as the forces, moments, and powers that change over the gait cycle. The forces are captured by the use of force plates embedded in a walkway or by using force sensing resistors embedded in shoe insoles, while the moments and powers are estimated from models of the body by using mechanics. Dynamic electromyography refers to the evaluation of muscle activity throughout the gait cycle. This is accomplished through the use of either surface or needle electrodes. Electromyography (EMG) techniques provide detection and monitoring of electrical muscle activity, however, it does not provide a direct measure of movement, and a substantial number of electrodes and huge amount of data 8 processing are required for studying complex movements such as gait (Den Otter et al., 2006; Frigo et al., 2000). Gait kinematics refers to the branch of biomechanics that deals with accelerations, velocities and most importantly angles, in particularly joint angular changes over the gait cycle. Gait analysis has become a widely used clinical tool, and increasing number of clinicians is choosing suitable treatments for their patients based on the information from kinematic and kinetic data. It also provides an effective tool for evaluating and quantifying the effects of surgical interventions. Measuring gait parameters is an important requirement in the orthopedic and rehabilitation fields. At the present time, gait analysis is primarily carried out in one of two ways (a) in a motion laboratory, with full analysis of the motion using highly accurate optical systems, and (b) in a doctor’s office with the physician making visual observation. A complete gait analysis system in the motion laboratory uses an optical motion setup for kinematic data combined with force platforms for assessment of ground reaction forces. Such system is expensive, requires a large space, and cannot be used outside the laboratory environment. The capture volume is also limited to several gait strides. Beside that, camera systems require dedicated working space, and are time demanding both for set up and data processing. Consequently, the use of motion analysis systems is mainly limited to research studies. A different approach is to use low-power portable recording systems carried by the subject for long-term ambulatory measurements. Such portable systems are also suitable for capturing gait information over larger distances and outside the laboratory environment. During the last two decades, with a progress made in micro-electromechanical systems (MEMS), body- mounted sensors consisting of accelerometers and/or rate gyroscopes have been used to obtain kinematic values, such as segments inclination angles, and joint angles (Favre et al., 2008; Williamson et al., 2001; Miyazaki, 1997; Sabatini et al., 2005; Tan et al., 2008). Force sensing resistors have been used for estimation of ground reaction forces (Morin et al., 2002). Although many devices exist for kinematic and kinetic assessment, there is a need for a low-cost system that could be used in routine practice. Such system should be reliable, ambulatory, and easy to use (Sekine et al., 2000; Culhane et al., 2005; 9 Aminian et al., 2004). Being motivated by practical needs we have developed a system that fulfils these requirements, which is described in this thesis. 10 Gait analysis systems As explained in the previous section, gait analysis is important for objective assessment of the effects of the rehabilitation intervention. The most accurate systems for gait analysis are camera-based systems with reflective markers (Furnee, 1997). Camera systems, together with force platforms, are considered as the gold standards among gait analysis systems. Optical gait analysis systems (Camera systems) are currently the most well known and most precise type of gait analysis systems. Two of the main suppliers to gait laboratories are Vicon (www.vicon.com) and Qualisys (www.qualisys.com). Although they are the most accurate system on the market, their high price (60- 150.000€), complexity of preparation for recordings and data analysis makes this system inadequate for clinical gait assessment (especially daily routine assessment of patients with severe gait impairments). The alternatives to camera-based systems are ultrasound systems (Kiss et al., 2004) and magnetic tracking systems (Kobayashi et al., 1997) which allow a complete 3-D kinematic analysis of human movements. Force plates measure the ground reaction forces generated by a body standing on or moving across them, to quantify balance, gait and other parameters of biomechanics. For studies of movements, such as gait analysis, force platforms with at least three pedestals and usually four are used to permit forces that migrate across the plate. Gait laboratories typically record from only one or two steps in the middle of the gait sequence. The platforms are 60 x 60 cm, so that aiming the platforms hinders subject’s natural gait pattern. Major manufacturers are Advanced Mechanical Technology, Inc. (AMTI), http://www.forceandmotion.com, Bertec, http://www.bertec.com, and Kistler International, http://www.kistler.com/us_en- us/KistlerCountryHome_KIC. Over the last decade, many systems using non-traditional methods for gait analysis have been developed. These systems, for example, use laser technology or measure near-body air flow (Bonomi et al., 2010; Palleja et al., 2009) in order to estimate kinematics and spatial gait parameters. Also, electronic carpet or wearable force sensors are used for estimation of ground reaction forces, centre of pressure, and temporal gait parameters (Yun, 2011; Liu et al., 2010). Since there is often a need for 11 gait recording in various environments, portable body-mounted systems are preferred (Allet et al., 2010; Zheng et al., 2005). The portable body-mounted systems allow data acquisition from many steps. The portable systems for kinematics data acquisition directly measure joint angles, or they can record accelerations or angular velocities of the body segments that carry the sensors. Measurement of joint angles can be done with various electrogoniometers (Zheng et al.,2005; Legnani et al., 2000; Shiratsu and Coury, 2003). Particularly convenient are flexible goniometers, which measure the relative angle between two small blocks that are fixed to the body segments (e.g., Biometrics flexible Penny & Giles sensors). The advantages of flexible goniometers are: their output is directly proportional to the angle and their mounting is simpler compared to some other measurement systems. However, they are not sufficiently robust for daily clinical usage. An alternative to goniometers, offered by the progress made in micro- electromechanical systems (MEMS), is the use of inertial measurement units (IMUs) comprised of accelerometers and gyroscopes, and sometimes magnetometers. The advantages of these sensors include their small size and robustness when compared with goniometers. However, the disadvantages of the accelerometers (and gyroscopes) are computational problems for determining the angles (Mayagotia et al., 2002; Aminian et al., 2002; Dejnabadi et al., 2005; Luczak et al., 2006). Most frequently used gait analysis technologies and main manufacturers for each type of system are shown in Fig. 1.1. 12 Figure 1.1. – Commercial technologies for gait analysis. Upper panel: gait analysis systems based on inertial sensors, lower panel: other gait analysis systems sorted according to requirements for recording space. Inertial sensors Researchers use IMU (Inertial Measurement Unit) sensor modules to analyze body movement in a range of different medical applications. Inertial sensors combine accelerometers, rate gyroscopes and sometimes earth’s magnetic field sensors and allow detection of any angular displacement within biomechanical bodies. Commercially available systems based on this technology have been used for a range of research projects into biomechanics, in which the user generally makes his/her own analysis of the data. Researchers can use these to develop their own algorithms, while for clinical applications, clinicians can use the company’s software solutions for specific measurements. XSens’s wireless 3-D motion tracker, MTw (www.xsens.com) is shown in Fig. 1.2. It provides drift-free 3-D orientation, acceleration, angular velocity, earth- magnetic field as well as static pressure that provides evaluation of height. 13 Figure 1.2. – Inertial measurement units from Xsens (MTw) Another commercial representative of the same technology is FAB System (Functional Assessment of Biomechanics, http://www.biosynsystems.net/f-a-b- system) developed by BIOSYN SYSTEMS INC. CA (Noraxon), shown in Fig. 1.3. Their software displays and calculates kinematic and kinetic data in real time and animates body motion with selectable graphical models. The small receiver box can also be used as a belt data logger. An onboard control panel allows setup of the system in the field and without PC and storage to SD card. The motion capture data are processed automatically and the system provides: angle data, force data, torque data, velocity, acceleration, power, pressure on the feet and weight. Figure 1.3. – FAB system: mounted on a subject (left), sensor units (middle), sensor insoles (right) KinetiSense (from CleveMed, Cleveland Medical Devices Inc., http://www.clevemed.com.) is a compact wireless device that integrates inertial motion sensing and electromyography for a wide variety of motion analysis research applications (Fig. 1.4.). KinetiSense units comprise 3-D accelerometers and rate gyroscopes, and support EMG recording of two channels per sensor unit. The KinetiSense system is not intended to measure absolute motion and the software does 14 not calculate joint angle. The software provides linear acceleration along and angular velocity about each axis. Figure 1.4. – KinetiSense sensors from CleveMed. Sensor units record kinematics and electromyography. Other gait analysis technologies Other existing solutions that are based on different technologies (Fig. 1.1, within dashed line) are also mentioned in this section. They provide complete gait analysis (camera based motion capture system with force plates, electromagnetic systems, ultrasound systems) or are specialized in either kinetics or kinematics, or just parts of them (GaitRite, stride analyzer, electrogoniometers etc.). In general, these partial gait analysis systems are less expensive than commercially available IMUs and very simple to use, hence, many clinics or rehabilitation centers will rather purchase some of those systems. Shoe Insole Sensor Systems - The Stride Analyzer (B&L Engineering, www.bleng.com) is a system designed to record foot-floor contact data from footswitches and calculates all the gait parameters obtainable from this data. Patient tests may be performed in any convenient walking area. This system doesn’t provide any information about movement of the leg segments (angles, trajectory etc.). (Electro)Goniometers – widely used (in clinics) for measuring joint angles. They are lightweight, portable, relatively easily applied, do not restrict movements nor interfere in patient activities and adapt well to different body segments. However, they are fragile and require calibration and often suffer form crosstalk between the axes. Most frequently used goniometers are strain gauge goniometers (or flexible goniometers) from Biometrics (http://www.biometricsltd.com/gonio.htm). 15 Optoelectronic (goniometers) systems offer good precision, but their calibration procedures and data analysis are time-consuming. Ultrasound system: 3-D gait analysis with the measuring system CMS-HS (ZEBRIS company, www.zebris.de). The measuring process is based on the technical ultrasound pulse time measurement. In this method small ultrasound markers as well as triple markers are attached to the sacrum, thigh, ankle and 1st metatarsals, on both sides of the body. System outputs the joint rotations in all three planes for the hip, knee and ankle joint. The individual gait phases, step length, cadence and average speed are also evaluated. The EMG and force data as well as the joint rotation angles are shown as normalized and averaged curves with standard deviation. Measuring can be done on a treadmill or a feeding unit can be used in order to extend the walking range on the floor to up to 4 meters. Electronic Walkways - flexible and portable walkways with embedded pressure sensitive sensors. One of the most frequently used electronic walkways is GaitRite (www.gaitrite.com). The advantage of a walkway is evident when transferring objective measures of gait to clinical practice since, for example, neurologically impaired patients tolerate poorly lengthy preparations and many attached cords for recording purposes. The pressure sensor system records the location of the foot activated sensors and the time of their activation or deactivation. It provides the spatial and temporal variables of gait along with a dynamic pressure mapping of each footprint during walking. 16 Our hardware: custom made wireless sensor system Hardware system for acquisition of sensor signals was developed in parallel with development of methods for gait analysis. We started with our first prototype, named BUDA, 24-channels sensor system which supported force sensing resistors and accelerometers, connected by wires to the acquisition device which sent data to remote PC by using Bluetooth communication. Although it worked very well, and we conducted many experiments with it, we discovered that clinical applications need faster setup of the system, and more comfortable feeling for the one who wears it. Therefore, switching to completely wireless sensor units was the next logical step. Creating the next prototype opened the door for many new studies, and with the further development of the MEMS-made gyroscopes and their decreasing price, we created the system it is today: 72-channel sensor system that supports force sensing resistors, 3-D accelerometers, gyroscopes, magnetometers or any other voltage controlled sensor. The current version (called SENSY) is simple for mounting on a subject, comfortable for subjects, simple for installing and recording, and presents reliable system which provides high quality sensor data. SENSY hardware comprises 6 peripheral inertial measurement units (IMUs), one per each leg segment of both legs, and central PC communication unit (connected to USB port of remote computer where signals are monitored and stored) (Fig. 1.5.) Foot IMUs are designed with connectors for force sensors which can be attached either incorporated in shoe insoles or as independent force sensors. IMUs comprise 3- D analog accelerometer, and 3-D analog gyro sensor. The IMUs are powered by rechargeable batteries. The system is also equipped with its appropriate battery charger for IMUs (each unit has an external charger port). 17 Figure 1.5 – Architecture of the SENSY system Peripheral nodes acquire signals from accelerometers, gyroscopes and FSRs, and send data using wireless link to the central node. The central node acts as a master in the network. It synchronizes data from different nodes and implements a scheduling algorithm assuring low packet data losses and high throughput. Every node can “hear” all other nodes enabling simple redirection techniques - if one sensor node looses link with the central node, it can send its data over any other node which is visible by both that node and the central node. That node is called redirection node. The criterion for choosing redirection node is based on the equalization of delays from all nodes – redirection node should be the one that has the shortest delay of its data, as explained in (Jovicic et al., 2012). Central node is connected with computer using USB connection. The technical data are summarized in Table 1.1. Table 1.1. Specifications of the system Sampling rate Fixed 100 Hz Communication (IMUs–central node) Proprietary based on IEEE802.15.4 standard Battery Li-Ion, continuous work for up to 1h, standby up to 5h Accelerometers 3-D analogue ADXL330, Analog Devices Gyroscopes 3-D analogue LPR530, LPY530, Analog Devices Max bit rate 250 kb/s (defined by IEEE5802.15.4 standard) Range 30 m Transmission power Class 1 (100 mW) Datastreams Synchronized Packet loss Proprietary mechanism for reduction of packet loss based on the use of internal buffering and retransmission protocol 18 This system was tested in different environments and for various clinical applications, and proved to have excellent performances in all investigated cases. This hardware is easy to mount, light, and does not hinder subject’s movements. The system is portable and convenient for storing, as shown in Fig. 1.6. Figure 1.6. – SENSY system; (a) sensor unit, (b) sensors mounted on a subject, (c) complete equipment for gait analysis in light, compact and portable package. 19 2. Processing data from inertial sensors As explained in the previous chapter, for many clinical applications there is a need for objective evaluation of gait parameters instead of conventional observational gait analysis. Using wearable sensor systems has numerous advantages over specialized gait laboratories (with camera based motion capture systems recording kinematics and force plates recording kinetics). Inertial sensors can be used in any environment, number of strides is not limited, setup is less time demanding and they are significantly cheaper. However, when it comes to accuracy, camera-based optical systems represent the golden standard offering 0.5 degrees and 3 mm accuracy. Inertial sensors and inertial measurement units The term IMU (Inertial Measurement Unit) is widely used to refer to a chip containing three-axial accelerometers and three-axial gyroscopes and optionally three- axial magnetometers, with axes placed in orthogonal pattern. Accelerometers detect total acceleration in respect to the fixed coordinate system which is a sum of Earth’s gravitational acceleration and acceleration of the body the sensor is attached to (inertial acceleration). Detected accelerations are given through three components of the local coordinate system of IMU. Gyroscopes work based on Coriolis force, and they detect rotational movements i.e., change of sensor orientation, also given through their three components called yaw, roll, and pitch. (These components represent rotation around axes of the local coordinate system.) Angular velocity measured by gyroscope is the same from both coordinates systems, global (room) and local (IMU) (Goldstein et al., 2000; Brizard, 2004). 20 Recently, more and more manufacturers also include three-axial magnetometers in IMUs. This allows better performance for dynamic orientation calculation in attitude and heading reference systems which base on IMUs. However, magnetometers are not included in this study. Inertial sensors suffer from offset which is present in the recorded signal, and by performing numerical integration (or other operations) on this signal, the results exhibit drift and they diverge in time. There are many contributors to the offset in inertial sensors - the sensor itself, the readout electronics, mechanical damping, and all electrical resistances. This is their main, if not the only, disadvantage and numerous researchers in science and industry struggle to find the optimal way to minimize or cancel it. IMUs for body motion sensing Accelerometers are used for long-term monitoring of human movements, for assessment of energy expenditure, physical activity, postural sway, fall detection, postural orientation, activity classification and estimation of temporal gait parameters (Yang et al., 2010; Mizuike et al., 2009; Eng et al., 1995, Takeda et al., 2009; Yang et al., 2011). However, only few papers report using only accelerometers for angle estimation, or position and orientation estimation due to their tendency to drift after integration. Therefore, most methods include additional types of sensors (gyroscopes, magnetometers, etc.) (Favre et al., 2008; O’Donovan et al., 2007; Williamson and Andrews, 2001; Ferrari et al., 2010). Combination of accelerometers and gyroscopes is the most frequent approach for angle and trajectory estimation. Based on different sensor configuration, it is possible to perform two-dimensional (2-D) or three- dimensional (3-D) gait analysis. Since the majority of movements are performed in sagittal plane, 2-D analysis describing segment inclinations or joint angle flexions/extensions in sagittal plane can be considered satisfactory. However, some pathologies and gait deformities exhibit movements outside of this plane, and they require 3-D movement analysis in order to quantify the performed movement. For obtaining full 3-D analysis, three-axial accelerometers and gyroscopes are not enough. Namely, we miss information about azimuth for each leg segment, but also for coordination among leg segments (the azimuth is not the same for each leg segment, it can change in time, and we need to have information about their 21 coordination). These problems can not be solved without magnetometers. However, implementing magnetometers and using their output for information about azimuth or correcting angles calculated from gyroscope is not completely reliable since they are prone to signal disturbances coming from the presence of ferromagnetic materials in the near by environment. Since our systems are designed to be used in any environment (clinical settings, at home, outdoors etc.), we chose not to use magnetometers until we develop reliable algorithm to detect signal distractions and neglect corrupted information. Angle estimation methods for lower limb kinematics When we consider gait measurements using IMUs, we can identify two Cartesian coordinate systems. The first system is stationary (fixed) system of the room. This is a global system. The second system is local. It is attached to a body sensor, so that it is movable, (and each leg segment has it own local system). General problem is to determine the relative position between these two coordinate systems, more specific, to determine position vector between centers of two systems, angles between local and global systems (segment angles), and angles between two local systems (joint angles). Since we don’t have to analyze the position vector for angle estimation, it is safe to assume that the centers of both coordinate systems (local and global) coincide (i.e. we translate the centers). One of the ways to define relationships between two systems is to use Euler angles. However, Euler angles can be problematic since there are cases there they are singular. The alternative to Euler angles is to use rotation matrices or quaternions. We use here the rotation matrices because they are most efficient for handling vector rotations and, hence, most appropriate for our algorithms. The application of transformation matrices will be explained in details in chapter 5. Typical requirement for the gait analysis system is to be able to provide angles between axes of two segments (joint angles) or between axes of the segment and fixed axis (segment angles). For clinical applications, the most frequent approach is to estimate the range of motion of the hip, knee and ankle in the sagittal plane by measurement of the flexion and extension at the joint angles However, both segment and joint angles can be 22 defined in all three body planes (sagittal, coronal, and transverse), as shown in Fig. 2.1. Figure 2.1 – Body planes: sagittal, coronal, and transverse. The most common definition of 3-D angles applied in medicine and biomedical engineering refers to 3-D angles of the joints defined by Grood and Suntay (1983). This definition applies only to joint angles and it assigns center of the joint to be the center of the coordinate system. However, there are also other definitions, such as the definition for the Helen Hayes model where segment angles are defined with respect to global coordinate system (by projecting the segments onto three orthogonal planes), and joint angles are defined based on the estimated segment angles (Woltring et al., 1991). Later approach was applied in all our angle estimation methods. Since the majority of gait movements is performed in the sagittal plane, gait is often evaluated by observing angles in this plane, while only few gait deformities require gait assessment in all three planes. Schematic of the system configuration with axes orientation for both coordinate systems is illustrated in Fig 2.2. The figure also shows definition of segment and joint angles in the sagittal plane. 23 x′ z′ y′ x′z′ y′ z y x thighθ shankθ footθ hipθ kneeθ ankleθ Figure 2.2 – Schematic of the system configuration with the coordinate systems. There are several ways to estimate segment or joint angles from inertial sensors. One method to calculate the angle is to double integrate the measured angular acceleration. However, the double integration leads to a pronounced drift (Aminian et al., 2002; Dejnabadi et al., 2005). Several techniques have been presented in literature for minimization of the drift. For example, Morris (1973) identified the beginning and the end of each walking cycle and made the signals at the beginning and the end of the cycle equal. Tong and Granat (1999) applied a low cut high pass filter on the shank and thigh inclination angle signals. However, these methods also removed the static and low-frequency information about the angles and they cannot be applied to real-time processing. The other method for estimation of angles from the measured accelerations is to use accelerometers as inclinometers. In this way, the segment angles are defined as the inclination angles between the segments (sensor) and the vertical, and joint angles by the subtraction of the angles for neighboring body segments. The results are acceptable only if the segment accelerations are small compared to the gravity (Dejnabadi et al., 2006). For lower limb kinematics, this condition if fulfilled only in some periods of gait cycle. 24 Adding Kalman filtering to integration procedure can decrease the drift and provides real-time application, but it requires calibration and data from other sensors (accelerometers, gyroscopes, and magnetic sensors in most cases) for error minimization, as well as noise statistics and good probabilistic models (Luinge and Veltink, 2004; 2005; Cooper et al., 2010). These algorithms can be applied in real- time and seem to give excellent accuracy for motions which exhibit lower accelerations than the leg segments, and which are not exposed to impacts like those at heel contacts. For the lower extremities, the performance of the Kalman filter is considerably reduced when measuring the orientation angles of segments that move fast (Dejnabadi et al., 2006). Inertial sensors that consist of accelerometers, gyroscopes, and magnetometers, along with Kalman filtering, allow a good accuracy for estimation of lower limb angles (Willemsen et al., 1990). However, good accuracy of angle estimation can also be achieved using fewer sensors and much simpler algorithms that are not sensitive to the presence of metals and ferromagnetic materials such as those that comprise magnetometers (O’Donovan et al., 2007). Willemsen et al. (1990) developed a technique to estimate joint angles without integration. This method considers 2-D motion of two neighboring leg segments (which have a hinge joint) and implements a pair of accelerometers on each segment. Based on the signals obtained from a pair of accelerometers on one segment, the acceleration of the joint is obtained in the form of a linear combination of the measured accelerations. This computation is repeated for the other pair of accelerometers located on the neighboring leg segment. By equating the two results for the joint acceleration, the angle between the limb segments is identified. The method requires adequate low-pass filtering, which introduces a delay and to a certain extent hinders the real-time applicability. Further, the accelerometer pairs need to be precisely oriented, so that their axes intersect at the joint, which is very difficult to achieve considering that the human joints are polycentric. Also, the distances between sensors and the joints are required for computation. (Dong et al.,2007). Using accelerometer pairs mounted in suggested configuration in sensor unit, we get reliable 2-D analysis. By replacing one accelerometer sensor with a gyroscope, we can estimate movements in 3-D. Gyroscopes are also prone to drift, but less than accelerometers, and in order to get movement displacement it takes one integration less from angular velocities than it takes from angular accelerations. 25 Trajectory estimation for lower limb kinematics To obtain the trajectory, we have to double integrate the corresponding accelerations measured by our sensors. Since accelerometers measure the accelerations in local coordinate system related to sensor i.e., leg segment, we first have to project accelerations to the axes of global coordinate system. For this we need angles between coordinate systems. After obtaining accelerations in global (fixed) coordinate system, by integrating them we get velocities and coordinates in the room and we can track movement trajectory in the room. As explained in previous chapters, performing numerical integration on acquired accelerations is susceptible to drift. One approach is to use Kalman filter to eliminate the drift (Luinge and Veltink, 2004). We used nonlinear optimization which is explained in chapter 4 which is based on drift removal by polynomials. 26 2-D Gait Analysis: Processing Accelerometer Pairs For this kind of sensor settings, we have developed two algorithms for angle estimation in sagittal plane. The first one is extremely simple to implement, but is limited to cyclic gait pattern, and gait velocities in the range of [0.2, 1.6] m/s. The second is computationally more complex, but has no limitation to movement velocity, and can be implemented for non-cyclic leg movements. 27 3. “Smart” filtering algorithm for 2-D angle estimation based on accelerometers∗ Introduction We have developed an accurate, yet simple method and instrumentation for estimation of absolute segment and joint angles during the gait (assuming kinematics in the sagittal plane) which minimizes the effects of drift. The proposed system is based only on accelerometer sensors, which is advantageous because their calibration is static and less complex than the dynamic calibration required for gyroscopes. Additional motivation for this research was the “bad reputation” of accelerometers due to the pronounced drift. We wanted to investigate if it is possible to use only accelerometers for angle estimations and evaluate the precision of the results. Experimental Section Sensor system The acquisition system that we developed for gait analysis is designed as distributed wireless sensor network. A set of battery powered sensor nodes is placed on subject, one sensor node for each leg segment of both legs. Sensor nodes establish ∗ This chapter is based on: Djurić-Jovičić M., Jovičić N., Popović D.B., “Kinematics of Gait: New Method for Angle Estimation Based on Accelerometers”, Sensors, Nov. 2011, Volume 11(11), pp. 10571-10585. 28 communication with the coordinator node through low power 2.4 GHz wireless communication link. The coordinator node is connected using USB interface to the computer. Wireless communication is bidirectional with coordinator node acting as a master, and sensor nodes as slaves. The coordinator node manages network traffic and USB connection with computer. Data streams from sensor nodes are synchronized and system operates with 100Hz sampling rate. Sensor node is realized as a sandwich structure of processor and sensor board with Li-Ion battery placed between the boards. The compact size design of sensor nodes, with dimensions 70x25x15 mm and 27 grams weight, enables comfortable wearing and does not hinder subject’s movements. Hardware design is based on Texas Instrument’s microcontroller CC2430, which integrates RF front end and 8051 core in the same case. Standard microcontroller peripherals enable interfacing to analog and digital sensors, and different sensor boards can be combined with the same processor board. In configuration used in this research, sensor board comprises two high performance 12-bit digital accelerometers LIS3LV02 (SGS-Thomson Microelectronics, USA). Range of sensors is either ±2g or ±6g, which can be selected in the acquisition software. Accelerometers are aligned to y axes with distance of 55 mm between centers. This configuration requires the clinician only to fix the sensor array along the body segment, approximately at the mid section of lateral side of leg (Fig. 3.1). 29 Figure 3.1. Setup of the sensor system. a) photo of the sensors mounted on the body during the gait analysis, b) schematic of the system configuration with the coordinate systems. Goniometers were attached to the leg segments by using double sided adhesive tape and secured with elastic bands with Velcro endings, mounted over the sensors and around the leg segment. Sensor nodes were placed in custom made tight sensor node-size elastic pockets placed on elastic bands with Velcro at their ends. The custom-designed software, created in CVI (LabWindows, National Instruments, USA), is used for online monitoring and storing of the acquired data. Algorithm The mechanics of importance for the analysis considers two sensors (denoted by S1 and S2), which are mounted on a rigid rod (Fig. 3.2). The distance between the sensors is l. The rod is freely moving with respect to the fixed global coordinate system (O'x'y'z'), shown in Fig. 3.1b and Fig. 3.2. The axis x' of the global coordinate system is walking direction, and the axis y' is vertical. The center of the rod (O) is determined by the position vector )(' 0 tOO r= . To analyze the movement in the sagittal plane, we consider the case when the rod moves in the O'x'y' plane (2-D model). We define the vector l, which connects the 30 centroids of the two sensors. The positions of the accelerometers are 2/01 lrr −= and 2/02 lrr += , respectively. O' r0y' x' z' r2 r1 O y x z S2 S1 l Figure 3.2. Rod with two accelerometers and the coordinate system for analysis of movement in the sagittal plane. Each accelerometer measures the two Cartesian components of the acceleration vector, with respect to the local coordinate system Oxy attached to the rod. The equivalent accelerations measured by the two sensors are glrgra −−=−= 2d d .. 0 .. 2 1 2 1 t (1) and glrgra −+=−= 2d d .. 0 .. 2 2 2 2 t . (2), where g is the gravity acceleration. The difference of the signals from these two sensors is proportional to the amplitude of the vector .. 21 laa −=− . In this way, we cancel out the influence of the movement of the rod centroid and of the gravity, and retain information only about the changes of the vector l. The second derivative of the vector l is yx llll iilul 2 0 2 0 ω−α=φ−φ= &&&&& , (3) where φ is the angle between the axes x and x', ω and α are the absolute angular velocity and angular acceleration of the rod, respectively, xi , yi , and zi are the unit 31 vectors of the x, y, and z axes of the rod, respectively, l=l , yl ill == /0 , and 00 liu ×= z (where zi is the unit vector of the axis of rotation). The difference of the outputs from the accelerometers in the direction along the rod axis ( yaΔ ) is proportional to the square of the angular velocity, and the difference between the outputs from the accelerometers in the perpendicular direction ( xaΔ ) is proportional to the angular acceleration of the segment. The proportionality coefficient is equal to the distance between the centers of the accelerometers. In this way we eliminated the gravity component from the signal, and eliminated the need for precise positioning of the rod on the body segment and calibration of the system. As explained in the introduction, one of the main problems with accelerometers is significant drift after integration, whether the integration is performed numerically or by means of analog integrators. A characteristic example of the drift, which resulted even with carefully calibrated accelerometers, is presented in Fig. 3.3. Figure 3.3. Joint angle (bottom panel) and angular velocity (middle panel) obtained by numerical integration of the measured acceleration of the segment (top panel). Dashed line envelope on the bottom panel is fitted through the points where the knee should be fully extended with zero degree joint angle. However, due to the integration drift, instead of remaining approximately constant, this line has a parabolic shape. We introduce a method for estimation of the joint angles based on digital filtering. In order to explain the method, we use the frequency domain. According to the Laplace transform, the integration in the time domain corresponds to multiplication by s/1 in the frequency domain, and the double integration corresponds to multiplication by 2/1 s , where s is the complex frequency. On the frequency axis (i.e., in the Fourier-transform domain, where ω= js ), this corresponds 32 to multiplication by 22 /1)j/(1 ω−=ω .Hence, we use a second-order low-pass filter, which mimics this multiplication. Further in the chapter, we shall write 2/1 ω− instead of 2/1 s , because we wish to emphasize the fact that this multiplicative term is purely real (although negative). Without the loss of generality, we can assume that the signal ltax /)(Δ is nearly periodic. Hence, it has pronounced spectral components at ,...2,1,/ == iTifi , where T is the stride period. All relevant spectral components of ltax /)(Δ should be in the roll-off region of the filter, where its transfer function is proportional to 2/1 ω− . For example, the transfer function of the second-order Butterworth filter is 12 1)( 0 2 0 2 +⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛ ω+⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛ ω = ss sB , (4) where 00 2 fπ=ω is the cutoff angular frequency. On the imaginary axis, when 0|| ω>>s , we get 2 2 0 2 2 0 2 )( ω ω−=ω≈ s sB . (5) In order to approximate the double integration, we pass the signal through the filter and divide the output by 20ω . The filter is, however, dispersive. Various spectral components have various delays and the filtered signal will only barely resemble the actual function )(tφ . A zero-delay filter can be obtained by bi-directional filtering, using the function filtfilt in Matlab. First, the signal is filtered in the forward direction. Then, the filtered sequence is reversed and run back through the filter. This procedure results in a real and positive transfer function (zero-phase distortion and zero group delay), whose order corresponds to double the filter order. For example, if we use the first-order Butterworth filter, whose transfer function is )1//(1)( 01 +ω= ssB , the result of the bidirectional filtering is the transfer function 220 2 0 /)1)//((1 ωω≈+ωω when 0ω>>ω . Hence, to obtain the 2/1 ω− transfer function, we use the first-order Butterworth filter with filtfilt function and divide the result by 20ω− . Fig. 3.4 shows 33 the normalized spectrum of the angular acceleration, along with the 220 / sω function, and the magnitude of the transfer characteristic of the low-pas filter (Bode plot), low- pass filter combined with a high-pass filter obtained by the filtfilt function. Figure 3.4. Normalized spectrum of the angular acceleration for knee angle (shown in Figure 3.3) and magnitudes of the function 22 0 / sω , transfer function of a second-order Butterworth low-pass filter (lpf), and function of this low-pass filter combined with a high-pass filter (lpf+hpf). The choice of 0f followed the heuristics (Fig. 3.5). If 0f is too low, joint angles exhibit drifting similar to the numerical integration. On the other hand, in order to keep the spectral components in the roll-off region of the filter, the condition Tff gc /10 =< should be fulfilled. If this condition is not respected and 0f is taken to be higher than gcf (where gcf is the gait cycle frequency), one or more spectral components are within the pass band of the filter, where the transfer function of the filter is approximately constant and close to 1. Increasing further the filter bandwidth, i.e., increasing 0f , these components are not affected by the filter. However, their magnitudes are divided by 20ω , so that the level of these components is reduced, and the result is distorted. On the contrary, the magnitudes of the spectral components that are in the roll-off region of the filter are insensitive to the modifications of 0f . Since filtering should replace the integration, all relevant spectral components of the gait should be in the roll-off region of the filter, i.e., well above the cutoff frequency of the filter. By comparing the filter amplitude characteristic, which is the modulus of (4), and 220 ωω , it can be verified that the error between these two functions is less than 34 1 dB for spectral components that are above two times the cutoff frequency of the filter. Similarly, the error is less than 0.5 dB for spectral components above three times the cutoff frequency. Hence, the cutoff frequency 0f of the filters is determined so that the lowest relevant spectral component of the signal is positioned between 02 f and 03 f . Figure 3.5. Influence of 0f on angle estimation: filtering for several cutoff frequencies compared with angle acquired from goniometers. The drift can be additionally reduced if a high-pass filter is used in conjunction with the low-pass filter. The cutoff angular frequency of the high-pass filter should be below 0ω , so that the major role of the low-pass filter is not affected. The order of the high-pass filter can be selected as an additional parameter to help keep the drift under control. As an example, Fig. 3.4 shows the magnitude of the transfer characteristic of the combined low-pas filter and a high-pass filter of 8th order, which is an actual filter used in computations in this algorithm. Since filtering distorts the DC level, we restore this information through the self-calibration in the following way. Before gait initiation, the subject needs to remain standing still (and sensors immobile) for at least two seconds. During this interval, the initial conditions are determined for each pair of accelerometers by using them as inclinometers. The procedure of approximating the double integration can be applied both for reconstruction of absolute angles (the angles between the rod axes and the vertical axis of the fixed coordinate system) and the reconstruction of the joint angles. For example, the knee angle is obtained directly from the difference 35 ltalta xx /)(/)( shankthigh Δ−Δ (6) by performing the bi-directional filtering, and divide the result by 20ω− . An analogue procedure is performed for the ankle angle. Experiments This algorithm was tested on 27 healthy subjects walking on the ground with their natural pace. In order to provide a more systematic validation, we additionally recorded 10 subjects (age: 26 ±1.5 mean±SD) walking with various velocities on treadmills (Life Fitness 9500hr and Panatta Advance Lux 1AD003), whose results are presented in this chapter. Four trials per subject were recorded. Besides walking, recording sequence also included standing still for at least 2 s before and after each walking sequence, which was used for self-calibration and checking. Subjects were walking with various velocities on a treadmill, starting from 0.15 m/s and incremented by 0.05 m/s until 2 m/s. As the reference system for this study, we used flexible goniometers SG110 and SG150 with the joint angle units for signal conditioning (Biometrics, Gwent, UK). Goniometers were mounted on the lateral side of the leg (at the ankle and knee joints) following the instructions of the manufacturer. Simultaneously, the sessions were recorded with a video camera for later analysis. Processing of measured data Based on the recorded accelerometer data, the joint angles were estimated by the proposed algorithm. Also, joint angles recorded by goniometers were computed. The accuracy of our algorithm was evaluated in terms of the root-mean-square error (RMSE) as well as the Pearson’s correlation coefficients (PCCs) between the goniometer results and angles provided by the proposed method. RMSE is expressed in degrees. PCC values range between −1 and 1, where 1 represents the best possible similarity between the two sets of angles (identical shapes). The first and the last stride were excluded from each trial, and comparison between goniometer signals and angles provided by our method was done on the remaining sequence. The data processing was done offline in Matlab 7.5 (Natick, MA, USA). 36 Results Two typical examples for the knee and ankle angles are shown in Fig. 3.6. The error was defined as the difference between the angles obtained by the proposed method and the angles obtained from goniometers. Based on the results obtained from treadmill recordings, the cutoff frequency for the knee angles should be in the range [ gcf /3, gcf /2], where gcf is the gait cycle frequency. The blue area in Fig. 3.6 shows the family of curves estimated by our method when filtered with various frequencies in the range [ gcf /3, gcf /2]. Figure 3.6. Knee and ankle joint angles measured with goniometers and estimated from accelerometers. Thick lines represent angles from goniometers (dashed line) and accelerometers (solid line) for optimal filter frequency. Blue areas show the range of angle values estimated from accelerometers when filtered with various frequencies in the range [ gcf /3, gcf /2]. Fig. 3.7 shows the optimal filtering frequency versus gait cycle frequency. The squares represent optimal points obtained by maximizing Pearson’s correlation coefficient and minimizing RMSE between our results and the angles obtained by goniometers for each walking trial (different gait velocity). The straight lines are obtained by fitting these data. It is obvious from Fig. 3.7 that signals for the estimation of the ankle angles should be filtered with about two times higher cutoff frequency than the signals used for the estimation of the knee angles. These findings are in agreement with the theoretical background from the previous section. 37 Figure 3.7. Optimal filtering frequencies for knee and ankle angles versus gait cycle frequency. The higher cutoff frequency of the filter for the ankle angle can be used because the spectrum of the ankle angle has a very pronounced second harmonic. This is convenient, because the influence of the drift is further suppressed. Using the proposed algorithm, for each subject and each trial, we evaluated the PCC and RMSE values between the angles estimated by our method and goniometer outputs. Fig. 3.8 shows the PCC and RMSE, respectively, as a function of normalized frequency, for various velocities. Although we recorded velocities from 0.15 to 0.2 m/s in the steps of 0.05 m/s, Fig. 3.8 shows results for a subset of the recorded velocities. 38 Figure 3.8. Pearson’s correlation coefficient (PCC) and RMSE between goniometers and estimated angles, for knee and ankle angles, vs. filtering frequency in the range [ gcf /3, gcf /2]. The normalization is with respect to the optimal frequency. Each curve corresponds to different gait velocity (in m/s). As shown in Fig. 3.8, all RMSE curves have broad minima, which are, on average about 1. These results are in accordance with PCC curves, confirming that the optimum normalized cutoff frequency is 1. Cumulative results for all walking trials and all subjects are presented in Figure 3.9. Figure 3.9. Boxplots presenting comparison between joint angles (for knee and ankle angles) calculated from accelerometers and goniometers. Left: Pearson’s correlation coefficient (PCC). Right: root-mean- square error (RMSE). 39 Discussion and Conclusions The presented results are based on a model which assumes that the lower limbs move in the sagittal plane. For healthy subjects, this 2-D model has proven to be sufficient, because the sagittal plane is the plane where the majority of movement takes place. Generally, for clinical applications, the proposed method provides an acceptable accuracy for angles and high correlation coefficients with the measurements obtained from goniometers. In particular, PCCs for the knee angle are higher than 0.97 and RMSE is within 6º for the angle values. Further, our results showed that a 5º RMSE is obtained for walking velocities in the frequency range gcf ∈ [0.35, 1.15]. This corresponds to all velocity curves in Fig. 3.8 except for v=0.15 m/s, 1.8 m/s, and 2.0 m/s (the slowest and fastest recorded walking). Regarding the ankle angles, PCCs are slightly lower and in the range from 0.85 to 0.97, while RMSE is from 2º to 4.7º. For both estimated joint angles (except for the extreme velocities for the knee angle), this error is in the range of 5º mean error limit accepted by the American Medical Association to consider the measurements reliable for the evaluation of movement impairments in a clinical context (Zheng et al., 2005). The accuracy of our simple system is comparable to the accuracy demonstrated in plots presented in (Ferrari et al., 2010), which were obtained by much more complex hardware and software. Joint angles were determined by subtracting the absolute angles of the neighboring leg segments. The error of our method was estimated based on joint angles and includes errors from the two segments. In this way, the total error of the joint angle estimation is different than the error of the absolute angles. Hence, comparison of the absolute angles with a camera system would be more appropriate for validation of the proposed method. However, such a comparison was not possible in our experiments because the treadmill would present a visual obstacle between cameras and markers. Since our main goal was to investigate how our method performs for various gait velocities, the treadmill was the important part of the experiments. Therefore, we selected goniometers as the reference system, which have ±2º accuracy and 1º repeatability (Biometics Datasheet). Although electrogoniometers are prone to errors due to potential misalignment with the femur and tibia segment in 40 the sagittal plane, this does not affect the validation of our method since we secured sensor units to be aligned with goniometer blocks. Skin motion artifacts cause errors to all body-fixed sensors. Sensors placed on thighs are more susceptible to skin and soft tissue related motion, because the majority of femur is concealed by a substantial amount of soft tissue. However, the errors that we report here are much smaller than errors due to rod misalignment in the procedure proposed by Willemsen et al. (1990). Another limit for this algorithm is the velocity of subject’s gait. As it can be seen from Figs. 3.7 and 3.8, if the gait is very slow, the quality of our method decreases. This is due to very low angular accelerations, whose major component comes from the impacts. This suggests that for a very slow walk, e.g., for subjects with high levels of disability, the quality of the angle estimation may not be acceptable. However, for very slow gait, accelerometers can be used as inclinometers and angles can be successfully estimated in this way (Mizuike et al., 2009). In our method, we do not need information about the distances to joint centers or distances between sensor rods placed on different segments, which is one of the benefits of this algorithm. The only request for mounting the sensors is that they follow the segment line (to be aligned with a line connecting adjacent joints, viz. hip and knee, knee and ankle, and along the foot and parallel to the ground). This algorithm could be used not only for level walking, but also for estimation of angles during slope walking, stair climbing, or any other rhythmical (periodic) leg movements. It can also be used for estimation of other segment and joint angles, as long as the movements are in 2-D. However, movements should be fast enough so that the angular acceleration signal is sufficiently above the noise floor. The proposed method is suitable for postprocessing of raw data. However, it can also be included into real-time algorithms to estimate the angles of the leg segments with a delay of one stride. The proposed method is simple and computationally efficient. We have demonstrated that it yields accurate shapes of the ankle and knee angles. The accuracy of the method is sufficient for quick diagnostics of gait, as well as for applications of gait control. 41 4. Nonlinear optimization for drift removal in estimation of gait kinematics based on accelerometers∗ Drift removal by polynomials Inertial sensors are prone to time varying offset and other problems (such as nonlinearities), which introduce an error in the measured signals. Upon integration of the signals, these errors cause drifting of the results (Fig. 3.3). As an example, we consider the translation in 3-D. We assume to have measured (by accelerometers) three components of the linear acceleration in the fixed coordinate system ( zyx aaa ,, ). By one integration with respect to time, we obtain the components of the linear velocity. By one more integration, we obtain the components of the displacement vector. Both the linear velocity and the displacement contain drift. Note that the first derivative of the drift in the velocity gives the error in the accelerometer signals. Similarly, the second derivative of the drift in the displacement, gives the error in the accelerometer signals. We present here a simple technique for reducing the drift. It is based on approximating the drift by a polynomial. Consequently, the error in the accelerometer ∗ This chapter is based on: Djurić-Jovičić M., Jovičić N., Popović D.B., Djordjević A.R., “Nonlinear Optimization for Drift Removal in Estimation of Gait Kinematics Based on Accelerometers”, Journal of Biomechanics, November 2012, Volume 45(16), pp. 2849-2854. 42 signals is also approximated by a polynomial, which is the second derivative of the polynomial that approximates the drift. We consider the time interval of one stride, between two anchoring points. At an anchoring point, the foot is assumed to be fully at rest, touching the ground. Hence, at an anchoring point, we know the following: 1. the linear velocity is zero ( 0=== zyx vvv ), 2. the linear acceleration is zero ( 0=== zyx aaa ), 3. the z coordinate of the foot is zero ( 0=z ). We assume that the drift in a signal can be approximated by a polynomial, which is added to the true signal. Based on the data at the anchoring points, we estimate the coefficients of this polynomial. Thereafter, we subtract the polynomial from the reconstructed signal. Thus we at least partly remove the drift and hence improve the reconstruction. For the x component and the y component of the acceleration, we have sufficient information to use a quadratic polynomial (with three coefficients). For the z component we use a cubic polynomial (with four coefficients). Drift removal from x and y components of acceleration We assume 0=t at the beginning of the stride (the first anchoring point) and Tt = at the end (the second anchoring point). By )(ta we denote the x component of the acceleration reconstructed from the algorithm up to this point. We consider only the x component. The algorithm for the y component is identical to this one. The actual acceleration (which is not known, but which we want to estimate better) is )()()( tetatA −= , where )(te is the error. We approximate the error by a quadratic polynomial, 22102 )()( tctcctpte ++=≈ , where 210 ,, ccc are unknown coefficients. Hence, we have approximately 2 210)()( tctcctatA −−−= . (1.16) At 0=t , 0)0( =A . From (1) we have 0)0(0 ca −= . (1.17) At Tt = , we also have 0)( =TA . Hence, 2 210)(0 TcTccTa −−−= . (1.18) 43 We now integrate (1) on the interval ),0( T . The integral of )(tA is the velocity at Tt = and it must be zero (because Tt = is an anchoring point). The integral of )(ta can be represented as aTtta T =∫ 0 d)( , where a is the average value. Hence, we have 32 0 3 2 2 10 TcTcTcaT −−−= . (1.19) From (2) we have )0(0 ac = , from (3) we have )0()(221 aTaTcTc −=+ , and from (2) and (4) we have )0( 32 2 21 aa TcTc −=+ . Now, we have a system of two equations in terms of 21, cc , whose solution is ( ))()0(222 TaaaTc −−= , ( )aTaa T c 2)()0(322 −+= . Having evaluated )2(2p , we now have a better approximation to the actual acceleration, )()()( 2 tptatA −≈ , and we use this approximation in subsequent computations. Drift removal from z component of acceleration The procedure is similar as for the x and y components. Now, 3 3 2 2103 )()( tctctcctpte +++=≈ , so that we have approximately 3 2 2 210)()( tctctcctatA −−−−= . (1.20) At 0=t , 0)0( =A so that )0(0 ac = . Although we could have evaluated the remaining coefficients analytically, we used a numeric approach. Based on the remaining anchoring conditions, we formed a system of linear equations in terms of 321 ,, ccc and solved it using matrix inversion. First, we integrated )0()( ata − numerically (using time stepping) to obtain the velocity, )(tv , and the displacement, )(tz , with the initial conditions 0)0( =v and 0)0( =z . Note that the acceleration, velocity, and displacement are plagued by the drift. From the condition 0)( =TA , we obtain the first equation, 44 )0()(33 2 21 aTaTcTcTc −=++ . (1.21) From the condition 0)( =TV , where )(tV is the actual velocity, i.e., the integral of )(tA , we obtain )( 432 4 3 3 2 2 1 Tv TcTcTc =++ . (1.22) Finally, from the condition 0)( =TZ , where )(tZ is the integral of )(tV , we obtain )( 20126 5 3 4 2 3 1 Tz TcTcTc =++ . (1.23) By solving the system (6)-(8), we obtain 321 ,, ccc , and, hence, )(3 tp becomes known. Application of nonlinear optimization The drift removal presented at the previous subsection has limited capabilities, because the correcting polynomial has only a few unknown coefficients. Hence, we have only a few degrees of freedom to estimate the error. Keeping the idea of approximating the error in signals by a polynomial, which is a relatively slowly varying function of time, we upgrade the technique from the previous subsection to enable a better approximation of the error, and thus more efficient drift removal, using a more sophisticated technique, based on nonlinear optimization using simplex algorithm (Nelder and Mead, 1965; Lagarias et al., 1998). Our approach assumes that we have two or more independent sensors that directly or indirectly measure the same kinematical quantity. The goal of the algorithm is to estimate the drift and cancel it. We demonstrate the method by using data from arrays of accelerometers. However, the method is applicable to data coming from gyroscopes or other sensors when integration of original signals is required. The method is applicable for real-time data inspection and postprocessing that is of interest for gait analysis. Method and Materials Sensor system with accelerometers Each sensor unit has the form of a short rod on which two accelerometers (A1 and A2) are mounted. The rods are placed along the direction of the long axes of the leg segments at the lateral side (Fig. 4.1). The model discussed here is related to the 45 planar analysis (2-D) in the sagittal plane. The local Cartesian coordinate system O'x'y'z' of a rod has the origin midway between the accelerometers, and two sensors aligned with the y' axis. The z' axis is perpendicular to the sagittal plane. We consider that the Oxy plane of the fixed (global) Cartesian coordinate system (Oxyz) coincides with the plane (O'x'y'). The position vector of O' is yx yxtOO iir 000 )(' +== , where xi and yi are the unit vectors of the Oxyz system. The vector ylil = connects the two sensors (l is the distance between the sensors). The position vectors of the accelerometers are 2/01 lrr −= and 2/02 lrr += . x′ z′ y′ x′z′ y′ thighψ shankψ footθ hipψ kneeψ ankleθ y′ x′ θ (a) (b) Figure 4.1 – (a) Experimental setup showing sensors, markers, and coordinate systems. (b) Sensor unit with two accelerometers. Each accelerometer measures the acceleration vector of the sensor with respect to the fixed system: glrgra −−=−= 2d d .. .. 02 1 2 1 t , glrgra −+=−= 2d d .. .. 02 2 2 2 t , (1) where g is the gravity acceleration. The acceleration vector is measured in terms of its Cartesian components in the system O'x'y'z'. The average signal from the two accelerometers gives ( ) 0..021 2 agraa =−=+ (this signal would be obtained from a single accelerometer located at O'). Following the derivation from Djuric-Jovicic et al. (2011), the 46 difference of the accelerometer signals is ' 2 ' .. 12 yx ll iilaa ω−α−==− , where tddθ=ω is the absolute angular velocity. The difference of the 'y components is proportional to the centripetal (radial) acceleration 2''1'2 ω−=Δ=− laaa yyy . The angular velocity vector is 'ω ziω= . The difference of the 'x components comes from the tangential components of acceleration, α−=Δ=− laaa xxx ''1'2 , where tddω=α is the absolute segment acceleration. Hence, we evaluate the squared angular velocity as l ay '2 Δ−=ω , (2) and the angular acceleration as l ax 'Δ−=α . (3) Algorithm The data determined by Eq. (2) and (3) are used to reconstruct the segment angles. The key novelty of the method is to model the drift by a slowly-varying function of time. We used a cubic spline polynomial, )(tp that is defined by a set of discrete points (knots). The knot time instants are predefined at intervals of 100 ms. We generate a new dataset by modifying the angular acceleration )(tα obtained from Eq. (3) by adding a polynomial, )(tp . We estimated the knot ordinates by nonlinear numerical optimization with the goal to reduce the drift. The knot ordinates are tuned using the simplex optimization procedure (Nelder and Mead, 1965; Lagarias et al., 1998) to minimize a cost function that reflects errors due to the drift in the angular velocity )(tω and in the angle )(tθ . We utilize two sets of data in order to estimate these errors. The first set consists of the squared angular velocity: we use the information about the centripetal acceleration to remove the drift from the angular velocity. The second set consists of “anchoring data”, which are based on the knowledge of the sensor array state at certain intervals (anchoring intervals) during each gait cycle, similarly to Tong et al. (1999) and Schepers et al. (2007). The foot is practically stationary during mid-stance phases, when both heel and metatarsal area are in contact with the ground, and the accelerations (linear and angular) and velocities are close to zero value. The foot elevation above the ground at this time and the foot angle can be assumed to be zero. 47 During the mid-stance phases, the foot angle is calculated by using the accelerometers as an inclinometer because accelerometers are affected only by the gravity since the velocity and accelerations are zero. This angle is approximately zero, the deviation being mainly caused by mounting errors. The foot elevation at these time intervals is assumed to be zero. In reality, the foot elevation at these intervals depends on the placement of the sensors on the foot, which depends on subject’s anatomy and shoe type. The sensors are mounted as low as possible on the medial side of the shoe, several centimeters above the ground. For all tested walking speeds the assumption that the elevation is zero showed to be valid, and it did not affect the accuracy. However, for applications like running and severely impaired gait where there are no foot-flat intervals, this assumption should be taken with care and requires some alterations in the model. The anchoring intervals are identified from the accelerometer signals as intervals, at least 0.1 s long, when all acceleration signals are sufficiently small (below the threshold). During these intervals, θ can be evaluated as ),(atan2 '' xy aa −=θ , (4) where atan2 is the four-quadrant inverse tangent function defined in Matlab program. At the anchoring intervals, for healthy gait patterns and proper sensor alignment, 0≈θ . We integrate the modified angular acceleration, )()()(~ tptt +α=α , to obtain the estimate of the angular velocity, )(~ tω . We integrate the result once again to obtain the estimate of the angle, )(~ tθ . Both integrations are performed from a time instant within an anchoring interval. Hence, the initial condition for )(~ tω is 0, and the initial condition for )(~ tθ is evaluated from Eq. (4). For off-line applications, we integrate up to the next anchoring interval, when we incorporate the next set of anchoring data. For real-time applications, we integrate up to the current time instant, until we run into the next anchoring interval. Thereafter, we restart the algorithm towards the next anchoring interval. The cost function is a weighted sum of two error terms, ( ) ( )∑∑ ωω =θ θ =ω ω θ−θ+ω−ω= N i ii N i ii N w N wc 1 2 1 222 ~~ , (5) 48 where ωw and θw are weights, ωN is the total number of samples from the integration, and θN is the total number of samples during the anchoring intervals. The first term in Eq. (5) reflects the difference between the estimated squared angular velocity, 2~ω , and the “true” 2ω (reference) obtained from Eq. (2). For illustration, Fig. 2 shows the measured and estimated squared angular velocities without the correcting polynomial and after obtaining the optimal polynomial, )(tp . The second term in Eq. (5) reflects the difference between the estimated angle )(~ tθ and the reference from Eq. (4). Figure 4.2 – The square of the angular velocity obtained directly from Eq. 2, and by integrating raw accelerometer data (Eq. 3) and using accelerometer data corrected by the optimized polynomial. The optimization starts from a randomly selected set of knot ordinates. During the optimization, the cost function is minimized by fine-tuning )(tp , thus minimizing the differences between the estimated and the reference values. The optimization yields an estimation of the absolute foot angle, )(~ tθ , which is virtually drift-free. To evaluate the trajectory of O', it is essential to have an accurate foot angle, which is used to transform the components of the linear acceleration from the moving coordinate system to the fixed system. We found heuristically that the optimal weights are 4)srad(1 −ω =w and 2rad100 −θ =w , which approximately equalize the two terms in Eq. (5). The real-time application of the proposed method for the foot angle and the importance of anchoring points are illustrated in Fig. 4.3. Square markers denote the beginnings of the anchoring intervals, and round markers denote the endings. The foot angle obtained from the optical system is shown as reference. 49 Three sets of results obtained from the accelerometers are shown in Fig. 4.3. The first set is obtained by the numerical integration of )(tα , without drift correction. The integration starts within an anchoring interval. The initial conditions are obtained based on the optimized angle from past strides. These curves show a strong drift, which increases as time progresses. Figure 4.3 – Comparison of the angle without drift correction, angle optimized up to the current time (real-time application), angle optimized until the next anchoring point (optimized full stride), and angle from reference optical system. The anchoring intervals are between square and round markers. For the second set, optimization starts within one anchoring interval and terminates within the next anchoring interval, thus spanning one complete stride. Such curves are typically obtained by off-line processing. They are also obtained in real- time applications for the strides preceding the current stride. The third set is obtained from optimizations that start within one anchoring interval and span less than a complete stride. These curves are shown for time increments of 0.1 s. For each curve, anchoring data are used at the beginning of the time interval involved in the optimization, whereas data for 2ω are used during the whole interval. When the termination of the time interval reaches the subsequent anchoring interval, the corresponding anchoring data are also used. Hence, the error of )(~ tθ increases as the time progresses until the anchoring data become available (at the end of the time interval). Thereafter, the error abruptly diminishes and the optimized curve for the current time coincides with the optimized curve for the full stride (Fig. 4.3). We proceed to the evaluation of the trajectory of O'. The components of the acceleration 0a in the moving coordinate system are transformed to the fixed coordinate system as )(cos)()(sin)()( '' ttattata yxx θ+θ= and )(sin)()(cos)()( '' ttattata yxy θ+θ−= , and the influence of the gravity is subtracted 50 from )(tay . We integrate )(tax and )(tay twice. The first integration gives the corresponding components of the linear velocity. The second integration gives the Cartesian coordinates, )(tx and )(ty . At the beginning of the analyzed gait sequence, we set both coordinates to zero. To remove the drift, we perform independent optimization procedures for each coordinate using similar correcting polynomials and identical anchoring intervals as for )(~ tθ . For )(~ tx , during each anchoring interval, the linear acceleration and the linear velocity are zero. For )(~ ty , we have an additional piece of information: the foot touches the ground and, hence, 0)(~ =ty . The cost functions for the coordinates are similar to the one given in Eq. (5). They include terms for the acceleration and velocity at the anchoring intervals, and for )(~ ty , there is another term for the elevation. We have now defined the complete motion of the foot. To “build” the remaining two segments of the leg (shank and thigh) on top of the foot segment, we need to calculate their segment angles defined in Fig. 4.1a. We perform similar optimizations as for the foot angle. The first reference for optimization is the squared angular velocity, as in Eq. (5). The second reference is the angle )(tψ obtained by starting from the foot angle and adding it to the estimation of ankle angle obtained as follows. The foot, shank, and thigh segments are regarded as a kinematic chain and the joint angles are evaluated after Willemsen et al. (1990, 1991) and using the generalization after Dong et al. (2007). The joint angle is calculated from the best-fit rotation matrix that matches the joint-center acceleration calculated from accelerometers on one segment (e.g., foot) with that calculated from the adjoining segment (e.g., shank). To that purpose, we need to know position of accelerometers relative to each joint The optimization yields the absolute shank angle (Fig. 4.4), whence we calculate the ankle angle (Fig. 4.5). The thigh optimization is performed in a similar manner, using the knee angle for anchoring, and resulting in the thigh angle (Fig. 4.4) and the knee angle (Fig. 4.5). This completes the reconstruction of the gait, resulting in the leg stick diagram shown in Fig. 4.6. 51 Figure 4.4 – Segment angles estimated from the accelerometer data and measured by the camera based system (reference). Figure 4.5 – Joint angles estimated from the accelerometer data and measured by the camera based system (reference). Figure 4.6 – Stick diagram of the leg and trajectories of leg joints estimated from the accelerometer data (dashed lines) and measured by the camera based system (solid lines). 52 Procedure for testing the method We recorded the gait pattern in 10 healthy young subjects (age: 31±4 mean±SD) who gave their informed consent. The procedure was approved by the local ethics committee. Subjects were asked to walk for 10 m on level ground at three speeds (slow, normal and fast). The data was collected from five trials at each gait speed. In parallel, the reference results (target trajectory) were measured by the EvArt system (6 cameras, 100 Hz sampling rate) in the Health Technologies Unit, Tecnalia, San Sebastian, Spain, and by the Qualisys system (8 cameras, 200 Hz sampling rate) in the Center for Sensory-Motor Interaction, Aalborg University, Aalborg, Denmark. For both set-ups, reflective markers were placed on the anatomical landmarks of the joints (hip, knee, ankle, heel, and 1st and 5th metatarsal points), and on the two ends of each sensor unit (Fig. 4.1). To ensure visibility, the foot markers were actually placed 2 cm above the SU ending points. The acquisition system used for data collection is a distributed wireless sensor network (Djuric-Jovicic et al., 2011, Jovicic et al., 2012). SU were placed as shown in Fig. 4.1a, one per each leg segment. Each SU comprises two 12-bit digital accelerometers LIS3LV02 (SGS-Thomson Microelectronics, USA) mounted on the print-circuit board at a distance of 55 mm, with the 'y axes aligned. The range of accelerations can be selected by software to g2± and g6± . Results For the estimation of angles, the accuracy of our algorithm was evaluated in terms of the root-mean-square error (RMSE, in degrees) and by the Pearson correlation coefficients (PCCs) between the camera results and angles provided by the proposed method (Fig. 4.7). Based on the results of the t-test for two independent samples, there is no significant difference between the angles obtained from the camera system and the proposed method )05.0( >p , which can also be seen from the boxplots in Fig. 4.7. Regarding the estimation of the stride length, the average error was 2% (min: 0.05%, max: 3.9%), i.e., around 2.5 cm per stride and 10 cm for the whole recorded sequence. 53 Thigh Shank Foot Knee Ankle 1 1.5 2 2.5 3 3.5 4 R M SE [° ] Thigh Shank Foot Knee Ankle 0.84 0.86 0.88 0.9 0.92 0.94 0.96 0.98 1 PC C Figure 4.7 –The RMS error and Pearson correlation coefficients between angles estimated from camera system and the proposed method for the ten subjects in the study. Discussion and Conclusion The proposed method gives accurate estimation of lower limb angles during gait and their trajectories in the sagittal plane. The results for all subjects are in agreement with those of the reference camera system, showing small RMS errors and high PCC (Fig. 4.7). The existing differences between results coming from our method and the reference system could be considered systematic, since their minima and maxima are periodical in the same time interval within the gait cycle. The average shank angle difference is 2.5°. It reaches 5° during mid-stances, falls below 2° during swing 54 phases, and decreases towards heel contact events (round grey markers in Fig. 4.4). Similarly, for the thigh angle, the maximal difference is during mid-swing and terminal swing phases, before heel contacts. The maximal differences occur at each turning point moment when the thigh segment goes forward, and then it starts going backwards in order to prepare the foot for contact with the ground. This change in movement direction is prone to soft tissue artifacts of the thigh, which provoke spikes in accelerometer signals. This could be corrected by the filtering accelerometer data during this gait sub-phase. However, it would provide additional delay and decrease the real-time applicability of the method. Since joint angles are evaluated from the differences of the segment angles for neighbouring segments, the errors for joint angles are directly influenced by the errors in the segment angles. However, for most gait analysis applications, the proposed method provides satisfactory accuracy and reliability. The range of errors is comparable with the errors coming from the movement of the sensor on the leg, and, even more important, comparable with stride- to-stride variability of the gait pattern. The simplex algorithm is robust and one of the most efficient derivative-free local optimizers. In our application, it converged towards acceptable solutions within several hundred iterations which can be implemented in real-time on a medium- performance computer. The optimization could be restarted from another randomly selected initial data, but this was never needed in the cases we examined. Other robust optimization techniques may be applied as well, such as the particle swarm optimization, simulated annealing etc. The drawback of the method is that it requires the positions of the accelerometers with respect to the joints as the reference for ankle and knee angle. However, the results are not very sensitive to the accuracy of these positions. Also, the reference angles are prone to errors due to tissue movement. However, our optimization function involves averaging, which smoothes-out these errors. The current model is limited to 2-D sagittal kinematical analysis. Although the 3-D gait analysis is important for some clinical applications, in many cases the information from 2-D kinematics gives sufficient data for pathologies related to knee and ankle (Sabatini et al., 2005). The optimization method for reduction of errors in estimated kinematical parameters is applicable in all cases where two sets of data describe the same 55 movement. The estimation of the 3-D movement based only on accelerometers, using the present approach, requires hardware with more accelerometers that are mounted in a 3-D arrangement. By using four accelerometers in a spatial arrangement, it is possible to assess the squared angular velocity for all three orthogonal axes; thereby, estimate the 3-D kinematics. The algorithm developed is applicable for real-time and off-line analysis of gait. The method does not need any adaptation with respect to gait velocity or individuality of gait. The method can also be applied for various gait modalities (subjects with different levels and types of gait disability). Since the appearance during the last few years of low-cost miniature gyroscopes on the market, many new methods for motion analysis implement fusion of different sensors trying to eliminate drift and minimize the estimation error. Hence, kinematical analysis based on accelerometers alone might be considered obsolete. However, regarding sensor technology, the state-of-art is such that the current consumption of MEMS accelerometers is three orders of magnitude lower than of rate gyroscopes, which could be a significant advantage of the proposed system compared to conventional IMU systems. Accelerometers are excellent for demonstrating the implementation of the optimization: our method compares favourably with other methods for estimation of segment and joint angles. Furthermore, the method can be easily used with IMUs consisting of one gyroscope and two accelerometers. Pairs of accelerometers can provide reference data for joint angles which can be used to reduce the drift. A single accelerometer can provide the reference angle only for a steady limb segment, which may be used for the foot at anchoring intervals, but not for the shank and thigh segments. However, accelerometer pairs can provide reference angles for all segments. The drift in the integration of the gyroscope signal could be removed using the information about these angles as in the proposed method, providing reliable 2-D and 3-D movement analysis. Extension to 3-D analysis Let us assume that four 3-D accelerometers (an accelerometer quad) are located in a spatial arrangement, as shown below. (The axes of the quad need not be aligned with the axes of the fixed coordinate system.) 56 x y z O d d d A B C The Cartesian coordinates of the accelerometers are )0,0,0(O , )0,0,(dA , )0,,0( dB , and ),0,0( dC . The difference of readings of x-components of accelerometers A and O yields the centripetal acceleration due to rotation about both axes y and z, i.e., d s d aa xzyOxAx = ω+ω=− 22 . (A proof may be needed that squared components of the angular velocity add up.) Similarly, d s d aa yxzOyBy =ω+ω=− 22 and d s d aa zyxOzCz = ω+ω=− 22 . From these three equations, we have 2 2 xzy x sss −+=ω , 2 2 yxz y sss −+=ω , and 2 2 zyx z sss −+=ω . We can now implement the method of drift reduction by optimization (DRO) to reconstruct the components of the angular velocity. In this way, the accelerometer quad can replace three gyroscopes. However, this technique can reduce the drift only in the angular velocities, and not the angles. More information is required for angles. For example, during stance phases, two angles for the foot can be obtained using accelerometers as inclinometers. However, the third angle is missing (same as we had for 3-D gait reconstruction, where we needed to include an angular correction by an educated guess). For the knee and ankle angles, it may be worth seeing if the Willemsen method can be generalized. 57 5. 3-D Gait Analysis from Accelerometers and Gyroscopes Introduction In Chapter 2, we introduced global and local coordinate systems and the need for defined relationships between them. In order to perform the gait analysis, we need to transform various vectors (e.g., position, linear and angular velocity and acceleration) from one coordinate system into the other system. Movements of the local system can be separated as translation of the coordinate origin and rotation of the system. We consider a position vector r, knowing it in one coordinate system (in terms of its Cartesian components), and we want to get its Cartesian components in the other coordinate system. For the purpose of introducing the rotation matrix, we assume that the origins of the two coordinate systems coincide. However, in the general case, the coordinate origins are distinct. Hence, the transformation from one coordinate system into the other one should include the translation between the two origins, in addition to the rotation. Transformation matrices In order to provide transformations between local coordinate system (sensors) and global coordinate system (room), we have to specify the orientation of one set of axes relative to another set. We assume that the two sets have a common origin. One typical approach is to state the direction cosines of one set of axes (primed) relative to the unprimed. In this way, the x axis could be specified by its three direction cosines 58 of angles between the unit vector of the axis under consideration and the respective axes of the unprimed system, cosines of the angles 131211 ,, θθθ in Fig. 5.1. x y z iz iy ix y x z iz ix iy Figure 5.1 – Angles between the axes of the primed and unprimed system. If zyx iii , , are three unit vectors along x, y, z, and zyx ' ,' ,' iii perform the same function in the primed system x’, y’, z’ (Fig. 1), then these direction cosines are defined as: xxxxxx iiiiii ⋅=⋅=∠=θ ''))',((coscos 11 , xyyxyx iiiiii ⋅=⋅=∠=θ ''))',(cos(cos 12 , (5.1) xzzxzx iiiiii ⋅=⋅=∠=θ ''))',(cos(cos 13 yxxyxy iiiiii ⋅=⋅=∠=θ ''))',(cos(cos 21 , and similarly for 312322 cos,cos,cos θθθ etc. The angle IJθ is defined so that the first index refers to the unprimed system and the second index to the primed system. These direction cosines can also be used to express the unit vector in the primed system in terms of the unit vectors of the unprimed system: 131211 cos'cos'cos' θ+θ+θ= zyxx iiii , 232221 cos'cos'cos' θ+θ+θ= zyxy iiii , (5.2) 333231 cos'cos'cos' θ+θ+θ= zyxz iiii . These sets of nine direction cosines ( IJθcos ) completely specify the orientation of the x, y, z axes relative to the x’, y’, z’ axes. In a similar way, we can express unit vectors zyx ' ,' ,' iii in terms of zyx iii , , . 59 Let us now assume that we know the position vector r in the global (unprimed) system. Using the above transformations, we can evaluate the components of the position vector in the primed coordinate system, as zyxzyx zyxzyx '''''' iiiiiir ++=++= (5.3) The coordinates of a point in a given reference frame are the components of the position vector, r, along the primed and unprimed axes of the system, respectively. The primed coordinates are given in terms of x, y, z, as shown in Eq. 5.4. 131211 cos'cos'cos' θ+θ+θ=⋅= zyxx xir 232221 cos'cos'cos' θ+θ+θ=⋅= zyxy yir (5.4) 333231 cos'cos'cos' θ+θ+θ=⋅= zyxz zir The direction cosines can be considered as coordinates describing the instantaneous orientation of the body, relative to a coordinate system fixed in space but with origin in common with the body system. These coordinates are not independent, from nine of them only three coordinates are needed to specify an orientation. The basis vectors in both coordinate systems are orthogonal to each other and have unit magnitude; therefore 0=⋅=⋅=⋅ xzzyyx iiiiii and 1=⋅=⋅=⋅ zzyyxx iiiiii , and similarly for zyx ' ,' ,' iii . ,0coscos 3 1 ' =θθ∑ = lm l lm 'mm ≠ , orthogonality property (5.5) ∑ = =θ3 1 2 1cos l lm , normality property (5.6) Equations 5.5 and 5.6 are together referred to as the orthonormality property. The transition from coordinates fixed in space to coordinates fixed to the rigid body (with a common origin) is accomplished by an orthogonal transformation. The array of transformation quantities (direction cosines) are called the matrix of transformation, which is defined as: ⎥⎥ ⎥ ⎦ ⎤ ⎢⎢ ⎢ ⎣ ⎡ ⎥⎥ ⎥ ⎦ ⎤ ⎢⎢ ⎢ ⎣ ⎡ θθθ θθθ θθθ = ⎥⎥ ⎥ ⎦ ⎤ ⎢⎢ ⎢ ⎣ ⎡ ' ' ' coscoscos coscoscos coscoscos 333231 232221 131211 z y x z y x (5.7) where 60 [ ] ⎥⎥ ⎥ ⎦ ⎤ ⎢⎢ ⎢ ⎣ ⎡ θθθ θθθ θθθ = ⎥⎥ ⎥ ⎦ ⎤ ⎢⎢ ⎢ ⎣ ⎡ = 333231 232221 131211 333231 232221 131211 coscoscos coscoscos coscoscos rrr rrr rrr R (5.8) is the rotation matrix. Following the same procedure, we can express x, y, z, in terms of x’, y’, z’. To obtain ⎥⎥ ⎥ ⎦ ⎤ ⎢⎢ ⎢ ⎣ ⎡ ⎥⎥ ⎥ ⎦ ⎤ ⎢⎢ ⎢ ⎣ ⎡ θθθ θθθ θθθ = ⎥⎥ ⎥ ⎦ ⎤ ⎢⎢ ⎢ ⎣ ⎡ z y x z y x 332313 322212 312111 coscoscos coscoscos coscoscos ' ' ' (5.9) Obviously, [ ] [ ]T1 RR =− where “T” denotes transposed transformation matrix. Therefore, the orthonormality is valid not only for the columns of the matrix [ ]R , but also for the rows of this matrix. 61 3-D gait analysis for SENSY system x′ z′ y′ x′z′ y′ thighθ shankθ footθ hipθ kneeθ ankleθ Figure 5.2 – Schematic of the system configuration with the coordinate systems. Initialization and initial angles Initialization is performed from the initial standing position, and initial angles are obtained from accelerometers by using them as inclinometers. We take 10 samples of the walking sequence at the (selected) beginning (when the subject is assumed to stand still), and average these samples. Hence, we obtain accelerometer readings from x’, y’, z’ which give cosines of the angles between the accelerometer axes and the vertical direction (the z axis of the fixed coordinate system). These angles are xθ , yθ , and zθ (marked angle0footx, angle0footy, and angle0footz, for foot segment), and they are evaluated by taking the inverse cosines. ,/acccos ',',',, gzyxzyx =θ (5.10) where g is the gravity acceleration and “acc” is the signal from accelerometer. Observing Fig. 5.2, ,thigh yθ=θ ,shank yθ=θ °+θ=θ 90foot y . xθ is defined as the angle between x’ and z axis, yθ is the angle between y’ and z axis, and zθ is the angle between z’ and z axis. Also, zyx θθθ ,, correspond to 333231 ,, θθθ from Fig. 5.1. Creating initial rotation matrix Assuming that the coordinate origins of the fixed and the rod coordinate system coincide, the rotation matrix maps the rod coordinates (x’, y’, and z’) into the fixed 62 coordinates (x, y, and z) according to ⎥⎥ ⎥ ⎦ ⎤ ⎢⎢ ⎢ ⎣ ⎡ ⎥⎥ ⎥ ⎦ ⎤ ⎢⎢ ⎢ ⎣ ⎡ = ⎥⎥ ⎥ ⎦ ⎤ ⎢⎢ ⎢ ⎣ ⎡ ' ' ' 333231 232221 131211 z y x rrr rrr rrr z y x , where [ ] ⎥⎥ ⎥ ⎦ ⎤ ⎢⎢ ⎢ ⎣ ⎡ = 333231 232221 131211 rrr rrr rrr R is the rotation matrix. We form initial rotation matrix for each leg segment (rot0shank, rot0thigh, rot0foot) in the following way: Foot segment Knowing xθ , yθ , and zθ , we can readily evaluate the elements of the last row of the rotation matrix, i.e., we can write [ ] ⎥⎥ ⎥ ⎦ ⎤ ⎢⎢ ⎢ ⎣ ⎡ θθθ = zyx rrr rrr coscoscos 232221 131211 footR . The axis along the rod coincides with the y’ axis. We assume that at the initial position, this axis is in the Oxz plane of the fixed coordinate system. In other words, the 'y axis of the foot points approximately in the direction of the x axis of the fixed system, i.e., in the walking direction. Hence, the projection of the 'y axis on the y axis of the fixed system is zero, i.e., 022 =r . This is not entirely correct, because the foot (or any leg segment) is not only in Oxz plane. We calculate the angle between rod and Oxz plane, but we assign it to be the angle between rod and x axis, which creates an error. Further, the mapping from the 'y axis to the x axis is simply yr θ= sin12 , so that at this point we have [ ] ⎥⎥ ⎥ ⎦ ⎤ ⎢⎢ ⎢ ⎣ ⎡ θθθ θ = zyx y rr rr coscoscos 0 sin 2321 1311 footR , where we don’t know the sign of yθsin , we just know it’s yθ−± 2cos1 . This sign will be discussed later. The rotation matrix is orthonormal. From the orthogonality of the first and the second column, we have 0coscos0sin 2111 =θθ+⋅+θ yxy rr . Hence, y yxr θ θθ−= sin coscos 11 . In a similar way, from the second and the third column we 63 obtain y zyr θ θθ−= sin coscos 13 . Now, the rotation matrix becomes [ ] ⎥⎥ ⎥⎥ ⎥⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎢⎢ ⎣ ⎡ θθθ θ θθ−θθ θθ− = zyx y zy y y yx rr coscoscos 0 sin coscos sin sin coscos 2321footR . Just to check, the norm of the second column is yy θ+θ 22 cossin . The norm of the first column is 1231 2 21 2 11 =++ rrr . Knowing 11r and 31r , we have 2 31 2 1121 1 rrr −−±= . Similarly, from the norm of the third column, we obtain 2 33 2 1323 1 rrr −−±= . The element 23r maps 'z into y. The position of the rod is such that 'z is directed towards right of the sagittal plane, whereas the y axis is towards the left of the sagittal plane. Hence, 023 r , the 'x coordinate contributes to y− . In that case, we select the lower sign for 21r . Otherwise, we select the upper sign. We can further rearrange the rotation matrix. z y zyr θ−⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛ θ θθ−±= 2 2 23 cossin coscos 1 y zyzyy θ θθ−θθ−θ±= sin cossincoscossin 22222 y zy θ θ−θ±= sin cossin 22 y zy θ θ−θ−±= sin coscos1 22 . Since 1coscoscos 222 =θ+θ+θ zyx , we have y x y x y xr θ θ±=θ θ±=θ θ±= sin cos sin cos sin cos2 23 because always 0sin >θ y . except when the rod is aligned with the z axis (when 0sin =θz ). However, we have established that 64 023 =π > =θ 0,0,atan 0,0,atan 0,0,0 0,0, 2 0,0, 2 0,atan ba a b ba a b ba ba ba a a b (5.21) One example of the calculated 3-D angles, compared with the corresponding angles obtained from motion capture system as the reference are shown in Fig. 5.4. The comparison is not entirely correct since the markers for camera system were placed on anatomical landmark (hip, knee, ankle, heel and metatarsals). Although 70 sensors were carefully placed in order to minimize the alignment error (as it can be seen from the results for sagittal plane and beginning intervals for all 3 planes), there is a different kind of error which couldn’t be decreased – soft tissue artifacts. These problems are not present on the places where reflective markers were attached while the placements of sensor nodes were highly exposed to them. More accurate comparison would include markers placed on both ends of each sensor node, where they would exhibit the same kind of movements as the SENSY sensors. 71 Fi gu re 5 .4 – C om pa ris on o f t he se gm en t a ng le e st im at io n fo r a ll an gl es in a ll th re e pl an es : s ag itt al , c or on al , a nd tr an sv er sa l. R ed li ne sh ow s S EN SY o ut pu ts , b la ck li ne sh ow s o ut pu ts o f t he c am er a sy st em . 72 Stride segmentation Stride segmentation is usually performed based on FSRs and gyroscopes from the foot unit. However, using FSRs for this was susceptible to errors in the analysis of some forms of pathological gait. Therefore, we switched to decision-making based on gyroscope sensors only. We take all three axes from gyro-sensor from the foot unit and calculate “gyro sum” in the following way: we eliminate DC component from each gyro axis on the foot, and then we summarize their absolute values which are powered on “p”, as shown in equation below. ∑ = −= zyxi ii gyromeangyrogyrosuma ,, p)( (5.22) The reason for making this kind of modification in the sum of gyro signals is to emphasize higher values and decrease lower values. The “p” value is selected to be 4, which was heuristically determined after inspecting recorded data from patients with gait impairment. After this, this “sum” is normalized to its maxima. It is shown in Fig. 4 with red line, and marked as “gyrosum a” (“a” is from analogue). Then we apply moving average filter (smoothing 30 points) and normalize to its maxima. This signal (“gyrosum f”) is shown with thick blue line in Fig. 4. Based on this filtered signal we make its binary form (“gyrosum b”) by applying heuristically determined threshold (we set this value to 0.002 but it could be changeable). “gyrosum b” is set to 0 when the foot is moving, and set to 1 when the movement is below the threshold value. This signal is shown with thin black line in Fig. 5.5. After this, we search for the moment when the foot movement (“gyrosum f”) is minimal within every interval where the “gyrosum b” is positive. These moments are shown with cyan square markers in Fig. 5.5. FSRs signals are also shown in Fig. 6. FSRs are drawn with green line, providing reference for recognition of gait phases (positive – stance phase, zero – swing phase). 73 5 10 15 20 25 30 0 0.2 0.4 0.6 0.8 1 t[s] Fo ot S en so rs gyrosuma gyrosumb gyrosumf fsrsum Figure 5.5 – Stride segmentation method. These cyan square markers represent moments within stride where the foot is still, and these moments are used for stride segmentation and resetting the calculated stride length after each stride. They are also used for fitting the polynomial function in order to eliminate drift from angles, as it will be described in the following section. Eliminating drifts from angles This correction was needed for longer files recorded at the Rehabilitation clinic, when the calibration was not appropriate. Additional drift minimization was performed for absolute angles only. Panels in Fig. 5.6 show, from top to bottom, segment angles for the foot, shank, and thigh) and vertical ground reaction force. Red line represents the original signal, and the green line shows the corrected signal after drift removal. Foot angles are corrected by fitting the polynomial through the points marked as the moments inside the stance interval where the inertial sensors on the foot register the minimal movement (black square markers in Fig. 5.6). 74 2 4 6 8 10 12 14 -60 -40 -20 0 20 F O O T [d eg ] Segment Angles - Drift Removal original removed drift 2 4 6 8 10 12 14 -40 -20 0 20 40 S H A N K [d eg ] original removed drift 2 4 6 8 10 12 14 -40 -20 0 20 T H IG H [d eg ] original semi-removed drift removed drift 2 4 6 8 10 12 14 0 0.5 1 G R F t [s] Figure 5.6 – Drift elimination for segment angles. Red line represents the original signal, and the green line shows the corrected signal. Black square markers show moments in the stance phase when the foot has minimal movements. Black square markers show the selected moment inside the stance interval where the inertial sensors on the foot register the minimal movement. These moments are used for fitting the polynomial through them and removal of the baseline from the foot angle. Drift removal from the shank angle is performed by polynomial fitting through the points between local maxima and minima for each stride (red square markers). Regarding the fitting through the thigh angles, drift removal has one additional step in comparison to the shank segment. The reason for this is that the thigh segment is very susceptible to soft tissue artifacts which impair the quality of angle estimation. Therefore, the first step in drift removal is fitting through the red square markers, which mark the same moments as black markers for the foot segment. This step is presented with yellow line in Fig. 5.6. After this, we perform fitting through the yellow square markers, using the same method as we did for the shank angle. 75 However, this kind of drift removal can damage the original baseline. This is prevented by calculating the baseline from the second stride (the reference baseline) from the beginning of the sequence, and shifting the estimated angles after the drift removal to this reference baseline. We selected the second stride as the reference baseline because we needed the stride that is both representative in its movement magnitudes and closer to the beginning of the recording i.e., the still period so that the drift would not be very much expressed. Since the first stride usually has smaller amplitudes i.e., smaller movements, flexions and extensions (as shown in Fig. 5.6), its baseline can be significantly different from the other baselines which is why it was discarded as the candidate for the “reference baseline” and the following stride was selected. Discussion and Conclusion There are several limitations that need to be considered. First, current placement of the sensors (lateral side of each leg segment) is extremely exposed to soft tissue artifacts, which makes segment reconstruction more difficult. It especially creates problems for reconstruction of the movement in coronal plane, which can be important for some forms of pathological gait. This agrees with Stagni et al.(2005) where it was showed that the biggest error provoked by soft tissue artifact comes from the movements performed in coronal plane i.e., for medio-lateral movements. One solution is to move sensors placement to the frontal side of leg segments, which is less susceptible to soft tissue artifacts then our original position. However, this was not applied for SENSY for two reasons: the first one is that most of our clinical studies included other types of sensors (EMG recordings) which had to be positioned in the frontal part and there was no place left for sensor nodes, and the second reason is that validation of angle estimation for pathological gait was performed with electrogoniometers that needed to be placed on lateral side (except for the foot segment which could be placed medially but due to diversity of patients footwear we couldn’t get really comparable results and real foot angles). Since the accuracy of angles in coronal plane is not satisfactory, we can not provide accurate information about the excursions for circumduction movements. Furthermore, current SENSY system does not have information about the azimuth which is needed for the estimation about the initial foot opening. Without that 76 there is a problem with trajectory reconstruction in 3-D which is artificially solved by entering approximate information in the program. This problem can be solved by using magnetometers. This may help reconstructing positions of two legs (though it would not be sufficient without using additional information, such as controlling the distance between the hips). However, the problem with magnetometers is that they are sensitive to the presence of ferromagnetic materials and metals and need to be considered with care and in controlled environment. These applications are not really suitable for clinical and home environment since it couldn’t be requested from patients to take so much care about the surrounding area. Furthermore, patients often have metal canes or walkarounds that could hinder sensor’s data. There is another limitation regarding the system hardware: if the sensors are not properly calibrated, the system is prone to higher errors. It means either to have reliable and more frequent calibration or better sensors. We planned to switch to digital inertial sensors of the new generation which would provide better signal quality, both its accuracy and resolution. However, this was not done. Drift removal was performed for signals from foot accelerometers and already calculated angles. However, it was not applied for trajectory estimation. This is not the problem if the sensors are well calibrated and there is a sequence in the beginning of the recorded file where the subject is standing still. These moments should be taken for estimation of initial conditions. Based on the information from our collaborators in two clinics, they prefer to observe angles and trajectory reconstruction is interesting for them only for illustrating if there is some kind of gait disturbance and the measure of that disturbance they like to read from estimated angles which is why we didn’t correct the 3-D trajectory estimation. There were two main protocols tested with SENSY: 1) 10 m long straight pathway for gait assessment of patients with stroke and 2) at least 50 m long complicated pathway with turns and turn overs for patients with Parkinson’s disease. Inside the SENSY algorithm, there is a block for performing integration reset after each stride (in the middle of the stance phase, when the foot is still). This block can be activated or deactivated, depending on the application we want to use. For straight- line pathway, it was shown that 3-D trajectory reconstruction works better if this reset is performed. However, for application in PD patients this reset disables correct 77 information about the turns and reconstruction of their complicated pathway which is crucial for interpreting what caused patient’s disturbances and freezing episodes. Besides these hardware modifications, some future work if we decide to get back to this could include additional drift removal of gyro data, for example by comparing average values of absolute angles obtained by gyroscopes and accelerometers (used as inclinometers) or filtering. 78 6. Processing of signals from rigid double pendulum: comparison of angle estimation methods In this chapter, we compare the angle estimation methods described in chapters 3, 4, and 5. Those methods were tested and validated in real working conditions, sensors were mounted on subjects and tested for different gait speeds. However, walking can never be completely in 2-D and there are always some segment rotations and skin and tissue motion artifacts. Precise estimation of segment and joint angles is always disturbed by these artifacts, and so is the validation of those methods. Therefore, we performed testing and comparison of these methods for angle estimation of a double pendulum composed of two rigid segments. The pendulum has a rotary joint to a fixed support and another rotary joint connecting the two segments (Figure 6.1). This configuration allows pendulum movements to be considered analog to a leg model in 2-D. Such double rigid pendulum allowed movements in one plane (that would correspond to the sagittal plane for human movements) without the risk of having crosstalk between the axes. In this chapter, we additionally point to some of the problems of the numerical analysis applied on data from inertial sensors, i.e., integration vs. derivation, by comparing angular accelerations (α), angular velocities (ω), and angles (θ) obtained from each other by calculus. Namely, errors in measurement or signal noise, however small, are accumulated from point to point. This leads to drift, or an ever-increasing difference between where the system thinks it is located, and the actual location. 79 tdd tdd ∫ ∫ Experiment setup Inertial sensors were carefully placed along the pendulum segments, connected by mechanical joints equipped with optical incremental rotary encoders, as illustrated in Fig. 6.1. The upper segment was attached to a fixed base by its upper joint (Encoder 1). The lower segment was attached to the upper segment by their common joint (Encoder 2). Both pendulum segments were able to move freely in 2-D. The mechanical properties, size, and weights of encoders and inertial sensors did not hinder the pendulum movements. Three sensor units with inertial sensors were placed along the upper segment, and two sensor units on the lower segment. The first (highest) sensor unit was placed especially carefully so that its accelerometer position matches with the encoder’s center of rotation. Round black markers show the actual location of accelerometer sensors within sensor unit, while square black markers show the location of gyroscopes. The axes orientations are as shown in Fig. 6.1. Signals measured from inertial sensors i.e., accelerometers and gyroscopes are marked as aij and gij, where i marks the pendulum segment (1-upper segment, 2-lower segment), and j marks the order of the sensor unit, from top to bottom. Joint angles are marked by iθ and segment angles by iφ (with respect to vertical). 80 11 φ=θ x z y 2φ 2θ 1111, ga 1212 , ga 1313 , ga 2121, ga 2222 , ga Figure 6.1– Double rigid pendulum with encoders placed as joints and inertial sensors placed along segments. Accelerometer positions in sensor units are denoted by round markers and gyroscope positions by square markers. The following recording protocol was performed with the double pendulum: Static phase with transients: 1. pendulum was vertical and still (neutral position), 2. pendulum’s lower segment was lifted to horizontal position and remained still, 3. pendulum’s upper segment was lifted so that the lower segment returns to neutral position and it was held still, Oscillating phase: 4. pendulum was released to swing. These phases can be seen from the angles measured by the encoders as shown in Fig. 6.2. 81 0 20 40 60 80 100 120 140 160 180 200 -1.5 -1 -0.5 0 0.5 1 1.5 t [s] A ng le s [r ad ] θ1 θ2 1. 2. 3. 4. Figure 6.2 – Joint angles based on encoders, complete recorded sequence. 82 0 20 40 60 80 100 120 140 160 180 200 -1 0 1 θ 1 [r ad ] 0 20 40 60 80 100 120 140 160 180 200 -1 0 1 θ 2 [r ad ] t [s] 0 20 40 60 80 100 120 140 160 180 200 -1 0 1 t [s] φ 2 [r ad ] Figure 6.3 – Joint and segment angles based on encoders, complete recorded sequence. The segment angle for the lower segment is obtained by adding the angles obtained from encoders 1 and 2. Results for the upper segment Basic methods for angle estimation The first angle (the angle of the first segment with respect to the vertical direction, θ1) was evaluated by the inclinometer based on acceleration components of the first IMU (a11) and by integrating the gyroscope (g11), starting from 2 s. The results are compared with the encoder, which is assumed to be the reference. 83 15 20 25 30 35 40 45 50 55 60 65 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 t [s] A ng le [ ra d] encoder θ1 inclinometer int(g11) Figure 6.4 – Angle of the upper segment obtained from encoder, inclinometer, and integrated gyroscope. A slow drift can be seen in the integrated signal over the whole recorded interval, as shown in Fig. 6.5. 0 20 40 60 80 100 120 140 160 180 200 -1 -0.5 0 0.5 1 t [s] A ng le [ ra d] encoder θ1 inclinometer int(g11) Figure 6.5 - Angle of the upper segment obtained from encoder, inclinometer, and integrated gyroscope (exhibiting drift), complete recorded sequence. The angle labeled “inclinometer” is obtained by using the accelerometers as inclinometers by applying “atan2”, as introduced in chapter 4. This is a good 84 assumption when the pendulum is still, but the results are very inaccurate when the pendulum is in motion. The upper panel of Fig. 6.6 shows the differences between the angle from encoder (reference) and the angle obtained from accelerometer which is used as an inclinometer. The differences between integrated gyroscope data and the angle from encoder can be seen in Fig. 6.6 (lower panel). For IMU#1, the position of accelerometer a11 coincides with rotation axis of the upper segment. For this reason, a11 provides good results when used as inclinometer, and neither centripetal nor tangential accelerations hinder its angles estimation. This accelerometer would be analog to accelerometer positioned on the ankle, for gait analysis. 0 20 40 60 80 100 120 140 160 180 200 -0.1 0 0.1 0.2 E nc od er -I nc li n. [ ra d] 0 20 40 60 80 100 120 140 160 180 200 -0.1 0 0.1 0.2 t [s] E nc od er -i nt (g 11 ). [ ra d] Figure 6.6 – Differences in the angles estimated by inclinometer and integrated gyroscope, compared to the encoder angle. Angular velocity and angular acceleration are shown in Fig. 6.7 and 6.8, obtained in different ways (directly from sensors or by differentiation, as appropriate): 85 15 20 25 30 35 40 45 50 55 60 65 -5 -4 -3 -2 -1 0 1 2 3 4 5 t [s] ω[r ad /s ] diff(θ1) diff(Inclin.) g11 Figure 6.7 – Angular velocity obtained as derivatives of encoder and inclinometer angle, and directly from gyroscope. 15 20 25 30 35 40 45 50 55 60 65 -400 -300 -200 -100 0 100 200 300 400 t [s] α [ ra d/ s2 ] diff2(θ1) diff2(inclin.) diff(g11) Figure 6.8 – Angular acceleration obtained as double derivatives of encoder and inclinometer angle, and as derivative from gyroscope. The differentiation 1produces a lot of noise in all cases, so that filtering is required to obtain meaningful results2. Clearly, the worst results are obtained by 1 The input signals are digital. Each signal is a train of samples, at a sampling rate of 100 Hz. Each signal is also quantized. The differentiation is performed here by taking first-order differences. When 86 double differentiation of the angle estimated by inclinometer. The reason for this is the quantization of the signals. Double differentiation of encoder angle also provides useless results. As explained in chapter 2, the integration introduces drift. Comparison of drifts after integration obtained from the three gyroscopes located on the first segment is shown in Fig. 6.9 and 6.10. 15 20 25 30 35 40 45 50 55 60 65 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 t [s] A ng le ( φ 1) [ ra d] int(g11) int(g12) int(g13) Figure 6.9 – Segment angles obtained by integrating signals from gyroscopes positioned along the segment, sequence at the beginning of recorded file (transition state) the signal is slowly varying, we subtract almost equal numbers. Hence, the quantization introduces large errors in the differences. 2 However, the filtering cannot make a total improvement of the results. To reduce the additive noise, the filtering is low-pass, but it not only decreases the rapidly-varying noise, but also distorts even the actual signal (changes the form and introduces a delay). 87 165 170 175 180 185 190 195 -0.4 -0.3 -0.2 -0.1 0 0.1 t [s] A ng le ( φ 1) [ ra d] int(g11) int(g12) int(g13) Figure 6.10 – Segment angles obtained by integrating signals from gyroscopes positioned along the segment, sequence at the end of recorded file. Different sensors exhibit different amount of drift. Engle estimation based on accelerometer pairs Using pairs of accelerometers, we obtained angular acceleration of the upper segment, as shown in Fig. 6.11. Angular acceleration is also calculated by differentiation of gyroscope signal, which exhibits noise. 15 20 25 30 35 40 45 50 55 60 65 -25 -20 -15 -10 -5 0 5 10 15 20 25 g p t [s] α [ ra d/ s2 ] diff(g11) a11 & a12 a11 & a13 a12 & a13 Figure 6.11 – Angular acceleration obtained from different accelerometer pairs from upper segment, and by differentiation of gyroscope signal. 88 As shown in Fig. 6.11, differentiation provides worse results for small values (static phase) than for oscillating phase where the signals are have larger magnitudes. Results from the oscillating phase show excellent agreement with the reference angular accelerations, while results in the static phase exhibit large noise. By integrating the angular acceleration obtained from the accelerometers, the results are obtained for the angular velocity and angle, exhibiting a large drift (Fig. 6.12 and 6.13). 15 20 25 30 35 40 45 50 55 60 65 -10 -5 0 5 t [s] ω [ ra d/ s] a11 & a12 a11 & a13 a12 & a13 g11 Figure 6.12 – Angular velocity obtained by integrating from different accelerometer pairs from upper segment, and directly from gyroscope signal. Different accelerometer combinations exhibit different amount of drift. 89 15 20 25 30 35 40 45 50 55 60 65 -200 -150 -100 -50 0 50 g p t [s] A ng le [ ra d] a11 & a12 a11 & a13 a12 & a13 g11 Figure 6.13 – Angles of the upper segment obtained by double integration of different accelerometer pairs and by integration of gyroscope signal. Different accelerometer combinations exhibit different amount of drift. Willemsen method The angle obtained by the Willemsen approach is shown in Fig. 6.14. 15 20 25 30 35 40 45 50 55 60 65 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 t [s] A ng le [ ra d] a12 & a13 θ1 Figure 6.14 – Angle obtained by Willemsen method from accelerometers a12 and a13, compared to reference encoder angle. 90 The Willemsen method provides estimation of the joint angles, but it can be applied to evaluate the absolute angle of the upper segment since this segment is attached to a stationary base (so that the absolute angle equals the joint angle). By differentiating the angle, the angular velocity and angular acceleration are obtained, as shown in Fig. 6.15 and 6.16. 15 20 25 30 35 40 45 50 55 60 65 -6 -4 -2 0 2 4 6 8 t [s] ω [ ra d/ s] a12 & a13 g11 Figure 6.15 – Angular velocity obtained by differentiating the angles obtained by Willemsen method, compared to gyroscope signal (reference). 91 15 20 25 30 35 40 45 50 55 60 65 -1000 -800 -600 -400 -200 0 200 400 600 800 1000 t [s] α [ ra d/ s] a12 & a13, Willemsen a12 & a13 Figure 6.16 – Angular acceleration obtained by double differentiating the angles obtained by Willemsen method, compared to angular acceleration directly from a12 and a13 (reference). Generally, Willemsen’s results for the angle are very noisy due to subtractions, and the subsequent differentiations make the results useless. Results of a differentiation could be improved by applying some more sophisticated form of numerical differentiation (e.g., interpolation by a polynomial and differentiation of the polynomial) or by applying low pass filtering. However, the filtering would suppress the noise, but distort the signals. “Smart” filtering method (integration in frequency domain) The following figures show results obtained by the “smart” filtering algorithm (SFA) described in chapter 3. Anchoring to the encoder results is at 65 s, where the signal is nearly periodic. This anchoring is necessary since filtering distorts the dc value and we have to anchor results of SFA to dc value of the angle, as explained in Chapter 3. 92 15 20 25 30 35 40 45 50 55 60 65 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 t [s] A ng le [ ra d] a12 & a13 θ1 Figure 6.17 – Angle obtained by smart filtering, applied on accelerometer pair, compared to encoder angle (black line), during the static sequence from the beginning of the recorded sequence. 175 180 185 190 195 -0.1 -0.05 0 0.05 0.1 0.15 0.2 g p ( g) t [s] A ng le [ ra d] a12 & a13 θ1 Figure 6.18 – Angle obtained by smart filtering from accelerometer pair, compared to encoder angle (black line), ending sequence from the recorded file. As discussed in chapter 3, where we introduced the SF method, this method is only applicable to movements that are quasi-periodic or periodic. Therefore, angles estimated by this method follow the reference angle in shape and amplitude for pendulum oscillations, while for transient state they show significant errors compared to encoder angles. 93 Nonlinear optimization for drift removal The optimization method, presented in chapter 4, was performed in windows of 20 s, with anchoring to the encoder reading at 10 s intervals, up to 70 s. Shown below is the squared angular velocity obtained from the accelerometers. The reading should be zero in the time intervals when the segment is still. This is true in the calibration intervals, but not when the pendulum is on the chair. The most probable cause is nonlinearities of the accelerometers. 15 20 25 30 35 40 45 50 55 60 65 -4 -2 0 2 4 6 8 10 12 14 t [s] ω2 [ (r ad /s )2 ] a12 & a13 preconditioned g11 2 Figure 6.19 – Quadratic angular velocity obtained from y axis of accelerometer pair (red line), preconditioned (all negative values are clipped with 5pt moving average filter), and directly from gyroscope The negative quadratic angular velocity phenomenon can be explained by mechanical vibrations (but we assume that the segment is rigid), or by slightly unsynchronized sensor sampling, etc. The biggest error is synchronized with pendulum transitions between different positions, around 25 s and 52 s. For the time interval 0-20 s, the input data for optimization is shown in Figure 6.20. 94 2 4 6 8 10 12 14 16 18 20 -10 -8 -6 -4 -2 0 2 In pu t da ta f or o pt im iz at io n t [s] ay ax α ω2 60 62 64 66 68 70 72 74 76 78 80 -20 -15 -10 -5 0 5 10 15 20 In pu t da ta f or o pt im iz at io n t [s] ay ax α ω2 Figure 6.20 – Input data for optimization, linear accelerations, angular acceleration, and quadratic angular velocity obtained from accelerometers. Upper panel: sequence at the beginning of file (static phase), lower panel: sequence from the oscillating phase The results obtained from the optimization are shown in Fig. 6.21. 95 0 2 4 6 8 10 12 14 16 18 20 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 t [s] ω2 [ (r ad /s )2 ] optimized from ay 60 62 64 66 68 70 72 74 76 78 80 0 2 4 6 8 10 12 t [s] ω2 [ (r ad /s )2 ] optimized from ay Figure 6.21 – Quadratic angular velocity: directly from the subtraction of x axes of the accelerometer pair (blue line), after optimization (red line) 96 0 2 4 6 8 10 12 14 16 18 20 -0.04 -0.03 -0.02 -0.01 0 0.01 0.02 0.03 0.04 t [s] ω [ ra d/ s] optimized g11 60 62 64 66 68 70 72 74 76 78 80 -4 -3 -2 -1 0 1 2 3 4 g y p t [s] ω [ ra d/ s] optimized g11 Figure 6.22 – Angular velocity after optimization (red line), compared to the gyroscope signal (blue line). 97 0 2 4 6 8 10 12 14 16 18 20 -8 -6 -4 -2 0 2 4 6 x 10 -3 Angle from optimization t [s] A ng le [ ra d] inclinometer optimized encoder (θ1) 60 62 64 66 68 70 72 74 76 78 80 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 t [s] A ng le [ ra d] optimized encoder (θ1) inclinometer Figure 6.23 – Angle after optimization, compared to the angle obtained by inclinometer and from encoder. For the whole sequence 0-80 s, the reconstructed angle is shown below: 98 0 10 20 30 40 50 60 70 80 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 t [s] A ng le [ ra d] optimized encoder (θ1) Figure 6.24 – Angle after optimization (red line), compared to the angle from encoder (blue line). The oscillations in the interval when the pendulum is in static phase, come from accelerometer nonlinearities3. The problem is because the accelerometers readings indicate that the pendulum is in motion, whereas it is actually completely still. Further fine tuning of various parameters used in the optimization could improve the reconstruction. Results for the lower segment Basic angle estimation methods The second angle (the angle of the second segment with respect to the vertical direction) was evaluated by integrating the gyroscopes (starting from 2 s). The results are compared with the encoders, which are assumed to be the reference. 3 The accelerometers are calibrated from the recorded sequence at 10 s and at 50 s, when the pendulum i still at two positions. The resulting oscillations at these two positions are relatively small, as expected. However, the oscillations at the intermediate still position (30-40 s) are much larger. If the accelerometers were totally linear, the two-point calibration would be sufficient for all positions. However, due to the nonlinearities, the accelerometers are not accurately calibrated in-between the calibration points. Due to different characteristics of the two accelerometers in a pair, the resulting error is much larger away from the calibration points, as is obvious at the intermediate still position. 99 15 20 25 30 35 40 45 50 55 60 65 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5 t [s] A ng le [ ra d] φ2 g21 g22 Figure 6.25 – Angle of the lower segment, obtained directly from encoders (red line) and by integrating g21 and g22 (green and blue line). Angle estimation from accelerometer pairs Using pairs of accelerometers, the following angular acceleration is obtained: 15 20 25 30 35 40 45 50 55 60 65 -50 0 50 t [s] α [ ra d/ s2 ] a21 & a22 g21 Figure 6.26 – Angular acceleration of the lower segment, measured by accelerometer pair a21 and a22, and by differentiation of gyroscope signal g21. By integrating the angular acceleration obtained from the accelerometers, the following results are obtained for the angular velocity and angle, again exhibiting a large drift: 100 15 20 25 30 35 40 45 50 55 60 65 -15 -10 -5 0 5 10 t [s] ω [ ra d/ s] int(a21 & a22) g21 diff(φ2) Figure 6.27 – Angular velocity obtained from accelerometer pair, from differentiation of encoder signal, and gyroscope (reference). 15 20 25 30 35 40 45 50 55 60 65 -80 -70 -60 -50 -40 -30 -20 -10 0 10 t [s] A ng le [ ra d] int2(a21 & a22) int(g21) Figure 6.28 – Angle of the lower segment obtained by double integration of angular acceleration obtained from accelerometer pair, and by integrating gyroscope signal. Willemsen method The following figure shows the angle obtained by the Willemsen approach: 101 15 20 25 30 35 40 45 50 55 60 65 -1.5 -1 -0.5 0 0.5 1 1.5 t [s] A ng le [ ra d] (a12 & a13)&(a21&a22) encoder (θ2) Figure 6.29 – Joint angle (angle between upper and lower pendulum segment), obtained from accelerometer pairs from upper and lower segment by Willemsen method, compared to encoder angle (reference). Smart filtering (integration in frequency domain) The following figures show results obtained by the filtfilt algorithm. Anchoring to encoder results is at 65 s, as explained for the upper segment. 102 15 20 25 30 35 40 45 50 55 60 65 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 t [s] A ng le [ ra d] a21&a22 φ2 80 85 90 95 -1 -0.5 0 0.5 1 t [s] A ng le [ ra d] a21&a22 φ2 Figure 6.30 – Angle of the lower segment obtained by smart filtering (red line) and from encoders (black line). Upper panel: sequence from the beginning of file – transition state. Lower panel: zoomed sequence from the middle of the signal (stationary state). Nonlinear optimization for drift removal Shown below is the squared angular velocity obtained from the accelerometers. Note that it is one order of magnitude larger than for the first segment. 103 15 20 25 30 35 40 45 50 55 60 65 0 20 40 60 80 100 120 t [s] ω2 [ (r ad /s )2 ] a21 & a22 preconditioned g21 Figure 6.31 – Quadratic angular velocity from directly from accelerometer pairs, from preconditioned accelerometer pairs, and from gyroscope signal. The reconstructed angle of the second segment is shown in Fig. 6.32. 0 10 20 30 40 50 60 70 80 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 t [s] A ng le [ ra d] optimized encoder (φ2) Figure 6.32 – Angle of the lower segment obtained from accelerometer pair by optimization method, and directly from encoder (reference). 104 Conclusion In this chapter, we compared all angle estimation methods investigated in previous chapters, and demonstrated their advantages, faults, and limitations on a pendulum model. This study confirmed our previous observations and recommendations from literature that accelerometers can be reliably used as an inclinometer only for static measurements or very slow ambulation. Angle estimation in the still periods showed the closest match to reference angles by using accelerometers as inclinometers, while for the oscillating sequence, it gave the worst results (if we neglect double integration of accelerometers leading to enormous drift). This founding implies that, although the inclinometer can be successfully implemented during anchoring intervals for zero- velocity updates and calculations for the foot segment, it can lead to errors if it is applied for shank and thigh segments. For patients who ambulate very slowly, this may be acceptable, but those gait velocities should be less than 0.2 m/s. Smart filtering method based on accelerometer pairs, as expected from the results in chapter 3, can successfully provide angle estimation for movements with repetitive nature, such as gait. However, for patients with severe gait deformities whose gait cannot be considered cyclic, this method distorts the signal and provides false angle values. Finally, the results of the optimization method, as expected from the results in chapter 4, showed that this method can provide reliable angle estimation for all gait types and velocities. However, the disadvantage of the optimization method is the requirement to know the distances from the sensors to the centers of joints. These values do not need to be very precise, but clinicians are required to take the measuring tape and measure the positions, which slightly prolongs the assessment process. Integrating signals from gyroscopes provides the best results for angle estimations, comparing all investigated techniques during oscillating period, if there can be anchoring intervals which could be used for integration reset. As explained in chapter 5, this can easily be applied for the foot segment, but not for shank and thigh. Gyroscope drift can be minimized by applying polynomial fitting, as explained in chapter 5, but in some cases there is a risk of distorting the signal. Furthermore, integration of angular velocities obtained from different sensors placed on same 105 segments provided different results. This is an important issue which proved that each sensor comprises different noise which results in different drift after integration. Regarding calculus problems, while transforming one variable to another by numerical integration or differentiation, the numerical differentiation turned out to be very noisy due to quantization errors, and thus almost useless. The integration, naturally, suffers from drift. Hence, the drift removal techniques are essential to obtain reliable results for angles, in particular if the observation interval is longer than several seconds. 106 7. Clinical Applications∗ Clinical gait assessment: no-tech, low-tech and high- tech methods Clinical gait assessment is being performed in many ways, from patient’s self- report measures and visual observation by therapist (no-tech assessment) to performance-based set of multiple tasks (either no-tech by filling up questionnaires based on visual observation, or low-tech by measuring distance and/or timing). High- tech assessments involve kinematic and kinetic analysis for assessment of balance and gait disorder. This kind on gait assessment can be performed by systems that are integrated in the clinical gait laboratory or by wearable systems which allow gait assessment in any environment and unlimited space. In self-report measures, respondent rank the presence or absence of a problem with walking or a walking-related task, with rankings ranging from no difficulty in ∗ This chapter is based on following publications: Popović M.B, Djurić-Jovičić M., Petrović I., Radovanović S., Kostić V., “A simple method to assess freezing of gait in Parkinson's disease patients”, Braz J Med Biol Res, September 2010, Volume 43(9), pp. 883-889. Djurić-Jovičić M., Jovičić N.S., Milovanovic I., Radovanović S., Kresojevic N., Popović M.B, “Classification of Walking Patterns in Parkinson’s Disease Patients Based on Inertial Sensor Data”, Proceedings from the 10th Symposium on Neural network Applications in Electrical Engineering, Neurel 2010, 23-25 September, pp. 3-6, Belgrade, Serbia. Kojović J., Djurić-Jovičić M., Došen S., Popović M.B., Popović D.B., “Sensor-driven four-channel stimulation of paretic leg: Functional electrical walking therapy”, Journal of Neuroscience Methods, June 2009, Volume 181(1), pp. 100-105. 107 task performance (independent), to unable to perform the task wither with or without human assistance (dependent) (Alexander et al., 2000). The performance-based measures include sets of functional gait and balance tasks (which includes gait-related tasks such as turning while standing) in order to detect and quantify abnormalities. These tasks are either timed or scored semi- quantitatively, usually based upon whether a subject is able to perform the task, and if able, how normal or abnormal the performance was. Typical gait parameter that is being used for gait assessment is gait speed, measured as part of a timed short distance, measuring the usual and maximal gait speed, within 6-minute walk test, or long-distance corridor walk (Alexander et al., 1996). Another assessment approach is set of multiple tasks. Dynamic Gait Index includes evaluating alterations in response to changing tasks demands, including gait speed, head turns, turning, clearing an obstacle, stair climbing, walking backwards, with eyes closed or on a narrowed support (Shunway-Cook et al., 2004). The Emory Functional Ambulation Profile (EFAP) measures the time to walk under five environmental circumstances, with and without the use of an assistive device in stroke patients: 1) 5-m walk on hard floor, 2) 5-m walk on short pile carpeted floor, 3) time up and go test, 4) step over a brick and turn around a trash can, 5) walk up four steps, turn around and return (Wolf et al., 1999). The Berg Balance scale evaluates static and dynamic balance abilities though 14 simple balance tasks ranging from standing up from sitting down position, to standing on one foot. The degree of success in achieving each task is given a score of zero (unable) to four (independent), and the final measure is the sum of all of the scores (Berg et al., 1989). The Brunnstrom scale separates the recovery process into six stages, from stage one (no movements of the affected limb), to stage six (healthy like movement) (Brunnstrom, 1970). The Fugl- Meyer (FM) sensorimotor test comprises groups of items measuring movement, coordination, reflexes, sensation and other motor behaviors (Fugl-Meyer et al., 1975). The FM test for lower extremities allocates 34 points for maximum performance. Functional Ambulation Classification (FAC) uses a five-point scale to rate the extent of human assistance (stand-by, intermittent touch, and continuous support) required to walk varying surfaces (level, nonlevel, stairs, and inclines) while using an assistive device if necessary (Holden et al., 1984). Performance-Oriented Mobility Assessment, 108 also known as the Tinetti Balance and Gait Scale, in one of the earliest and most widely used methods to assess balance, gait, and fall risk in older adults. This test includes an evaluation of balance under perturbed conditions (such as while rising from a chair, after a nudge, with eyes closed, and while turning) as well as an evaluation of gait characteristics (such as gait initiation, step height, step length, continuity and symmetry, trunk sway, and path deviation) (Tinetti et al., 1988). Timed Up and Go Test (TUG) is a measure of time taken to stand u from a chair with armrests, walk 3 m, turn, walk back to the chair and sit down. Difficulty and/or unsteadiness in TUG performance is recognized as an important part of fall risk assessment (American Geriatrics Society, 2001). Dual task performance has been linked to an increased risk of falls based on walking performance while performing a simultaneous cognitive (dual) task. However, changes due to therapy or rehabilitation are often slow and susceptible to subjective evaluation of the clinician. Besides that, there are also issues with subjectivity which is based on the observer’s experience and individual bias, the restriction to a specific pathology, and low sensitivity to changes (Konig et al., 1997, Flansbjer et al., 2005). Hence, in order to create real clinical image of the patient, assessment of the therapy effects require implementation of hi-tech measurements which offer objective evaluation of the results and which enhance the clinician’s ability to assess the overall outcome. Among measuring techniques used for rehabilitation outcome evaluation, gait analysis is well standardized and currently most used (Langhorne et al., 2011). A clinical gait laboratory usually record patients walking with motion capture system (based on stereophotogrammetry) which tracks patient movements (usually used in parallel with standard video system), force plates that measure ground reaction force (described in chapter 1). Electromyography systems are also often used to record muscle activity while walking. Some other gait analysis systems that could be found in clinical settings are treadmills with integrated force plates, or force-sensitive walking paths (such as GaitRite, described in chapter 1). However, all these measuring tools are accessible only in a few specialized laboratories, as they are complex, expensive, require a lot of room space, and they can hinder patient’s usual gait pattern. 109 Another approach is to use wearable sensor system and record patient’s gait in any environment (in clinical settings, outdoors, at home). This provides more precise assessment of patient’s gait patterns because the space is not limited and there is enough time for patient to “forget” he is being recoded. Flexible shoe insoles with integrated force sensors are frequently used to evaluate load distribution and temporal gait characteristics (such as Stride analyzer, Pedar, TecScan). They provide lower resolution and precision than force plates but allow unlimited number of steps to be made, which is their main advantage. Besides force distribution, wearable systems can comprise inertial sensors which can record gait kinematics (accelerations and velocities) and after signal processing, position and orientation estimations resulting in joint angle estimations, the most significant parameter for many clinical assessments (Perry et al., 1992). Although these systems have been developed over the last decade and they are increasing their accuracy and decreasing their prices, they are still not standardized as clinical assessment tool. IMU-based systems have been proven through numerous research publications (Aminian et al., 2004; Ferrari et al., 2010; Schwesig et al., 2011, van den Noort et al., 2012) to be valuable clinical assessment tool, yet, they are still considered to be more adequate for researchers that for clinical staff (physicians and therapist). In this thesis, we proposed an ambulatory system based on kinematic sensors attached on the lower limbs in order to overcome the limitations of the previously mentioned techniques. The main features of the measurement device were to be portable and wireless, easy to use, accurate, robust, and capable of continuously recording data in long-term without hindrance to subject’s gait pattern. Moreover, new algorithms were proposed to accurately measure joints and segments angles in the sagittal plane. These data were then used to develop a gait analysis system providing spatio-temporal parameters, kinematic curves, and skeleton visualization. In this chapter, we applied the gait analysis methods, proposed in chapter 5, in real clinical applications. System outputs were combined to provide clinicians with an objective image about patient’s state and recovery process. 110 Gait analysis for patients with Parkinson’s disease Postural and gait disorders are the most disabling cardinal signs found in people with Parkinson's disease (PD). These patients have been reported to have postural impairments including a reduction of limits of stability, impaired postural adjustment, and poor responses to perturbation. Gait disturbance includes slow speed, shorter step length, and increased variability of step time. Some patients also suffer from episodic features such as freezing of gait (FOG). Gait and balance deficits predispose people with PD to falls. In a 20-year follow-up study, it has been reported that 87% of individuals with PD experienced one or more falls and 35% sustained injuries resulting from falls during walking. Falls can lead to physical injuries and psychological traumas such as fear of falling. This results in functional mobility restriction, loss of independence, social isolation, decrease in quality of life with increased risk of institutionalization, and increased mortality rate. Therefore, the objective gait assessment and monitoring progress of the disease can give clinicians important information about changes in the gait pattern and potential gait deviations, especially for patients who exhibit FOG episodes and concomitant falls. In this study, we used the proposed system (SENSY) for assessment of gait pattern for patients with PD. We used SENSY gait parameters analysis to provide follow-up of the treatments effects or progress of the disease. In people with Parkinson’s disease, the relationship between stride length and stride frequency exhibits abnormalities (Morris et al., 1998). In some cases, medications can stabilize a person’s gait by increasing the stride length and/or decreasing stride-to-stride variability. The benefit of this treatment can be quantified by monitoring the stride characteristic before and after treatment. Also, this kind of analysis can provide a warning about changed gait pattern and possibility of subsequent falls. As an upgrade of this application, we also developed a method for identification of freezing episodes for patients with PD, which is a special focus of the study. Freezing of gait phenomenon About one third of Parkinson’s disease patients experience sudden, transient block of movement performance, the phenomenon known as freezing of gait (FOG) (Giladi et al., 1992). A FOG episode is defined as the state when the patient is not responding within 1 s to the instruction to walk, or if it appears as he/she is trying 111 unsuccessfully to initiate or continue locomotion, as well as break in gait performance for no apparent reason. These episodes typically last from few seconds to 30 s., but in more cases they may last even for several minutes, and they are often resistant to therapy (Giladi et al., 2001). FOG can be manifested in several ways. Patients can come to a complete motor blocks (akinesia), which patients often describe as a feeling of the feet being glued to the floor. This “glued feet” feeling can also be accompanied by the trembling in place, i.e., the alternated tremor of the legs (knees) usually at a frequency 3-8 Hz (Moore et al., 2008; Hausdorf et al., 2003). Another FOG feature is walking with very short steps, which is a consequence of increase of cadence with a decrease in step length. Besides short steps, patients can also experience shuffling forward, defined as very short shuffling steps, several mm to several cm long, or walking in place where the foot or toe does not leave the ground or only barely clears the support surface (Hausdorf et al., 2003). Investigators commonly distinguish five typical provocative factors for FOG scenarios: start hesitation, turn hesitation, hesitation in tight quarters (FOG through narrow space), destination hesitation, and open-space hesitation (Schafsma et al., 2003). Turns and narrow spaces are especially important for provoking FOG episodes. During turns, gait control mechanisms might produce different motor programs to the axial (inner) and pivotal (outer) leg, and this mismatch challenges coordination (Plotnik et al., 2012). Step length reduction of the inner leg challenges step scaling control and may trigger the sequence effect. On the other hand, walking through narrow spaces may lead to slowness of gait and step length reduction while choosing a leading leg for passing through the passage (attention and coordination demanding) (Plotnik et al., 2012). Circumstances that can cause FOG episodes include approaching doorways, dual-tasking, distractions, crowded places, and being under mental pressure or stress. Furthermore, FOG can be provoked by emotions, excitement, auditory cueing at the proper pace, targets for stepping, and climbing stairs. Patients with FOG exhibit increased variability of stride time and length, disordered bilateral coordination, reduction of stride amplitude, and increased cadence to abnormally high rates during a turn compared to PD without FOG and to healthy subjects (Hausdorf et al., 2003; Plotnik et al., 2008; Chee et al., 2009). 112 FOG abnormalities are more pronounced in patients in the off state of the disease (i.e., without dopamine therapy). Episodic abnormalities of gait accompanying FOGs are: incremental decrease in stride length, highly reduced joint angles ranges (hip, knee, and ankle), disordered temporal control of gait cycles, and high-frequency alternate trembling-like leg movements with power peaks in frequency bands of 3-8 Hz*. Since FOG can be asymmetric and affect one side or be provoked more by turning on one side, movements of both legs should be assessed in order to obtain an objective image about patient’s gait pattern and appearance and frequency of the FOG episodes (Nutt et al, 2011). Besides being very uncomfortable, and hindering patient’s daily activities, FOG represents a common cause of falls and consequent injuries in PD, which stresses the importance of clinical assessment of FOG (Moore et al., 2008). Nowadays, the assessment has been mostly based on subjective patient reports and questionnaires. Furthermore, FOG episodes are difficult to be elicited in a routine clinical examination, requiring performance of complex walking patterns with turns and obstacles. Still, there is no objective method to identify the FOG phenomenon or the type, duration, and intensity of disorder episodes. 113 Gait assessment: Recognizing gait disturbances and freezing episodes Inertial sensors (accelerometers and gyroscopes) are often used for gait pattern classification and activity recognition (Najafi et al., 2003; Han et al., 2003; Liu et al., 2007). In order to classify the gait or detect freezing, most algorithms use spectral analysis of signals from triaxial accelerometers placed on various parts of the lower body (Nyan et al., 2006; Wand et al., 2007; Zabaleta et al., 2009). Some investigations suggest the analysis of electromyographic profiles as a method to predict FOG (Nieuwboer et al., 2004). However, these studies do not offer proper FOG type classification. In this section, we describe three methods for recognition of FOG episodes. The fist one is based on FSR sensors. It is simple, yet efficient algorithm for recognition of gait disturbances provoked by FOG episodes. The second and third methods are based on inertial sensors. The second method is the most sophisticated. Besides recognition of gait disturbances, it also provides their classification (by discriminating regular walking, small steps, shuffling, intentional standing, festinations, or akinesia). The third method provides spectrum analysis of the recorded gait sequence, presented in a form of spectrogram. It is not designed as an expert system, but clinician can interpret change of leg segments frequencies and interpret the data according to them. 1. Gait disturbances and FOG recognition based on FSR sensors We have developed a simple method for the detection of gait disturbances provoked by FOG episodes (Popovic et al., 2009b; Djuric-Jovicic et al., 2009b; Djuric-Jovicic et al., 2010a). This method is based on Pearson’s correlation between a selected regular (“normal”) stride and the entire gait sequence. Our choice to apply the correlation analysis as a method to detect gait disturbances was based on the theory considered by Rodgers and Nicewander (1988), and Aldrich et al. (1995). The Pearson’s correlation coefficient (PCC) is a measure of the linear dependence between two signals. If the two signals have the same shape, the coefficient is 1. For the gait, given its quasiperiodic nature, the correlation between a “normal” stride and the rest of the strides should decrease during gait disturbances typical for patients with PD, such as smaller steps, shuffling, festinations, and complete akinesia. 114 Consequently, our hypothesis was that FOG may be detected as time intervals when the PCC diminishes from values that are typical for normal locomotion. Experiment setup and recording protocol By using force sensing resistors placed under heel, metatarsals, and toes area of both feet, we recorded force profiles from nine PD patients (7 males and 2 females, age: 70.4 ± 7.9). All patients reported a clinical history of FOG episodes. The experiments were performed at the Institute of Neurology, Clinical Centre of Serbia. The study was performed in accordance with the ethical standards of the Declaration of Helsinki. Institutional Ethics Committee approval was obtained and participants gave informed written consent prior to the inclusion in the experiment. Patients were asked to stand up from the chair placed in the corridor, walk toward the room, pass a doorway, turn 180° to the left (U-turn), and walk the same route back, ending with a turn to sit back in the chair, as shown in Figure 7.1. The complex path included gait initiation, doorway passes, a U-turn, and a destination. The distance to walk was approximately 13 meters in each direction. Patients were asked to complete two trials separated by a rest period of at least 10 minutes. All experiments were recorded with video camera. Clinicians used videos to identify gait disturbances and FOG episodes (type and duration), and these data were further used for validation of our method. Figure 7.1. - Walking protocol. Patients were asked to stand up from the chair placed in the corridor, walk toward the room, pass a doorway, turn 180° to the left (U-turn), and walk the same route back, ending with a turn to sit back in the chair. 115 Method We estimated ground reaction force profiles based on signals recorded form force sensitive resistors (FSR) placed under heel, metatarsal, and toe areas. Ground reaction force (GRF) was normalized with respect to its maximal value within each gait record. For each recorded gait sequence, one stride from the sequence of “normal” locomotion was selected, and we used this stride to compute the Pearson’s correlation coefficients with the entire gait record. Strides from the GRF records yielding Pearson’s peaks that substantially deviated from ±1 were considered to be “irregular”. We hypothesized that the decrease of these peaks was the result of gait disturbances, such as FOG. Video recordings were used for the validation of our method. The computed PCC was plotted in parallel to the video signal and GRF record. Visually identified FOGs from the video were compared and confirmed with those determined via PCC (GRF). We compared the duration of FOG intervals determined by the two methods, and the relative difference in duration was estimated as videopccvideo ttt )( −=Δ , and reported as a percentage. An example of freezing of gait (FOG) detection using Pearson’s correlation coefficient (PCC) is shown in Fig. 7.2. The top bar shows locomotion examined from video with recognized FOG episode (yellow dotted bar). The middle panel shows the ground reaction force calculated from FSRs. One “normal” stride from the sequence is highlighted by the red line. The result of computed PCC between the highlighted stride and the entire recorded sequence is shown in the bottom panel. A linear envelope (dashed line) was created from PCC peaks with values under ±1. 116 Figure 7.2 –An example of FOG detection using PCC analysis. Middle panel shows the ground reaction force of the entire sequence and one “normal” stride (red line). Bottom panel: result of computed PCC between the highlighted stride and the entire recorded sequence. Results We analyzed data from 24 episodes of FOG in 9 patients with advanced PD. An example of the PCC-based evaluation of freezing episodes in relation to several different obstacles is presented in Figure 7.3. This procedure involved the selection of a “normal/regular” step (thick red line) from a plot of ground reaction forces (upper panel) followed by the computation of the PCC (lower panel). During normal walking, PCC values oscillated between ±1, and just before the obstacles ±1 peaks were observed. This result occurred six times during this record, as shown in Figure 7.3. These times were compared to the times of FOG episodes estimated from the video (Figure 7.3, yellow bars at bottom). The relative difference for 6 FOG episodes was, at most, 0.6%. Figure 7.3 also illustrates various irregular freezing patterns. Akinetic freezing was observed before the door; “short strides” were observed in four instances (both times before the line, approaching the door after U- turn. and approaching the chair); and asynchronous complex movement, followed by “short strides”, was observed during turning. 117 Figure 7.3. Illustration of the PCC method in different walking situations. Upper panel: ground reaction force with one “normal/regular” stride selected (thick red line). Lower panel: results of computed PCC between selected “normal” stride and the entire gait sequence. Bottom bar: FOG episodes identified from video recordings. 2. FOG recognition and classification based on inertial sensors The SENSY algorithm for trajectory reconstruction, described in chapter 5, was upgraded to recognize and classify gait deformities and freezing episodes: festinations, shuffling, short steps, and akinesia (Djuric-Jovicic et al., 2010b). The upgrade includes applying rule-based thresholding classification of the stride length, stride time, speed, and frequency of the stride. Finding the optimal way to describe the stride spectrum by one numerical value was a challenging task. Considering the nature of FOG episodes, which could be manifested on one leg segment and absent on the other, we decided to follow all three segments (thigh, shank, and foot). As signals of interest, we selected the axes describing vertical accelerations. After calculating the power spectrum of the accelerometer signals for each stride, we needed to describe it by one numerical value. Chung et al. (2008) performed activity classification (standing, walking, running) based on the vertical acceleration component by analyzing the area under the signal magnitude and median frequency of the spectrum. The median frequency is defined as the frequency which divides the area of the power spectrum into two equal parts. The area under the signal magnitude was used to discriminate standing from moving, and this parameter combined with median frequency was used for discrimination of walking and running. We compared the performance of the median frequency with the mean frequency (frequency at which the average power within 118 epoch is reached), the peak frequency (frequency at which the maximum power is reached), and the mean power frequency (the average power of power spectrum within epoch), as shown in Fig. 7.4. The median frequency seems to provide the best results for this application. Therefore, we applied a similar procedure as in Chung et al. (2008) to differentiate the akinesia from FOG with leg trembling. 0 0.2 0.4 0.6 0.8 1.0 1.2 -10 0 10 20 t[s] ac c [m /s 2 ] 0 5 10 15 20 0 0.5 1 f[Hz] |F FT (a cc )|2 N power spectrum median f peak f MPF mean f Figure 7.4. Choosing characteristic frequency of the spectrum. Upper panel: shank accelerometer signal (vertical axis) in time domain, stride #37. Lower panel: its power spectrum (normalized) and characteristic frequencies – median frequency, peak frequency, mean power frequency, and mean frequency The classification algorithm starts by performing the stride segmentation explained in chapter 5. Each stride is then analyzed by using the stride length (sl), stride time (st), stride speed (ss), and median frequency of vertical accelerations from all three leg segments (mft, mfs, and mff). Gait pattern types are defined as following: normal stride, short stride, shuffling stride, FOG with involuntary leg movements (leg trembling), and FOG with complete akinesia. Each stride is classified as one of these five types, by applying rule-based thresholding classification. Each stride is presumed as “normal”, and then it is tested for four types of gait disturbances. For all four pathological states (gait disturbances) we have heuristically selected thresholds for the following parameters: stride length, stride duration, stride velocity, and median frequency of the stride spectrum. Normal (typical) values of 119 these parameters are shown in Table 7.1, as well as applied thresholds for rule-based classification of FOG episodes. Thresholds for stride length are selected according to criteria from literature (Nutt et al., 2011, Plotnik et al., 2012), describing values that discriminate shuffling, small strides, and normal strides. Thresholds for other three parameters were selected heuristically, based on inspecting recorded data from patients with FOG. Table 7.1. Typical values, lower and higher thresholds for gait parameters for rule- based classification Gait parameter Normal (typical) THlow THhigh Stride length (sl) [m] 0.7-1.3 0.15 (0.2) 0.5 Stride time (st) [s] 0.8-1.2 0.5 2 Stride speed (ss) [m/s] 0.5-1.2 0.5 2 Median frequency (mf) [Hz] /* 5 10 *depends on a leg segment and the observed axis. The decisions about FOG types are made according to rules form Table 7.2. Table 7.2. Rule-based classification of gait disturbances for PD patients. Gait disturbance Rule (Condition) Small stride sl(stride)> THSL,low and sl(stride) THST,high and ss(stride)< THSS,low and mft(stride)< THMF,high or mfs(stride) THST,high and ss(stride)< THSS,low and (mft(stride) < THMF,high or mfs(stride) < THMF,high or mff(stride)> THMF,high) Since the stride segmentation is performed from one foot contact with the ground to another (next) contact of the same foot with the ground, this means that both FOG types (FOG with involuntary movements of the leg segments and FOG with complete akinesia) will be captured within a complete stride that will be followed by eventual lifting of the foot and swing phase until the next foot contact. In this way, the stride length of this stride may belong to any category, depending on the way a patient exits the FOG episode. For this reason, the FOG detection and classification is based on the analyzing median frequency of the leg segments during the stance phase, which provides more reliable results than analyzing the entire stride. Higher frequencies indicate festinations, while lower frequencies indicate motor blocks, i.e., akinesia. 120 Experiment setup and protocol Ten patients with PD and clinical history of FOG were recorded for this study (7 males, 3 females, age: 64.3 ± 7.4), performed at the Neurology clinic, Clinical Centre of Serbia. All patients reported a clinical history of FOG. The study was performed in accordance with the ethical standards of the Declaration of Helsinki and all participants gave informed written consent prior to participate in the study. The sensor system was mounted on the patient (one IMU per each leg segment of both legs and FSR clusters integrated in shoe insoles). Patients were asked to walk along the created complex pathway, as illustrated in Fig. 7.5. Each patient was asked to start from the red chair where he/she was sitting, to stand up upon a voice command and start walking straight towards door 1 (normal door), to pass the door, turn left towards door 2 (very narrow), to make U-turn, return through door 2, then to go straight along the corridor, where he would pass door 3 (wide door), several strides later to make another U-turn, pass through door 3 again, turn left, pass through door 1, return to the chair, and sit down. This recording protocol (walking sequences) was repeated four times per subject. All walking sequences were recorded with video camera. Figure 7.5. Floor plan and the pathway along which patients are walking. 121 Results An example of a recorded sequence from a patient with FOG episodes is presented in Fig. 7.6, showing ground reaction forces of both legs obtained from FSR signals, which were summed and normalized to their maximal value within the entire sequence. Figure 7.6 also shows stride segmentation by red square markers. The horizontal bars on the bottom show FOG episodes recognized and labeled from the video recording by an experinced clinician, while the arrows mark moments when the patient was reaching each expected obstacle in the pathway. 0 0.2 0.4 0.6 0.8 1 G R F le ft 0 10 20 30 40 50 60 70 80 0 0.2 0.4 0.6 0.8 1 G R F ri gh t t [s] Figure 7.6 – Stride segmentation shown on FSR signals, upper panel – left leg, lower panel – right leg. Each red square marker represents one recognized stride. Horizontal bars on the bottom show FOG episodes recognized by clinician from video recording. Based on performed segmentation, and calculated gait parameters, we can observe stride-to-stride variability and apply previously defined thresholds and classification rules (Fig. 7.7). 122 Fi gu re 7 .7 . V ar ia bi lit y of st rid e le ng th , s tri de ti m e, st rid e sp ee d, a nd m ed ia n fre qu en ci es a lo ng th e re co rd ed se qu en ce , le ft le g an d ri gh t l eg . H or iz on ta l a xe s s ho w th e st rid e nu m be r w ith in th e re co rd ed se qu en ce . U pp er a nd lo w er th re sh ol d fo r c la ss ifi ca tio n ru le s a re sh ow n w ith b la ck d as he d lin es . 5 10 15 20 25 30 35 40 45 50 55 60 0 0. 51 1. 5 Stride L. [m] 5 10 15 20 25 30 35 40 45 50 55 60 024 Stride T. [s] 5 10 15 20 25 30 35 40 45 50 55 60 0 0. 51 1. 5 Speed [m/s] 5 10 15 20 25 30 35 40 45 50 55 60 051015 f median [Hz] st ri de # th ig h sh an k fo ot 5 10 15 20 25 30 35 40 45 50 55 60 0 0. 51 1. 5 Stride L. [m] 5 10 15 20 25 30 35 40 45 50 55 60 024 Stride T. [s] 5 10 15 20 25 30 35 40 45 50 55 60 0 0. 51 1. 5 Speed [m/s] 5 10 15 20 25 30 35 40 45 50 55 60 051015 f median [Hz] st ri de # th ig h sh an k fo ot 123 The results of the performed classification are shown in Fig. 7.8, by applying the following color code for each stride: black marker – normal stride, blue marker – small stride, green marker – shuffling stride, red marker – FOG with trembling, pink marker – akinesia. 0 10 20 30 40 50 60 70 80 0 0.5 1 G R F L 0 10 20 30 40 50 60 70 80 0 0.5 1 G R F R t [s] Figure 7.8. Ground reaction forces for both legs and color-coded stride segmentation which shows the following classification: normal strides (black markers), small strides (blue markers), shuffling strides (green markers), FOG with trembling (red markers) and FOG with motor blocks (pink markers). Clinical tool The described stride classification is merged with the algorithm for the stride reconstruction (explained in chapter 5), which resulted in an illustrative clinical tool that allows inspection of the subject’s movement and trajectory reconstruction with highlighted gait disturbances. This tool plots the subject’s 3-D trajectory as a stick diagram (Fig. 7.9), one “leg” per each stride, and the gait disturbances are color coded in the same way as in Fig. 7.8. 124 As explained, for this application, i.e., for detecting gait disturbances, we had to record the walking along complicated pathways and to analyze patients’ movements in 3-D. For this reason, we had to cancel elimination of the drift from the algorithm for 3-D trajectory reconstruction in order to keep all the details of these movements and festinations, but also keep all movements provoked by this complex pathway. In this way, we traded the precision for trajectory assessment, but this can provide clinicians with valuable information about the situation where freezing episode occurred. 125 -5 0 5 10 -10 -5 0 5 10 0 5 x [m]y [m] z [m ] (a) -5 0 5 10 15 -10 -5 0 5 10 0 5 x [m]y [m] z [m ] (b) Figure 7.9. Illustration of patient’s 3-D trajectory with recognized and classified FOG episodes, shown by color coding. Comparison of (a) Left leg and (b) right leg show asymmetry of the present gait disturbances. The algorithm marks freezing episodes by different colors and indicates gait sequences that should be studied and analyzed more carefully. The changes in the gait pattern were present at the beginning of the sequence (start hesitation), which is shown at (0,0,0) coordinate. The patient stopped after the 126 first step. Further on, he experienced gait disturbances before reaching door 1 (shuffling steps). After passing though the door 1, patient had another FOG episode where he made sudden stop (motor block). After start hesitation, we proceeded towards door 2, where he had another FOG episode (motor block). Passing through doorway was followed by very short strides and then he made a U-turn. He proceeded through the hallway and through the door 3, and made another U-turn where he stopped. After that, he started walking towards door 1, which manifested as regular gait until re reached the door where he exhibited FOG with leg trembling and nearly lost his balance at the door. After this episode, proceeded forward and returned to the chair. These FOG episodes are also identified by our method, as it can be seen from Fig. 7.9., proving this method to be valuable FOG assessment tool. 3. Gait assessment by spectrogram As explained in introduction of this chapter, PD patients may exhibit changes in frequencies of the movements. Normal walking is typically characterized with frequencies from 0.5 Hz to 3 Hz (vertical shank acceleration), FOG with alternate leg trembling typically manifests with tremor in the range from 3 Hz to 8 Hz, while patient can also experience motor blocks with no movement at all. For this reason, it is especially convenient to show spectrum analysis as a function of time, where clinician can observe patient’s changes of stride frequency, correlate them to existing obstacles along the path, and measure their duration and intensity. Very intuitive visual tool which provides this information is spectrogram, calculated from the time signal using the short-time Fourier transform preformed for each stride independently. For PD gait assessment application, spectrogram is used as a graph with the horizontal axis showing time or number of strides, and the vertical axis frequency. There is also a third dimension indicating the amplitude of a particular frequency at a particular time, and it is represented by the intensity or color of each point in the image. We selected Jet colormap in Matlab, where the low amplitudes are represented with cold color tones (starting from navy blue for the lowest amplitude) and heading towards warmer colors with amplitude increase (finishing with dark red color for the highest amplitude). 127 Figure 7.10. Example of FOG interpretation from spectrogram: FOG episode began as a motor block (stride#2) and was followed by festinations episode with frequencies between 5 and 7 Hz (stride#3) 128 Fi gu re 7 .1 1. S pe ct ro gr am s o f p ow er sp ec tru m s f ro m v er tic al a cc el er at io n ax es fo r a ll th re e le g se gm en ts . ( a) le ft le g, (b ) r ig ht le g. F or e ac h sp ec tro gr am , t he h or iz on ta l ax is sh ow s t he st rid e nu m be r, w hi le th e ve rti ca l a xi s s ho w s s ho rt- te rm p ow er sp ec tru m o f t he si gn al fo r e ac h st rid e. 129 Discussion and conclusions This study demonstrated three different tools for recognition of gait disturbances in PD. The first one, based on signals from FSR, performed computing Pearson’s correlation coefficient between one “regular” stride and the entire record of ground reaction forces from the patient’s foot. It was sensitive to various freezing patterns, such as sudden akinetic posture with or without leg trembling, as well as to small strides or shuffling. The method can not distinguish between these two motions. However, the advantage of this method is that it only requires FSR sensors which can be integrated in shoe insoles, and this would be sufficient to give a clinician an objective evaluation of the frequency and duration of FOG episodes. These FSR sensors with accompanying acquisition electronics can be hidden under clothes. Therefore, the patient could wear it as a holter monitor during his/her daily routines and it would not be visible to others or hinder his/her daily routines. The results of this method are published in (Popovic et al., 2010). The second method, based on inertial sensors, performed recognition and classification of gait disturbances by applying rule-based classification based on stride length, stride time, and limb frequencies. This method provides full information about FOG episodes and changes in gait patterns which preceded these episodes. However, the decisions about the classification of gait disturbances are sensitive to threshold values and prone to errors. Namely, some episode might be missed due to having strict thresholds, the observed parameter can have the value very close to the assigned threshold, but this would not satisfy the mathematical condition for classification to corresponding gait disturbance. Therefore, this final visualization tool (Fig. 7.9) should be considered in combination with stride-to-stride variability curves (Fig. 7.7) with shown thresholds, and according to them clinician should change the threshold of the critical gait parameter. Although this could be easily applied in practice, significant improvement of the method would include automatic threshold adaptation. Nevertheless, this study is still going on and these results should be confirmed when we collect significant number of patients for statistics. The third method, based on inertial sensors, performed gait pattern analysis by spectrogram. This illustrative tool is easy to interpret and provides clinician necessary 130 information about the temporal changes of the frequencies of the lower limbs. This method discriminates walking, FOG with leg trembling, and FOG with motor block. This provides clinician enough information about walking unsteadiness. Results of this method were shown at Symposium of Clinical Neurology (Djuric-Jovicic et al., 2010). When it comes to patients with PD, gait analysis performed in clinical environment frequently doesn’t actually capture a patient’s real state i.e., gait pattern. Since the patients are aware their gait is being observed, they often (consciously or subconsciously) change their pattern, by trying to walk better, or to emphasize their movement disorders. Therefore, having gait assessment system which could be used in home environment, as a holter monitor, would provide clear and objective image about patient’s state, as well as frequency and duration of experienced gait disturbances. This could be arranged by simple hardware adaptation, allowing sensor units to store data internally, instead of sending it to remote PC. SENSY software can be used as a software platform for this application, with a performed upgrade related to classification of FOG episodes. 131 Gait analysis for patients with hemiplegia The hemiplegia, which commonly follows neurological conditions such as a stroke, traumatic brain injury, rupture of an aneurysm or a vascular malformation, reduces patients' ability to walk. The primary goals of patients with hemiplegia include being able to walk independently and to manage to perform their daily activities from (Dittuno et al., 2005). The restoration of walking continues to be a major goal of rehabilitation for patients with hemiparesis. The success of rehabilitation process depends on many factors, including the severity of motor impairment and gait dysfunction, as well as intensity, specificity, and duration of the provided interventions. SENSY system was used for assessment of the recovery of gait by providing objective evaluation and quantification of gait parameters and gait kinematics (Djuric- Jovicic et al., 2010c). By using the proposed system, it is possible to have comparison of patient’s state before and after the therapy, and assessment of walking with different assistive devices (e.g., cane and Walkaround). SENSY was also used in synergy with EMG systems, and with electrical stimulators for functional electrical therapy (for sensor driven stimulation and assessment of FES effects). In the same way, SENSY was used for patients with various injuries that resulted in impaired walking abilities. The goal of this study was to provide gait parameters as a new objective method to assess rehabilitation effects. Recording protocol The proposed gait analysis system was used for gait assessment for patients with hemiplegia who were hospitalized at the Rehabilitation clinic “Dr. Miroslav Zotovic”, Belgrade, Serbia. The sensor setup comprised six IMUs (one per each leg segment of both legs) and FSRs integrated in shoe insoles, placed under heel, metatarsal and toe areas (Fig. 7.12). 132 Figure 7.12. SENSY clinical setup, assessment of the hemiplegic patients. We recorded signals from 22 patients with hemiplegia (age: 55±11 years, Fugl- Meyer score (FM) for lower extremities: 22±5 out of 34). The inclusion criteria for recruiting patients were the ability to stand on the paretic leg with some support, to be cognitively ready to understand the instructions, and without other neurological deficits. All subjects signed the informed consent approved by the local ethics committee. Subjects were asked to walk with their natural pace and to use the assistance they regularly use (therapist, cane, Walkaround). We recorded four walking sequences per assistive device for each subject. Walking was recorded in a 10 m long and a 3 m wide hallway, without any visual obstacle, providing clear space for continuous walking. All experiments were recorded with a video camera in parallel. Quantifying rehabilitation progress/success Assessment protocol for this part of study included gait recording at the beginning of rehabilitation process, follow-up three weeks after the beginning of therapy, and then six months after therapy (Fig. 7.13). 133 0 0.5 1 G R F -20 0 20 H ip 0 20 40 60 K ne e 1 2 3 4 5 6 7 -20 0 20 A nk le t[s] before therapy 3 weeks later 6 months later Figure 7.13. Therapy follow-up, comparison of ground reaction forces and joint angles before therapy, three weeks later and six months since the beginning of therapy. 0 0.5 1 Left -20 0 20 0 20 40 60 0.5 1 1.5 2 2.5 3 3.5 -20 0 20 t[s] 0 0.5 1 Right (Paretic) -20 0 20 0 20 40 60 0.5 1 1.5 2 2.5 3 3.5 -20 0 20 t[s] before after Figure 7.14. Therapy follow-up, comparison of ground reaction forces and joint angles before and after therapy (6 months later). Horizontal axes show that time interval of one stride before therapy (grey areas) corresponds to two and a half strides after therapy (red lines). 134 Table 7.3. Therapy follow-up, comparison of gait patterns before, after three weeks of therapy and after the therapy (6 months later). Before therapy (average±std.dev) Follow-up (average±std.dev) After therapy (average±std.dev) Gait parameters Left/Total Right* Left/Total Right* Left/Total Right* Cadence [stride/min] 27,10 / 58.44 / 76,19 / Stride Length [m] 0,90±0,09 / 0.82±0.11 / 0,99±0,05 / Step Time Dif [s] -0,10±0,23 / 0.09±0.04 / -0,06±0,04 / Cycle Time Dif [s] -0,08±0,19 / -0.01±0.06 / 0,01±0,03 / Speed [m/s] 0,26±0,01 / 0.45±0.06 / 0,71±0,03 / Step time [s] 1.78±0.20 1.67±0.05 0.88±0.04 0.97±0.03 0.73±0.02 0,67±0.02 Stride time [s] 3.45±0.18 3.37±0.13 1.85±0.06 1.85±0.05 1.40±0.03 1.4±0.03 Swing cycle [%] 18.57±2.23 18.31±1.25 32.25±1.46 31.31±2.04 33.72±1.53 37.94±1.32 Stance cycle [%] 81.43±2.23 81.69±1.25 67.75±1.46 68.69±2.04 66.28±1.53 62.06±1.32 Single supp. cyc.[%] 17.65±0.9 19±2.17 31.42±1.94 32.49±1.97 38.25±1.56 33.91±1.17 Doub. supp. cyc.[%] 63.78±2.94 62.69±3.12 36.33±2.62 36.2±2.73 28.03±1 28.15±2.05 * paretic leg. Table 7.4. Therapy follow-up, comparison of outputs for three clinical scales, before therapy, after three weeks of therapy and after the therapy (6 months later). Before therapy Follow-up After therapy Clinical scale BI FM BB BI FM BB BI FM BB Score 94 30 46 100 32 52 100 33 55 * Bartel index (BI), Fugl-Meyer score (FM), Berg Balance scale (BB) The results indicated that the gait functions of the patient were improved considerably after rehabilitation therapy. The numerical values of the gait parameters are shown in Table 7.3. As it can be seen from the Table, the patient increased speed by increasing stride length and cadence, and shortening stride time. Furthermore, the rehabilitation increased symmetry between the healthy and the paretic side, which can be seen from the step and stride time differences. Significant improvement can be noticed in the first follow-up, three weeks after the beginning, where the swing and stance ratio improved and became closer to the ratio typical for healthy gait pattern (Langhorne et al., 2011). Improvement of the gait functions can also be notices from the decrease of variability (expressed as standard deviations) of the step and stride time differentials, stride time, and swing/stance phases. However, perhaps the most significant parameters which testify the improvement of the stability and gait functions are the single-support and double-support cycle time, expressed as percentages of the gait cycle. The double support cycle decreased after three weeks rehabilitation and then 135 slowly continues the trend, which is confirmed in the second follow-up, six months after the beginning of therapy. Significant improvements are also present in the range of joint angle flexions/extensions, as shown in Fig. 7.14. As the result, the functional improvement was associated with an increase in walking performance (higher speed, higher range of rotations) and an improvement of walking regularity (lower gait variability). We provided an objective gait assessment, using ambulatory gait analysis, for assessing functional recovery for patients after stroke. These results cannot be obtained through other clinical evaluations, and complement the clinical scores by a useful objective evaluation. 136 Overall investigated clinical applicability of the system SENSY was used and tested in many experiments including healthy individuals and patients with gait disorders. Until now, we have recorded more than 150 patients in total. The description of the experiments, patients’ diagnoses, and the institution where the experiments were held are shown in Table 7.5. Table 7.5 – List of clinical experiments performed with SENSY system Institution Patients (Diagnosis) Experiment Patients # Neurology clinic, KCS, Belgrade, Serbia PD Recording gait patterns, 2xACC SENSY 11 Children’s Health Care Institute, Novi Sad, Serbia Cerebral palsy Recoding gait patterns, 2xACC SENSY 8 Rehabilitation clinic “Dr. Miroslav Zotovic”, Belgrade, Serbia Stroke Recoding gait patterns, 2xACC SENSY 7 Rehabilitation clinic “Dr. Miroslav Zotovic”, Belgrade, Serbia Stroke Comparison of gait patterns when walking with different assistive devices, 2xACC SENSY 10 Rehabilitation clinic “Dr. Miroslav Zotovic”, Belgrade, Serbia Stroke Recoding gait patterns – collecting sensor data for development of sensory driven electrical stimulation for control of movements 2xACC SENSY 12 Rehabilitation clinic “Dr. Miroslav Zotovic”, Belgrade, Serbia Stroke Implementation of sensory-driven electrical stimulation and follow-up of the effects of FES 2xACC SENSY 9 Clinics for physical medicine, rehabilitation and protetics, Clinical Center, Nis, Serbia. Stroke Implementation of sensory-driven electrical stimulation and follow-up of the effects of FES 2xACC SENSY 11 Neurology clinic, KCS, Belgrade, Serbia PD Healthy SENSY/ Gait Rite Comparison 2xACC SENSY 7 pat. + 5 h. Rehabilitation clinic “Dr. Miroslav Zotovic”, Belgrade, Serbia Stroke SENSY/ Goniometers Validation of angle estimation algorithms for pathological gait 2xACC SENSY 13 Rehabilitation clinic “Dr. Miroslav Zotovic”, Belgrade, Serbia Stroke Comparison of walking with and without Walkaround (before and after therapy) 2xACC SENSY 11 Neurology clinic, KCS, Belgrade, Serbia PD Gait assessment ACC+GYRO SENSY 22 Rehabilitation clinic “Dr. Miroslav Zotovic”, Belgrade, Serbia Stroke Kinematics and polymyography study (assistance with cane and Walkaround) ACC+GYRO SENSY 21 Rehabilitation clinic “Dr. Miroslav Zotovic”, Belgrade, Serbia Stroke SENSY/ Goniometers Validation of angle estimation algorithms for pathological gait ACC+GYRO SENSY 13 Neurology clinic, KCS, Belgrade, Serbia PD Monitoring therapy effects during the day ACC+GYRO SENSY 7 Neurology clinic, KCS, Belgrade, Serbia PD Recording FOG ACC+GYRO SENSY 8 137 The table content describes various clinical applications for which the sensor system (hardware) was used. Table 7.5 presents our clinical experience. Software was being developed in parallel with these experiments, and some of the described experiments were used for validation of the algorithms, as shown in Table 7.5. As it can be seen from the Table 7.5, SENSY was included in many clinical experiments, where system performances were tested and validated. In parallel with these recordings, and in collaboration with clinicians and patients, SENSY was upgraded and improved to present state. Based on this experience, we designed hardware and software in order to help clinicians and increase the quality of gait assessment in clinical setting, as well as to decrease the requested time for performing analysis. Nonetheless, for clinical measurement studies of movement impaired patients, SENSY is a highly efficient, cost-effective assessment tool. 138 8. Graphic user interface for gait analysis software Based on the methods explained in chapter 5, we have developed gait analysis software that provides both quantitative and qualitative evaluation of gait pattern. Software outputs include kinematic and kinetics diagrams, animated graphic images simulating the patients’ gait at various conditions, spectrum analysis, and spatio- temporal parameters of gait. SENSY software for gait analysis (algorithms and graphic interface) is developed in Matlab (MathWorks, USA). The features of this software correspond to other advanced commercial gait analysis systems, and its accuracy and performances and comparable to other gait analysis systems available on the market. Figure 8.1 – Initialization of graphic user interface Upon software initialization (Fig. 8.1), user is offered to select and import a recorded file of interest and it plots recorded raw calibrated data from all sensors. User can then choose either to observe complete recorded file or to select the 139 sequence of interest. From this panel user can observe raw data in more details, and select or deselect channels from sensor nodes (Fig. 8.2). Also, user can change the sequence he wants to observe or zoom in/out the sensor data. After the desired sequence is selected, data processing algorithms are performed on those data, and user can continue to the next panels that provide the analysis of kinematics, kinetics, and calculations of spatio-temporal parameters. Figure 8.2 – Initial software panel showing original sensor data. Kinematics panel comprises analysis of joint angles, segment angles and reconstruction of leg trajectories, as shown in Fig. 8.3, 8.4 and 8.5. The kinematic diagrams provide information about segment and joint angles, as well as the reconstruction of movement (gait) trajectory. These graphs help clinicians qualitatively assess time evolution of lower limb movements, variability at different phases of gait, symmetry, and ranges of movements. 140 Figure 8.3 –Kinematics panel showing joint angles of the hip, knee, and ankle for both legs. Figure 8.4 –Kinematics panel showing joint angles of the hip, knee, and ankle for both legs, strides plotted over each other showing stride-to-stride variability in time and range of motion. 141 Figure 8.5 –Kinematics panel 2-D reconstruction of leg trajectories, for left and right leg. Besides angle estimation, software also provides trajectory estimation through a stick diagram which mimics lower limb skeleton performing the same actions as the recorded subject. This tool provides visually appealing and easy to interpret information about subject’s gait pattern. The animated skeleton (stick diagram) can also be seen in 3-D and it can be viewed from arbitrary angles, thereby further helping the clinician to interpret the results. Kinetics panel shows force profiles, estimation of ground reaction forces (horizontal and vertical), as well as spatial visualization of force profiles, as shown in Figure 8.5. By clicking the Top View button (or 3-D model button) user can see how the forces change in time. Velocity of this presentation can be increased or decreased by sliding the vertical bar next to the image. 142 Figure 8.6 – Kinetics panel showing force profiles, ground reaction forces and graphical visualization of force profiles. In spectrum analysis panel, user can select any accelerometer or gyroscope signal (any axis from any segment) and watch its spectrum (FFT) and spectrogram. The vertical axis of the spectrogram shows frequency and the horizontal shows number of samples (100 samples=1 s) Figure 8.7 – Spectrum analysis, selected signal in time and frequency domain (left side of panel) and its spectrogram (right side of panel) 143 From the statistics panel, user can observe spatio-temporal gait parameters, their average value within selected sequence and standard deviation, as presented in Fig. 8.8. The spatio-temporal parameters provide a tool for objective outcome measures to quantify the expected gait improvement during rehabilitation of patients, monitoring therapy effects, or progress of a disease. Temporal gait parameters are by default calculated based on the recognized gait phases from FSRs. However, they can also be calculated from inertial sensors if there was a problem with FSRs or FSRs were not used for recording. Software also offers the possibility to store patient’s data related to the file, including name of the patient, date of birth, height/weight, diagnosis, walking assistance, etc. After examination, the processed data can be exported to Excel file which can be stored in database or attached to patient’s medical chart. Figure 8.8 – Statistics panel showing average values and standard deviations of spatio-temporal gait parameters (left side of the panel). Right side of panel provides export of the processed data to Excel file. 144 9. Conclusion and future research Contribution of dissertation The primary objective of this thesis was to design and validate new methods for gait analysis based on signals acquired from inertial sensors which can be used for clinical applications (starting hypothesis H1 and H2). Through investigated clinical applications the emphasis was put on kinematic parameters affected by locomotion disorders in order to provide an objective outcome evaluation based on these parameters. The main results and contributions of this thesis can be summarized as follows: Ambulatory recording system: A specific wireless sensor system based on inertial sensors and force sensing resistors was designed. The system is a portable ambulatory device with the following design criteria: wearable sensors are lightweight, donning and doffing is easy and fast, mounting does not require precise positioning, and sensors do not hinder the gait pattern. New methods for estimation of segment and joint angles based on accelerometer pairs were introduced, different in their complexity, applicability and required sensor configuration. The first two methods are based on sensor units with accelerometer pairs and their applicability is limited to 2-D (angles in saggital plane), while the third method uses combination of accelerometers and gyroscopes in order to estimate 3-D kinematics of a subject. In the first method, segment and joint angles were estimated by applying band- pass filtering with cut-off frequency related to subject’s gait pace, i.e., without need for integration. In this way, the obtained angles are free from drift. The method was validated on 12 healthy subjects walking at various speeds from 0.2 m/s to 2.0 m/s, 145 and 15 subjects with hemiplegia with different levels of gait impairment who were walking at their own comfortable pace. The obtained angles were compared to goniometers outputs which were used as a reference system. Comparison of the outputs of these two systems showed that accelerometer pairs can be used in this way with satisfying reliability and precision. However, this method is based on the rhytmicity of the gait, it does not proved good results for non-cyclic ambulatory movements which are present in more severely impaired gait patterns. The second method uses simplex optimization to model the offset by a slowly varying function of time (a cubic spline polynomial) and evaluate the polynomial coefficients by nonlinear numerical simplex optimization with the goal to reduce the drift in processed signals (angles and movement displacements). Besides segment and joint angles, this method also provides reconstruction of subject’s trajectory. These two methods represent positive answer to thesis hypothesis H1. The third method uses transformation matrices for estimation of kinematic parameters (angles and trajectories) based on accelerometer and gyroscope sensor combination and polynomial fitting to eliminate the drift of the estimated outputs. This method is the answer to starting hypothesis H2. Based on signals acquired from this sensor system and the proposed methods, a specific gait analysis software was designed. The goal of this software is to allow quantitative gait analysis for clinicians that would provide them objective evaluation of therapy effect or progress of a disease. Clinical applications: based on developed algorithms described in chapter 5, two independent clinical applications were designed: module for gait assessment of the stroke patients with hemiplegia, and module for gait assessment of the patients with Parkinson’s disease. These modules allow clinicians to have objective evaluation of a patient’s gait pattern, rehabilitation and therapy effects. In this way, clinicians can react in time, and change the therapy if needed, and prevent possible near-fall situation. These are the confirmations for starting hypothesis H1 and H2. For patients with PD, objective measure of gait disturbances directly provides a very important tool not only for following progress of the disease but also respond to treatments (medicines) and prescribed therapies. Also, since the proposed system could be used as a holter monitor, required time for patient to stay in hospital is significantly shortened. Patient can be examined during the day and returned home 146 with SENSY standalone holter unit. In this way, clinicians can get real picture about patient’s gait disturbances (types and timings) and save money for hospital days. These results are affirmative answer to starting hypothesis H3. For stroke patients, the proposed system can be used for assessment of walking with various assistive devices and orthoses. In this way, the proposed system can help clinicians to objectively select optimal walking assistance and provides a proof or a measurement of the assistance’s effect. Also, the proposed system can be used for assessment of the effects of functional electrical stimulation (FES), which is also time consuming part of the rehabilitation process. Providing gait assessment before, during and after walking with FES, short-term and long-term following of the therapy’s effects helps better planning of the therapy (better planning of stimulation protocol, intensity of stimulation, correction of stimulated muscles and also duration of the therapy). This application proved our starting hypothesis H4 to be true. SENSY should also be used for elderly population, in their annual (routine) examination. Following some specific gait characteristics (stride time, stride length, double support etc.), SENSY provides more objective view about aging effect and possible risks of falling. Responding to this kind of change in time and preventing falls (which in elderly population usually end with fractures) and prescribing adequate assisting device or similar, objective gait assessment saves patients and money that would otherwise be spent on hospital costs for healing fractures. In total, SENSY offers quantitative assessment of gait pattern and rehabilitation success. In comparison to traditional methods of gait analysis (observational and with scales, e.g. Fugl-Meyer), it provides more objective monitoring of patient’s state and better quality of service. Suitability of the proposed system Clinicians in the National Health Service (NHS) in the UK identified the need for a widely accessible and valid objective motion analysis tool. A 2003 survey of 1.826 physiotherapists in the UK revealed that only 23.1% of people with gait impairments were measured in gait laboratories. The main reasons cited were that it’s too time consuming (41.8%), too expensive (38.8%), and there is a lack of technical knowledge (27%). These results were confirmed in a 2009 survey conducted by Oxford Brookes University. 147 There are various systems for gait analysis present on the market today. They range in accuracy, price, outputs characteristics, applied technology, and suitability of clinical applications. The physiotherapist or the hospital clinician would greatly benefit from a faster, easy-to-use solution that is also small, lightweight, and preferably wireless. Inertial measurement units (IMUs) provide the basis for such a clinical tool. Our proposed system fulfills major clinical demands. The donning takes in average around 3 minutes, and there is a short self-calibration (1s) which does not require additional tools or special action. The doffing of the system is around 1min. The positioning of the sensors is well described, and the recordings do not depend on the precise positioning or orientation with respect to the body segment. The recordings can start immediately after donning. Wireless design eliminates all cables. The system also comprises instrumented shoe insoles. The mass of units that are at the body segment is about 50 g. The units are battery powered (rechargeable) and provide up to 3 hours of continuous recording. The distance from the base commute can be up to 30 m, and the data transfer is secured and basically noise free. The software supports online visualization of data recorded and storing of the data. Based on our clinical experience and great number of recorded patients, we recognized a necessity for objective quantification of gait impairment and recovery (or progress of the disease). For most of the diagnoses, there are functional scoring scales that usually detect if some kind of impairment is or isn’t present but there is a large gap between those levels. For this reason, it is impossible for clinicians to precisely quantify patient's recovery based on these clinical scales, because recovery is often a slow process and improvements can be very small and need to be observed and recorded in order to optimize therapy. Perspective and future research The thesis can be extended to the following directions: Extending 2-D to 3-D The proposed algorithms for kinematics analysis based on accelerometer pairs, described in chapters 3 and 4, offer high reliability and accuracy of the estimated segment and joint angles. However, these algorithms are limited to 2-D analysis of 148 movements in the sagittal plane. By adding more accelerometers in spatial arrangement, these methods could be extended to 3-D analysis, which could be beneficiary for some gait deformities which require information about other angles besides flexion/extension. Movement analysis The sensor system could be enhanced by adding 3-D magnetometers to existing sensor units. In this way, with the missing azimuth orientation, we could have complete 3-D positioning and orientation estimation that would allow more precise estimation during turns and better trajectory reconstruction of the complete lower limb trajectory (both legs together). Also, signals from magnetometers can be used for drift correction, but with careful restrictions regarding the surrounding area i.e., presence of ferromagnetic materials that create signal disturbances in the magnetometers. The proposed sensor system is, as explained, designed as an open platform with possibility to include different type of sensors. Besides magnetometers, system could also support ultra-sound sensors that could provide clinicians one more important gait parameter which is not possible to be assessed by inertial sensors – stride width. This parameter can provide insight about patient stability, especially significant in early stroke rehabilitation stages. Gait analysis The gait analysis program, explained in chapter 9, can be used for assessment of other gait pathologies. Future extension of the algorithm should focus to automatically detect different types of walking and activities of a patient at various conditions and focus on long-term outcomes. In general, the proposed clinical protocol can be used for other assessment and rehabilitation programs related to lower extremity treatments. Movement control Some parts of developed methods which are applicable in real time can be used for gait event detection and movement control by creating rules for functional electrical stimulation. This could have many applications, in stroke patient rehabilitation, stimulation of the PD patient when he exhibits FOG episode. It could 149 also be extended to other gait pathologies, such as cerebral palsy etc. Furthermore, movement control based on these sensors could be used for sports and fitness to improve exercise performance and help burning calories, which is an application that is becoming more popular each day. Holter monitor When it comes to patients with PD. gait analysis performed in clinical environment frequently doesn’t actually capture a patient’s real state i.e.,gait pattern. Since the patients are aware their gait is being observed, they often (consciously or subconsciously) change their pattern, by trying to walk better, or to emphasize their movement disorders. Therefore, having gait assessment system which could be used in home environment, as a holter monitor, would provide clear and objective image about patient’s state, as well as frequency and duration of experienced gait disturbances. This could be arranged by simple hardware adaptation, allowing sensor units to store data internally, instead of sending it to remote PC. SENSY software can be used as a software platform for this application, with a performed upgrade related to classification of FOG episodes. Hardware adaptations: miniaturization and minimization Although the proposed system has quite satisfying hardware regarding their size, weight, and comfort for their users, there is ever-lasting struggle to decrease their size up to invisibility. It is in human nature to want to fit in, i.e., not to look sick or unable, or helpless. Therefore, patients often prefer to have assessment hidden as much as possible, which directs engineers towards miniaturization and applicability of other already existing solutions, such as smart phones with integrated inertial sensors and strong platform that could support many computational tasks. Besides miniaturization, another goal is minimization of the assessment system. Attaching one sensor unit per each leg segment certainly provides detailed and reliable information about gait pattern and possible gait deformities. Although the proposed system does not hinder patient’s movements, and it is simple and not time consuming for donning and doffing, it would be better to use less sensors and to get reliable description of the gait pattern. Reduction in number of sensors would allow more elegant usage of the sensor system, and decrease number of things clinicians have to take care of while mounting and during assessment protocol. 150 Further upgrade of the proposed system combines sensors with actuator system for functional electrical stimulation. In this way, sensor data can be used for motor control and, in the same time it can perform the assessment of the movement made with it. This system is already developed and its architecture and preliminary testing are described in journal publication (Jovicic et al., 2012). 151 References Aldrich J., 1995. Correlations genuine and spurious in Pearson and Yule, Stat. Sci, 10, 364-376 Allet, L., Knols, R., Shirato, K., Bruin, E., 2010. Wearable Systems for Monitoring Mobility-Related Activities in Chronic Disease: A Systematic Review. Sensors 10(10), 9026-9052. Allexander N.B, Guire K.E., Thelen D.G., 2000. Self-reported walking ability predicts functional physical performance in frail older adults. Journal of American Geriatrics Society, 48, 1408-1413. American Geriatrics Society, 2001. Guideline for prevention of falls in older persons. Journal of American Geriatrics Society, 49, 664-672. Aminian, K., Najafi, B., 2004. Capturing human motion using body-fixed sensors: outdoor measurement and clinical applications. Computer animation and virtual worlds, 15, 79-94. Aminian, K., Najafi, B., Bnla, C., Leyvraz, P.-F., Robert, P., 2002. Spatio-temporal parameters of gait measured by an ambulatory system using miniature gyroscopes. Journal of Biomechanics 35(5), 689-699. Aminian, K., Trevisan, C., Najafi, B., Dejnabadi, H., Frigo, C. et al., 2004. Evaluation of an ambulatory system for gait analysis in hip osteoarthritis and after total hip replacement. Gait and Posture, 20, 102-107. Arden, N., Nevitt, M., 2006. Osteoarthritis: epidemiology. Best Pract Res Clin Rheumatol, 20(1), 3-25. Bachmann, E.R.., Duman, I., Usta, U.Y., McGhee, R. B., Yun, X.P. Zyda M.J., 1999. Orientation tracking for humans and robots using inertial sensors. In Proceedings of IEEE International Symposium on Computational Intelligence in Robotics and Automation CIRA '99, 187-194. Berg, K., Wood-Dauphinee, S., Williams, J.I., Gayton, D., 1989. Measuring balance in the elderly: preliminary development of an instrument. Physiotherapy Canada, 41,304-11. Bobath, B., 1990. Adult hemiplegia: evaluation and treatment. Butterworth-Heinemann Medical. Bonomi, A., Salati, S., 2010. Assessment of Human Ambulatory Speed by Measuring Near-Body Air Flow. Sensors, 10(9), 8705-8718. Brizard, A.J. Motion in a non-inertial frame, http://www.dartmouth.edu/~phys44/lectures/Chap_6.pdf, Acceleration in noninertial frames.pdf, 2004. Brunnstrom S. 1970. Movement theory in Hemiplegia: A Neurophysiological approach. New York: Harper and Row. 152 Burce, E.N., 2001. Biomedical signal processing and signal modeling. New York: Hogn Wiley and sons. Chung, W.Y., Purwar, A., Sharma, A., 2008. Frequency domain approach for activity classification using accelerometer. Proceedings of the EMBS 2008, 1120-1123. Cooper, G., Sheret, I., McMillian, L., Siliverdis, K., Sha, N., Hodgins, D. et al., 2010. Inertial sensor- based knee flexion/extension angle estimation. Journal of Biomechanics 42, 2678-2685. Culhane, K.M., O’Connor, M., Lyons, D., Lyons, G.M., 2005. Accelerometers in rehabilitation medicine for older adults.Age and Ageing, 34, 556-560. Dejnabadi, H., Jolles, B.M., Aminian, K., 2005. A new approach to accurate measurement of uniaxial joint angles based on a combination of accelerometers and gyroscopes. IEEE Trans. Biomed. Eng. 52, 1478-1484. Dejnabadi, H., Jolles, B.M., Casanova, E., Fua, P., Aminian, K., 2006. Estimation and visualization of sagittal kinematics of lower limbs orientation using body-fixed sensors. IEEE Transactions on Biomedical Engineering 53(7), 1385-1393. Den Otter, A.R., Geurts, A.C.H., Mulder, T., Duysens, J., 2006. Gait recovery is not associated with changes in the temporal patterning of muscle activity during treadmill walking in patients with post-stroke hemiparesis. Clinical Neurophysiology 117, 4-15. Ditunno, P.L., Patrick, M., Stineman, M., Morganti, B., Townson, A.F., Ditunno, J.F., 2005. Cross- cultural differences in preference for recovery of mobility among spinal cord injury rehabilitation professionals. Spinal cord, 44,567-575. Djurić-Jovičić, M., Milovanović, I., Jovičić, M., Radovanović, S., 2009a. Gait analysis: BUDA vs. GAITRITE. Proceedings of the 53rd ETRAN Conference, June 15-18, 2009, Serbia, ME 2.2., http://etran.etf.rs/etran2009/sekcije.htm Djurić-Jovičić, M., Milovanović, I.P., Jovičić, N.S., Popović, D.B., 2009b. Walkaround assisted walking of stroke patients. Conference Medical Physics and Biomedical Engineering, Sept. 7-12, 2009 Munich, Germany, 299-301. Djuric-Jovičić M., Popović M., Petrović I., Radovanović S., Kostić V., 2010a. Freezing vs no freezing steps in Parkinsons disease patients. Abstracts of the XVIII Congress of the International Society of Electrophysiology and Kinesiology, ISEK 2010, 16-19 June 2010, Aalborg, Denmark [CD- ROM]. ed. / Deborah Falla; Dario Farina. Aalborg: Department of Health Science and Technology. Aalborg University. Djurić-Jovičić, M., Jovičić, N.S., Milovanović, I., Radovanović, S., Kresojevic, N., Popović, M.B., 2010b. Classification of Walking Patterns in Parkinson’s Disease Patients Based on Inertial Sensor Data, Proceedings from the 10th Symposium on Neural network Applications in Electrical Engineering, Neurel 2010, 23-25 September, pp 3-6, Belgrade, Serbia Djurić-Jovičić, M., Milovanović, I., Janković, M., Dragin, A., 2010c. Synergies of a gait - Clinical measurements of kinematics and polymiography. Abstract in Clinical Neurophysiology, 122(7), e6. Djuric-Jovicic, M., Jovicic, N., Popovic, D.B., 2011. Kinematics of Gait: New Method for Angle Estimation Based on Accelerometers. Sensors 11(11), 10571-10585. 153 Djurić-Jovičić M., Jovičić N., Popović D.B., Djordjević A.R., 2012. Nonlinear Optimization for Drift Removal in Estimation of Gait Kinematics Based on Accelerometers. Journal of Biomechanics 45(16), pp. 2849-2854. Dorsey, E.R., Constantinescu, R., Thompson, J.P., Biglan, K.M. et al., 2007. Projected number of people with Parkinson disease in the most populous nations, 2005 through 2030. Neurology, 384- 386. Eng, J.J., Winter, D.A., 1995. Kinetic analysis of the lower limbs during walking: What information can be gained from a three-dimensional model? Journal of Biomechanics, 28, 753-758. Esser, P., Improving the Identification and Treatment of Walking Impairments. Oxford Brookes University, Clinical study European Brain Council's Cost of Brain Disorders in Europe, study Favre, J, Jolles, B.M., Aissaoui, R., Aminian, K., 2008. Ambulatory measurement of 3D knee joint angle. Journal of Biomechanics, 41, 1029-1035. Favre, J., Jolles, B. M., Aissaoui, R., Aminian, K., 2008. Ambulatory measurement of 3D knee joint angle. Journal of Biomechanics, 41(5), 1029-1035. Ferrari, A., Cutti, A.G., Garofalo, P., Raggi, M., Heijboer, M. et al., 2010. First in vivo assessment of “Outwalk”: a novel protocol for clinical gait analysis based on inertial and magnetic sensors. Medical and Biological Engineering and Computing. 48, 1–15. Flansbjer, U.B., Holmbäck, A.M., Downham, D., Patten, C., Lexell, J., 2005. Reliability of gait performance tests in men and women with hemiparesis after stroke. Journal of Rehabilitation Medicine, 37, 75-82. Frigo, C., Ferrarin, M., Frasson, W., Pavan, E., Thorsen, R., 2000. EMG signals detection and processing for on-line control of functional electrical stimulation. Journal of Electromyography and Kinesiology 10, 351-360. Frigo, C., Bardare, M., Corona, F., Casnaghi, D., Cimaz, R., Naj Fovino, P.L., 1996. Gait alterations in patients with juvenile chronic arthirtis: a computerised analysis. J Orthop Rheumatol, 9, 82-90. Fugl-Meyer A.R., Jaasko L., Leyman I., Olsson S., Steglind S., 1975. The post-stroke hemiplegic patient. 1. a method for evaluation of physical performance. Scand J Rehabil Med, 7,13-31. Furnee, H., 1997. Real-time motion capture systems in Three-Dimensional Analysis of Human Locomotion, Allard P., Cappozzo A., Lundberg A., Vaughan C.L. Eds. New York: Wiley, 85-108. Giladi, N., Treves, T.A., Simon, E.S., Shabtai, H., Orlov, Y., et al., 2001. Freezing of gait in patients with advanced Parkinson's disease”, Journal of Neural Transmission, 108(1), 53-61. Goldstein, H., Poole, C., Safko J., 2000. Classical Mechanics, Addison Wesley, San Francisco. Grood, E.S., Suntay, W.J., 1983. A Joint Coordinate System for the Clinical Description of Three- Dimensional Motions: Application to the Knee. Transactions to the ASME, 105,136-144. Han, J.H., Lee, W.J., Ahn, T.B., Jeon, B.S., Park, K.S., 2003. Gait analysis for freezing detection in patients with movement disorder using three dimensional acceleration system. Proceedings of the EMBS 2003, 2, 1863-1865. Hausdorff, J.M., Balash, Y., Giladi, N., 2003. Time series analysis of leg movements during freezing of gait in Parkinson’s disease: rhyme or reason?. Physica A: Stat Mechanics and Appl., 321, 565-570. 154 Hodgins, D., 2008. The Importance of Measuring Human Gait. Medical Device Technology, 19(5), 42- 47. Holden, M.K., Gill, K.M., Magliozzi, M.R., Nathan, J., Piehl-Baker, L., 1984. Clinical gait assessment in the neurologically impaired. Reliability and meaningfulness. Phys Ther, 64, 35-40. http://www.biometricsltd.com/gonio.htm Jovanov, E., Wang, E., Verhagen, L., Fredrickson, M., Fratangelo, R., 2009. deFOG – a Real Time System for detection and unfreezing of Gait of Parkinson’s Patients. Proceedings of the EMBS 2009, USA, 5151-5154. Jovicic, N.S, Saranovac L.V, Popovic, D.B, 2012. Wireless Distributed Functional Electrical Stimulation System. Journal of NeuroEngineering and Rehabilitation. Jurman, D., Jankovec, M., Kamnik, R., Topič, M., 2007. Calibration and data fusion solution for the miniature attitude and heading reference system. Sensors and Actuators A: Physical 138(2), 411– 420. Kiss, R.M., Kocsis, L., Knoll, Z., 2004. Joint kinematics and spatial-temporal parameters of gait measured by an ultrasound-based system. Med. Eng. Phys., 26, 611-620. Kobayashi, K., Gransberg, L., Knutsson, E., Nolén, P., 1997. A new system for three-dimensional gait recording using electromagnetic tracking. Gait and Posture, 6, 63-75. Kojović, J., Djurić-Jovičić, M., Došen, S., Popović, M.B., Popović, D.B., 2009. Sensor-driven four- channel stimulation of paretic leg: Functional electrical walking therapy. Journal of Neuroscience Methods, 181(1), 100-105. Konig, A., Scheidler, M., Rader, C., Eulert, J., 1997. The need for a dual rating system in total knee arthroplasty. Clinical Orthopaedics and Related Research, 161-167. Lagarias, J.C., Reeds, J. A., Wright, M. H., Wright, P. E., 1998. Convergence Properties of the Nelder- Mead Simplex Method in Low Dimensions. SIAM Journal of Optimization 9(1), 112-147. Langhorne, P., Bernhardt, J., Kwakkel, G., 2011. Stroke rehabilitation. The Lancet, 377,1693-1702. Legnani, G., Zappa, B., Casolo, F., Adamini R., Magnani P.L., 2000. A model of an electro-goniometer and its calibration for biomechanical applications. Med. Eng. Phys., 22, 711-722. Liu, R., Zhou, J., Liu, M., Hou, X., 2007. A wearable acceleration sensor system for gait recognition.Proceedings of 2nd IEEE Conference on Industrial Electronics and Applications, ICIEA 2007, Harbin, 2654-2659. Liu, T., Inoue, Y., Shibata, K., 2010. A Wearable Ground Reaction Force Sensor System and Its Application to the Measurement of Extrinsic Gait Variability. Sensors, 10(11), 10240-10255. Luczak, S., Oleksiuk, W., Bodnicki, M., 2006. Sensing Tilt With MEMS Accelerometers. IEEE Sensors Journal, 6(6), 1669-1675. Luinge, H.J., Veltink, P.H., 2004, Inclination measurement of human movement using a 3-D accelerometer with autocalibration. IEEE Transactions on Neural Systems and Rehabilitation Engineering, 12(1), 112-121. Luinge, H.J., Veltink, P.H., 2005. Measuring orientation of human body segments using miniature gyroscopes and accelerometers. Medical and Biological Engineering and Computing, 43, 273–282. 155 Madgwick, S.O.H., Harrison, A.J.L., Vaidyanathan, R., 2011. Estimation of IMU and MARG orientation using a gradient descent algorithm. Proceedings on IEEE International Conference on Rehabilitation Robotics (ICORR), 1-7. Mahony, R., Hamel, T., Pimlin J.-M., 2008. Nonlinear complementary filters on the special orthogonal group. IEEE Transactions on Automatic Control, 53(5), 1203 -1218. Mak, M.K.Y., Hui-Chan, C.W.Y., 2004. Audiovisual cues can enhance sit-to-stand in patients with Parkinson's disease, Movement Disorders, 19(9), 1012-1029. Mayagoitia, R.E., Nene, A.V., Veltink, P.H., 2002. Accelerometer and rate gyroscope measurement of kinematics: an inexpensive alternative to optical motion analysis systems. Journal of Biomechanics, 35(4), 537-542. Miyazaki, S., 1997. Long-Term Unrestrained Measurement of Stride Length and Walking Velocity Utilizing a Piezoelectric Gyroscope. IEEE Trans. On Biomedical Eng, 44(8), 753-759. Mizuike, C., Ohgi, S., Morita, S., 2009. Analysis of stroke patient walking dynamics using a tri-axial accelerometer. Gait and Posture, 30, 60-64. Moore, S.T., MacDougall, H.G., Ondo, W.G., 2008. Ambulatory monitoring of freezing of gait in Parkinson's disease. Journal of Neuroscience Methods, 167(2), 340-348. Morin, E., Reid, S., Eklund, J.M., Lay, H., Lu, Y. et al., 2002. Comparison of Ground Reaction Forces Measured with a Force Plate, F-Scan and Multiple Individual Force Sensors, Queen's University, Kingston, Ontario, Canada, Unpublished. Morris, J.R.W., 1973. Accelerometry – A technique for the measurement of human body movements. Journal of Biomechanics, 6, 729-736. Morris, M., 1998. Abnormalities in the Stride Length-Cadence Relation in Parkinsonian Gait. Movement Disorders, 13(1), 61–69. Nait-Ali A., 2009. Advanced biosignal processing. Berlin, Heidelberg: Springer. Najafi, B., Aminian, K., Paraschiv-Ionescu, A., Loew, F., Bula, C.J., Robert, P., 2003. Ambulatory system for human motion analysis using a kinematic sensor: monitoring of daily physical activity in the elderly. Biomedical Engineering, IEEE Transactions on, 50(6), 711-723. Nelder, J., Mead, R., 1965. A simplex method for function minimization, Computer Journal, 7, 308- 313. Niewboer, A., Dom, R., Weerdt, W.D., Janssens, L., Stijn, V., 2004. Electromyographic profiles of gait prior to onset of freezing episodes in patients with Parkinson’s disease, Brain, 127, 1650-1660. Nieuwboer, A., Fabienne, C., Anne-Marie, W., Kaat, D., 2007. Does freezing in Parkinson’s disease change limb coordination? A kinematic analysis. J Neurol, 254, 1268–77. Nutt, J.G., Bloem, B.R., Giladi, N., Hallett, M., Horak, F.B., Nieuwboer A., 2011. Freezing of gait: moving forward on a mysterious clinical phenomenon. The Lancet Neurology, 10(8), 734–744. Nyan, M.N., Tay, F.E.H., Seah, K.H.W., Sitoh Y.Y., 2006. Classification of gait patterns in the time- frequency domain. Journal of Biomechanics, 39(14), 2647-2656. O'Donovan, K.J., Kamnik, R., O'Keeffe, D.T., Lyons, G.M., 2007. An inertial and magnetic sensor based technique for joint angle measurement, Journal of Biomechanics 40, 2604-2611. Olney, S.J., Richards, C., 1996. Hemiparetic gait following stroke. Part I: Characteristics. Gait and Posture, 4, 136-148. 156 Pallejà, T., Teixidó, M., Tresanchez, M., Palacín, J., 2009. Measuring Gait Using a Ground Laser Range Sensor. Sensors, 9(11), 9133-9146. Pereda, E, Quian Quiroga, R, Bhattacharya, J., 2005. Nonlinear multivariate analysis of neurophysiological signals. Prog Neurobiol., 77, 1-37 Perry, J,P., 1992. Gait Analysis: Normal and Pathological Function. Thorofare, NJ: Slack, 1992. Peurala, S.H., Titianova, E.B., Mateev, P., Pitkanen, K., Sivenius, J., Tarkka, I.M., 2005. Gait characteristics after gait-oriented rehabilitation in chronic stroke. Restor Neurol Neurosci, 23,57- 65. Plotnik, M., Giladi, N., Dagan, Y., Hausdorff, J.M.., 2012. Is freezing of gait in Parkinson's disease a result of multiple gait impairments? Implications for treatment. Parkinson’s Disease, 2012, Article ID: 459321. Plotnik, M., Hausdorff, J.M., 2008. The role of gait rhythmicity and bilateral coordination of stepping in the pathophysiology of freezing of gait in Parkinson's disease, Movement Disorders, 23(2), supplement 2, S444–S450. Popovic, D.B., Sinkjær, T., 2003. Control of movement for the physically disabled: control for rehabilitation technology. 2 ed. Aalborg: Center for Sensory-Motor Interaction (SMI), Department of Health Science and Technology, Aalborg University, 488 p. Popović, D.B., Sinkjær, T., Popović M.B., 2009a. Electrical stimulation as a means for achieving recovery of function in stroke patients, Journal of NeuroRehabilitation, 25: 45-58. Popović M., Djuric-Jovičić M., Petrović I., Radovanović S., Kostić V., 2009b. Freezing of gait (FOG) in Parkinson’s disease patients: time analysis. Book of Abstracts, 9th Congress of Clinical Neurophysiology with International Participation, October 15-17, 2009, Belgrade, Serbia. Abstract published in Clinical Neurophysiology, Vol. 121, No. 4, 2010, p. e16, No. 48. Popović, M.B., Djurić-Jovičić, M., Petrović, I., Radovanović, S., Kostić, V., 2010. A simple method to assess freezing of gait in Parkinson's disease patients. Braz J Med Biol Res, 43(9), 883-889. Popović, D.B., Popović, M.P., 2011. Advances in the use of electrical stimulation for the recovery of motor function. Prog Brain Res, 194, 215-225 Reinsch, S.,MacRae, P., Lachenbruch, P.A., Tobis, J.S., 1992. Attempts to Prevent Falls and Injury: A Prospective Community Study. The Gerontologist, 32(4), 450-456. Rodgers, J.L., Nicewander, W.A., 1988. Thirteen ways to look at the correlation coefficient. Am Stat, 42, 59-66. Roetenberg, D., Luinge, H. J., Baten, C.T.M., Veltink, P.H., 2005. Compensation of magnetic disturbances improves inertial and magnetic sensing of human body segment orientation, IEEE Transactions on Neural Systems and Rehabilitation Engineering, 13(3), 395-405. Romanò, C.L., Frigo, C., Randelli, G., Pedotti, A., 1996. Analysis of the gait of adults who had residua of congenital dysplasia of the hip. J Bone Joint Surg, 78A, 1468–1479. Sabatini, A.M., 2006. Quaternion-Based Extended Kalman Filter for Determining Orientation by Inertial and magnetic Sensing, IEEE Transactions on Biomedical Engineering, 53(7), 1346-56. Sabatini, A.M., Martelloni, C., Scapellato, S., Cavallo, F., 2005. Assessment of Walking Features From Foot Inertial Sensing. IEEE Trans. On Biomedical Eng, 52(3),486-494. 157 Saup, P., 1995. Beat-to-beat variation of heart rate reflect modulation of cardiac autonomic outflow. News Physiol Sci., 5, 32-37 Schaafsma, J.D., Balash, Y., Gurevich, T., Bartels, A.L., Hausdorff, J.M., Giladi N., 2003. Characterization of freezing of gait subtypes and the response of each to levodopa in Parkinson’s disease. European Journal of Neurology, 10, 391–398. Schwesig R., Leuchte S., Fischer D., Ullmann R., Kluttig A., 2011. Inertial sensor based reference gait data for healthy subjects. Gait and Posture, 33(4), 673-678. Sekine, M., Abe, Y., Sekimoto, M., Higashi, Y., Fujimoto, T., et al., 2000. Assessment of gait parameter in hemiplegic by accelerometry. Proceedings of the 22nd Annual EMBS International Conference, July 23-28, 2000, Chicago IL. Shiratsu, A., Coury, H.J.C.G., 2003. Reliability and accuracy of different sensors of a flexible electrogoniometers. Clinical Biomechanics, 18, 682-684. Shumway-Cook A., Gruber W., Baldwin M., 1997. The effect of multidimensional excercises on balance, mobility and fall-risk in community-dwelling older adults. Physical Therapy, 77, 46-57. Sinkjær, T., Popovic, D.B., 2009. Neurorehabilitation technologies : present and future possibilities. NeuroRehabilitation, 25(1), p. 1-3. Stagni, R., Fantozzi, S., Cappello, A., Leardini A., 2005. Quantification of soft tissue artefact in motion analysis by combining 3D fluoroscopy and stereophotogrammetry: a study on two subjects. Clinical Biomechanics, 20, 320-329. Takeda, R., Tadano, S., Natorigawa, A., Todoh, M., Yoshinari, S., 2009. Gait posture estimation using wearable acceleration and gyro sensors. Journal of Biomechanics, 42, 2486-2494. Takeda, R., Tadano, S., Todoh, M., Morikawa, M., Nakayasu, M., Yoshinari, S., 2009. Gait analysis using gravitational acceleration measured by wearable sensors. Journal of Biomechanics, 42(3), 223-233. Tan, H., Wilson, A.M., Lowe, J., 2008. Measurement of stride parameters using a wearable GPS and inertial measurement unit. Journal of Biomechanics, 41, 1398-1406. Tinetti M.E., Williams T.F., Mayewski R., 1986. Fall risk index for elderly patients based on number of chronic disabilities. American Journal of Medicine, 80, 429-434. Tong, K., Granat, M.H., 1999. A practical gait analysis system using gyroscopes. Medical Engineering & Physics, 21, 87-94. Truelsen, T., Ekman, M., Boysen, G., 2005. Cost of stroke in Europe, European Journal of Neurology, Vol. 12, Issue Supplement s1, 78–84. van den Noort J.C., van der Esch M., Steultjens M.P., Dekker J., Schepers H.M., Veltink P.H., Harlaar J., 2012. The knee adduction moment measured with an instrumented force shoe in patients with knee osteoarthritis. Journal of Biomechanics, 45(2), 281-288. Vaughan, C.L, Davis, B.L, O'Coonor, J.C., 1992. Dynamics of human gait. Champaign, IL:Human Kinetics Publishers. Veg, A., Popovic, D.B., 2008. Walkaround: Mobile balance support for therapy of walking. IEEE Trans Neural Syst Rehabil Eng, 16, 264-269. 158 Wang, N., Ambikairajah, E., Lovell, N.H., Celler, B.G., 2007. Accelerometry based classification of walking patterns using time-frequency analysis. Proceedings of 29th Annual International Conference of IEEE-EMBS 2007. Lyon., 4899-4902. Willems, A.M., Nieuwboer, A., Chavret, F., Desloovere, K., Dom, R. et al., 2006. The use of rhythmic auditory cues to influence gait in patients with Parkinson’s disease, the differential effect for freezers and nonfreezers, an explorative study. Disability and Rehabilitation, 28(11), 721-728. Willemsen, A.T., van Alste, J.A., Boom, H.B.K., 1990. Real-time gait assessment utilizing a new way of accelerometry. Journal of Biomechanics, 23, 859-863. Willemsen, A.T.M., Frigo, C., Boom, H.B.K. 1991. Lower extremity angle measurement with accelerometers-error and sensitivity analysis. IEEE Transactions on Biomedical Engineering 38(12), 1186-1193. Williamson, R., Andrews, B.J., 2001. Detecting absolute human knee angle and angular velocity using accelerometers and rate gyroscopes. Medical and Biological Engineering and Computing, 39(3), 294-302. Winter, D.A.,1984. Kinematic kinetic patterns in human gait: variability and compensating effects. Human. Movement Science, 3, 51-76 Wolf S.L., Catlin P.A., Gage K., et al., 1999. Establishing the reliability and validity of measurements of walking time using the Emory Functional Ambulation Profile. Physical Therapy, 79, 1122- 1133. Woltring H.J., 1991. Representation and calculation of 3-D joint movement, Human Movement Science, 10(5), 603-616. Xu, W., Chang, C., Hung, Y.S., Kwan, S.K., Chin Wan Fung, P., 2007. Order statistics correlation coefficient as a novel association measurement with application to biosignal analysis. IEEE Trans. Signal Processiing, 55, 5552-5563. Yanagisawa, N., Hayashi, R., Mitoma, H.,2001 Pathophysiology of frozen gait in Parkinsonism. Adv Neurol, 87, 199–07. Yang, C.C., Hsu, Y.L., 2007. Algorithm design for real-time physical activity identification with accelerometry measurement. Proceedings of 33rd Annual Conference of the IEEE Industrial Electronics Society, IECON 2007, Taipei, 2996-3000. Yang, C.C., Hsu, Y.L. 2010. A Review of Accelerometry-Based Wearable Motion Detectors for Physical Activity Monitoring. Sensors, 10, 7772-7788. Yang, C., Hsu, Y., Shih, K., Lu, J. 2011. Real-Time Gait Cycle Parameter Recognition Using a Wearable Accelerometry System. Sensors, 11(8), 7314-7326. Yun, J., 2011. User Identification Using Gait Patterns on UbiFloorII. Sensors, 11(3), 2611-2639. Zabaleta, H., Keller, T., Marti Masso, J.F., 2009. Power spectral distribution for detection of freezing of gait in patients with Parkinson’s disease. IFMBE Proceedings, 22(16), 2089-2092. Zheng, H., Black, N.D., Harris, N.D., 2005. Position-sensing technologies for movement analysis in stroke rehabilitation. Med Biol. Eng. Comp., 43, 413-420. Zorowitz, R.D., 1999. Neurorehabilitation of the stroke survivor. Textbook of Neural Repair and Rehabilitation by Michael Selzer, 2, 579-592. 159 Appendix 160 Gait terminology Walking uses a repetitious sequence of limb motion to move the body forward while simultaneously maintaining stance stability. Each gait cycle is divided to two periods, stance and swing. Stance is the term used to designate the entire period during which the foot is on the ground. Stance begins with initial contact and ends with toes off which marks the beginning of the next phase- swing. The word swing applies to the time the foot is in the air for limb advancement. It begins as the foot is lifted from the floor (toe-off) and ends with initial contact. Stance is subdivided to three intervals according to the sequence of floor contact by two feet. Both the start and end of stance involve a period of bilateral foot contact with the floor (double support), while the middle portion of stance has one foot contact (single support). Initial double support begins the gait cycle. It is the time when both feet are on the floor after initial contact. Single support begins when the opposite foot is lifted for swing. Terminal double support is the third subdivision. It begins with floor contact by the other foot (contralateral initial contact) and continues until the original stance limb is lifted for swing (ipsilateral toe off). The gross normal distribution of the floor contacts is 60% for stance and 40% for swing. Timing for the phases of stance is 10% for each double support interval and 40% for single support. Single support of one limb equals swing of the other, as they occur in the same time. Gait phases In order to provide the basic functions required for walking, each stride involves an ever-changing alignment between the body and the supporting foot during stance and selective advancement of the limb segments in the swing. These reactions result in a series of motion patterns performed by the hip, knee and ankle. Each stride contains eight functional patterns called sub phases of the gait cycle. 161 Analysis of a person’s walking pattern by phases more directly identifies the functional significance of the different motions occurring at the individual joints. The phases of gait also provide a means for correlating the simultaneous actions of the individual joint into patterns of total limb function. This is particularly important for approach for interpreting the functional effects of disability. The relative significance of one joint’s motion compared to the other’s varies among the gait phases. Each of the eight gait sub phases has a functional objective and a critical pattern of selective synergistic motion to accomplish this goal. The sequential combination of the phases also enables the limb to accomplish three basic tasks: weight acceptance (WA), single limb support (SLS) and limb advancement (LA). Weigh acceptance begins the stance period and uses the first two sub phases, initial contact and loading response. Single limb support continues stance with the next two gait sub phases, mid stance and terminal stance. Limb advancement begins in the stance sub phase (pre-swing) and then continues through the three sub phases of swing, initial swing, mid swing, and terminal swing. Figure A.1 - Gait cycle with sub phases of gait Weight acceptance For this task, three functional patterns are needed: shock absorption, initial limb stability, and the preservation of progression. The challenge is the abrupt transfer of 162 the body weight onto a limb that has just finished singing forward and has an unstable alignment. 1. Initial contact Interval: 0-2% of GC This phase includes the moment when the foot just touches the ground. The joints postures present at this time determine the limb’s loading response pattern. Gait pattern (for normal gait): The hip is flexed, the knee is extended, the ankle is dorsiflexed to neutral. Floor contact is made with the heel. The other limb is at the end of terminal stance. Goal: The limb is positioned to start stance with a heel rocker. 2. Loading response Interval: 0-10% of GC This is the initial double support period. The phase begins with initial floor contact and continues until the other foot is lifted for swing. Gait pattern (for normal gait): Using the heel as a rocker, the knee is flexed to shock absorption. Ankle plantar flexion limits the heel rocker by forefoot contact with the floor. The opposite limb is in the pre-swing phase. Goals: Shock absorption, weight-bearing stability, preservation of progression Single limb support Lifting the other foot for swing begins the single limb support interval for the stance limb. This continues until the opposite foot again contacts the floor. During the resulting interval, one limb has the total responsibility for supporting body weight in both the sagittal and coronal planes while progression must be continued, Two phases are involved in single limb support: mid stance and terminal stance. 3. Mid stance Interval: 10-30% of GC This is the first half of the single limb support interval. It begins as the other foot is lifted and continues until body weight is aligned over the forefoot. Gait pattern (for normal gait): In the first half of single limb support, the limb advances over the stationary foot by ankle dorsiflexion (ankle rocker) while the knee and hip extend. The opposite limb is advancing in its mid swing phase. Goals: progression over the stationary foot and limb and trunk stability. 4. Terminal stance 163 Interval: 30-50% of GC This phase completes single limb support. It begins with heel rise and continues until the other foot strikes the ground. Throughout this sub phase body weight moves ahead of the forefoot. Gait pattern (for normal gait): During the second half of single limb support, the heel rises and the limb advances over the forefoot rocker. The knee increases its extension and then just begins to flex slightly. Increased hip extension puts the limb in a more trailing position. The other limb is in terminal swing. Goal: progression of the body beyond the supporting foot. Limb advancement In order to perform the limb advancement, preparatory posturing begins in stance. Then the limb swings through three postures as it lifts itself, advances and prepares for the next stance interval. Four sub phases of the gait cycle are involved in this: pre- swing (end of stance), initial swing, mid swing, and terminal swing. 5. Pre-swing Interval: 50-60% of GC This final phase of stance is the second (terminal) double support interval in the gait cycle. It begins with initial contact of the opposite limb and ends with ipsilateral toe-off. Other names for this sub phase are weight release and weight transfer. The limb is unloaded and the limb uses its freedom to prepare for the rapid demands of swing. Gait pattern (for normal gait): Floor contact by the other limb has started terminal double support. The reference limb responds with increased ankle plantar flexion, greater knee flexion and loss of hip extension. The opposite limb is in loading response. Goal: to position the limb for swing. 6. Initial swing Interval: 60-73% of GC This first sub phase of the swing is approximately one third of the swing period. It begins with lift of the foot from the floor and ends when the swinging foot is opposite the stance foot. 164 Gait pattern (for normal gait): The foot is lifted and limb advanced by hip flexion and increased knee flexion. The ankle only partially dorsiflexes. The other limb is in early mid stance. Goal: Foot clearance from the floor and advancement of the limb from its trailing position. 7. Mid swing Interval: 73-87% of GC This second phase of the swing period begins as the swinging limb is opposite the stance limb. The phase ends when the singing limb is forward and the tibia is vertical (i.e., hip and knee flexions postures are equal). Gait pattern (for normal gait): Advancement of the limb anterior to the body weight line is gained by further hip flexion. The knee is allowed to extend in response to gravity while the ankle continues dorsiflexing to neutral. The other limb is in late mid stance. Goals: Limb advancement and foot clearance from the floor. 8. Terminal swing Interval: 87-100% of GC This final phase of the swing period begins with a vertical tibia and ends when the foot strikes the floor. Limb advancement is completed as the leg (shank) moves ahead of the thigh. Gait pattern (for normal gait): Limb advancement is completed by knee extension. The hip maintains its earlier flexion, and ankle remains dorsiflexed to neutral. The other limb is in terminal stance. Goals: Complete limb advancement and to prepare the limb for stance. Spatio-temporal gait parameters The gait cycle is also identified by the term stride which describes actions of one limb. The duration of a stride is the interval between two sequential initial floor contacts by the same limb (i.e., right IC and the nest right IC). Step refers to the timing between the two limbs. There are two steps within each stride (gait cycle). At the midpoint of one stride the other foot contacts the ground to begin its next stance period. The interval between an initial contact by each foot is step (i.e., left and then right). 165 Stride length: the distance between the heel points of two consecutive strides (left to left, right to right). (The unit of measure is meters.) Step length: the distance from the heel center of the current footprint to the heel center of the previous footprint on the opposite foot. Stride time: the time elapsed between the first contacts of two consecutive footfalls of the same foot. Also called gait cycle time. Step time: the time elapsed from first contact of one foot to first contact of the opposite foot. Ambulation Time: the time elapsed between first contact of the first and the last footfalls. Cadence: number of strides per minute. Base width: the distance from heel center of one footprint to the line of progression formed by two strides of the opposite foot. This parameter is not possible to be estimated by inertial sensors. Toe in/Toe out: It is the angle between the line of progression and the midline of the foot. Angle is zero if the geometric midline of the foot is parallel to the line of progression; positive, toe out, when the midline of the footprint is outside the line of progression and negative, toe in, when inside the line of progression. (The unit of measure is degrees.) This parameter can not be estimated without magnetometers. Single support: the time elapsed between the last contact of the current footfall to the first contact of the next footfall of the same foot. It is measured in seconds and expressed as a percent of the Gait Cycle time of the same foot. Double Support: The two periods when both feet are on the floor, are called initial double support and terminal double support. Initial double support occurs from heel contact of one footfall to toe-off of the opposite footfall. Terminal double support occurs from opposite footfall heel strike to support footfall toe-off. Total double support is the sum of the initial double support and the terminal double support. Kinematic gait parameters Kinematic variables describe the extent, speed, and direction of movement of joints or body segments. They can be divided to translatory movements (displacement, velocity, and acceleration) and rotation movements (angular 166 displacement, angular velocity, and angular acceleration). Definition of these angles and various methods for their assessment are described in chapter 2. Gait event detections In order to estimate gait parameters, we need to have reliable and accurate estimation of the characteristic gait events such as heel contact and toe off which delimit swing and stance phases. There are numerous ways to do it, but one of the easiest is by applying thresholds to heel switch or FSR data (Djuric et al., 2008). The alternative is to use inertial sensors, usually from IMU placed on foot or shank segment, and to apply different processing methods for detection of these characteristic gait events (Sabatini et al., 2005; Kotiadis et al., 2010). Gait event detection from force sensors represents simpler and more reliable technique, especially for applications for pathological gait, where is difficult to adapt rules or thresholds for IMU-based gait event detection methods. For our applications, FSRs were used to detect gait phases by summing their signals from one leg, normalizing to sum’s maxima and applying threshold (THfsr) at 5% of the normalized sum. Intervals where the signal is beyond THfsr are considered to belong to the stance phase, while the other intervals (normalized sum bellow THfsr) are considered to be the swing phase. 850 900 950 1000 1050 1100 1150 1200 0 0.2 0.4 0.6 0.8 1 F S R s fsrsum mtt1 mtt5 heel Figure A.2. Gait event detection: threshold is applied to normalized FSR signals in order to detect heel contact (triangular marker pointing downwards), heel off (square marker), and toe off* moments (triangular marker pointing upwards). By applying previously given definitions of gait phases and temporal gait parameters to the detected gait events, SENSY calculates temporal gait parameters for each stride within recorded sequence. 167 Gait Phases Nomenclature Traditional nomenclature: the events taking place during the phase are named, for the most part, according to the events that take place at the foot, for example, heel strike. Ranchos Los Amigos nomenclature: the events taking place during the phases are named, for the most part, according to the purpose of the phase, for example, initial contact. Terms used to describe gait for observational analysis Traditional Rancho Los Amigos STANCE PHASE Heel strike: The beginning of the stance phase when the heel contacts the ground Foot flat: The portion of the stance phase that occurs immediately after the heel strike, when the sole of the foot contacts the ground Midstance: The point at which the body passes directly over the reference extremity Heel off: The point following midstance at which time the heel of the reference extremity leaves the ground Toe-off: The point after heel-off when only the toe of the reference extremity is in contact with the ground. Initial contact: The beginning of the stance phase when the heel or another part of the foot contacts the ground Loading response: The portion of the stance phase from immediately after initial contact until the contralateral extremity leaves the ground. Midstance: The portion of the stance phase that begins when the contralateral extremity leaves the ground and ends when the body is directly over the supporting limb. Terminal stance: The portion of the stance phase from midstance to a point just prior to initial contact of the contralateral extremity. Preswing: The portion of the stance from initial contact of the contralateral extremity to just prior to the liftoff of the reference extremity. This portion includes toe-off. Terms used to describe gait for observational analysis Traditional Rancho Los Amigos SWING PHASE Acceleration: The portion of beginning swing from the moment the toe of the reference extremity leaves the ground to the point when the reference extremity is directly under the body. Midswing: Portion of the swing pause when the reference extremity passes directly below the body. Midswing extends from the end of acceleration to the beginning of deceleration. Initial swing: The portion of the swing from the point when the reference extremity leaves the ground to maximum knee flexion of the same extremity. Midswing: Portion of the swing phase from maximum knee flexion of the reference extremity to a vertical tibial position. 168 Deceleration: The swing portion of the swing phase when the reference extremity is decelerating in preparation for heel strike. Terminal swing: The portion of the swing phase from a vertical position of the tibia of the reference extremity to just prior to initial contact. From O’Sullivan, SB and Schmitz, TJ: Physical Rehabilitation: Assessment and Treatment, ed 3. FA Davis, Philadelphia, 1994, p.169. 169 Kinematics and Kinetics curves: Normal gait pattern Figure A.3 Joint angles of hip, knee, and ankle joint.Typical kinematic curves for healthy gait Figure A.4 Ground reaction forces: load force (vertical reaction force) and friction force (horizontal reaction force). Typical force profiles for healthy gait 170 SENSY – Data acquisition system Software for signal acquisition is developed in CVI LabWindows (National Instruments, USA). This software is compatible with Windows operative system, easy to install and use, and does not require any additional installation. It provides online monitoring and automatic storing of data in ASCII format. Although the extension of files is *.zgz, they can be opened with Notepad or Wordpad, or imported to Matlab or Excel. Graphic interface of the acquisition software is user-friendly and comprises 8 graphic panels showing signals from sensors during acquisition, as shown in Fig. A.5. The upper panels display signals collected from right leg and lower panels those from left leg. Left panels (upper and lower) show data collected from force sensors placed in the shoe insoles (or attached to feet or shoes). The other panels, from left to right, display data from sensor nodes attached to thigh, shank and foot segments, respectively. Figure A.5 – Running SENSY Data acquisition software. The horizontal axes are showing the number of samples and the vertical axes the value of the samples in a 12 bit resolution. Accelerometers integrated in sensor units can be used either in 2g or 6g resolution, which can be selected for each IMU/leg segment individually by sliding 2g/6g bar in the upper right corner of the corresponding panel. For subjects who walk with stronger foot impacts, it is suggested to use 6g for foot sensors so the sensors wouldn’t go to saturation. For upper segments, it is recommended to use 2g (i.e., 171 higher resolution). Although the resolution is changed, the values displayed on graphs and saved in the file are in the same range because they are automatically scaled. An important feature of the software is data delay status bar, placed under each panel. This bar shows status of data transmission for each IMU. The intensity of the delay is shown in colors - green, yellow, and red. The intensity of the bar graph is directly proportional to the data delay. If the delay is increased, the color changes from green through yellow to red. If the delay reaches its maxima, data will be lost and bar graph will be circled with red line (Fig. A.6). Figure A.6 – Data delay status bar: green status showing there is no data delay , yellow status showing there is some transmission delay and the buffer is accumulating unsent data, red status indicates that the buffer is nearly full, and circled red status showing the buffer is full and all future data will be lost. This software is compatible with Windows operative system, easy to install and use, and does not require any additional installation. Most importantly, it is not space or memory demanding, so it can be use from PC with any configuration. Since it provides online monitoring and automatic naming (“date_time”) and storing in its own folder in Program files, there is no risk of loosing data and users can be devoted to observing subject and/or acquired signals, without focusing on managing the software. 172 Curriculum Vitae Name and surname: Milica Djurić-Jovičić Date and place of birth: 19.02.1982., Belgrade, Serbia E-mail: milica.djuric@etf.rs EDUCATION: 2007-2012 PhD studies in Biomedical Engineering (Department for System Control and Signal Processing), School of Electrical Engineering, University of Belgrade, Serbia. 2000-2007 School of Electrical Engineering, University of Belgrade, Serbia, Telecommunications department. 1996-2000 V Belgrade high school, mathematical- scientific curriculum. 1988-1996 Elementary school “Vuk Karadžić”, Belgrade. DIPLOMA WORK: April 2007 “Analysis of ground reaction forces while walking in shoes with different heel heights” within the course Principles of biomedical engineering (supervisor: prof. Dejan Popović) WORK EXPERIENCE: 2007- Research assistant at Signal and systems, School of Electrical Engineering, University of Belgrade, Serbia. 2008-2012 Researcher in the field of motor control and rehabilitation in “Tecnalia Serbia“ (previously called “Fatronik Serbia”), Belgrade, Serbia. Project manager for Wireless sensor system for gait analysis. 173 PROJECTS: 2008-2010 #145041, MNTR, “Functional Electrical Therapy (FET) for forming motor patterns after cerebrovascular insult”, National project funded by Serbian Ministry of Science and Education. PI: Professor Mirjana Popović. 2009-2011 SENSY (Wireless sensor system for gait analysis) in Tecnalia Serbia, Belgrade, Serbia, PI: Milica Djurić- Jovičić. 2009-2011 TEMPUS JP 144537-2008, ”Curricula Reformation and Harmonisation in the Field of Biomedical Engineering”, International project funded by EU, PI: Professor Dejan Popović. 2011-2014 #175016, MNS, “The effects of assistive systems in neurorehabilitation of sensory- motor systems”, National project funded by Serbian Ministry of Science and Education applicated, PI: Professor Mirjana Popović. PUBLICATIONS: International journals: *Djordjević A.R., Djurić M.D., Tošić D.V., Sarkar T.K., “On compact printed-circuit transmission lines”, Microwave and Optical Technology Letters, Nov. 2007, Vol. 49, pp. 2706-2709, ISSN 0895-2477 (M23). Kojović J., Djurić-Jovičić M., Došen S., Popović M.B., Popović D.B., “Sensor-driven four-channel stimulation of paretic leg: Functional electrical walking therapy”, Journal of Neuroscience Methods, June 2009, Volume 181(1), pp. 100-105 (M23). Popović M.B, Djurić-Jovičić M., Petrović I., Radovanović S., Kostić V., “A simple method to assess freezing of gait in Parkinson's disease patients”, Braz J Med Biol Res, September 2010, Volume 43(9), pp. 883-889 (M23). Djurić-Jovičić M., Jovičić N., Popović D.B., “Kinematics of Gait: New Method for Angle Estimation Based on Accelerometers”, Sensors, Nov. 2011, Volume 11(11), pp. 10571-10585 (M21). Djurić-Jovičić M., Jovičić N., Popović D.B., Djordjević A.R., “Nonlinear Optimization for Drift Removal in Estimation of Gait Kinematics Based on Accelerometers”, Journal of Biomechanics, , November 2012, Volume 45(16), pp. 2849-2854. National journals: *Djordjević, A.R., Djurić, M.D., "Printed one-dimensional periodic structures for guiding electromagnetic waves", 8th meeting of the Division of Technical Sciences, Serbian Academy of Sciences and Arts, June 2007, Glas of Serbian Academy of Sciences and Arts, Division of Technical Sciences, vol. CDXI, no. 36, pp. 13-33. 174 International conferences: *Djordjević A.R., Djurić M.D., Niccolai L., “Multiband modem antenna for cellular networks”, Applied Computational Electromagnetics Society (ACES) Conference, Syracuse, N.Y., April 19- 23, 2004., S06P02 (M33). Djurić M., “Automatic Recognition of Gait Phases from Accelerations of Leg Segments”, Proceedings from the 9th Symposium on Neural network Applications in Electrical Engineering, Neurel 2008, September 25-27, 2008, Belgrade, Serbia, ISBN 878-1-4244-2904-2, IEEE Catalog Number: CFP08481-PRT, pp. 121-124 (M33). Djurić-Jovičić M., Milovanović I.P., Jovičić N.S., Popović D.B., “Reproducibility of BUDA Multisensor System for Gait Analysis”, Proceedings of IEEE EUROCON Conference, EUROCON 2009, May 18-23, 2009, St.Petersburg, Russia, pp. 108-11 (M33). Djurić-Jovičić M., Milovanović I.P., Jovičić, N.S., Popović D.B., “Walkaround assisted walking of stroke patients”, Proceedings from the Medical Physics and Biomedical Engineering Conference, September 7-12, 2009 Munich, Germany, pp. 299-301 (M33). Djuric-Jovičić M., Popović M., Petrović I., Radovanović S., Kostić V., “Freezing vs no freezing steps in Parkinsons disease patients”, Abstracts of the XVIII Congress of the International Society of Electrophysiology and Kinesiology, ISEK 2010, 16-19 June 2010, Aalborg, Denmark [CD-ROM]. ed. / Deborah Falla; Dario Farina. Aalborg: Department of Health Science and Technology. Aalborg University. , 2010. Milovanović I., Djurić-Jovičić M., “Polymyography during hemiplegic walking: Implications for control of FES”, Proceedings from the 15th International Functional Electrical Stimulation Society Conference, IFESS 2010, September 8-12, 2010, Vienna, Austria, pp. 203-205 (M22). Djurić-Jovičić M., Jovičić N.S., Milovanovic I., Radovanović S., Kresojevic N., Popović M.B, “Classification of Walking Patterns in Parkinson’s Disease Patients Based on Inertial Sensor Data”, Proceedings from the 10th Symposium on Neural network Applications in Electrical Engineering, Neurel 2010, September 23-25, 2010, Belgrade, Serbia, pp. 3-6 (M33). Djurić-Jovičić M., Jovičić N., Popović D.B., “The Influence of Heel Height on Gait Pattern”, Proceedings from the 5th European Conference of the International Federation for Medical and Biological Engineering, IFMBE 2011, vol. 37, isbn: 978-3-642-23507-8, September 14-18, 2011 Budapest, Hungary, pp. 872-875 (M33). Djurić-Jovičić M., Miler-Jerković V., “Intra-subject stride-to-stride variability: selecting subject’s representative gait pattern”, Proceedings of the 19th Telecommunication forum, TELFOR 2011, November 22-24, 2011, Belgrade, Serbia, pp. 51-54 (M34). Miler-Jerković V., Djurić-Jovičić M., Popović M., “PCA Sensitivity: The Role of Representative and Outlier Strides in Gait Sequence“,. Proceedings from the 11th Symposium on Neural network Applications in Electrical Engineering, Neurel 2012, September 20-22, 2010, Belgrade, Serbia, pp. 123-126 (M33). 175 National conferences: *Djordjević, A.R., Niccolai, L., Olćan, D.I., Djurić, M.D., “Miniature GPS antenna with ring resonator”, Proceedings of XII Telfor, 2004, Belgrade, Serbia. Jovičić N., Djurić M., Popović D.: “Portable Data Acquisition System for Gait Analysis Based on Bluetooth Communication”, Proceedings from the 15th Telecommunication forum, TELFOR 2007, Belgrade, Serbia, pp. 484-487 (M63). Djurić M.D., Došen S., Popović M.B., Popović D.B., “Sensors for Control of Assistive System for Restoring of the Walking”, 52nd Conference for Electronics, Telecommunications, Computers, Automation, and Nuclear Engineering, ETRAN 2008, ME1.1, June 9-12, 2008, Palić, Serbia, http://etran.etf.rs/etran2008/sekcije.htm (M63). Djurić-Jovičić M., Milovanović I., Jovičić N., Radovanović S., “Gait analysis: BUDA vs. GAITRITE”, Proceedings from the 53rd Conference for Electronics, Telecommunications, Computers, Automation, and Nuclear Engineering, ETRAN 2009, ME 2.2., June 15-18, 2009, Vrnjačka Banja, Serbia, http://etran.etf.rs/etran2009/sekcije.htm (M63). Popović M., Djuric-Jovičić M., Petrović I., Radovanović S., Kostić V., “Freezing of gait (FOG) in Parkinson’s disease patients : time analysis”, Book of Abstracts, 9th Congress of Clinical Neurophysiology with International Participation, October 15-17, 2009, Belgrade, Serbia. Naucna KMD, 2009. p. 108-109. Abstract published in Clinical Neurophysiology, Vol. 121, No. 4, 2010, p. e16, No. 48. Djurić-Jovičić M., Milovanović I., Janković M., Dragin A., “Synergies of a gait - Clinical measurements of kinematics and polymiography”, Book of Abstracts from Symposium of clinical neurophysiology with international participation, October 8-9, 2010, Belgrade, Serbia, pp. 49 (M63). Djurić-Jovičić M., Milovanović I., Popović D.B., “Statistical method for gait classification based on data recorded with inertial sensors”, Book of Abstracts from Symposium of clinical neurophysiology with international participation, November 4, 2010, Belgrade, Serbia, pp. 33 (M63). Publications marked with * are not from the area of thesis. LANGUAGES: • English (fluent) • Russian (good) • Italian (good) • Spanish (moderate)