Estimation procedures affect the center of pressure frequency analysis

Braz J Med Biol Res, July 2009, Volume 42(7) 665-673

Estimation procedures affect the center of pressure frequency analysis

T.M.M. Vieira^1,2, L.F. Oliveira² and J. Nadal³

¹Laboratorio di Ingegneria del Sistema Neuromuscolare, Politecnico di Torino, Torino, TO, Italia

²Laboratório de Biomecânica, ³Programa de Engenharia Biomédica, Universidade Federal do Rio de Janeiro, Rio de Janeiro, RJ, Brasil

References ^{Correspondence and Footnotes}

Abstract

Even though frequency analysis of body sway is widely applied in clinical studies, the lack of standardized procedures concerning power spectrum estimation may provide unreliable descriptors. Stabilometric tests were applied to 35 subjects (20-51 years, 54-95 kg, 1.6-1.9 m) and the power spectral density function was estimated for the anterior-posterior center of pressure time series. The median frequency was compared between power spectra estimated according to signal partitioning, sampling rate, test duration, and detrending methods. The median frequency reliability for different test durations was assessed using the intraclass correlation coefficient. When increasing number of segments, shortening test duration or applying linear detrending, the median frequency values increased significantly up to 137%. Even the shortest test duration provided reliable estimates as observed with the intraclass coefficient (0.74-0.89 confidence interval for a single 20-s test). Clinical assessment of balance may benefit from a standardized protocol for center of pressure spectral analysis that provides an adequate relationship between resolution and variance. An algorithm to estimate center of pressure power density spectrum is also proposed.

Key words: Center of pressure; Spectral analysis; Quiet standing; Estimators; Median frequency

Introduction

Despite being a simple task, maintaining a quiet standing posture requires the interaction of different stimuli, which are coded by specialized sensory receptors, and the coordination of muscle responses to stabilize center of gravity spontaneous sways (1). Since the center of pressure (COP) under the feet and center of gravity displacements are equivalent at low frequencies, the ability to control body sway is usually assessed using a force plate (2). However, the choice of suitable stabilometric parameters is a crucial matter for pathologic classification and diagnosis and for monitoring rehabilitation.

Several studies focused on the identification of reliable and non-redundant parameters to describe COP time series (3-5). Chiari et al. (6) outlined the influence of five biomechanical factors on force plate measurements, while Lafond et al. (7) determined the number of trials to increase the intraclass correlation coefficient of averaged stabilometric parameters up to 0.90. According to these investigators (7), more than 20 trials are recommended to reduce the within-subject variability of mean and median frequency of COP, although they did not mention the procedures applied to estimate the COP power spectrum.

Frequency domain analysis seems to be useful for evaluating postural disturbances caused by specific diseases or conditions (8-11). Patients affected by a chronic tension-type headache exhibit increased COP frequencies compared with a group of healthy subjects (9). On the other hand, children with myelomeningocele show an opposite trend, characterized by reduced median frequency values (10).

Although the COP time series result from a stochastic process, most clinical studies reporting frequency analysis of COP sway generally neglect estimation errors of the computed power spectral density (PSD) function. Some experimental designs lack a description of the spectral estimation procedures (12), while others do not take the estimate bias and variance into account (3).

The objective of this study was to evaluate the influence of four different procedures on the estimated COP power density spectrum (ŜC_COP) during quiet standing tasks. Two procedures relate to COP signal acquisition, the sampling rate (SR) and test duration, while the number of segments (K) and the use of detrending apply to COP time series. Some guidelines are proposed and a script is provided in ^{Appendix A} to estimate the actual COP power density spectrum (S_COP).

Material and Methods

Thirty-five healthy subjects (20 males and 15 females), with a mean age of 23.7 years (range, 20-51 years), weight 57.62 kg (range, 54-95 kg) and height 1.64 m (range, 1.6-1.9 m), participated in the experiment. The Ethics Committee for medical research of the Federal University of Rio de Janeiro approved the experimental protocol, and all subjects gave written informed consent before enrollment in the study. Subjects presenting neurological or musculoskeletal impairments that affect equilibrium were excluded from the study.

Prior to the 60 s stabilometric test, all subjects were instructed to stand in the upright position with heels 6 cm apart, 10° of feet abduction, arms relaxed along the body and facing forward. Figure 1 briefly depicts the protocol for positioning the foot. To assess reliability according to different test durations, subjects repeated this task five times with a 3-min rest period. COP time series was measured using a force plate (AccusWay^PLUSS, AMTI, USA) with a built-in 16 bits A/D converter and 0-100 Hz bandwidth, and thus sampled at 100 Hz and low-pass filtered using a bidirectional second order Butterworth filter with 12.5 Hz cut-off frequency. Spectral analysis was carried out using custom algorithms developed in Matlab 6.5 (The Mathworks, USA).

Applying the fast Fourier transform (FFT) algorithm to a single and short segment of COP time series erroneously presumes an underlying deterministic process, so that the respective spectrum estimates do not converge due to inconsistency and bias effects (13). Some elementary concepts to provide an unbiased and consistent estimator of the PSD of a stabilometric record will be outlined in the following subsections.

Partitioning COP spectrum

Since body sway results from a stochastic process, spectral analysis is affected by statistical properties of the spectral estimator. The estimator based on the discrete Fourier transform (ŜC_COP) applied to a single segment of COP record is inconsistent, since its variance is equal to the quadratic magnitude of S_COP (14). Dividing the COP signal into K non-overlapped segments and averaging the resulting spectra is an alternative to reduce the PSD estimator variance known as periodogram or Bartlett method (13):

(Equation 1)

where Ŷ_k (ƒ) is the PSD estimated for each segment of the signal:

(Equation 2)

I_L(ƒ) is the Fourier transform of each epoch of COP time series, multiplied by a Hamming window with L samples. The window size (L) is the nearest integer defined by the total number of samples (N), L = round(N/K). As the use of a Hamming window attenuates the boundaries of each segment, the normalization factor U assures that Ŷ_k (ƒ) total energy is unitary:

(Equation 3)

where z(l) represents the applied window function. Note in Equation 1 that for a single segment, ŜC_COP becomes equal to Ŷ_k.

For each frequency the representation of I_L(ƒ) in the complex Cartesian plane is composed of a real and imaginary part, except for the real valued 0 Hz (mean value) and SR/2 frequencies, and each sample of ŜC_COP is given by the addition of both I_L(ƒ) squared components. As the real and imaginary components are Gaussian random variables and the segments are independent, which was assured by non-overlapping, each value of the estimated COP power spectrum adheres to a χ² distribution function with 2K degrees of freedom (13). Thus, a confidence interval for S_COP is estimated:

(Equation 4)

where LL and UL are, respectively, the lower and upper limits for the interval, obtained directly from the table of χ² distribution for 2K degrees of freedom and a significance level α. In this study, the confidence interval to S_COP was estimated for two, six and twelve degrees of freedom, according to COP data partitioning into 1, 3, and 6 segments, using α = 0.05 and a linear fit to remove each segment trend before applying FFT.

Sampling strategy

A single COP time series results from the multiplication between the continuous signal by an impulse train with 1/SR period. As in frequency domain multiplication corresponds to the convolution of both signals, S_COP becomes periodic at SR and only frequencies within ±SR/2 are represented. The actual signal energy residing at frequencies beyond these limits artificially increases the estimated power of higher frequencies. This phenomenon, known as "aliasing", is prevented when the Nyquist theorem is applied, assuring that SR ≥2 f_max, where f_max is the maximum frequency of the signal (14). As stabilometric signals are real valued, the estimated power spectra are symmetrical and only the components between 0 Hz and +SR/2 are considered, enclosing a total of N/2 components.

Test duration

Since the discrete PSD is defined at frequencies given by M(SR/N), where M is an integer between 0 and N/2, the length of the COP time series affects ŜC_COP spectral resolution and the amplitude of harmonics unrepresented in the estimated spectrum. Spectral resolution can be improved by increasing the number of samples while keeping sampling frequency constant, which means a longer stabilometric record.

In frequency domain, a stabilometric record corresponds to spectral convolution between the actual COP power spectrum and a sine(x)/x function, which produces a biased spectrum estimate due to leakage phenomenon. This leakage is proportional to the width of sine(x)/x central lobe, which shortens as the record length increases (see Ref. 14 for a full description). Therefore, a stabilometric test with increased duration improves spectral resolution and reduces bias error of ŜC_COP.

Detrending flag

Partitioning COP time series into epochs of smaller duration may produce non-stationary epochs, mostly due to the low frequencies in the actual record (i.e., a 0.1-Hz harmonic would be observed as a trend in epochs of 5 s). This effect is suppressed by removing each segment trend (14), in general given by the fitted first order polynomial. A linear detrending also reduces discontinuities at segment boundaries caused by truncation, eliminating spurious components at higher frequencies, at the cost of reduced power in very low frequency components.

Signal processing

Three estimators were applied to evaluate the number of segments and sampling rate effect on ŜC_COP spectrum, according to the signal division in K segments (1, 3, and 6) and to the length of the digital COP signal (preserved and decimated to 50 and 20 samples/s), respectively. The influence of the test duration was investigated using the initial duration of 60 s (6000 samples), which was then truncated to 40 s (4000 samples), and 20 s (2000 samples), keeping the same 100 Hz sampling rate. Figure 2 shows a brief overview of these procedures. The mean value and the first order polynomial fit were considered for assessing the detrending effects on COP estimated spectra.

To assess K and SR effects, COP spectra were estimated using the Bartlett method described in the previous section, with a fixed number of 6 segments and linear detrending. Considering test duration and trend removal, ŜC_COP was estimated using the Welch periodogram method (13), dividing the signal into 7 segments with 50% of overlapping and applying a linear detrend on each ^segment1. Although overlapping increases the variance of the averaged power spectrum estimator, it allows an increased frequency resolution due to longer epoch duration.

The median frequency of anterior-posterior COP displacements (F50), which divides COP spectrum into two portions of equal power, was the spectral descriptor used to evaluate the influence of acquisition (sampling rate and test duration) and estimation (data segmentation and detrending) parameters on ŜC_COP.

Statistical analysis

One-way ANOVA was used to test F50 equality between different estimators (15). By means of the t-test for paired samples, the linear detrending effect on F50 was evaluated. As a consistent estimate is expected from ŜC_COP estimated with the Welch method, F50 reliability was assessed according to the different test durations of 20, 40, and 60 s, using the intraclass correlation coefficient proposed by Shrout and Fleiss (16). This coefficient is the ratio between the variance due to subjects and total variance evaluated for all five trials.

Results

The F50 values across subjects decreased significantly when reducing the number of COP signal divisions, increasing test duration or applying mean detrending (Table 1). As the COP spectrum is mainly composed of low frequencies, F50 was invariant with sampling rate. Longer test durations exhibited high intraclass correlation for F50 values (Table 2), which means increased ŜC_COP reliability.

ŜC_COP shape and magnitude did not change with sampling rate strategy (Figure 3A), maintaining F50 value almost constant (P = 0.997). To emphasize this effect, all spectra presented in Figure 3A were depicted within 0-1 Hz, although ŜC_COP upper limit increased as a function of SR. Note that the first harmonic magnitude and position remains the same independent of the SR used.

The confidence interval of S_COP estimates reduced with COP time series fragmentation using an increased number of segments (Figure 3B). A reduction of first harmonic power is also evident, leading to a significant increase of median frequency (P < 0.001). When K = 1 segment, the estimated COP spectrum (thicker line) was not observed due to a large confidence interval, which corresponded to 2 ŜC_COP/UL (= 2040 ŜC_COP).

Shorter record lengths resulted in power attenuation and peak frequency displacement (Figure 4), which corresponds to the first harmonic of each estimated spectrum. As expected, median frequency presented a statistically significant decrease with test duration (P < 0.001), although increasing intraclass correlation and improving F50 reliability (Table 2).

The linear trend removal greatly reduced ŜC_COP power below 0.3 Hz (Figure 5), providing median frequency values significantly higher compared with mean subtraction (P < 0.001). Moreover, the trend was not removed from the windowed COP segments when estimating its power spectrum using the "psd" function supplied by Matlab 6.5, which provides a biased COP spectrum. By applying t-tests to compare the F50 values estimated with the Matlab "psd" function and the proposed algorithm, a slight increase (P = 0.16) was observed for the median frequency values estimated with the algorithm and using the mean value for detrending. The F50 significantly increased (P < 0.001) for the spectra estimated with the proposed algorithm and removing the linear trend of each segment (Table 3).

Figure 1.
Schematics of foot positioning on the force plate and an example of a short epoch of a typical stabilometric record. The trace on the force plate surface corresponds to the centre of pressure (COP) trajectory on the surface, starting at 0 s (square) and lasting 10 s (circle). The bottom panel illustrates the time series of the same COP trace along the anterior-posterior axis in arbitrary units (a.u.).

Figure 2.
Schematic representation of the analysis procedures applied to investigate sampling rate (SR), number of segments (K) and test duration (T) effect on the center of pressure (COP) estimated spectrum. A, Sampling each record at 100, 50 and 20 Hz; B, Signal division according to grid lines as 1, 3, and 6 segments; C, COP time series truncation at 20, 40, and 60 s.

Figure 3.
Center of pressure (COP) power spectrum estimated for 20, 50, and 100 Hz sampling rates (SR) by partitioning the record into six non-overlapped epochs (A). B, COP power spectrum (thick line) and its 95% confidence interval (thin line) estimated according to signal partition into 1, 3, and 6 segments (K).

Figure 4.
Spectral estimates for a single COP record with different durations (full line = 60 s, long dashed line = 40 s, and short dashed line = 20 s). COP = center of pressure.

Figure 5.
Effect of the mean (dashed line) and linear (solid line) detrending flag on the center of pressure (COP) power density spectrum averaged for all estimates. Note the difference between COP spectrum estimated with Matlab "psd" function (top) and with the algorithm proposed in ^{Appendix A} (bottom), which removes zero frequency value and provides the density spectrum.

Table 1.
F50 values for each procedure employed for the COP spectral estimate.

Table 2.
Mean values of F50 for each trial and test duration and the intraclass correlation coefficient (ICC) calculated across trials.

Table 3.
Mean ± SD values of F50 estimated with the proposed algorithm (Custom) and the Matlab "psd" function, removing the mean of linear trend of each segment.

Discussion

This study focused on a technical approach to investigate the validity of COP frequency analysis widely applied to assess postural stability in clinical settings. When using direct methods (periodogram) for estimating the PSD of COP displacements (S_COP), the observed median frequency values were dependent on the number of segments, test duration and the trend removal approach.

Spectral invariability with sampling rate

Corroborating other studies (3,6), COP sway during quiet standing was restricted to low frequencies (<1 Hz) independent of the spectral estimation procedure carried out, likely explaining F50 and ŜC_COP invariability with the time resolution applied. However, contribution of frequencies up to 2 Hz during quiet standing is reported in the literature (17-19). Isolating the lower limb proprioceptive reflex by strapping subjects to a vertical support and applying a random disturbance of small amplitude, Fitzpatrick et al. (17) observed reflex responses within 1-2 Hz frequency band, suggesting an operating range for lower limb proprioceptive receptors. Furthermore, for COP spectrum estimated using 20 Hz sampling rate and 12.5 Hz cut-off frequency, any energy within 10.0-12.5 Hz would emerge as an "alias" into the 7.5-10.0 Hz band and increase F50 values, which was not confirmed. Considering the lack of evidence for contributions above 2 Hz, a conservative minimum sampling rate of 5 Hz should be considered, not only to assure COP time series properties, but also the representation of ŜC_COP higher frequencies. In addition, an anti-aliasing (low-pass, 2 Hz cut-off) filtering of COP signal should precede any A/D conversion or decimation.

Consistent spectrum estimate and segmentation effects

Although the duration effect of a single segment (15, 30, 60, and 120 s) on ŜC_COP was previously addressed (3), the estimate variability was not. From the theory of periodograms, it is well known that spectrum consistency increases with signal segmentation, as the confidence interval of S_COP becomes proportionally narrow (13,14). On the other hand, when using fixed signal duration, the number of segments and their length are inversely related so that partitioning decreases spectral resolution. Power attenuation of lower COP frequencies and leakage to neighboring frequencies caused by signal partition may overestimate the ŜC_COP higher frequencies contribution.

Since very low frequencies may represent a linear trend for a short COP epoch, shortening segments duration emphasizes the linear trend removal effect, explaining the significant increase of F50 values when increasing K, in addition to the low resolution produced by an increase in frequency spacing. Resolution and consistence must be considered before test application, since both are inversely dependent for the same test duration.

Resolution, detrending and F50 reliability

As test duration increased COP spectrum compressed toward lower frequencies, so that the F50 calculated with 2000 samples segment was 137% greater than the F50 estimated using non-truncated data. The significant increase of F50 obtained for shorter test duration results merely from spectral estimation method related to resolution and detrending effects. While lower spectral resolution caused a larger spacing between harmonics, affecting the median determination, the linear trend removal applied to assure stationarity before DFT calculation produced a major effect on F50 estimations. When reducing test duration, detrending effects were critical, attenuating the power of ŜC_COP frequencies below 0.3 Hz (Figure 4).

As expected, detrend played a significant role in the median frequency increase. Instead of an energy leakage as observed for shorter segment estimations in Figure 3B, linear detrending affected only lower frequencies, although enough to increase the relative contribution of higher frequencies to F50. On the other hand, the mean removal eliminates only the 0 Hz component preserving ŜC_COP shape (Figure 5).

The periodograms estimated with 50% overlapping allowed a satisfactory relationship between spectral resolution and consistence, with the first harmonic below 0.1 Hz for a 60-s epoch and small variance as well. The lower intraclass coefficient values of F50 mentioned by Carpenter et al. (3), even to longer test duration (0.45-s to a 60-s test duration), may reflect the higher variance of COP spectrum estimated from a single segment. In the present study, intraclass correlation coefficient values ranged from 0.82 to 0.96 for F50 calculated with signal durations ranging from 20 to 60 s, respectively (Table 2). The estimation method based in overlapping segments confers high within-subject repeatability even for estimates of short duration (20 s), which is adequate in clinical surveys that avoid longer test durations (i.e., fall risk assessment, rehabilitation monitoring).

Standardization

Even though median frequency is usually evaluated in clinical studies, the lack of standardization and description of the procedures for spectral estimation constrains the comparison of outcomes from different studies. Comparing studies that apply spectral analysis in similar samples of 13 young healthy subjects, the values reported for anterior-posterior median frequency in each study were, respectively, 0.06, 0.6, and 0.18 Hz (8,11,12). Such differences likely emerge from anthropometry and base of support features or spectral estimation procedures. However, Chiari et al. (6) observed that F50 is almost independent of individual biomechanical factors and support base configuration, attributing the between-studies divergence in F50 to estimation procedures. The lower F50 values mentioned by Nardone et al. (12) may be derived from the longer test duration (51 s) compared with the other two studies (30 s). Unfortunately, except for Davidson et al. (8) who estimated COP spectra using sinusoidal tapers, the others did not provide enough information concerning estimation procedures.

Meyer et al. (20) observed an increase of F50 values when the tactile sensitivity of forefoot was suppressed but not with complete inhibition of plantar region receptors, suggesting that remaining sensorial mechanisms did not compensate for the forefoot sensitivity loss. Although providing information about S_COP estimation, the authors did not account for ŜC_COP variance. COP power spectrum was estimated partitioning a 30-s stabilometric record into 2 segments with 50% overlapping, leading to unreliable F50 values.

The standardization of a protocol that uses spectral analysis to evaluate postural control must take into account the relationship between spectral resolution and estimator variance. Stabilometric tests of 60 s are recommended to reduce ŜC_COP variance without constraining low frequency analysis by applying overlapped segmentation. When longer tests must be avoided, 20-s duration provides reliable spectra with a low 0.2-Hz spectral resolution. In such cases, the use of modern estimators based on autoregressive modeling could be advantageous, since these methods are adequate for shorter signals (14).

The procedures carried out for spectral analysis, including signal partition, test duration and trend removal, produced major alterations in amplitude distribution of COP spectrum. Anterior-posterior median frequency increased with COP time series segmentation, shorter test durations and linear detrending. These factors must be considered before applying spectral analysis. For general clinical purposes, applying stabilometric tests of 60 s, using a sampling rate of at least 20 Hz and partitioning the whole record into 6 or 7 segments seems to provide consistent estimates for COP spectrum. Tests of 20 s should be restricted to studies focusing on COP higher frequencies, or when longer tests are not advisable, and providing that frequency resolution is not so relevant.

_____________________________________

1Matlab "psd" function allows for trend removal of each segment; however, instead of a density spectrum the values supplied by "psd" correspond to the bilateral and absolute spectrum and the trend is removed before applying the desired window, biasing the estimate toward 0 Hz. The algorithm described in ^{Appendix A} corrects for the original Matlab "psd" function.

Address for correspondence: L.F. Oliveira, Departamento de Biociências, EEFD/UFRJ, Av. Carlos Chagas Filho, 540, Cidade Universitária, 21949-900 Rio de Janeiro, RJ, Brasil. E-mail: [e-mail]

Research partially supported by CNPq, Fundação Universitária José Bonifácio (FUJB), and FAPERJ. T.M.M. Vieira was the recipient of MSc scholarship from FAPERJ and CNPq. Received July 18, 2008. Accepted March 3, 2009.

• 1. Windhorst U. Muscle proprioceptive feedback and spinal networks. Brain Res Bull 2007; 73: 155-202.
• 2. Vieira TMM, Oliveira LF, Nadal J. An overview of age-related changes in postural control during quiet standing tasks using classical and modern stabilometric descriptors. J Electromyogr Kinesiol 2009 (doi: 10.1016/j.jelekin.2008.10.007).
• 3. Carpenter MG, Frank JS, Winter DA, Peysar GW. Sampling duration effects on centre of pressure summary measures. Gait Posture 2001; 13: 35-40.
• 4. Karlsson A, Frykberg G. Correlations between force plate measures for assessment of balance. Clin Biomech 2000; 15: 365-369.
• 5. Rocchi L, Chiari L, Cappello A. Feature selection of stabilometric parameters based on principal component analysis. Med Biol Eng Comput 2004; 42: 71-79.
• 6. Chiari L, Rocchi L, Cappello A. Stabilometric parameters are affected by anthropometry and foot placement. Clin Biomech 2002; 17: 666-677.
• 7. Lafond D, Corriveau H, Hebert R, Prince F. Intrasession reliability of center of pressure measures of postural steadiness in healthy elderly people. Arch Phys Med Rehabil 2004; 85: 896-901.
• 8. Davidson BS, Madigan ML, Nussbaum MA. Effects of lumbar extensor fatigue and fatigue rate on postural sway. Eur J Appl Physiol 2004; 93: 183-189.
• 9. Giacomini PG, Alessandrini M, Evangelista M, Napolitano B, Lanciani R, Camaioni D. Impaired postural control in patients affected by tension-type headache. Eur J Pain 2004; 8: 579-583.
• 10. Karlsson A, Norrlin S, Silander HC, Dahl M, Lanshammar H. Amplitude and frequency analysis of force plate data in sitting children with and without MMC. Clin Biomech 2000; 15: 541-545.
• 11. Vuillerme N, Vincent H. How performing a mental arithmetic task modify the regulation of centre of foot pressure displacements during bipedal quiet standing. Exp Brain Res 2006; 169: 130-134.
• 12. Nardone A, Tarantola J, Giordano A, Schieppati M. Fatigue effects on body balance. Electroencephalogr Clin Neurophysiol 1997; 105: 309-320.
• 13. Shiavi R. Introduction to applied statistical signal analysis London: Academic Press; 1999.
• 14. Bendat JS, Piersol AG. Random data analysis and measurement procedures New York: John Wiley & Sons; 2000.
• 15. Dawson-Saunders B, Trapp RG. Basic and clinical biostatistics Norwalk: Appleton & Lange; 1994.
• 16. Shrout PE, Fleiss JL. Intraclass correlations: Uses in assessing rater reliability. Psych Bull 1979; 86: 420-428.
• 17. Fitzpatrick RC, Gorman RB, Burke D, Gandevia SC. Postural proprioceptive reflexes in standing human subjects: bandwidth of response and transmission characteristics. J Physiol 1992; 458: 69-83.
• 18. Baratto L, Morasso PG, Re C, Spada G. A new look at posturographic analysis in the clinical context: sway-density versus other parameterization techniques. Motor Control 2002; 6: 246-270.
• 19. Winter DA, Patla AE, Prince F, Ishac M, Gielo-Perczak K. Stiffness control of balance in quiet standing. J Neurophysiol 1998; 80: 1211-1221.
• 20. Meyer PF, Oddsson LI, De Luca CJ. The role of plantar cutaneous sensation in unperturbed stance. Exp Brain Res 2004; 156: 505-512.

Appendix A

Correspondence and Footnotes

Publication Dates

• Publication in this collection
  26 June 2009
• Date of issue
  July 2009

History

• Accepted
  03 Mar 2009
• Received
  18 July 2008