Simple Tools for Motion and Gait Analysis

There is a great satisfaction in building good tools for other people to use.
                                                                                                                             Freeman Dyson
Calculating angles in 2D from marker coordinates 2Dangles.xls Zong-Ming Li, PhD
Musculoskeletal Research Center 
Department of Orthopaedic Surgery 
University of Pittsburgh 
See Discussion
The tan function is usually used to derive angle from 2D coordinate data. Unfortunately, it goes negative in the second quadrant (>90º). This can be avoided by using the atan2 function. However, this too goes negative in the fourth quadrant, so a test must be made to detect this and correct the result by adding 360º. So the complete algorithm is as follows:

ATAN2(x,y) returns an angle between -p and p
If angle < 0, angle = angle + 2p

A solution using atan is also given in case the atan2 function is not available.

Javascript digitizer Cycling example #1, #2, #3, #4, #5, #6 Chris Kirtley Video needs to be saved as separate JPG image files using, e.g. QuickTime export function, and numbered 01, 02, 03... etc. Wait for each image to be loaded.
Automatic marker tracking MaxTRAQ Innovision Systems Free 15 day trial
Finite Impulse Response (FIR) Critically-damped Low-pass Filter
4th order Butterworth
FIRfilter.xls Jim Martin
See Discussion
The Hamming coefficients are calculated first:

etc. where i = 1 to N, the number of taps; freq is the sampling

...then the tap coefficients themselves:
a(k)=H(k)*(0.54 + 0.46*cos(k*pi/5)) where k = 0 to N
Correction (for unity gain) = sum(a(0):a(5))+sum(a(1):a(5))

The filter equation is then:
y =(a(0)*x(t) + a(1)*x(t-1) + a(1)*x(t+1) + a(2)*x(t-2) + a(2)*x(t+2) +
a(3)*x(t-3) + a(3)*x(t+3) + a(4)*x(t-4) + a(4)*x(t+4) + a(5)*x(t-5) +

Infinite Impulse Response (IIR) Critically-damped Low-pass Filter 4th order (two-pass) Butterworth IIRfilter.xls Chris Kirtley (from David Winter's "Biomechanics and
Motor Control of Human Movement")
The equation is of the form:
yt = a.x + b.xt-1  + c.xt-2  +  +
where y is the filtered data point at the present time (time t ), and x the raw data at this time: t-1 means the last sample of this variable (e.g. ankle angle) etc. Notice that this is a recursive filter, because the output is fed back to the input. This can lead to problems (instability). The coefficients (a, b, c, d, e) can be calculated as follows:
   wc = tan(fc * p/N) where fc = cutoff frequency required and N is the sampling rate (e.g. frame rate of video for kinematics). Divide cutoff by 0.802 to take account of two passes (see below): e.g. for 6 Hz cutoff, set frequency to 6/0.802 = 7.48 Hz

   K1 = Ö2 * wc , where wc = 2pfc
   K2 = wc2
   a = K2 / (1 + K1 + K2)
   b = 2 * a
   c = a
   K3 = b / K2
   d = -2 * a + K3
   e = 1 - 2 * a - K3

After the first (forward) pass of the filter, phase distortion is introduced (as shown in the demo). This is removed by a second (backward) pass (note the improved correlation), of equation:

zt = ax + bxt+1  + cxt+2  + dyt+1  + eyt+2
where z is the final output. Notice that a couple of samples are needed to prime the filter in each direction - these can be set to the original (raw or first-pass) samples.
High-pass Filter
4th order (two-pass) Butterworth
HighPass.xls Peter Sinclair
School of Exercise and Sport Science
University of Sydney
A double pass filter like the butterworth filter provided by Chris
Kirtley.  Filter equations were published by Murphy and Robertson, JAB 10(4):  374-381, 1994. Like Chris's filter, there is potential for end-point distortion because of the lack of prior filter values. I have chosen to automatically use the raw values here, knowing that there can be some distortion and caution should be exercised near the end point.

For use of the filter, simply drop your data into column A and
enter appropriate frequencies in cells $A5 and $D$5. Formulae from columns B and C can be copied down as appropriate.

Resampling data to enable matching of two data sets recorded at different sampling rates Resample.xls Li Li, PhD (modified from a program
by D. Gordon Robertson, University of Ottawa) 
Department of Kinesiology
Louisiana State University
Linear interpolation using Excel Macro. Can be used for both down and up sampling. Since it uses a macro, you'll need to adjust your security settings to make it operational. Select Tools-Macros-Edit to see the script.
Fourier analysis to calculate Median Power Frequency MedianFreq.xls Peter Sinclair
School of Exercise and Sport Science
University of Sydney
Uses the Fourier Analysis tool within Excel to analyse a signal
and report Median Power Frequency. Fourier Analysis is an option from the Analysis Toolpack Add-In that comes with Excel. It only works with a binary number of data points (1024, 2048, 4096, etc). This particular sheet is designed for a sample frequency of 1000 Hz and analysis of 4096 points, but can be 
modified for other values.

Data is dropped into column A. Select Tools - Data Analysis - Fourier Analysis. Input range must be cells $A$2:$A$4097, output range must be cell $H$3.

Fourier Analysis outputs an array of complex numbers specifying
the amplitude for a range of frequencies. Column G is set up to give the frequency for each amplitude in column H. This must be modified if you use other than 1000 Hz collection frequency or 4096 data points. Column I gives a real number corresponding to the amplitude of each frequency. Column J 
calculates signal Power as the square of amplitude.

Median Frequency is calculated as the frequency corresponding to half the area of the Power - Frequency curve (Bartlett, 1997, p245). Column L calculates a cumulative sum of the powers and the Median Power Frequency is displayed in cell $0$6.

Residual Analysis to calculate suitable filter cutoff frequencies Resid.xls Peter Sinclair
School of Exercise and Sport Science
University of Sydney
This sheet contains a macro to repeatedly apply a 4th order
Butterworth filter to analyse the pattern of residual error vs smoothing frequency (Winter, 1990).

Data is dropped into column A and the formulae from columns B-D
filled down appropriately. Cell $H$7 contains the sum of squares of the residuals. In order to avoid end-point errors, the range of cells analysed in $H$7 should stay as far from the end points as possible.

Cells G10-G20 contain frequencies to be applied in the residual
analysis. These can be changed to whatever you want. Hitting the "Calculate" button runs the macro to smooth data using these frequencies.

The resulting graph has a line ready for you to move yourself by
freehand. Extreme caution should be used at this stage in choosing a suitable smoothing frequency. My experience is that the closer you look to find a suitable frequency, the more problems you find (the curve above the cut-off frequency will not be exactly linear). This method is rather approximate. If smoothing frequency is critical (eg if calculating acceleration) then more modern "Optimal" methods should be used (e.g. Giakas and Baltzopoulos).

Center of Mass from Force Platform Integration com.xls Chris Kirtley Integration of Ground Reaction Force over three successive steps, to derive acceleration, velocity and displacement of body CoM. CoM from total body kinematics shown for comparison. Note importance of accurate body mass and initial velocity estimates. In the vertical direction, the starting velocity can be found by integrating (using initial velocity = 0) to find the mean velocity, vo: - vo can then be used as the initial velocity to correct the integration. Note this only works for the vertical direction, and only if the mean vertical velocity is truly zero.
Savitsky-Golay Smoothing SG.xls D.Acosta, Univ. Florida Savitzky-Golay seeks to preserve shapes of peaks.
Statistics WinStat - free 30 day trial R. Fitch Broad range of statistical funtions, including ANOVA, cross-correlation, polynomial regressions, etc.

Please send your submissions or suggestions to Chris Kirtley

BackBack to Clinical Gait Analysis home page