- Email: [email protected]
Application of singular spectrum analysis to the smoothing of raw kinematic signals
ARTICLE IN PRESS
Journal of Biomechanics 38 (2005) 1085–1092
Application of singular spectrum analysis to the smoothing of raw kinematic signals F.J. Alonsoa,*, J.M. Del Castilloa, P. Pintadob a
Department of Electronics and Electromechanical Engineering, University of Extremadura, Avda. de Elvas s/n, 06071, Badajoz, Spain b Department of Applied Mechanics, University of Castilla La Mancha, Avda. Camilo Jose´ Cela s/n, 13071, Ciudad Real, Spain Accepted 24 May 2004
Abstract Motion capture systems currently used in biomechanical analysis introduce systematic measurement errors that appear in the form of noise in recorded displacement signals. The noise is unacceptably amplified when differentiating displacements to obtain velocities and accelerations. To avoid this phenomenon, it is necessary to smooth the displacement signal prior to differentiation in order to eliminate the noise introduced by the experimental system. The use of singular spectrum analysis (SSA) is presented in this paper as an alternative to traditional digital filtering methods. SSA is a novel non-parametric technique based on principles of multivariate statistics. The original time series is decomposed into a number of additive time series, each of which can be easily identified as being part of the modulated signal, or as being part of the random noise. Several examples that demonstrate the superiority of this technique over other methods used in biomechanical analysis are presented in this paper. r 2004 Elsevier Ltd. All rights reserved. Keywords: Singular spectrum analysis; Signal differentiation; Smoothing
1. Introduction Biomechanical analysis is prone to numerous sources of error that affect its results and reduce its usefulness. The most important sources of error have recently been pointed out by Hatze (2002). Among them is the inaccuracy of velocities and accelerations derived from experimentally measured displacements. The high-frequency noise in the displacement signal is dramatically amplified when obtaining velocities and accelerations. It is absolutely necessary to filter or to smooth the displacement signal prior to differentiation in order to increase the accuracy of inverse dynamic analysis (van den Bogert and de Koning, 1996). The problem of filtering displacement time series to obtain noiseless velocities and accelerations with minimum loss of information has been widely studied. Digital Butterworth filters, splines, and filters based on spectral analysis are among traditional smoothing methods (D’Amico and Ferrigno, 1990,1992; Giakas *Corresponding author. Tel.: +34-924289600; fax: +34-924289601. E-mail address: [email protected] (F.J. Alonso). 0021-9290/$ - see front matter r 2004 Elsevier Ltd. All rights reserved. doi:10.1016/j.jbiomech.2004.05.031
and Baltzopoulos, 1997a, b; Vaughan, 1982). Comparison studies have been carried out to test the capabilities of these methods (Giakas and Baltzopoulos, 1997b; Walker, 1998), and results show that spline-based methods are usually more stable for different signal-tonoise ratios. Nonetheless, traditional filtering methods are not suited for smoothing non-stationary signals. This drawback is particularly problematic in biomechanical analysis since physical activity involves impact-like floor reaction forces (Woltring, 1995; Giakas et al., 2000; Georgiakis et al., 2002). Discrete wavelet transforms (Adham and Shihab, 1999) and the Wigner function (Giakas et al., 2000) are used in advanced filtering methods that may be used for non-stationary signals. Nonetheless, the drawback in these cases is the complexity of devising an automatic and systematic procedure. A mother wavelet function must be selected when using discrete wavelet transforms, and the filtering function parameters must be chosen when using the Wigner function. The goal of this paper is to demonstrate the advantages of smoothing methods based on singular spectrum analysis (SSA) techniques. One of these
ARTICLE IN PRESS 1086
F.J. Alonso et al. / Journal of Biomechanics 38 (2005) 1085–1092
advantages is the ease with which an automatic and systematic procedure may be devised. The method will be applied to four different signals: two of them taken from the literature, and two others acquired with our own motion capture system.
2. Methods 2.1. Background: SSA SSA is a novel non-parametric technique used in the analysis of time series and based on principles of multivariate statistics. Its usefulness has been proven in the analysis of climatic, meteorological and geophysical time series (Golyandina et al., 2001). In fact, SSA was first applied to extract tendencies and harmonic components in meteorological and geophysical time series (Vautard and Ghil, 1989; Vautard et al., 1992), as well as to identify periodic motion in complex dynamic systems (Pasmanter, 1995). A concise description of the method will be given in this section, whereas Golyandina et al. (2001) have presented a complete derivation. The method starts by producing a Hankel matrix from the time series itself by sliding a window that is shorter in length than the original series. This step is referred to as ‘‘embedding’’. The columns of the matrix correspond to the terms inside the window for every position of the said window. This matrix is then decomposed into a number of elementary matrices of decreasing norm. This step is called singular value decomposition (SVD). Truncating the summation of elementary matrices yields an approximation of the original matrix. The approximation is one in which those elementary matrices that hardly contribute to the norm of the original matrix have been neglected. This step is called ‘‘grouping’’. Thus, the result is no longer a Hankel matrix, but an approximated time series may be recovered by taking the average of the diagonals. This new signal is the smoothed approximation of the original. This step is the ‘‘reconstruction’’ or the ‘‘diagonal averaging’’. The above description may be stated in formal terms as follows: Step 1: Embedding Let F ¼ ðf0 ; f1 ; y; fN1 Þ be the length N time series representing the noisy signal. Let L be the window length, with 1oLoN and L an integer. Each column Xj of the Hankel matrix corresponds to the ‘‘snapshot’’ taken by the sliding window: X j ¼ ðfj1 ; fj ; y; fjþL2 ÞT ; j ¼ 1; 2; y; K; where K=NL+1 is the number of columns, that is the number of different possible positions of said window. The matrix X=[X1, X2,y,XK] is a Hankel matrix since all elements in diagonal i+j=constant are equal. This matrix is sometimes referred to as the trajectory matrix.
Step 2: SVD of the trajectory matrix It can be proven that the trajectory matrix (or any matrix, for that matter) may be expressed as the summation of d rank one elementary matrices X ¼ E 1 þ ? þ E d ; where d is the number of non-zero eigenvalues of the L L matrix S ¼pXffiffiffiffi X T : The elementary matrices are given by E i ¼ li U i V Ti ði ¼ 1; y; dÞ; where l1 ; y; ld are the non-zero eigenvalues of S, in decreasing order, U 1 ; y; U d are the corresponding eigenvectors, pffiffiffiffi and vectors Vi are obtained from V i ¼ X T U i = li ; i ¼ 1; y; d: pffiffiffiffi The norm of elementary matrix Ei equals li : Therefore, the contribution of the first matrices to the norm of X is much higher than the contribution of the last matrices, and it is likely that these last matrices represent noise in the signal. The plot of the eigenvalues in decreasing order is called the singular spectrum and is essential in deciding the index from where to truncate the summation. Step 3: Grouping This step is very simple when the method is used for smoothing a time series. It consists of approximating matrix X by the summation of the first r elementary matrices. More complex types of grouping may be necessary when one intends to extract tendencies from the time series behaviour. In our case, matrix X is approximated by XEE 1 þ E 2 þ ? þ E r : Step 4: Reconstruction (diagonal averaging) The approximated matrix described above is no longer a Hankel matrix, but an approximated time series may be recovered by taking the average of the diagonals. Nevertheless, it may be more practical to carry out this averaging for each elementary matrix independently in order to obtain time series that represent the different components of the behaviour of the original time series. These ‘‘elementary’’ time series are referred to as ‘‘principal components’’. Let Y be any of the elementary matrices Ei, the elements of which are yij; 1pipL; 1pjpK: The time series g0,y,gN1 (principal component) corresponding to this elementary matrix is given by 8 P 1 kþ1 > > ym;kmþ2 for 0pkoL 1; > > > k þ 1 m¼1 > > > > < L 1 P ym;kmþ2 for L 1pkoK ; gk ¼ L > m¼1 > > > > > NK P þ1 1 > > > ym;kmþ2 for K pkoN; :N k m¼kK þ2 where L=min(L,K), and K=max(L,K). The smoothed time series is obtained by adding the first r principal components. It is worth pointing out that the application of the basic SSA algorithm requires selecting the values of just two parameters: the window length L, and the number r
ARTICLE IN PRESS F.J. Alonso et al. / Journal of Biomechanics 38 (2005) 1085–1092
of principal components to be retained in the reconstruction.
1087
150 100
There are no general rules for selecting the values of parameters L and r that participate in the SSA algorithm. Their values depend on the type of signal to be analysed, and on the type of information to be extracted from the signal: tendencies, harmonics, noise. A few general rules will be given in this section, whereas, again, Golyandina et al. (2001) have presented a complete description of parameter selection. The main goal is to select a window length that produces separable and independent principal components. The number of these components to be retained is then easily determined, not only from the appearance of the principal component itself, but also from the relative value of the corresponding eigenvalue. A few remarks about the window length will help in determining its value: (1) The SVD performed on matrices obtained with a window length L is equivalent to that performed on matrices obtained with the complementary window length K=NL+1. Therefore, increasing the window length above half the time series length would reproduce results already tested with shorter window lengths. (2) The longer the window, the more detailed the decomposition. This statement, along with the previous one, indicates that the best detailed decomposition is obtained with LEN/2. (3) Despite the two previous statements, a large window length may mix noise and certain components in the case of time series with a complex structure. In this case (an example of which will be presented in Section 3), it is worth testing different window lengths, starting with LEN/2, and proceeding downwards until the appropriate separation is achieved. 2.3. Test data Four different signals will be smoothed using the method described in the previous section. A double differentiation is performed on the signals in order to quantify the effect of smoothing. First-order finite differences are used to calculate the higher derivatives (Burden et al., 1981). Signal 1: Stationary sinusoidal signal A reflective marker was attached to a uniformly rotating crank. The marker motion was captured for 8 s at a frequency of 100 Hz using three infrared cameras (Qualisys Medical AB). The vertical component of this motion (Fig. 1) was smoothed with an SSA algorithm implemented in Matlab (The Mathworks Inc.). Signal 2: Biomechanical signal (Vaughan et al., 1992; Giakas and Baltzopoulos, 1997a)
Displacement (mm)
2.2. Selection of SSA parameters 50 0 -50 -100 -150 0
1
2
3
4 Time (s)
5
6
7
8
Fig. 1. Noisy displacement (Signal 1) recorded using a motion capture system.
The second signal to be processed has been taken from the literature. It corresponds to the ‘woman’ file from GAITLAB (Vaughan et al., 1992). The signal measures the lateral displacement of a marker attached on the right tibial tubercle. Data were acquired for 0.94 s with a frequency of 50 Hz. The resulting time series length is 48 elements. The reference signal is taken as that obtained after filtering the raw data at 6.25 Hz, but different amplitude white noises are added in order to test smoothing methods. An amplitude 1, time generated, white noise was used for comparisons (Giakas and Baltzopoulos, 1997a). Signal 3: Signal with impact (Dowling, 1985; Giakas et al., 2000) This third signal has also been taken from the literature. The signal is a measure of the angular coordinate of a pendulum impacting against a compliant wall (Dowling, 1985; Giakas et al., 2000). The angular acceleration obtained from the motion capture system is compared to that obtained directly (after dividing by pendulum length) from accelerometers. Three accelerometers were used in order to average their measurements to reduce noise. The average signal, logged at a sampling rate of 512 Hz is used as the acceleration reference signal. Signal 4: Slider motion A vertical slider was moved by hand to obtain a nonstationary mono-dimensional motion of biomechanical origin. A subject was asked to move the slider randomly with fast upward and downward movements. The vertical position (Fig. 2) was acquired with the use of a marker attached to the slider and the Qualisys camera system (Qualisys Medical AB). The second derivative of this record, or of the record obtained after smoothing, is
ARTICLE IN PRESS 1088
F.J. Alonso et al. / Journal of Biomechanics 38 (2005) 1085–1092 109 108
Eigenvalue
107 106 105 104 103 102
Fig. 2. Original noisy displacement (Signal 4) recorded using a motion capture system.
101
compared to the acceleration obtained directly from an accelerometer attached to the slider. Displacement and acceleration were sampled at 200 Hz during 5.90 s, obtaining records of 1182 elements each.
0
5
10
15 20 25 Eigenvalue number
30
35
40
Fig. 3. Singular spectrum (in logarithmic scale) of Signal 1, obtained with a window length L = 40.
100 0
3. Results
-100 0
1
2
3
4
5
6
7
8
-100 0 0.5
1
2
3
4
5
6
7
8
1
2
3
4
5
6
7
8
1
2
3
4
5
6
7
8
1
2
3
4 Time (sec)
5
6
7
8
100
3.1. Signal 1: Stationary sinusoidal signal Displacement (mm)
The first time series has a total of 800 elements, and the trajectory matrix was generated by sliding a window of 40 elements in length. The plot of the eigenvalues of the L L matrix S ¼ X X T (X being the trajectory matrix) is given in Fig. 3. It has been mentioned that the norm of the elementary matrices into which matrix X may be decomposed equals the square root of the eigenvalues of S. The eigenvalues in Fig. 3 are plotted in logarithmic scale. Therefore, the contribution of the first few eigenvalues to the norm of X is much higher than the contribution of the rest. In fact, 99.9992% of the norm of X is contributed by the first two eigenvalues, whereas the rest are much lower and vary slightly. This fact usually indicates that these last eigenvalues represent noise, or at least a nonperiodic component without any latent structure (Golyandina et al., 2001). It has been pointed out that the signals reconstructed using one elementary matrix at a time are called principal components. The first five principal components for the case under analysis are plotted in Fig. 4. It is clear from the figure that the two leading components represent the main trend (note that the vertical axis is scaled differently for each component), whereas the next three components represent noise. The approximated reconstruction will then be carried out using the first two principal components. The reconstructed signal is plotted in Fig. 5a along with the original signal. No differences can be made out at the scale of the plot. Slight divergence is revealed
0
0 -0.5 0 0.1 0 -0.1 0 0.1 0 -0.1 0
Fig. 4. Individual reconstruction of the five leading components (plotted in different scales) obtained from the SVD of Signal 1, using a window length L = 40.
when the difference between the two signals (Fig. 5b), or the signals themselves (Fig. 5c), is plotted on a much greater scale. The signal-to-noise ratio is very large in this example. Nonetheless, this small noise heavily contaminates the numerically obtained acceleration signal, and it is here, in the acceleration, where the efficiency of the method can be tested. The acceleration obtained using the original raw signal is plotted in Fig. 6 along with the acceleration obtained using the approximated reconstructed signal. It is clear that smoothing or filtering is absolutely mandatory in order to obtain meaningful acceleration signals. Moreover, SSA smoothing is superior to Butterworth filtering. This statement is supported in the example
ARTICLE IN PRESS F.J. Alonso et al. / Journal of Biomechanics 38 (2005) 1085–1092
1089
150
50 0 -50 -100 -150
0
2
4 Time (s)
Displacement (mm)
(a)
(b)
6
8
Displacement (mm)
Displacement (mm)
100
2
(c)
120 119 118 117 1.72 1.74 1.76 1.78 Time (s)
1.8
1.82
0
-2
0
1
2
3
4
5
6
7
8
Time (s)
Fig. 5. (a) Noisy displacement (dotted line) and SSA-smoothed displacement (continuous line) corresponding to Signal 1. (b) Residual signal defined as the difference between the original and smoothed signals. (c) Zoom in the neighbourhood of t = 1.77 s to appreciate the slight difference between the original and SSA-smoothed signals.
Fig. 6. (a) Acceleration calculated from Signal 1 original displacement data. (b) Acceleration calculated from Signal 1 after having been smoothed using SSA. (c) Acceleration calculated from Signal 1 after having been passed through a 4 Hz cut frequency Butterworth filter.
also shown in Fig. 6. The original signal was passed through a second-order Butterworth filter with a 4 Hz cut frequency, and the result was differentiated twice to obtain acceleration. The superiority is particularly relevant in the elimination of the so-called end-point errors. Methods based on signal extension have been proposed to reduce end-point error (Giakas and Batzopoulos, 1998). No extension is necessary in the
case of SSA smoothing. Errors due to this phenomenon are negligible (Fig. 6). 3.2. Signal 2: Biomechanical signal (Vaughan et al., 1992; Giakas and Baltzopoulos, 1997a) Decomposition was carried out on a trajectory matrix obtained with a window length L=N/2=24. The first 9
ARTICLE IN PRESS F.J. Alonso et al. / Journal of Biomechanics 38 (2005) 1085–1092
1090
200
4000 3000
100
1000
Acceleration (rad/s2)
Acceleration (mm/s2)
2000
0 -1000 -2000 -3000 -4000
0
-100
-200
-300
-5000
-400
-6000 0
0.1
0.2
0.3
0.4 0.5 Time (s)
0.6
0.7
0.8
0.9
0
0.2
0.4
0.6 Time (s)
0.8
1
1.2
Fig. 7. Acceleration obtained from reference signal 2 (continuous line) and acceleration calculated from SSA-smoothed displacement data (dotted line).
Fig. 8. Reference angular acceleration corresponding to Signal 3 (continuous line) and acceleration calculated from angular displacement data smoothed with a sequential SSA decomposition (dotted line).
principal components were used for reconstruction. The tenth component is a low-amplitude high-frequency motion which is most likely due to noise. Comparison with traditional filtering techniques is quantified by taking the root mean square of the error signal (RMSE) of displacement, velocity and acceleration. The error signal is the difference between the reference signal (or its derivatives) and the signal obtained after smoothing or filtering the raw data (or its derivatives). Our method achieved RMSE=0.11 mm for displacements, versus RMSE=0.4 mm obtained by Giakas and Baltzopoulos (1997a) using power spectrum assessment (PSA). Similar improvements are achieved in velocity: RMSE=4.35 mm/s versus RMSE=10 mm/s; and in acceleration: RMSE=256.89 mm/s2 versus RMSE=400 mm/s2. Accuracy in acceleration may also be appreciated in Fig. 7, where the acceleration obtained from the reference signal is plotted along with the acceleration calculated from SSA-smoothed displacement data.
form a sequential SSA decomposition (Golyandina et al., 2001), which consists in applying SSA smoothing to an already smoothed record. A second SSA decomposition was applied to the previously smoothed signal. Different window lengths and grouping strategies were tested. The best results were obtained using the first two new components for reconstruction with L2 = 20. A value of RMSE = 23.04 rad/s2 (see Fig. 8) was reached in this case. The result is similar to the value RMSE = 23.60 rad/s2 obtained by Giakas et al. (2000) with the help of the Wigner distribution. Our reconstruction right at the instant of impact is a little worse, but our signal in the vicinity of this time instant is more accurate. The slight loss of accuracy in the SSA method is compensated by the ease with which the method is applied: fewer parameters to be chosen, and no need to extend the ends of the record in order to eliminate end-point errors.
3.3. Signal 3: Signal with impact (Dowling, 1985; Giakas et al., 2000)
3.4. Signal 4: Slider motion
Giakas et al. (2000) have recently shown that secondorder Butterworth filters do not work properly in cases like this due to the non-stationary nature of the signal. Thus, this is a good example for testing SSA smoothing versus the use of the Wigner distribution. The record of 600 elements from motion capture was processed using a window length L=100. A value of RMSE = 64.98 rad/s2 is obtained when reconstructing with the first nine components. This high error is due to the fact that part of the noise has been included in the reconstruction. To minimize this effect, one may per-
Decomposition was attained with a window length L = 100. Two different reconstructed signals will be used for comparisons, one obtained from the first 16 principal components, and another obtained from the first five components alone. These signals were numerically differentiated twice to produce acceleration records to be compared to the accelerometer output. These comparisons are made in Fig. 9, where the acceleration obtained from the raw data, from the two reconstructed signals, and directly from the accelerometer, are shown. The errors in terms of RMSE are 2.04 and 1.47 m/s2, respectively.
ARTICLE IN PRESS F.J. Alonso et al. / Journal of Biomechanics 38 (2005) 1085–1092
1091
Fig. 9. (a) Acceleration calculated from Signal 4 original displacement data (dashed line) and acceleration measured by the accelerometer (continuous line). (b) Acceleration calculated from displacement reconstructed using the sixteen leading components (dotted line) and acceleration measured by the accelerometer (continuous line). (c) Acceleration calculated from displacement reconstructed using the five leading components (dotted line) and acceleration measured by the accelerometer (continuous line).
4. Discussion The results of our study show the superiority of SSA smoothing over filtering techniques found in the literature. The SSA algorithm decomposes the original signal into independent additive components of decreasing weight. This fact is key in the success of the method when it comes to extracting the latent trend in the signal from the random noise inherent to the motion capture system. One of the main advantages to the method is the fact that the algorithm requires selection of just two parameters, namely, the window length and the number of components to use for reconstruction. The method is also very intuitive in the sense that the appearance of the different components, along with the corresponding eigenvalue, is a very good indicator for distinguishing signal from noise. It has also been shown that the method works properly on both stationary and nonstationary signals. Moreover, the method can be very easily programmed as a stand-alone automatic algorithm. Complex records may need sequential decomposition to better separate noise from signal. Nevertheless, this sequential process can also be made automatic with little effort. Convergence may be measured by means of the difference between the current and previous values of the acceleration RMS. The algorithm stops when this difference is sufficiently small.
As drawbacks, one may mention the fact that there are no fixed, objective rules for selecting the window length. Nevertheless, it has been shown in the examples that the results are not very sensitive to window length. Moreover, the grouping strategy is usually clearly indicated by the singular spectrum. However, in some cases, as in signal 3, a choice of a large window length L produces a poor separation between signal trend and noise. On the other hand, for small L, only the overall trend would be extracted, relevant information being lost in the original signal (Golyandina et al., 2001). One way to overcome the problem is to apply sequential SSA. Some components of the initial series are extracted using standard SSA, and then the components of interest are extracted by applying SSA once more. In order to automate such a procedure, one could check the power spectral density (PSD) of the reconstructed signal after each decomposition to detect the presence of noise and to decide whether or not additional decompositions are necessary. Future studies will focus on automating the process. Techniques will be devised to automatically choose the window length and the grouping strategy, and a procedure will be implemented to automatically make the decision on whether or not sequential SSA is recommended for the given signal. In conclusion, we believe that the biomechanics community will benefit from this new smoothing
ARTICLE IN PRESS 1092
F.J. Alonso et al. / Journal of Biomechanics 38 (2005) 1085–1092
technique that has proven its effectiveness with complex signals. Future studies will need to focus on the possibilities of producing an automatic algorithm, and on the possibilities of embedding the algorithm in commercial biomechanical analysis packages. The application of the method to signals with rapid abrupt jumps, such as those used by Adham and Shihab (1999), will also need to be studied in greater detail. It will probably be necessary to devise an extension of SSA to detect structural changes in the time series, and use different parameters (window length and grouping strategy) before, during, and after the structural change.
References Adham, R.I., Shihab, S.A., 1999. Discrete wavelet transform: a tool in smoothing kinematic data. Journal of Biomechanics 32, 317–321. Burden, R.L., Faires, J.D., Reynolds, A.C., 1981. Numerical Analysis. PWS Publishers, Boston. D’Amico, M., Ferrigno, G., 1990. Technique for the evaluation of derivatives from noisy biomechanical displacement data using a model-based bandwith selection procedure. Medical and Biological Engineering and Computing 28, 407–415. D’Amico, M., Ferrigno, G., 1992. Comparison between the more recent techniques for smoothing and derivative assessment in biomechanics. Medical and Biological Engineering and Computing 30, 193–204. Dowling, J., 1985. A modelling strategy for the smoothing of biomechanical data. In: Johnsson, B. (Ed.), Biomechanics, Vol. XB. Human Kinetics, Champaign, IL, pp. 1163–1167. Georgiakis, A., Stergioulas, L.K., Giakas, G., 2002. Wigner filtering with smooth roll-off boundary for differentiation of noisy nonstationary signals. Signal Processing 82, 1411–1415. Giakas, G., Baltzopoulos, V., 1997. A comparison of automatic filtering techniques applied to biomechanical walking data. Journal of Biomechanics 30, 847–850.
Giakas, G., Baltzopoulos, V., 1997a. Optimal digital filtering requires a different cut-off frequency strategy for the determination of the higher derivatives. Journal of Biomechanics 30, 851–855. Giakas, G., Baltzopoulos, V., 1998b. Improved extrapolation techniques in recursive digital filtering: a comparison of least squares and prediction. Journal of Biomechanics 31, 87–91. Giakas, G., Stergioulas, L.K., Vourdas, A., 2000. Time-frequency analysis and filtering of kinematic signals with impacts using the Wigner function: accurate estimation of the second derivative. Journal of Biomechanics 33, 567–574. Golyandina, N., Nekrutkin, V., Zhigljavsky, A., 2001. Analysis of Time Series Structure—SSA and Related Techniques. Chapman & Hall/CRC, Boca Raton, FL, pp. 13-78. Hatze, H., 2002. The fundamental problem of myoskeletal inverse dynamics and its implications. Journal of Biomechanics 35, 109–115. Pasmanter, R.A., 1995. Variational search of periodic motions in complex dynamical systems. Chaos, Solitons and Fractals 6, 447–454. van den Bogert, A.J., de Koning, J.J., 1996. On optimal filtering for inverse dynamics analysis. Proceedings of the Nineth Canadian Society of Biomechanics Congress, Burnaby, BC, pp. 214–215. Vaughan, C.L., 1982. Smoothing and differentiation of displacementtime data: an application of splines and digital filtering. International Journal of Bio-Medical Computing 13, 375–386. Vaughan, C., Davis, B., O’Connor, J., 1992. Gait Analysis Laboratory. Human Kinetics, Champaign, IL. Vautard, R., Ghil, M., 1989. Singular spectrum analysis in nonlinear dynamics, with applications to paleoclimatic time series. Physica D: Nonlinear Phenomena 35, 395–424. Vautard, R., Yiou, P., Ghil, M., 1992. Singular-spectrum analysis: a toolkit for short, noisy chaotic signals. Physica D: Nonlinear Phenomena 58, 95–126. Walker, J.A., 1998. Estimating velocities and accelerations of animal locomotion: a simulation experiment comparing numerical differentiation algorithms. Journal of Experimental Biology 201, 981–995. Woltring, H.J., 1995. Smoothing and differentiation techniques applied to 3-D data. In: Allard, P., Stokes, I.A.F., Blanchi, J.P. (Eds.), Three Dimensional Analysis of Human Movement. Human Kinetics. Champaign, IL, pp. 79–99.
Journal of Biomechanics 38 (2005) 1085–1092
Application of singular spectrum analysis to the smoothing of raw kinematic signals F.J. Alonsoa,*, J.M. Del Castilloa, P. Pintadob a
Department of Electronics and Electromechanical Engineering, University of Extremadura, Avda. de Elvas s/n, 06071, Badajoz, Spain b Department of Applied Mechanics, University of Castilla La Mancha, Avda. Camilo Jose´ Cela s/n, 13071, Ciudad Real, Spain Accepted 24 May 2004
Abstract Motion capture systems currently used in biomechanical analysis introduce systematic measurement errors that appear in the form of noise in recorded displacement signals. The noise is unacceptably amplified when differentiating displacements to obtain velocities and accelerations. To avoid this phenomenon, it is necessary to smooth the displacement signal prior to differentiation in order to eliminate the noise introduced by the experimental system. The use of singular spectrum analysis (SSA) is presented in this paper as an alternative to traditional digital filtering methods. SSA is a novel non-parametric technique based on principles of multivariate statistics. The original time series is decomposed into a number of additive time series, each of which can be easily identified as being part of the modulated signal, or as being part of the random noise. Several examples that demonstrate the superiority of this technique over other methods used in biomechanical analysis are presented in this paper. r 2004 Elsevier Ltd. All rights reserved. Keywords: Singular spectrum analysis; Signal differentiation; Smoothing
1. Introduction Biomechanical analysis is prone to numerous sources of error that affect its results and reduce its usefulness. The most important sources of error have recently been pointed out by Hatze (2002). Among them is the inaccuracy of velocities and accelerations derived from experimentally measured displacements. The high-frequency noise in the displacement signal is dramatically amplified when obtaining velocities and accelerations. It is absolutely necessary to filter or to smooth the displacement signal prior to differentiation in order to increase the accuracy of inverse dynamic analysis (van den Bogert and de Koning, 1996). The problem of filtering displacement time series to obtain noiseless velocities and accelerations with minimum loss of information has been widely studied. Digital Butterworth filters, splines, and filters based on spectral analysis are among traditional smoothing methods (D’Amico and Ferrigno, 1990,1992; Giakas *Corresponding author. Tel.: +34-924289600; fax: +34-924289601. E-mail address: [email protected] (F.J. Alonso). 0021-9290/$ - see front matter r 2004 Elsevier Ltd. All rights reserved. doi:10.1016/j.jbiomech.2004.05.031
and Baltzopoulos, 1997a, b; Vaughan, 1982). Comparison studies have been carried out to test the capabilities of these methods (Giakas and Baltzopoulos, 1997b; Walker, 1998), and results show that spline-based methods are usually more stable for different signal-tonoise ratios. Nonetheless, traditional filtering methods are not suited for smoothing non-stationary signals. This drawback is particularly problematic in biomechanical analysis since physical activity involves impact-like floor reaction forces (Woltring, 1995; Giakas et al., 2000; Georgiakis et al., 2002). Discrete wavelet transforms (Adham and Shihab, 1999) and the Wigner function (Giakas et al., 2000) are used in advanced filtering methods that may be used for non-stationary signals. Nonetheless, the drawback in these cases is the complexity of devising an automatic and systematic procedure. A mother wavelet function must be selected when using discrete wavelet transforms, and the filtering function parameters must be chosen when using the Wigner function. The goal of this paper is to demonstrate the advantages of smoothing methods based on singular spectrum analysis (SSA) techniques. One of these
ARTICLE IN PRESS 1086
F.J. Alonso et al. / Journal of Biomechanics 38 (2005) 1085–1092
advantages is the ease with which an automatic and systematic procedure may be devised. The method will be applied to four different signals: two of them taken from the literature, and two others acquired with our own motion capture system.
2. Methods 2.1. Background: SSA SSA is a novel non-parametric technique used in the analysis of time series and based on principles of multivariate statistics. Its usefulness has been proven in the analysis of climatic, meteorological and geophysical time series (Golyandina et al., 2001). In fact, SSA was first applied to extract tendencies and harmonic components in meteorological and geophysical time series (Vautard and Ghil, 1989; Vautard et al., 1992), as well as to identify periodic motion in complex dynamic systems (Pasmanter, 1995). A concise description of the method will be given in this section, whereas Golyandina et al. (2001) have presented a complete derivation. The method starts by producing a Hankel matrix from the time series itself by sliding a window that is shorter in length than the original series. This step is referred to as ‘‘embedding’’. The columns of the matrix correspond to the terms inside the window for every position of the said window. This matrix is then decomposed into a number of elementary matrices of decreasing norm. This step is called singular value decomposition (SVD). Truncating the summation of elementary matrices yields an approximation of the original matrix. The approximation is one in which those elementary matrices that hardly contribute to the norm of the original matrix have been neglected. This step is called ‘‘grouping’’. Thus, the result is no longer a Hankel matrix, but an approximated time series may be recovered by taking the average of the diagonals. This new signal is the smoothed approximation of the original. This step is the ‘‘reconstruction’’ or the ‘‘diagonal averaging’’. The above description may be stated in formal terms as follows: Step 1: Embedding Let F ¼ ðf0 ; f1 ; y; fN1 Þ be the length N time series representing the noisy signal. Let L be the window length, with 1oLoN and L an integer. Each column Xj of the Hankel matrix corresponds to the ‘‘snapshot’’ taken by the sliding window: X j ¼ ðfj1 ; fj ; y; fjþL2 ÞT ; j ¼ 1; 2; y; K; where K=NL+1 is the number of columns, that is the number of different possible positions of said window. The matrix X=[X1, X2,y,XK] is a Hankel matrix since all elements in diagonal i+j=constant are equal. This matrix is sometimes referred to as the trajectory matrix.
Step 2: SVD of the trajectory matrix It can be proven that the trajectory matrix (or any matrix, for that matter) may be expressed as the summation of d rank one elementary matrices X ¼ E 1 þ ? þ E d ; where d is the number of non-zero eigenvalues of the L L matrix S ¼pXffiffiffiffi X T : The elementary matrices are given by E i ¼ li U i V Ti ði ¼ 1; y; dÞ; where l1 ; y; ld are the non-zero eigenvalues of S, in decreasing order, U 1 ; y; U d are the corresponding eigenvectors, pffiffiffiffi and vectors Vi are obtained from V i ¼ X T U i = li ; i ¼ 1; y; d: pffiffiffiffi The norm of elementary matrix Ei equals li : Therefore, the contribution of the first matrices to the norm of X is much higher than the contribution of the last matrices, and it is likely that these last matrices represent noise in the signal. The plot of the eigenvalues in decreasing order is called the singular spectrum and is essential in deciding the index from where to truncate the summation. Step 3: Grouping This step is very simple when the method is used for smoothing a time series. It consists of approximating matrix X by the summation of the first r elementary matrices. More complex types of grouping may be necessary when one intends to extract tendencies from the time series behaviour. In our case, matrix X is approximated by XEE 1 þ E 2 þ ? þ E r : Step 4: Reconstruction (diagonal averaging) The approximated matrix described above is no longer a Hankel matrix, but an approximated time series may be recovered by taking the average of the diagonals. Nevertheless, it may be more practical to carry out this averaging for each elementary matrix independently in order to obtain time series that represent the different components of the behaviour of the original time series. These ‘‘elementary’’ time series are referred to as ‘‘principal components’’. Let Y be any of the elementary matrices Ei, the elements of which are yij; 1pipL; 1pjpK: The time series g0,y,gN1 (principal component) corresponding to this elementary matrix is given by 8 P 1 kþ1 > > ym;kmþ2 for 0pkoL 1; > > > k þ 1 m¼1 > > > > < L 1 P ym;kmþ2 for L 1pkoK ; gk ¼ L > m¼1 > > > > > NK P þ1 1 > > > ym;kmþ2 for K pkoN; :N k m¼kK þ2 where L=min(L,K), and K=max(L,K). The smoothed time series is obtained by adding the first r principal components. It is worth pointing out that the application of the basic SSA algorithm requires selecting the values of just two parameters: the window length L, and the number r
ARTICLE IN PRESS F.J. Alonso et al. / Journal of Biomechanics 38 (2005) 1085–1092
of principal components to be retained in the reconstruction.
1087
150 100
There are no general rules for selecting the values of parameters L and r that participate in the SSA algorithm. Their values depend on the type of signal to be analysed, and on the type of information to be extracted from the signal: tendencies, harmonics, noise. A few general rules will be given in this section, whereas, again, Golyandina et al. (2001) have presented a complete description of parameter selection. The main goal is to select a window length that produces separable and independent principal components. The number of these components to be retained is then easily determined, not only from the appearance of the principal component itself, but also from the relative value of the corresponding eigenvalue. A few remarks about the window length will help in determining its value: (1) The SVD performed on matrices obtained with a window length L is equivalent to that performed on matrices obtained with the complementary window length K=NL+1. Therefore, increasing the window length above half the time series length would reproduce results already tested with shorter window lengths. (2) The longer the window, the more detailed the decomposition. This statement, along with the previous one, indicates that the best detailed decomposition is obtained with LEN/2. (3) Despite the two previous statements, a large window length may mix noise and certain components in the case of time series with a complex structure. In this case (an example of which will be presented in Section 3), it is worth testing different window lengths, starting with LEN/2, and proceeding downwards until the appropriate separation is achieved. 2.3. Test data Four different signals will be smoothed using the method described in the previous section. A double differentiation is performed on the signals in order to quantify the effect of smoothing. First-order finite differences are used to calculate the higher derivatives (Burden et al., 1981). Signal 1: Stationary sinusoidal signal A reflective marker was attached to a uniformly rotating crank. The marker motion was captured for 8 s at a frequency of 100 Hz using three infrared cameras (Qualisys Medical AB). The vertical component of this motion (Fig. 1) was smoothed with an SSA algorithm implemented in Matlab (The Mathworks Inc.). Signal 2: Biomechanical signal (Vaughan et al., 1992; Giakas and Baltzopoulos, 1997a)
Displacement (mm)
2.2. Selection of SSA parameters 50 0 -50 -100 -150 0
1
2
3
4 Time (s)
5
6
7
8
Fig. 1. Noisy displacement (Signal 1) recorded using a motion capture system.
The second signal to be processed has been taken from the literature. It corresponds to the ‘woman’ file from GAITLAB (Vaughan et al., 1992). The signal measures the lateral displacement of a marker attached on the right tibial tubercle. Data were acquired for 0.94 s with a frequency of 50 Hz. The resulting time series length is 48 elements. The reference signal is taken as that obtained after filtering the raw data at 6.25 Hz, but different amplitude white noises are added in order to test smoothing methods. An amplitude 1, time generated, white noise was used for comparisons (Giakas and Baltzopoulos, 1997a). Signal 3: Signal with impact (Dowling, 1985; Giakas et al., 2000) This third signal has also been taken from the literature. The signal is a measure of the angular coordinate of a pendulum impacting against a compliant wall (Dowling, 1985; Giakas et al., 2000). The angular acceleration obtained from the motion capture system is compared to that obtained directly (after dividing by pendulum length) from accelerometers. Three accelerometers were used in order to average their measurements to reduce noise. The average signal, logged at a sampling rate of 512 Hz is used as the acceleration reference signal. Signal 4: Slider motion A vertical slider was moved by hand to obtain a nonstationary mono-dimensional motion of biomechanical origin. A subject was asked to move the slider randomly with fast upward and downward movements. The vertical position (Fig. 2) was acquired with the use of a marker attached to the slider and the Qualisys camera system (Qualisys Medical AB). The second derivative of this record, or of the record obtained after smoothing, is
ARTICLE IN PRESS 1088
F.J. Alonso et al. / Journal of Biomechanics 38 (2005) 1085–1092 109 108
Eigenvalue
107 106 105 104 103 102
Fig. 2. Original noisy displacement (Signal 4) recorded using a motion capture system.
101
compared to the acceleration obtained directly from an accelerometer attached to the slider. Displacement and acceleration were sampled at 200 Hz during 5.90 s, obtaining records of 1182 elements each.
0
5
10
15 20 25 Eigenvalue number
30
35
40
Fig. 3. Singular spectrum (in logarithmic scale) of Signal 1, obtained with a window length L = 40.
100 0
3. Results
-100 0
1
2
3
4
5
6
7
8
-100 0 0.5
1
2
3
4
5
6
7
8
1
2
3
4
5
6
7
8
1
2
3
4
5
6
7
8
1
2
3
4 Time (sec)
5
6
7
8
100
3.1. Signal 1: Stationary sinusoidal signal Displacement (mm)
The first time series has a total of 800 elements, and the trajectory matrix was generated by sliding a window of 40 elements in length. The plot of the eigenvalues of the L L matrix S ¼ X X T (X being the trajectory matrix) is given in Fig. 3. It has been mentioned that the norm of the elementary matrices into which matrix X may be decomposed equals the square root of the eigenvalues of S. The eigenvalues in Fig. 3 are plotted in logarithmic scale. Therefore, the contribution of the first few eigenvalues to the norm of X is much higher than the contribution of the rest. In fact, 99.9992% of the norm of X is contributed by the first two eigenvalues, whereas the rest are much lower and vary slightly. This fact usually indicates that these last eigenvalues represent noise, or at least a nonperiodic component without any latent structure (Golyandina et al., 2001). It has been pointed out that the signals reconstructed using one elementary matrix at a time are called principal components. The first five principal components for the case under analysis are plotted in Fig. 4. It is clear from the figure that the two leading components represent the main trend (note that the vertical axis is scaled differently for each component), whereas the next three components represent noise. The approximated reconstruction will then be carried out using the first two principal components. The reconstructed signal is plotted in Fig. 5a along with the original signal. No differences can be made out at the scale of the plot. Slight divergence is revealed
0
0 -0.5 0 0.1 0 -0.1 0 0.1 0 -0.1 0
Fig. 4. Individual reconstruction of the five leading components (plotted in different scales) obtained from the SVD of Signal 1, using a window length L = 40.
when the difference between the two signals (Fig. 5b), or the signals themselves (Fig. 5c), is plotted on a much greater scale. The signal-to-noise ratio is very large in this example. Nonetheless, this small noise heavily contaminates the numerically obtained acceleration signal, and it is here, in the acceleration, where the efficiency of the method can be tested. The acceleration obtained using the original raw signal is plotted in Fig. 6 along with the acceleration obtained using the approximated reconstructed signal. It is clear that smoothing or filtering is absolutely mandatory in order to obtain meaningful acceleration signals. Moreover, SSA smoothing is superior to Butterworth filtering. This statement is supported in the example
ARTICLE IN PRESS F.J. Alonso et al. / Journal of Biomechanics 38 (2005) 1085–1092
1089
150
50 0 -50 -100 -150
0
2
4 Time (s)
Displacement (mm)
(a)
(b)
6
8
Displacement (mm)
Displacement (mm)
100
2
(c)
120 119 118 117 1.72 1.74 1.76 1.78 Time (s)
1.8
1.82
0
-2
0
1
2
3
4
5
6
7
8
Time (s)
Fig. 5. (a) Noisy displacement (dotted line) and SSA-smoothed displacement (continuous line) corresponding to Signal 1. (b) Residual signal defined as the difference between the original and smoothed signals. (c) Zoom in the neighbourhood of t = 1.77 s to appreciate the slight difference between the original and SSA-smoothed signals.
Fig. 6. (a) Acceleration calculated from Signal 1 original displacement data. (b) Acceleration calculated from Signal 1 after having been smoothed using SSA. (c) Acceleration calculated from Signal 1 after having been passed through a 4 Hz cut frequency Butterworth filter.
also shown in Fig. 6. The original signal was passed through a second-order Butterworth filter with a 4 Hz cut frequency, and the result was differentiated twice to obtain acceleration. The superiority is particularly relevant in the elimination of the so-called end-point errors. Methods based on signal extension have been proposed to reduce end-point error (Giakas and Batzopoulos, 1998). No extension is necessary in the
case of SSA smoothing. Errors due to this phenomenon are negligible (Fig. 6). 3.2. Signal 2: Biomechanical signal (Vaughan et al., 1992; Giakas and Baltzopoulos, 1997a) Decomposition was carried out on a trajectory matrix obtained with a window length L=N/2=24. The first 9
ARTICLE IN PRESS F.J. Alonso et al. / Journal of Biomechanics 38 (2005) 1085–1092
1090
200
4000 3000
100
1000
Acceleration (rad/s2)
Acceleration (mm/s2)
2000
0 -1000 -2000 -3000 -4000
0
-100
-200
-300
-5000
-400
-6000 0
0.1
0.2
0.3
0.4 0.5 Time (s)
0.6
0.7
0.8
0.9
0
0.2
0.4
0.6 Time (s)
0.8
1
1.2
Fig. 7. Acceleration obtained from reference signal 2 (continuous line) and acceleration calculated from SSA-smoothed displacement data (dotted line).
Fig. 8. Reference angular acceleration corresponding to Signal 3 (continuous line) and acceleration calculated from angular displacement data smoothed with a sequential SSA decomposition (dotted line).
principal components were used for reconstruction. The tenth component is a low-amplitude high-frequency motion which is most likely due to noise. Comparison with traditional filtering techniques is quantified by taking the root mean square of the error signal (RMSE) of displacement, velocity and acceleration. The error signal is the difference between the reference signal (or its derivatives) and the signal obtained after smoothing or filtering the raw data (or its derivatives). Our method achieved RMSE=0.11 mm for displacements, versus RMSE=0.4 mm obtained by Giakas and Baltzopoulos (1997a) using power spectrum assessment (PSA). Similar improvements are achieved in velocity: RMSE=4.35 mm/s versus RMSE=10 mm/s; and in acceleration: RMSE=256.89 mm/s2 versus RMSE=400 mm/s2. Accuracy in acceleration may also be appreciated in Fig. 7, where the acceleration obtained from the reference signal is plotted along with the acceleration calculated from SSA-smoothed displacement data.
form a sequential SSA decomposition (Golyandina et al., 2001), which consists in applying SSA smoothing to an already smoothed record. A second SSA decomposition was applied to the previously smoothed signal. Different window lengths and grouping strategies were tested. The best results were obtained using the first two new components for reconstruction with L2 = 20. A value of RMSE = 23.04 rad/s2 (see Fig. 8) was reached in this case. The result is similar to the value RMSE = 23.60 rad/s2 obtained by Giakas et al. (2000) with the help of the Wigner distribution. Our reconstruction right at the instant of impact is a little worse, but our signal in the vicinity of this time instant is more accurate. The slight loss of accuracy in the SSA method is compensated by the ease with which the method is applied: fewer parameters to be chosen, and no need to extend the ends of the record in order to eliminate end-point errors.
3.3. Signal 3: Signal with impact (Dowling, 1985; Giakas et al., 2000)
3.4. Signal 4: Slider motion
Giakas et al. (2000) have recently shown that secondorder Butterworth filters do not work properly in cases like this due to the non-stationary nature of the signal. Thus, this is a good example for testing SSA smoothing versus the use of the Wigner distribution. The record of 600 elements from motion capture was processed using a window length L=100. A value of RMSE = 64.98 rad/s2 is obtained when reconstructing with the first nine components. This high error is due to the fact that part of the noise has been included in the reconstruction. To minimize this effect, one may per-
Decomposition was attained with a window length L = 100. Two different reconstructed signals will be used for comparisons, one obtained from the first 16 principal components, and another obtained from the first five components alone. These signals were numerically differentiated twice to produce acceleration records to be compared to the accelerometer output. These comparisons are made in Fig. 9, where the acceleration obtained from the raw data, from the two reconstructed signals, and directly from the accelerometer, are shown. The errors in terms of RMSE are 2.04 and 1.47 m/s2, respectively.
ARTICLE IN PRESS F.J. Alonso et al. / Journal of Biomechanics 38 (2005) 1085–1092
1091
Fig. 9. (a) Acceleration calculated from Signal 4 original displacement data (dashed line) and acceleration measured by the accelerometer (continuous line). (b) Acceleration calculated from displacement reconstructed using the sixteen leading components (dotted line) and acceleration measured by the accelerometer (continuous line). (c) Acceleration calculated from displacement reconstructed using the five leading components (dotted line) and acceleration measured by the accelerometer (continuous line).
4. Discussion The results of our study show the superiority of SSA smoothing over filtering techniques found in the literature. The SSA algorithm decomposes the original signal into independent additive components of decreasing weight. This fact is key in the success of the method when it comes to extracting the latent trend in the signal from the random noise inherent to the motion capture system. One of the main advantages to the method is the fact that the algorithm requires selection of just two parameters, namely, the window length and the number of components to use for reconstruction. The method is also very intuitive in the sense that the appearance of the different components, along with the corresponding eigenvalue, is a very good indicator for distinguishing signal from noise. It has also been shown that the method works properly on both stationary and nonstationary signals. Moreover, the method can be very easily programmed as a stand-alone automatic algorithm. Complex records may need sequential decomposition to better separate noise from signal. Nevertheless, this sequential process can also be made automatic with little effort. Convergence may be measured by means of the difference between the current and previous values of the acceleration RMS. The algorithm stops when this difference is sufficiently small.
As drawbacks, one may mention the fact that there are no fixed, objective rules for selecting the window length. Nevertheless, it has been shown in the examples that the results are not very sensitive to window length. Moreover, the grouping strategy is usually clearly indicated by the singular spectrum. However, in some cases, as in signal 3, a choice of a large window length L produces a poor separation between signal trend and noise. On the other hand, for small L, only the overall trend would be extracted, relevant information being lost in the original signal (Golyandina et al., 2001). One way to overcome the problem is to apply sequential SSA. Some components of the initial series are extracted using standard SSA, and then the components of interest are extracted by applying SSA once more. In order to automate such a procedure, one could check the power spectral density (PSD) of the reconstructed signal after each decomposition to detect the presence of noise and to decide whether or not additional decompositions are necessary. Future studies will focus on automating the process. Techniques will be devised to automatically choose the window length and the grouping strategy, and a procedure will be implemented to automatically make the decision on whether or not sequential SSA is recommended for the given signal. In conclusion, we believe that the biomechanics community will benefit from this new smoothing
ARTICLE IN PRESS 1092
F.J. Alonso et al. / Journal of Biomechanics 38 (2005) 1085–1092
technique that has proven its effectiveness with complex signals. Future studies will need to focus on the possibilities of producing an automatic algorithm, and on the possibilities of embedding the algorithm in commercial biomechanical analysis packages. The application of the method to signals with rapid abrupt jumps, such as those used by Adham and Shihab (1999), will also need to be studied in greater detail. It will probably be necessary to devise an extension of SSA to detect structural changes in the time series, and use different parameters (window length and grouping strategy) before, during, and after the structural change.
References Adham, R.I., Shihab, S.A., 1999. Discrete wavelet transform: a tool in smoothing kinematic data. Journal of Biomechanics 32, 317–321. Burden, R.L., Faires, J.D., Reynolds, A.C., 1981. Numerical Analysis. PWS Publishers, Boston. D’Amico, M., Ferrigno, G., 1990. Technique for the evaluation of derivatives from noisy biomechanical displacement data using a model-based bandwith selection procedure. Medical and Biological Engineering and Computing 28, 407–415. D’Amico, M., Ferrigno, G., 1992. Comparison between the more recent techniques for smoothing and derivative assessment in biomechanics. Medical and Biological Engineering and Computing 30, 193–204. Dowling, J., 1985. A modelling strategy for the smoothing of biomechanical data. In: Johnsson, B. (Ed.), Biomechanics, Vol. XB. Human Kinetics, Champaign, IL, pp. 1163–1167. Georgiakis, A., Stergioulas, L.K., Giakas, G., 2002. Wigner filtering with smooth roll-off boundary for differentiation of noisy nonstationary signals. Signal Processing 82, 1411–1415. Giakas, G., Baltzopoulos, V., 1997. A comparison of automatic filtering techniques applied to biomechanical walking data. Journal of Biomechanics 30, 847–850.
Giakas, G., Baltzopoulos, V., 1997a. Optimal digital filtering requires a different cut-off frequency strategy for the determination of the higher derivatives. Journal of Biomechanics 30, 851–855. Giakas, G., Baltzopoulos, V., 1998b. Improved extrapolation techniques in recursive digital filtering: a comparison of least squares and prediction. Journal of Biomechanics 31, 87–91. Giakas, G., Stergioulas, L.K., Vourdas, A., 2000. Time-frequency analysis and filtering of kinematic signals with impacts using the Wigner function: accurate estimation of the second derivative. Journal of Biomechanics 33, 567–574. Golyandina, N., Nekrutkin, V., Zhigljavsky, A., 2001. Analysis of Time Series Structure—SSA and Related Techniques. Chapman & Hall/CRC, Boca Raton, FL, pp. 13-78. Hatze, H., 2002. The fundamental problem of myoskeletal inverse dynamics and its implications. Journal of Biomechanics 35, 109–115. Pasmanter, R.A., 1995. Variational search of periodic motions in complex dynamical systems. Chaos, Solitons and Fractals 6, 447–454. van den Bogert, A.J., de Koning, J.J., 1996. On optimal filtering for inverse dynamics analysis. Proceedings of the Nineth Canadian Society of Biomechanics Congress, Burnaby, BC, pp. 214–215. Vaughan, C.L., 1982. Smoothing and differentiation of displacementtime data: an application of splines and digital filtering. International Journal of Bio-Medical Computing 13, 375–386. Vaughan, C., Davis, B., O’Connor, J., 1992. Gait Analysis Laboratory. Human Kinetics, Champaign, IL. Vautard, R., Ghil, M., 1989. Singular spectrum analysis in nonlinear dynamics, with applications to paleoclimatic time series. Physica D: Nonlinear Phenomena 35, 395–424. Vautard, R., Yiou, P., Ghil, M., 1992. Singular-spectrum analysis: a toolkit for short, noisy chaotic signals. Physica D: Nonlinear Phenomena 58, 95–126. Walker, J.A., 1998. Estimating velocities and accelerations of animal locomotion: a simulation experiment comparing numerical differentiation algorithms. Journal of Experimental Biology 201, 981–995. Woltring, H.J., 1995. Smoothing and differentiation techniques applied to 3-D data. In: Allard, P., Stokes, I.A.F., Blanchi, J.P. (Eds.), Three Dimensional Analysis of Human Movement. Human Kinetics. Champaign, IL, pp. 79–99.