Fourier Series Approximations and Low Pass Filtering: Facilitating Learning of Digital Signal Processing for Biomechanics Students Steven J Elmer and James C Martin Sportscience 13, 1-8, 2009
(sportsci.org/2009/sjejcm.htm) |

Filtering raw biomechanical data (e.g., position, ground reaction force) to remove noise is a key first step that must be performed prior to further analysis such as determining velocity and acceleration or performing inverse dynamic calculations of joint torque. Most techniques for filtering raw biomechanical data are designed to remove noise above (low pass filtering) or below (high pass filtering) a specified cutoff frequency. Although low frequency noise and high pass filtering are used in some biomechanical applications (e.g., EMG) this paper and the related spreadsheet will deal only with low pass filtering and will use the Butterworth low-pass filter (Winter, 2005). Regardless of the type of filter used, the student must grasp the concept that data contain components that occur at various frequencies. Several authors of biomechanics textbooks (e.g., Enoka, 2008;
Griffiths, 2006; Robertson et al., 2004; Winter, 2005) introduce the notion
of frequency content in data within the context of a Fourier series
approximation. This technique can be
used to approximate data representing any movement cycle that occurs over a
known time interval (T) or at a known frequency ( We have used this Excel spreadsheet to teach graduate and honors undergraduate biomechanics courses and have found that it facilitates hands-on learning of a topic that can be quite abstract. We believe that this paper and related Excel spreadsheet will be useful as a teaching tool in higher education and will be helpful for anyone interested in analyzing biomechanical data who does not have a strong background in digital signal processing. The Excel spreadsheet may also be helpful to many other spreadsheet users because it demonstrates a general method for determining any type of regression coefficients (including non-linear) using the first-principles approach of minimizing sum of squared error. Please note that this spreadsheet and the instructions were derived using Excel from Office 2003 and newer versions of Excel may function differently. Also, help with the solver function is available from this link at the Microsoft website. ## Fourier Series Approximation## Zero OrderThe most basic approximation to a signal (e.g., position, force) is the mean of that signal over the entire time interval. In the Fourier series approximation the mean is referred to as the zero order approximation and given the coefficient a0. A zero order approximation to a signal would be written, for example, as X(t) = a0, where t represents time. The mean will give some information about the signal but it will provide no detail about the variation of the signal within the movement (at any specific point in time). Thus, the a0 coefficient is not generally used by itself but rather in combination with higher order coefficients. ## First Order and Basic IdeasMany movement patterns are generally sinusoidal in nature (e.g.,
limb kinematics during gait) and thus a sine function may provide a
reasonable approximation to the movement: X(t) = sin(t). In this case, the time variable “t” refers
to the instantaneous time at a point within the movement. That time value must be scaled to the
overall time (T) or frequency ( Further, the input to a sine function must be an angle, so the time variable must be expressed as a portion of a complete angular cycle (2p) and the variable becomes 2pt/T so that X(t) = sin(2pt/T). In this way, a complete cycle of a sinusoid is generated as “t” goes from 0 to T. The signal might not oscillate equally above and below zero but rather above and below the mean value (a0) and thus the approximation will be improved by adding the mean to the sine function: X(t) = a0 + sin(2pt/T). The sine function must also be modified to take into account the amplitude (a1) of the movement (how far above and below the mean the movement oscillates; Figure 1) so the approximation becomes: X(t) = a0 + a1sin(2pt/T). Finally, the movement may not follow a true sine function in that it may not be zero at the initial point of the cycle. To correct for such an offset, the angle within the sine function should be corrected by adding some offset value q1 (Figure 1). Thus, the equation for a 1st order Fourier series approximation is: X(t) = a0 + a1sin(2pt/T+q1). This is termed a 1st order approximation, because it only includes the frequency of the overall movement. The meanings and functions of variables T, a0, a1, and q1 are illustrated in Figure 1.
In the Excel spreadsheet, the tab labeled 1st Order illustrates
this approximation. The coefficients
(a0, a1, and q1) were
determined iteratively using the Solver function in Excel. Note that the a0 coefficient determined
by the Solver function is equal to the arithmetic mean of the raw data (shown
in cell G9 of the 1st Order tab). The
Solver function must be installed as an Add-in (Figure 2; select Tools menu;
click Add-in, check Solver Add-in, click OK).
The reader can examine the workings of this function by entering some
other values in cells G5-G7 (e.g., 0’s) and then using Solver to achieve a
minimum value of the sum of squared error term (SSE; cell E2; highlighted in
red). Solver will do this by changing
coefficient values in cells G5-G7 (highlighted in yellow) to minimize the sum
of squared error term. From the tools
menu select Solver and the Solver dialog box will appear (Figure 3). Once the correct cells and “Min” criterion
are selected, click the solve button and the coefficients will be determined
iteratively. Note that this first
order approximation captures a great deal of the character of the PRF data
with an r
## Higher Order ApproximationsFourier approximations are based on harmonics, or multiples of the
central frequency of the movement. In this
particular data set, the central frequency was 2 Hz. Thus each harmonic is 2
Hz and each multiple of that harmonic represents frequencies within the
movement. That is, patterns that oscillate once within the movement occur at
the central frequency, patterns that oscillate twice within the movement
occur at twice the central frequency, and so on. In the spreadsheet tabs labeled 2 The additional frequencies allow the higher order approximations
to follow additional nuances of the original signal. As series order is increased the approximations
tend to converge upon the raw PRF data (Figure 4). Indeed, the 4th order approximation
captures most of the nuances of the raw PRF data (r
The spreadsheet can be used to examine the functions of the Fourier coefficients by performing the following exercise. In the 1st Order tab, reset the a0 coefficient to a vale of zero. When this change is made, the approximated data will oscillate about zero, rather than about the mean demonstrating that the a0 coefficient forces the approximation to vary about the mean. As mentioned above, the reader may note that the a0 coefficient is the same as the mean for the raw PRF data (cell G9). Once the effect has been observed, the reader may use the “undo” function (select edit; click on undo) to restore the coefficients or may use solver to reestablish them. Second, the reader may reset the a1 coefficient to zero and will note that the approximation becomes a straight line equal to the mean. This demonstrates that the a1 coefficient represents the amplitude of the excursion above and below the mean. Next, the reader may set the q1 coefficient to zero and will note that the approximated data is shifted horizontally relative to the raw data. This will demonstrate that the function of the q1 coefficient is to adjust the relative position of the sine function. These effects are also illustrated in Figure 1. Additionally, in the 2nd Order tab, setting the a1 coefficient to zero will reveal the effect of the a2 coefficient which causes the approximation to oscillate at twice the central frequency (twice with the movement cycle). This is only a brief exercise and the reader may wish to explore other coefficients and tabs. ## General MethodWe have used Solver to determine coefficients that minimize the sum of squared error term between the original data and the modeled approximation. This technique can be used to determine coefficients for any linear or non-linear approximation. To do this, one simply writes an equation for the model of choice which includes coefficients and independent variables (column “C” in our approximation spreadsheets). That equation should be “filled down” a column of an Excel spreadsheet so that an approximation is calculated for each data point. Further, a set of cells should be used as a location for the coefficients (column “G” in our spreadsheet). Any value can be selected for the initial coefficients (e.g., 0’s or 1’s). In the next column, one should form a squared error term as the difference in the raw value and the approximation squared (column “D” in our spreadsheet). All of these squared error terms are then summed to form the sum of squared error term (cell “E2” in our spreadsheet). Solver performs an iterative process to determine coefficient values that minimize this sum of squared error term. We have used this technique to determine coefficients for functions as basic as quadratic equations and as specialized as a Hill type force-velocity equation (rectangular hyperbola). The interested reader will likely find many applications for this technique. ## Butterworth FilterThe fourth-order zero-lag Butterworth filter is often used in biomechanical analyses and is described in detail by Winter (2005). Briefly, this filter produces a weighted average of data from several time points and the weight on each time point (the a0, a1, a2, b1, b2 coefficients) determines the cutoff frequency (the frequency above which noise is removed). Note that these Butterworth coefficients are not related to the Fourier approximation coefficients of similar names. In our example filter, the weighted average includes the point of interest as well as the two preceding time points. This process of averaging time points prior to the time of interest causes the filtered data to “lag” behind the raw with respect to time. To correct for this lag, data are filtered once in forward time order, then again in reverse time order to produce filtered data that are properly aligned in time. In addition to correcting for lag, the second filtering in the reverse direction creates a sharper cutoff and this two-pass filter is referred to as a fourth-order zero-phase shift (or zero-lag) filter. We encourage the reader to consult a textbook (e.g., Winter 2005) for a thorough treatment of Butterworth filter. In the Excel spreadsheet tab labeled Butterworth, we have implemented just such a filter. We will not present the details of this filtering technique in this paper but rather provide it to illustrate the effects of filtering the PRF data at various cutoff frequencies and to allow the reader to explore the similarities and differences between Fourier series approximations and filtered data. The Butterworth coefficients (K1, K2, K3, a0, a1, a2, b1, b2) are calculated in the spreadsheet as described by Winter (2005, page 47). The a and b coefficients are the weights applied to the data. The K coefficients are used to calculate the a and b coefficients based on the desired cutoff frequency and the known sampling frequency. The reader can change the cutoff frequency simply by changing the value in cell D2 (highlighted with a yellow background). As an exercise, we encourage the reader to observe data filtered at a range of cutoff frequencies from 1 to 16 Hz. As the cutoff frequency increases, the filtered data will more closely approximate the raw data, including the inherent noise. Some intermediate cutoff frequency (e.g., 8 Hz) will “pass” the signal through while filtering the higher frequency noise. Determining the proper cutoff frequency for any specific data set requires knowledge of the signal and noise and goes beyond the intended scope of this paper. The interested reader may see one technique for determining cutoff frequency in Winter (2005). ## Velocity and AccelerationIt is common practice to perform finite differentiation of position data to obtain values for velocity and acceleration of limb segments (velocity equals change in position divided by change in time, and acceleration equals change in velocity divided by change in time). Therefore, we have also implemented a fourth-order zero-lag Butterworth filter in the tab labeled Velocity and Acceleration. For this tab we used kinematic data for foot angle during cycling (movement frequency 2 Hz, sampled at 240 Hz). The reader may note that even the unfiltered data appear to be quite smooth (smoother than the raw PRF data) and one might believe that the raw kinematic data could be used for subsequent analyses. However, as demonstrated in Figure 5 and in the Excel spreadsheet, the effects of filtering on angular velocity (w) and acceleration (a) are much more dramatic and these comparisons serve to underscore the importance of properly filtering data prior to differentiation. Note that these angular velocity and acceleration data were obtained by finite differentiation of the filtered position data, not by filtering the differentiated raw data. Filtering the raw position data removes a smaller magnitude of noise and should provide more accurate derivatives. ## Fourier vs ButterworthIn the tab labeled Comparison, we have implemented a fourth-order zero-lag Butterworth filter and linked to the Fourier series approximations (every other point has been omitted to reduce overlap in the figure). As with the other filters, the cutoff frequency can be adjusted by simply typing a value in cell I2 (highlighted with a yellow background). This will facilitate an exploration of frequency content by comparing various values of cutoff frequency with the calculated Fourier series approximations. Note that in our data the movement frequency was 2 Hz and so the 4th order approximation generally includes data up to 8 Hz. When the cutoff frequency for the Butterworth filter is set 8 Hz, the data essentially overlap those of the 4th order approximation. Although the filtered data and approximation data are quite similar they are not identical. The differences arise mainly for two reasons. First, the Butterworth filter allows some noise above the cutoff frequency and eliminates some signal below the cutoff frequency (Winter, 2005). For additional information on the “sharpness” of a cutoff frequency, we direct the interested reader to Figure 1.26 in Enoka (2008), and for a comparison of 2nd vs. 4th order Butterworth filter, to Figure 4.16 in Bartlett’ (2008). Secondly, the Fourier approximation only includes discrete frequencies that are integer multiples of the central frequency, whereas the actual signal as well as the Butterworth filtered signal include a continuous spectrum of frequencies. Thus, while filtering at a cutoff frequency of 8 Hz produces similar data to those obtained with a Fourier approximation they are not identical and should not be confused. Also, in some cases a Fourier Transform procedure can be used as a filtering technique but that technique is a different from what we have done in our spreadsheet and we do not intend to discuss that in this manuscript.
## ConclusionWe believe this brief tutorial and Excel spreadsheet will facilitate learning in graduate and upper division undergraduate biomechanics courses. How it will ultimately be used will depend on the curiosity and background of the reader. We do not intend this paper and Excel spreadsheet to be all encompassing but rather to be one helpful aid to understanding of digital signal process within a biomechanics course. An added benefit is that the general method used to determine the Fourier coefficients can be used to determine the coefficients for any linear or non-linear model within an Excel spreadsheet. ## ReferencesBartlett R (2007). Enoka RM (2008). Griffiths IW (2006). Robertson DG, Caldwell G, Hamill J, Kamen G, Whitlesey SN (2004). Winter DA (2005). Published March 2009 |