IN
DEGREE PROJECT MEDICAL ENGINEERING, SECOND CYCLE, 30 CREDITS
,
STOCKHOLM SWEDEN 2020
Evaluation of the validity of IMU
sensors measuring wrist angular
velocity by comparison with an
optical motion tracking system
MOGOS TSELETU TESFALDET
KTH ROYAL INSTITUTE OF TECHNOLOGY
Evaluation of the validity of IMU sensors
measuring wrist angular velocity by comparison
with an optical motion tracking system
MOGOS TSELETU TESFALDET
Master of Science Thesis in Medical Engineering Date: October 6, 2020
Supervisor: Farhad Abtahi Reviewer: Mikael Forsman Examiner: Matilda Larsson School of engineering science in Chemistry, Biotechnology, and Health KTH Royal Institute of Technology Swedish Title: Utvärdering av validiteten hos IMU- mätningar av handledshastighet
Abstract
There is a need for objective methods for wrist angular velocity measurements for accurate risk assessments because there is a high frequency of musculoskeletal disorder in workers. The goal of this project was to validate the accuracy of inertial measurement unit sensors to measure the angular velocity. More specifically, the purpose of this master thesis project was to apply an alternative algorithm to compute the markers velocity, other than the one from the optical system that Jenny Wingqvist, and Josephine Lantz used. The project used an experimental data of 10 participants from the previous project done by Jenny Wingqvist and Josephine Lantz. To validate the accuracy, the data of angular velocity of the sensors was compared with the data of angular velocity of markers. The lowest mean value of the root mean square differences value was 23.5 degrees/s during flexion and deviation standard movements at 40 BPM (Beats Per Minute) and the maximum value was 110.5 degrees/s at 140 BPM. The mean value of the correlation coefficients between markers and sensors angular velocities in standard movements of flexion and deviation were 0.85, 0.88, and 0.89 at 40 BPM, 90 BPM, and 140 BPM, respectively. The smallest and the largest mean value of the absolute difference in 50th percentile was found in 40 BPM (19.4±11.3), and 140 BPM (51.2±28.5) respectively. The decorrelation coefficient between the subjects 50th percentile of the angular velocity was 0.91 for the standard movements. The upper limit of agreement for the standard movements was 78.36 degrees/s, while the lower limit of agreement was -13.76 degrees/s. The results show that the error was too large, so there is a need of further research to measure the wrist angular velocity using IMU sensors.
Keywords: Angular velocity, Optical motion tracking, IMU sensors, percentile, correlation
Sammanfattning
Det finns ett behov av objektiva metoder för mätningar av vinkelhastighet från handleder för att noggranna riskbedömningar, på grunda av att det råder en hög frekvens av skelettmuskelproblem hos arbetare. Målet med detta projekt var att bekräfta noggrannheten hos interna mätsensorer för att mäta vinkelhastigheten. Mer specifikt var syftet med detta examensarbete att använda en alternativ algoritm för att beräkna markörhastigheten, annan än den från det optiska systemet som Jenny Wingqvist och Josephine Lantz använde. Projektet använde en experimentell data från tio deltagare från det tidigare projektet gjort av Jenny Wingqvist och Josephine Lantz.
För att bekräfta noggrannheten, jämfördes vinkelhastighetsdata från sensorer med data från vinkelhastighetsmarkörer. Det lägsta medelvärdet av root mean square avvikelserna var 23.5 grader/s under flexion och deviation standardrörelser om 40 slag/minut. Det maximala värdet var 110.5 grader/s vid 140 slag/minut. Medelvärdet av korrelationskoefficienterna mellan markörer och sensorer i standardrörelser av flexion och deviation var 0,85, 0,88 och 0,89 vid 40, 90 och 140 slag/minut respektive. Det minsta och det största medelvärdet av den absoluta differensen i den 50:e percentilen hittades vid 40 slag/minut (19.4±11.3) och 140 slag/minut (51.2±28.5) respektive. Dekorrelationskoefficienten mellan undersökandes 50:e percentil av vinterhastigheten var 0.91 för standardrörelser. Den övre gränsen av överenskommelse för standardrörelser var 78.36 grader/s, medan den lägre gränsen av överenskommelse var -13.76 grader/s. Resultaten visar att felet var mycket stort, varför det finns ett behov av framtida forskning gällande mätning av vinkelhastigheten hos handleder med IMU sensorer.
Nyckelord: Vinkelhastighet, Optiskt rörelsespårning, IMU sensorer, percentil,
Acknowledgments
List of figures
Figure 1: Location of the Markers and Sensors Figure 2: Location of the four markers Figure 3: Bland - Altman plot between markers velocity and sensors velocity
Figure 4: Angular velocity from the Markers in standard movements of flexion
and deviation during 40 BPM
Figure 5: Angular velocity from the Markers in standard movements
of flexion and deviation During 90 BPM
Figure 6: Angular velocity from the Markers in standard movements of flexion and deviation During 140 BPM Figure 7: Bland- Atman plot explained the difference between the markers and sensors 50th percentile of the angular velocity under standard movement of flexion
and deviation in 40 BPM , 90 BPM, and 140 BPM Figure 8: Correlation between angular velocity of the sensors and markers In Standard
movements of flexion and deviation in 40 BPM ,90BPM, and140 BPM
Figure 9: IMU Sensors
List of Tables
Table 1: Data base Table 2: Correlation coeffient between markers and sensors angular velocity
Table 3: RMSD between Sensors and markers velocity of standard movement at 40 BPM, 90BPM, and 140 BPM. Table 4: The mean and standard deviation of the participants in 10th, 50th, and 90th
Table of contents
Chapter 1 ... 1
1. Introduction ... 1
1.1 Aim ... 2
Chapter 2 ... 3
2. Materials and methods ... 3
2.1 Participants ... 3
2.2 Wrist Motion Measurement systems ... 3
2.2.1. Inertial measurement units ... 3
2.2.2 Optical Motion Tracking System ... 4
2.3 Wrist Angular Velocity Calculation ... 4
2.4 Data Processing ... 6
2.5 Statistical Analysis ... 7
2.5.1 Bland-Altman plot ... 7
2.5.2 Correlation Plot ... 8
2.5.3 Low Pass Butterworth Filter ... 8
Chapter 3 ... 9
3 Results ... 9
3.1 Angular velocity of the markers and sensors ... 9
3.2 Results from Bland-Altman Plot ... 10
3.3 Results from Correlation plot ... 11
3.4 Results of mean and standard deviation during 10th, 50th, and 90th percentiles ... 12
3.5 Results from Root Mean Squar Deviation (RMSD) ... 13
Chapter 4 ... 14
4.1 Discussion of results ... 14
Chapter 5 ... 15
5.1 Conclusions and future work ... 15
Appendix A ... 16
State of the Art ... 16
A.1 Anatomy of the hand and wrist ... 16
A. 2 Wrist Movement ... 16
A.2.1 Flexion and Extension ... 16
A.2.2 Radial and ulnar deviation ... 16
A.3 Angular velocity ... 17
A.4 Motion measurements ... 17
A.4.2 Mechanical Gyroscopes ... 19
A.4.2 Micro electromechanical system (MEMS) gyroscopes ... 19
A.5 Optical tracking system ... 21
A.6 Measurement methods in Ergonomics ... 21
1
Chapter 1
1. Introduction
Musculoskeletal disorder is a work-related physical disorder that occurs in different parts of the companies around the world. Today in western industrial countries, it affects thousands of workers and cost employers millions of dollars [1]. These problems are due to different reasons; generally, the leading causes are high demand workloads, rapid movements, and stressful work positions. These physical risk factors are due to various work activities related to the development of work-related musculoskeletal disorders [1],[2].
In addition to work-related risk factors, there are individual risk factors such as increasing age, smoking, and high body mass index that cause workplace problems. Besides, other activities outside the work like family take care, and opportunity for relaxation, recovery, and exercise may have a significant effect in all the individual risk factors directly or indirectly connected with the physical workloads [1].
Measurement methods are needed to find a solution for the above problems and design a sustainable work environment. These measurement methods can measure the wrist angles and angular velocity easily, accurately, and efficiently. The advantage of such measurement methods is to improve the optimization of the workplace design and to take decisions on recommended maximum values that would not exceed.
During the past few years, an ergonomic movement was measured by an instrument called Biometrics Goniometer. These measurements can provide high accuracy, but it is costly and needs individual staff for data analysis. Therefore, a new method should be set up, inexpensive, offers excellent efficiency, and measures the muscular activities or movement.
However, such problems should be solved using mechanical sensors on the human body, these sensors are called inertial measurement units (IMU)[1]. These sensors are containing a 3D-accelerometer, 3D gyroscope, and magnetometer. These devices are small, low cost, compact, and low processing power [3].
2
1.1
Aim
This thesis work was aimed to estimate the wrist angular velocity from IMU sensors. More specifically, the purpose of this master thesis project was to apply an alternative algorithm to compute the markers velocity, other than the one from the optical system that Jenny Wingqvist, Josephine Lantz used [3].
Then the goal was to minimize the difference between these calculated markers velocity and the IMU sensors. Finally, this thesis work will help develop methods so that practitioners, non-researchers, can do a valid risk assessment on workplaces and provide a solution.
These experimental works were done by Jenny Wingqvist and Josephine Lantz in uMove, an optical motion tracking lab at Karolinska institute.
3
Chapter 2
2. Materials and methods
2.1 Participants
The database consists of twelve participants, five men, and seven women, who participated in the experimental study, with the mean age of 25 years (range 19 -52), height 172.9 cm, and weight 65kg. The mean age, weight, and height for the men were 29 years, 75 kg, and 179.5 cm, respectively, and for the women were 23 years, 56.4 kg, and 167.3 cm, respectively, see table1. Before the experiment, all the information needed was given to the main participants and even the specific goal of the validation experiment's thesis project, and they had signed informed consent. The Regional Ethics Committee approved the experimental study in Stockholm (DNR 2018 / 1756-32) [3]. Table 1. Data base
Participants Number of Participants
Mean age Weight Height
Men Women 5 7 29 years 23 years 75 kg 56.4 kg 179.5 cm 167.3 cm
2.2 Wrist Motion Measurement systems
For testing the accuracy, measurements were taken from Inertial measurement units (IMU) on the hand, and on the forearm, and optical tracking system, because both these devices have high precision and accuracy.
Finally, the angular velocity from the IMU sensors compared with the data of wrist angular velocity from the markers.
2.2.1. Inertial measurement units
4
2.2.2 Optical Motion Tracking System
The devices placed on the study participant were markers and sensors. The optical motion tracking system used in the experiment was Vicon Nexus 2.8.1 [15]. These devices recorded the coordinates of each marker in the X, Y, and Z-axis. The markers on the middle finger, elbow, and the two markers between the wrist were used to calculate the angle of flexion and deviation. Then the wrist angular velocity from the markers data was calculated by change of flexion and deviation angle per change of time.
2.3 Wrist Angular Velocity Calculation
There were eleven markers used in the project, and they were located on the clavicle, sternum, right shoulder, right elbow: ulnar, and radial, right wrist, right middle finger, seventh cervical, and tenth thoracic. The right shoulder blade marker was added as a reference marker. The Sensors were placed on the middle arm, forearm, upper arm, and thoracic.
In this master thesis project, only the coordinates data from the markers on the middle finger, elbow, and the two markers on the wrist were interested to calculate markers velocity. For sensors velocity, the sensors on the hand and forearm were used, See figure 1 below. As shown in figure 1, the markers ware with black and white colours, and the sensors were bigger than the markers with blue colour.
Figure 1. Location of the Markers and Sensors [3]
5
Figure 2. Location of the four markers
From figure 2, the following information were obtained and used to calculate for flexion and deviation angle of the wrist.
• M1 indicates Marker 1 on the wrist with coordinates (x1, y1, z1) • M2 represents Marker 2 on the wrist with coordinates (x2, y2, z2) • M3 is Marker 3 on the elbow with coordinates (x3, y3, z3) and
• M4 indicates Marker 4 on the middle finger with coordinates (x4, y4, z4)
The point Q (d, e, f) is the midpoint between the wrist marker (M1) and wrist marker (M2), and it was calculated by applying the middle point formula between these two points in a line segment
.
𝑄(𝑑, 𝑒, 𝑓) = (𝑋1+𝑋2 2 ) , ( 𝑌1+𝑌2 2 ) , ( 𝑍1+𝑍2 2 ) (2.1)
Then, the planes Q (d, e, f), M4 (x4,y4,z4) and M3 (x3,y3,z3) formed a triangle with sides a, b, and c as showed in figure 2. The three sides were calculated using the distance formula between the two points as follows. 𝑎 = √((𝑋4− 𝑋3)2+ (𝑌4− 𝑌3)2+ (𝑍4− 𝑍3)2) (2.2) 𝑏 = √((𝑋3− 𝑑)2+ (𝑌 3− 𝑒)2 + (𝑍3− 𝑓)2) (2.3) 𝑐 = √((𝑋4− 𝑑)2+ (𝑌 4− 𝑒)2+ (𝑍4− 𝑓)2) (2.4)
Using the cosine law, angle of flexion and deviation was calculated by applying the equation 2.5, and the three sides (a, b, and c) of a triangle were found in equations 2.2, 2.3, and 2.4, respectively. Angle of flexion and deviation , 𝜃 = 𝐶𝑜𝑠−1(𝑐2+𝑏2−𝑎2
6
The angle ( 𝜃) was calculated in radians, so it should be converted to degrees and the conversion from radians to degrees was giving in equation 2.6.
1 rad = 1800 ÷ 𝜋 = 57.2957795130 (2.6)
From equation 2.7below, the markers angular velocity of the wrist ( 𝑉 ) is calculated by the change of angle of flexion and deviation per change of time,and it is giving in degrees/second.
𝑉 = (Δ𝜃/Δ𝑡) (2.7) In the experiment, the sampling frequency used was 100 Hz, and frequency (f) is defined as a number of samples per unit second. Therefore, the sampling time (t) was given by 1/100, which is 0.01 seconds.
The angular velocity from the sensors was calculated from the gyroscopes data on the hand and forearm, by taking the difference between velocity of the middle hand sensor and velocity of the forearm sensor. These gyroscope data include data of angular velocity along x, y, and z directions. Then, total angular velocity of sensors could be calculated by taking the absolute amount of the velocities along the three axes, see equation 2.8.
Total velocity from sensors (VT) = √𝑥2 + 𝑦2+ 𝑧2 (2.8)
Where (X) is the angular velocity difference between the gyroscope data of the hand and forearm along x-axis, and (Y) is the value of gyroscope data of the hand minus the gyroscope data of the forearm along y-axis. Similarly (Z) is calculated by taking the difference between the angular velocity of the hand and forearm along z-axis.
2.4 Data Processing
The data from the sensors (only sensors on the middle hand and forearm data) and data from the markers (only the data of coordinates of elbow marker, middle finger marker and two wrist markers), were processed in MATLAB R2019b.
In Matlab, the data file from the middle hand sensors and forearm sensors were loaded, any data lost exists in the file were removed using NaN values function in MATLAB. Since the file size of the two gyroscopes data were in different sizes, a rough interval was selected for the two jumps for the hand data and forearm data. Then the gyroscope data were interpolated, the sampling frequency were applied for both data hand and forearm. When the process of interpolation was finished, the data from the middle hand sensor and forearm sensor were cross correlated for their data to be in phase.
Then, another interval was selected for calculating the velocity of sensors and markers, because they have different file sizes. This interval corresponds to an exercise in standard movements. Finally, the angular velocity of the sensors was calculated using equation 2.8.
7
that, the velocity of the markers was calculated, which was done by diving the change of an angle to delta time, as shown in equation 2.7.
When the angular velocity of both markers and sensors were calculated, low pass filtered, 6th order
Butterworth filter with a cut-off frequency of 5 Hz. For better results, the values of velocity for each system exceeds 550 degrees / second, were converted to 550 degrees/second. Next, the velocity of both systems was correlated, and the result of correlation was shown in figure 4 ,5 and 6.
Finally, different statistical analysis method was evaluated: correlation coefficient see table 2, RMSD see table 3, and 10th, 50th, and 90th percentile of standard movements of flexion and deviation at 40
BPM, 90 BPM, and 140 BPM. Besides, the mean of markers and mean difference between markers and sensors were calculated.
2.5 Statistical Analysis
To compare the angular velocity obtained from the sensors and markers under static condition, Bland-Altman plot were applied to test the difference between the two data and correlation coefficients were calculated. First, the signals from the markers and sensors were synchronized, to check the difference between markers and sensors velocity. Secondly, the mean and minimum correlation coefficient, and mean and maximum root mean square differences were calculated.
2.5.1 Bland-Altman plot
Bland-Altman plot specifies the similarities between two different quantitative measurements [32]. Its graph is plotted by taking the difference between the two measurements along the Y axis, and the mean of the two measurements along the X axis . In this project, the two measurements were IMU sensors and optical tracking motion systems. Then the result from the two measurements, that is the markers angular velocity (Vm), and sensors angular velocity (Vs) were considered to describe the agreement between the two instrument measurements. These difference between the wrist angular velocity of the markers and velocity of the sensors were plotted in the Y axis, and the average of the marker’s velocity, and sensors velocity was drawn in X axis, see figure 3 [32].
In Bland -Altman plot, the mean difference value for the two measurements was located on the middle line of the graph, and the limits of an agreements lines was on other side. This limit of an agreement calculated by taking the mean and standard deviation (SD) of the difference between sensors angular velocity and markers velocity [32]. In this project, the 50th percentiles of the sensors
and 50th percentile of the markers used to plot the Bland Altman plot between sensors and markers
velocity.
8
For plotting Bland-Alman plot ,the following data from the markers and sensors were used :-
• 50th percentile value of the markers total velocity during standard movements of flexion and deviation at a rate of 40 BPM
• 50th percentile difference between the sensors and markers total velocity at a rate of 40 BPM
• 50th percentile of the markers total velocity during standard movements of flexion and deviation at 90 BPM
• 50th percentile difference between the sensors and markers total velocity at 90 BPM
• 50th percentile difference between the sensors and markers total velocity at140 BPM
• 50th percentile of the markers total velocity during standard movements of flexion and deviation at 140 BPM
2.5.2 Correlation Plot
Correlation plot is one of the statistical analysis methods, describes how the two measurements are related to each other. In this project, the sensors velocity was correlated with the wrist angular velocity of the markers, and as a result the correlation coefficient determined the strong linear relationship between these two measurements. The value of correlation coefficient (r) varies from -1 to +-1, and when the value of the correlation coefficient is closer to +-1, shows that there is great linear relationship between markers velocity and sensors velocity) [32].
In this project the data needed to plot the correlation plot were listed below :- • 50th percentile total markers data at a rate of 40 BPM
• 50th percentile total markers data at a rate of 90 BPM
• 50th percentile total markers data at a rate of 140 BPM
• 50th percentile total sensors data at a rate of 40 BPM
• 50th percentile total sensors data at a rate of 90 BPM
• 50th percentile total sensors data at a rate of 140 BPM
2.5.3 Low Pass Butterworth Filter
9
Chapter 3
3 Results
3.1
Angular velocity of the markers and sensors
Figures 4, 5, and 6 described the angular velocity of the markers and sensors after being correlated with each other during standard movement of flexion and deviation at a rate of 40 BPM, 90 BPM, and 140 BPM. Data from the left of the graph and data to the right of the graph described the standard movement of flexion and deviation respectivelly. As shown in the figure below, the amplitude of markers velocity is smaller than the amplitude of sensors velocity during flexion movement.
Figure 4. Angular velocity of the markers and sensors during standard movement of flexion and deviation in 40 BPM.
10
Figure 6. Angular velocity of the markers and sensors during standard movement of flexion and deviation in 140 BPM.
3.2 Results from Bland-Altman Plot
Figure 7 shows the Bland - Atman plot for the standard movements of flexion and deviation at 40 BPM, 90 BPM, and 140 BPM.The mean difference was 32.3 degrees/s, and the upper and lower limits of angreement were 78.36 degrees/s and -13.76 degrees/s respectively.
11
3.3 Results from Correlation plot
The correlation coefficients was 0.9127 between the 50th percentile of the sensors and markers, as
shown in figure 8 below.
12
From table 2 below, the minimum mean correlation coeffients is 0.85 degrees/s at 40 BPM, and the maximum mean is 0.89 degrees/s at 140 BPM. So, as the speed increase, the mean value also increased.
Table 2. Correlation coeffient between markers and sensors angular velocity
Number of Participants 40 BPM 90 BPM 140 BPM FP4 0.82166 0.7944 0.9317 FP5 0.9170 0.9261 0.7746 FP6 0.9137 0.9549 0.9393 FP7 0.9218 0.9253 0.9369 FP8 0.7435 0.7831 0.8380 FP9 0.7328 0.8189 0.9179 FP10 0.8509 0.90757 0.8840 FP11 0.8484 0.9388 0.9522 FP12 0.8260 0.8837 0.9614 FP13 0.9015 0.8368 0.7702 Mean 0.85 0.88 0.89 Minimum Value 0.7328 0.7831 0.7702
3.4 Results of mean and standard deviation during 10
th, 50
th, and 90
thpercentiles
From table 3 below, the mean value varied between 3.6 degrees/s and 162.6 degrees/s, while the standard deviation varied between 1.76 degrees/s and 38.4 degrees/s. The largest difference was appeared in 90th percentile at a rate of 140 BPM of the standard movements, in contrast the smallestdifference was found in 10th percentile at a rate of 40 BPM.
13
3.5 Results from Root Mean Squar Deviation (RMSD)
Table 4 below described the RMSD value between Sensors and markers velocity of standard movements at a rate of 40 BPM, 90 BPM, and 140 BPM. As shown in figure 4, the maximum value of RMSD was 110.46 degrees/s at 140 BPM and the maximum mean was 99.2 degrees/s at 140 BPM.
Therefore, as the speed increased the mean value increased and the maximum value also increased.
14
Chapter 4
4.1 Discussion of results
In figure 4, figure 5, and figure 6, which describes the angular velocity signal of the markers and sensors in standard movement of flexion and deviation at 40 BPM, 90 BPM, and 140 BPM. As shown in the figure, the curves correlated along the time axis, during flexion movement the markers' amplitude was lower than the sensors amplitude, and the difference (error) is too large. These is one of the problems in the project, and it is difficult to mention the cause of the problem. These is the future work needed in this project. Secondly as we compare the speeds between both systems, sensors have a higher rate than the markers, the reason is may be one of the markers were placed on the elbow, which is far from the wrist.
In figure 7, the Bland-Altman plots showed the difference between the sensor’s velocity and markers velocity. The result of the Bland-Altman plot showed that the errors were large and positive, the reason is because there is an amplitude difference between sensors velocity and markers velocity. The amplitude of the sensor’s velocity is double as markers velocity. As shown in the figure 7, the mean difference between the 50th percentiles were 32.3 degrees/s, which was very large and not
good.
From figure 8, the correlation coefficient between the 50th percentile of the markers and sensors was
0.91273, which is closer to 1, and showed the strong linear relationship between both systems. In table 2, the mean correlation coefficient values are 0.85 degrees/s, 0.88 degrees/s, and 0.89 degrees/s at 40 BPM, 90 BPM, and 140 BPM, respectively. As we compared the results, they were almost the same.
In table 3, which explained the maximum, minimum, and mean value of RMSD for the standard movements at 40 BPM, 90 BPM, and 140 BPM. The highest mean value is 99.2 obtained at 140 BPM and the minimum value is 34.9 at a rate of 40 BPM.
In previous project done by Jenny Wingqvist and Josephine Lantz[3], they evaluated the accuracy of the sensors , and computed the angular velocity of the markers by taking the absolute value of the gyroscope velocity along x, y, and z axis from the optical tracking system. In this project, I used the coordinates of the markers for calculating the angular velocity of the markers.
From the previous project, the correlation coefficient between the 50th percentile between the
markers and sensors was 0.9565, which showed that the two measurements have linear relationship between them. But in this project, the correlation coefficient was 0. 9127, which is less than 0. 9565. These two values are closer to 1, which showed that there is linear relationship between the two measurements.
15
Chapter 5
5.1 Conclusions and future work
In this project, I compared the sensors velocity with markers velocity of the wrist in terms of recognition of wrist movements (flexion /extension and deviation). From experiments with 10 participants, three standard measurements were taken in 40 BPM, 90 BPM, and 140 BPM for testing the accuracy of IMU sensors. However, based on the statistical analysis, there were more significant difference were detected between the sensor’s velocity and velocity of the markers during flexion movements of the wrist. These showed that the errors were large during flexion movement.
16
Appendix A
State of the Art
A.1 Anatomy of the hand and wrist
The human hand structure is very complex due to more numbers of different functional joints and can perform excellent mobility and flexibility in the digits [4]. It helps us to grip, grasp, and form some good movement. In human anatomy, the wrist is a complex joint that bridges the forearm to the hand. It includes joints, multiple bones, ligaments, and muscles [4][5]. The wrist can flex and extend the fingers, and this movement is due to the internal and external muscle groups.
A place more than one bones meet is known as a joint. The radiocarpal joint articulates the convex proximal row of carpal bones and the concave radius and disk. Between the distal rows and proximal carpals are midcarpal joints [5]. Both the radiocarpal and midcarpal joints are responsible for the wrist movements: flexion, extension, radial deviation, and ulnar deviation—the carpometacarpal located between distal carpals and metaphalangeal [4]. The bones found on the fingers are metacarpophalangeal and interphalangeal joints.
Trapezoid metacarpal joint found between the trapezium and the first metacarpals. It is described as a critical joint of the hand because of its exceptional mobility. As a result, there is a pure movement of flexion, extension, abduction, and adduction [4].
The two muscle groups responsible for producing wrist movement are known as extrinsic and intrinsic muscles [7]. External and inherent tissue is found on the forearm and hand, respectively. The main difference between these two muscles is that extrinsic muscle is considerable and gives more strength, but intrinsic muscles are small and provide quick movement for the hand and fingers[7].
A. 2 Wrist Movement
To make more accessible the hand function, the upper limp must be mobile and stable. If the upper limp is not portable, the hand is not functioning and positions in the right place. The wrist is very mobile, move and rotate in different directions. The stability of the wrist depends on the interaction between the different articular surfaces and soft tissue parts [8]. The wrist is capable of three rotational movements: flexion, extension, and deviation [8]. Wrist rotation in the sagittal plane (Y-axis) corresponds to flexion and extension. The radial and ulnar deviation along the coronal plane (X-axis). These wrist movements are due to passive external and internal muscle groups, [7] [8].
A.2.1 Flexion and Extension
Wrist flexion is a movement bending the palm down towards the wrist, while an extension of the wrist is raising the hand. The range of motion of the wrist extension is between 70-80 0, while for the
flexion are 75-85 0 [4]. This range of motion varies among different groups. The finger joints' position
can influence it due to the constant length of the extrinsic finger of the flexion and extension muscle. When the fingers expanded, the flexion range is higher than the extension range [28].
A.2.2 Radial and ulnar deviation
flexion-17
extension. The radial deviation is described as rotation when the thump is very closer to the forearm's radial aspects. The average range of the active radial-ulnar deviation is 210 and 28 0
respectively [ 8]. The range of motion for the radial deviation is 15-20 0, while for the ulnar deviation
is 30-40 0 [4].
A.3 Angular velocity
Angular velocity is defined as a time rate at which an object rotates, revolves around an axis, or angular displacement between the two bodies changes. It describes how fast the angular position and orientation of an object changes with time, see equation A.1 below. The symbol for the Angular velocity is a small omega (ω) and measured in angles per unit time. The SI unit of angular velocity is radian per second [4]. The angular velocity and the rotation angle related through these equations.
𝜔 = (Δ𝜃/Δ𝑡) (A.1) wrist angular velocity is calculated because musculoskeletal disorders in hand and wrist regions
increased among the workers in different companies and industries. There is a direct relation between the angular velocities and the work injuries. Finally, the angular velocity of the sensors would compare with the marker’s angular velocity[1].
A.4 Motion measurements
History of IMU started in the 1930s, where it was used in aircraft navigation and large devices. Because of its constraints in size, cost, and power consumption, IMU usage is restricted to bulk applications and unpopular to smaller size devices [9]. In recent years, IMU sensors was developed with a good feature of low cost, compact, and low processing power. As a result, its application areas increased from industrial to ergonomic and from biomechanical to life science [10]. In this thesis project, the LP-Research Motion Sensor Bluetooth version 2(LPMS-B2) series is nine-axis Bluetooth IMU sensors and used for direct measurement during wrist movement [13]. These are very versatile, performing an accurate and high-speed orientation, acceleration, and relative displacement measurements [13][14].
IMU sensors are MEMS device which was used to measure the angular velocity, orientation, acceleration, and the gravitation force. It consists of a three-axial accelerometer, three axial gyroscopes, and a magnetometer [14]. Gyroscope sensors are used to measure the angular velocity
, while the magnetometer is used to measure the magnetic field strength. The accelerometer measured the sensor's linear acceleration by finding small vibrations when an
18
Figure 9. IMU SENSORS IMU sensors wanted to be calibrated from day to day, because they should be ensuring that it is working properly. Magnetometers and accelerometer are calibrated using the magnetic field of the earth and gravity of the earth respectively, by measuring some data points and then take the mean of the data points. Besides, the gyroscopes are calibrated by rotating a known angular velocity. As show in figure 9, these devices form three different axes and a total of nine degree of freedom. The degree of freedom is vital for determining of how many independent parameters can exist in a system. The merit of magnetometer is to measure the yaw angle; thus, it can be calibrated with the gyroscope data to improve the big drift issue [25].
A.4.1 Gyroscopes
The world gyroscope derived from the ancient Greek language, which means the precession motion [21] [26]. It serves as for measuring the rotational rate (angular velocity) around the three axes: pitch (x-axis), roll (y-axis), and yaw (z-axis).These MEMS gyroscopes reads the angular rate when the frame is rotating in any direction, because it is placed around the frame. It has great advantage than other devices for different applications, the reason is that they can be used alone or combine with other complex systems, such as Inertial Measurement Unit [23], Inertial Navigation System, Gyrocompass [22] and Attitude heading reference system [26].
19
A.4.2 Mechanical Gyroscopes
A mechanical gyroscope is known as a displacement gyroscope and rate gyroscope. It consists of spinning mass rotates around its axis, rigid frame with rotating bearing, and sensing systems. According to Physicist Léon Foucault's studies (1852), when the mass is rotating around its axis, it remains parallel to itself to oppose any attempt to change its orientation [26]. Mechanical gyroscopes show several physical phenomena, for example, a change in the orientation of the rotational axis of the freely moving object, and this process is called precession [25]. The main principle application of the mechanical gyroscopes is, when the spinning mass tends to keep its angular velocity concerning an inertial reference frame and when the external torque applied to the mass, its rotational axis undergoes a precision motion at a constant angular speed, which travels along the direction of applied torque. In contrast, the rotating part of the centre of the mass is not affected by the external force [26].
The physical phenomenon of a mechanical gyroscope explained in equations A.2, and A.3 below ∁𝛾 = −Ι𝛀𝝎z (A.2) ∁𝑧 =𝙸𝛺𝟂𝒚 (A.3)
Where I is the mass moment of inertia of the spinning mass, 𝟂y and 𝟂z are the precision speed in y
and z axis respectively, Cy and Cz are the torque acting along y and z axis respectively and angular velocity of the spinning mass along the rotational axis, figure 10.
Figure 10. Elements of Mechanical gyroscopes and main parameters, © 2017 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/), [25]
A.4.2 Micro electromechanical system (MEMS) gyroscopes
MEMS is defined as motion sensors that can detect and measure the angular motion of an object. These classes of gyroscope use a vibrating mechanical element as a sensing element for detecting the angular velocity. MEMS gyroscopes measure the angular velocity in one axis, two-axis, and three axes. MEMS gyroscopes is used in different applications, such as in medical applications, military applications, electronics, and industrial applications, [26]. These gyroscopes are based on the transfer of energy between two vibration modes caused by Coriolis's acceleration [26].
20
Figure 11. A. Coriolis acceleration acting on a moving particle b). Mass spring model of MEMS
gyroscope,[26]. © 2017 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/)
MEMS gyroscopes depend on the Coriolis force to measure the angular velocity. The angular velocity and Coriolis acceleration direct proportional to each other. From figure 11 a, Coriolis force's effect should be described through a rotating system that both the object and observer are in the same fixed reference system. If the object is rotating with the angular velocity (Ω) along the x-axis, the observer stands on the z-axis and sees the purpose of turning along the z-axis. As a result of these effects, the mathematical relation between the Coriolis acceleration(ac) and angular velocity is given below.
ac = 2 𝑣Ω (A.4)
In figure11 b, the primary excitation mode is excited along driving axis (Y) by applying a force ( F drive ), while the secondary mode is excited along X (sense axis) by applying the Coriolis force [25]. The general equation for the mass spring model of MEMES gyroscopes is giving below.
m(d2y/dx2 ) = - Κ yy - ∁y 𝑑𝑦 𝑑𝑥 + Fdrive (A.5) m𝑑𝑦 𝑑𝑥 = - Κzz – ∁z 𝑑𝑧 𝑑𝑥 + Ϝcz (A.6) FZ = 2m Ω v (A.7)
Cy and CZ are the damping coefficients, while
Κ
y andΚ
z are elastic stickiness parameters, Ω is theangular velocity,
F
z is Coriolis force, andF
drive are applying force. The displacement along z axis is21
A.5 Optical tracking system
Gold standard data is a method related to something that has already been checked, known standard, and gives good results and accuracy [14]. During motion detection, optical tracking systems serve as a gold standard. It is a method where the cameras were used to detect the wrist movement of the human hand. In this process, the cameras can record the markers' position, angle, and angular velocity along with the three coordinate systems (X, Y, Z) concerning the time axis [15]. The computer software programming Vicon Nexus 2.8.1 and the camera Vicon MX40 is essential in the motion detection method.
A.6 Measurement methods in Ergonomics
According to the International Ergonomics Association, the concept of Ergonomics defined as follows: Ergonomics is described as a scientific discipline concerned with understanding the interaction between human being, and the profession that applies theory, [33] principle, data, and methods to design to optimize human well-being and overall system performance [18].
Three methods were used to make measurements: self-reports, observational, and direct measurements to reduce the risk of strain injuries in load ergonomics. Self-reporting is a method of collecting data, where the worker provides information about the work, easy to use, and cheap. The observation method describes how another person observes the angles and posture of the worker. It is difficult to measure the small and fast movement in both self-reporting and observation techniques because they have very low accuracy. In contrast to the above two methods, the direct measurement method provides an accurate result, but costly and additionally, it needs advanced technology for the data analysing [19].
22
B. Bibliography
[1] Balogn,l., Arvidsson, I., Björk, J. Et al. “Work-related neck and upper limb disorders- quantitative exposure-response relationships adjusted for personal characteristics and psychosocial conditions Musculoskeletal Disorder (2019). https://doi.org/10.1186/s12891-019-2491-6.
[2] Nancy N. Byl, Mary F. Barbe, Carolyn Byl Dolan, Grant Glass, Chapter 27 – RepetitiveStress Pathology, Soft Tissue. Pathology and Intervention in MusculoskeletalRehabilitation (Second Edition), W.B. Saunders, (2016), ISBN 9780323310727 ,https://doi.org/10.1016/B978-0-323-31072-7.00027-0.
[3] Jenny Wingqvist, Josephine Lantz, et al. Evaluation of the accuracy of IMU sensors measuring angular velocity compared with optical motion tracking (2019).
[4] Wadsworth, C. (1983). Clinical anatomy and mechanics of the wrist and hand. The Journal of Orthopaedic and Sports Physical Therapy, 4(4), 206-216.
[5] Kijima, Y., & Viegas, S. Wrist Anatomy and Biomechanics. Journal of Hand Surgery, (2009). [6] An, K., Berger, R., & Cooney, W. Biomechanics of the Wrist Joint (1991).
[7] Andris Freivalds ,Biomechanics of the upper limp: Mechanics, Modelling, and Musculoskeletal injuries (2004).
[8] Apergis, E. Fracture-Dislocations of the Wrist (1st ed. 2013. ed.) (2013).
[9] Kai Daniel oberländer, Dr, rer, nat, Inertial Measurement Unit (IMU) Technology (September 2015).
[10] Norhafizan M. Khari, Reviews on various Inertial Measurement Unit (IMU) sensorApplication (2 December 2013).
[11] Aditya Chaudhry, Chris shin, Alex skillin, Derek witlpalek, Inertial measurement units, (12 November 2018).
[12] Life performance research,” Introducing IMU core,” [Online].
Available: https://lpresearch.com/imucore-sensor-fusion/. (15 -03- 2020).
[13] The LP-RESEARCH,” LPMS-B2: 9-Axis Inertial Measurement Unit (IMU) / AHRS with Bluetooth Classic and BLE Connectivity,” 2019. [Online]. Available: https://lpresearch.com/lpms-b2/. (16 03 2020).
[14] Cardoso, Jefferson, Pereira, Ligia, Iversen, Maura, Ramos Adilson, what is gold standard and what is ground truth? (2014-11-09).
[15] Vicon motion, Plug-in Gait Reference Guide, (10-08- 2017). [16] Ergonomic products and service, wrist movement,
23
[17] F. Schuind et al, Advanced in the biomechanics of the human hand, (1994).
[18] Chartered Institute of Ergonomics & Human Factors, what is ERGONOMICS, (2017). [19] Liyun Yang, Wilhelmus J.A. Grooten, Mikael Forsman, An iPhone application for upper arm posture and movement measurements, Elsevier Ltd, (12 February 2017).
[20] BIOPAC Systems, Goniometers and Torsiometers, ( 10 21 2014), www.biopac.com. [21] NASA Available online: https://solarsystem.nasa.gov/news/2005/03/14/ brief history of gyroscopes(accessed on 24 June 2017).
[22] Elliott-Laboratories. The Anschutz Gyrocompass and Gyroscope Engineering; Wexford College Press: Kiel, Germany, (2003); pp. 7–24. ISBN 978-1-929148-12-7.
[23] Honeywell. Available online: https://aerospace.honeywell.com/en/products/navigation-and-
sensors/ inertial-measurement-units (accessed on 24 June 2017).
[24] Inertial Labs. Available online: https://inertiallabs.com/ahrs.html (accessed on 24 June 2017). [25] Passaro, Vittorio M. N, Cuccovillo, Antonello, Vaiani, Lorenzo, DeCarlo, Martino,& Campanella, Carlo Edoardo. Gyroscope Technology and Applications: A Review in the Industrial Perspective. Sensors (Basel, Switzerland), 17(10), 2284. (2017).
[26] Kempe, Volker. Inertial MEMS - Principles and Practice - 8.1 Some Basic Principles, (2011) Cambridge University Press. Retrieved
https://app.knovel.com/hotlink/pdf/id:kt008MUHJC/inertial-mems-principles/some basic-principles. [27] Seel, Thomas, Raisch, Jörg, & Schauer, Thomas. IMU-Based Joint Angle Measurement for Gait Analysis. Sensors (Basel, Switzerland), 14(4), 6891-6909. (2014).
[28] Rechy-Ramirez, Ericka Janet, Marin-Hernandez, Antonio, & Rios-Figueroa, Homero Vladimir. A human–computer interface for wrist rehabilitation: A pilot study using commercial sensors to detect wrist movements (2019). The Visual Computer, 35(1), 41-55.
[29] Greiff, P ,Boxenhom, B.; Niles, L. Silicon Monolithic Micromechanical Gyroscope. In
Proceedings of the International Conference on Solid-State Sensors and Actuators, San Francisco, CA, USA, (24–27 June 1991); pp. 966–968.
[30] Yazdi, N , Ayazi, F , Najafi, K. Micromachined inertial sensors. (Proc. IEEE 1998),86, 1640–1659. [31] Sensors 2013 sensors - MDPI. https://www.mdpi.com/1424-8220/13/8/11114/pdf
doi:10.3390/s130811114.
[32] David Giavarina , Understandign Bland Altman Analysis ,(2015) , DOI: https://doi.org/10.11613/BM.2015.015 .
24
APENDIXI B
MATLAB Code
%% Data used for calculating Sensors angular velocity
close all clear all clc
Read Data, Data from the sensors on hand and forearm
% data_forearm = readtable('FP1_190405_092508_RightForearm_node_data.csv'); % data_hand = readtable('FP1_190405_092508_Right Hand_node_data.csv'); % data_forearm = readtable('FP4_190405_113306_RightForearm_node_data.csv'); % data_hand = readtable('FP4_190405_113306_Right Hand_node_data.csv'); % data_forearm = readtable('FP5_190411_085855_RightForearm_node_data.csv'); % data_hand = readtable ('FP5_190411_085855_Right Hand_node_data.csv'); % data_forearm = readtable('FP6_190411_092901_RightForearm_node_data.csv'); % data_hand = readtable ('FP6_190411_092901_Right Hand_node_data.csv'); % data_forearm = readtable('FP7_190411_095742_RightForearm_node_data.csv'); % data_hand = readtable ('FP7_190411_095742_Right Hand_node_data.csv'); % data_forearm = readtable('FP8_190411_104010_RightForearm_node_data.csv'); % data_hand = readtable ('FP8_190411_104010_Right Hand_node_data.csv'); % data_forearm = readtable('FP9_190411_111236_RightForearm_node_data.csv'); % data_hand = readtable('FP9_190411_111236_Right Hand_node_data.csv');
% data_forearm = readtable('FP10_190411_114047_RightForearm_node_data.csv'); % data_hand = readtable('FP10_190411_114047_Right Hand_node_data.csv'); % data_forearm = readtable('FP11_190412_093259_RightForearm_node_data.csv'); % data_hand = readtable('FP11_190412_093259_Right Hand_node_data.csv'); % data_forearm = readtable('FP12_190412_100114_RightForearm_node_data.csv'); % data_hand = readtable('FP12_190412_100114_Right Hand_node_data.csv'); % data_forearm = readtable('FP13_190412_104307_RightForearm_node_data.csv'); % data_hand = readtable('FP13_190412_104307_Right Hand_node_data.csv');
% Data used for calculating the Markers angular velocity
25 %marker_data = readtable('FP12 40BPM.csv'); %marker_data = readtable('FP12 90BPM.csv'); %marker_data = readtable('FP12 140BPM.csv'); %marker_data = readtable('FP13 40BPM.csv'); %marker_data = readtable('FP13 90BPM.csv'); %marker_data = readtable('FP13 140BPM.csv'); % Convert CHAR to DOUBLE
data_forearm = str2double(table2cell(data_forearm(2:end,:))); data_hand = str2double(table2cell(data_hand(2:end,:)));
marker_data = str2double(table2cell(marker_data(7:end,:)));
% Remove NaN rows
data_forearm(any(isnan(data_forearm),2),:) = []; data_hand(any(isnan(data_hand),2),:) = [];
marker_data(any(isnan(marker_data),2),:) = [];
% Take out data from matrix to specific vectors
Ax_forearm = data_forearm(1:end,3); Ay_forearm = data_forearm(1:end,4); Az_forearm = data_forearm(1:end,5); Ax_hand = data_hand(1:end,3); Ay_hand = data_hand(1:end,4); Az_hand = data_hand(1:end,5);
disp(['(forearm) Choose a rough interval that includes two jumps']); jump_1_forearm = input('First jump: ');
jump_2_forearm = input('Second jump: ');
disp(['(hand) Choose a rough interval that includes two jumps']); jump_1_hand = input('First jump: ');
jump_2_hand = input('Second jump: ');
%% Update data with new interval
Gx_forearm = data_forearm(jump_1_forearm:jump_2_forearm,6); Gy_forearm = data_forearm(jump_1_forearm:jump_2_forearm,7); Gz_forearm = data_forearm(jump_1_forearm:jump_2_forearm,8); Gx_hand = data_hand(jump_1_hand:jump_2_hand,6); Gy_hand = data_hand(jump_1_hand:jump_2_hand,7); Gz_hand = data_hand(jump_1_hand:jump_2_hand,8); Ax_forearm = data_forearm(jump_1_forearm:jump_2_forearm,3); Ay_forearm = data_forearm(jump_1_forearm:jump_2_forearm,4); Az_forearm = data_forearm(jump_1_forearm:jump_2_forearm,5); Ax_hand = data_hand(jump_1_hand:jump_2_hand,3); Ay_hand = data_hand(jump_1_hand:jump_2_hand,4); Az_hand = data_hand(jump_1_hand:jump_2_hand,5); t_forearm = data_forearm(jump_1_forearm:jump_2_forearm,1); t_hand = data_hand(jump_1_hand:jump_2_hand,1); %% Interpolation
% Make a time vector with the frequency of 100 Hz
t2_forearm = t_forearm(1):1/100:t_forearm(end); t2_hand = t_hand(1):1/100:t_hand(end);
% forearm G
Gx_forearm_interp = interp1(t_forearm, Gx_forearm, t2_forearm); Gy_forearm_interp = interp1(t_forearm, Gy_forearm, t2_forearm); Gz_forearm_interp = interp1(t_forearm, Gz_forearm, t2_forearm);
% forearm A
26
% hand G
Gx_hand_interp = interp1(t_hand, Gx_hand, t2_hand); Gy_hand_interp = interp1(t_hand, Gy_hand, t2_hand); Gz_hand_interp = interp1(t_hand, Gz_hand, t2_hand);
% hand A
Ay_hand_interp = interp1(t_hand, Ay_hand, t2_hand);
%% Cross correlation sequence between sensor forearm and hand
[corr_coeff, lag] = xcorr(Ay_hand_interp, Ay_forearm_interp); [~, I] = max(corr_coeff);
shift = lag(I);
%% Align the data using the delay info if shift < 0 Ay_forearm_corr = Ay_forearm_interp(abs(shift):end); Ay_hand_corr = Ay_hand_interp; Gx_forearm_corr = Gx_forearm_interp(abs(shift):end); Gy_forearm_corr = Gy_forearm_interp(abs(shift):end); Gz_forearm_corr = Gz_forearm_interp(abs(shift):end); Gx_hand_corr = Gx_hand_interp; Gy_hand_corr = Gy_hand_interp; Gz_hand_corr = Gz_hand_interp; else Ay_hand_corr = Ay_hand_interp(shift:end); Ay_forearm_corr = Ay_forearm_interp; Gx_hand_corr = Gx_hand_interp(abs(shift):end); Gy_hand_corr = Gy_hand_interp(abs(shift):end); Gz_hand_corr = Gz_hand_interp(abs(shift):end); Gx_forearm_corr = Gx_forearm_interp; Gy_forearm_corr = Gy_forearm_interp; Gz_forearm_corr = Gz_forearm_interp; end
%% Calculating angular velocities (av) at a certain interval % Choosing interval
disp(['Choose an interval for which you want to calculate the angular velocity']) int_start = input('Start: ');
int_end = input('End: ');
% Calculating angular velocity
x = Gx_hand_corr(int_start:int_end) - Gx_forearm_corr(int_start:int_end); y = Gy_hand_corr(int_start:int_end) - Gy_forearm_corr(int_start:int_end); z = Gz_hand_corr(int_start:int_end) - Gz_forearm_corr(int_start:int_end); vel_tot_sensors = sqrt (x.^2 + y.^2 + z.^2);
%% calculating angular velocity from the markers using coordinates.
x1 = marker_data (1:end, 24); %% x1 is coordinates of the marker in the right wrist in x axis
x2 = marker_data (1:end, 27); %% x2 is coordinates of the marker in the left wrist in x axis
x3 = marker_data (1:end, 30); %% x3 is coordinates of the marker on the middle in x axis
x4 = marker_data (1:end, 21); %% x4 is coordinates of the marker on the elbow right wrist in x axis
27
y2 = marker_data (1:end, 28); %% y2 is coordinates of the marker in the left wrist in y axis
y3 = marker_data (1:end, 31); %% y3 is coordinates of the marker on the middle in y axis
y4 = marker_data(1:end,22); %% y4 is coordinates of the marker on the elbow right wrist in y axis
z1 = marker_data (1:end, 26); %% z1 is coordinates of the marker in the right wrist in z axis
z2 = marker_data(1:end, 29); %% z2 is coordinates of the marker in the left wrist in z axis
z3 = marker_data (1:end,32); %% z3 is coordinates of the marker on the middle in z axis
z4 = marker_data (1:end, 23); %% z4 is coordinates of the marker on the elbow right wrist in z axis
% Choosing interval
disp(['Choose an interval for which you want to calculate the angular velocity']) int_start = input('Start: ');
int_end = input('End: '); x1 = marker_data (int_start:int_end, 4); x2 = marker_data (int_start:int_end, 7); x3 = marker_data (int_start:int_end, 10); x4 = marker_data (int_start:int_end, 1); y1 = marker_data (int_start:int_end, 5); y2 = marker_data (int_start:int_end, 8); y3 = marker_data (int_start:int_end, 11); y4 = marker_data (int_start:int_end, 2); z1 = marker_data(int_start:int_end, 6); z2 = marker_data (int_start:int_end, 9); z3 = marker_data (int_start:int_end, 12); z4 = marker_data (int_start:int_end, 3);
xq = (x1+x2)./2; % xq ,yq and zq are coordinates of the middle points between two wrist markers
yq = (y1+y2)./2; zq = (z1+z2)./2;
% m , n and k are the sides of the triangle
m = sqrt((x4-x3).^2+(y4-y3).^2+(z4-z3).^2); n = sqrt((x3-xq).^2+(y3-yq).^2+(z3-zq).^2); k = sqrt((x4-xq).^2+(y4-yq).^2+(z4-zq).^2);
angle = acos ((k.^2+ n.^2-m.^2)./(2.*(k.*n))); % angle of flexion and deviation
B = 57.295779513 ; % conversion from radian to degree
anglenew = angle .*B ; % angle in degrees
ts=0.01;
for i=1: p % p is file size of the angle in degree
t(i)=ts; ts=ts + 0.01;
28
% velocity of the marker
vel_tot_markers=zeros(length(anglenew)-1,1); for i =1:length(anglenew)-1 vel_tot_markers(i)= (abs(anglenew(i+1)-anglenew(i)))./(t(i+1)-t(i)); end %% Filter [b,a] = butter(6,0.05); % 5/100 Hz vel_tot_sensors_filter = abs(filter(b,a,vel_tot_sensors)); vel_tot_markers_filter = abs(filter(b,a,vel_tot_markers));
% Take away "extreme" values
if length(vel_tot_markers_filter) < length(vel_tot_sensors_filter) stop = length(vel_tot_markers_filter); else stop = length(vel_tot_sensors_filter); end for i = 1:stop if vel_tot_markers_filter(i) > 550 vel_tot_markers_filter(i) = 550; elseif vel_tot_sensors_filter(i) > 550 vel_tot_sensors_filter(i) = 550; end end
%% Cross correlation between marker velocity and sensor velocity
[corr_coeff, lag] = xcorr(vel_tot_markers_filter, vel_tot_sensors_filter); [~, I] = max(corr_coeff);
shift = lag(I);
%% Align the data using the delay information if shift > 0 vel_tot_markers_filter_corr = vel_tot_markers_filter(1+abs(shift):end); vel_tot_sensors_filter_corr = vel_tot_sensors_filter; else vel_tot_sensors_filter_corr = vel_tot_sensors_filter(1+abs(shift):end); vel_tot_markers_filter_corr = vel_tot_markers_filter; end
%% Control alignment and visually
if length(vel_tot_markers_filter_corr) < length(vel_tot_sensors_filter_corr) int_end = length(vel_tot_markers_filter_corr); else int_end = length(vel_tot_sensors_filter_corr); end figure(7)
subplot(311), plot([1:int_end]/100, vel_tot_markers(1:int_end),'b'
,[1:int_end]/100, vel_tot_sensors(1:int_end), 'r'), title('Before filter
and correlation. Vel tot markers (blue) and vel tot sensors (red)');
subplot(312), plot([1:int_end]/100, vel_tot_markers_filter(1:int_end),'b' ,[1:int_end]/100, vel_tot_sensors_filter(1:int_end), 'r'), title('After
filter. Vel tot markers (blue) and vel tot sensors filter (red)');
29
subplot(313), plot([1:int_end]/100,
vel_tot_markers_filter_corr(1:int_end),'b' ,[1:int_end]/100,
vel_tot_sensors_filter_corr(1:int_end), 'r'), title('After correlation. Vel
tot markers corr (blue) and vel tot sensors filter corr (red)');
xlabel('Time [s]')
%% Statistical Analysis
% Calculate difference between markers velocity and sensors velocity
vel_diff = vel_tot_markers_filter_corr(1:int_end) - vel_tot_sensors_filter_corr(1:int_end)';
% percentile to table, Bland-Altman, and Correlation Plot
percentile_10_sensors = prctile(vel_tot_sensors_filter_corr(1:int_end),10); percentile_50_sensors = prctile(vel_tot_sensors_filter_corr(1:int_end),50)
% Correlation Plot, y-axis
percentile_90_sensors = prctile(vel_tot_sensors_filter_corr(1:int_end),90)
% Correlation Plot, y-axis
percentile_10_markers = prctile(vel_tot_markers_filter_corr(1:int_end),10); percentile_50_markers = prctile(vel_tot_markers_filter_corr(1:int_end),50)
% Correlation Plot & Bland Altman, x-axis
percentile_90_markers = prctile(vel_tot_markers_filter_corr(1:int_end),90)
% Correlation Plot & Bland Altman, x-axis
percentile_10_diff = percentile_10_sensors - percentile_10_markers
percentile_50_diff = percentile_50_sensors - percentile_50_markers % Bland
Altman y-axis
percentile_90_diff = percentile_90_sensors - percentile_90_markers % Bland
Altman y-axis % To Tables corrCoef = corr2(vel_tot_markers_filter_corr(1:int_end),vel_tot_sensors_filter_corr(1: int_end)') RMSD = sqrt(mean(vel_diff.^2))
MATLAB code for Bland -Altman Plot
%50th percentile value of the markers total velocity_40BPM
%50th percentile difference between the sensors_40BPM and markers_40BPM
%50th percentile of the markers total velocity_90BPM
%50th percentile difference between the sensors_90BPM and markers_90BPM
%50th percentile difference between the sensors_140BPM and markers_140BPM
%50th percentile of the markers total velocity_140BPM
function [dataX1, dataY1, dataX2, dataY2, dataX3, dataY3] =
BlandAltman_func_bpm(dataX1,dataY1,dataX2,dataY2,dataX3,dataY3)
%dataX1= 50th percentile value of the markers total velocity_40BPM
%dataY1 = 50th percentile difference between the sensors_40BPM and markers_40BPM
%dataX2 =50th percentile of the markers total velocity_90BPM
30
%dataX3 = 50th percentile difference between the sensors_140BPM and markers_140BPM
%dataY3 = 50th percentile of the markers total velocity_140BPM [m,n] = size(dataX1); [i,j] = size(dataX2); [k,l] = size(dataX3); if(n>m) dataY1 = dataY1'; end if (j>i) dataY2 = dataY2'; end if (l>k) dataY3 = dataY3'; end if(size(dataX1)~=size(dataY1))
error('Data matrices must be the same size')
end
if(size(dataX2)~=size(dataY2))
error('Data matrices must be the same size')
end
if(size(dataX3)~=size(dataY3))
error('Data matrices must be the same size')
end
Y_axis = [dataY1', dataY2', dataY3']; md_bpm = mean(Y_axis) % Mean of difference
between instruments
sd_bpm = std(Y_axis); % Std dev of difference
between instruments sd196_bpm = 1.96*sd_bpm figure(1) % Bland Altman plot
plot(dataX1,dataY1,'.','MarkerSize',20)
hold on,plot(dataX2,dataY2,'.','MarkerSize',20) hold on,plot(dataX3,dataY3,'.','MarkerSize',20)
plot(md_bpm*ones(300),'-k') % Mean difference line
plot(md_bpm+sd196_bpm*ones(300),'-k') % Mean plus 2*SD
line
plot(md_bpm-sd196_bpm*ones(300),'-k') % Mean minus
2*SD line % %% Correlation plots % % Correlation plot - 40BPM, 90BPM, 140BPM percentile_50_markers_tot_bpm = [percentile_50_markers_tot_40BPM, percentile_50_markers_tot_90BPM, percentile_50_markers_tot_140BPM]; percentile_50_sensors_tot_bpm = [percentile_50_sensors_tot_40BPM, percentile_50_sensors_tot_90BPM, percentile_50_sensors_tot_140BPM]; figure(3) plot(percentile_50_sensors_tot_bpm,percentile_50_markers_tot_bpm,'o') corrcoef = corr2(percentile_50_sensors_tot_bpm,percentile_50_markers_tot_bpm) title('Correlation between markers´ and sensors´ velocity at 40BPM, 90BPM and 140 BPM');
xlabel('Markers velocity'); ylabel('Sensors velocity');
str=[' r = ',num2str(corrcoef)]
T = text(min(get(gca, 'xlim')), max(get(gca, 'ylim')), str);
TRITA TRITA-CBH-GRU-2020:235
