Data Smoothing - Summary of Responses

[Robert Newton]

Dear networkers,

There has been many reponses to my posting requesting help on data
smoothing and differentiation. As well a dozen or so have requested a summary
of the responses. Following is a general posting of the replies that I
received. I would like to thank all who took the time to respond. You have
been a great help and I am sure others will benefit from your input. It is
this sharing of knowledge which makes Biomch-L such a valuable resource.

Cheers
--
Robert Newton Internet: rnewton@loki.une.oz.au
Lecturer in Biomechanics ACSnet: rnewton@loki.une.oz
University of NewEngland,
Northern Rivers Phone: +61 (066) 203762
P.O. Box 157
Lismore NSW Australia

************************************************** ***************************
>From roger@armus.me.utexas.edu Thu Feb 4 03:08:35 1993

Look at this reference for the basics of Recursive Butterworth Filters.

Barr, R.E., and Chan, E.K.Y., (1986) Design and implementation of digital
filters for biomedical signal processing. Journal of Electrophysiology.
Tech. 13: 73-93.

Roger Gonzalez
Mechanical/Biomedical Engineering
Engineering Teaching Center, 3.104
Austin, Texas 78712
(512) 471-5718 / 471-5615
FAX: (512) 471-7683
E-mail: roger@armus.me.utexas.edu
------------------------------------------------------------------------

>From KIRSCH@bmeucl.medcor.mcgill.ca Thu Feb 4 07:44:35 1993

cccccccccccccccccccccccccccccccccccccccccccccccccc cccccccccccccccccccccc
c
c
c DIFF: a subroutine for taking numerical derivatives of data.
c
c This algorithm is based on taking the derivative of a parabola
c which is determined from (ipts-1)/2 points on either side of the
c point of interest. This least-squares polynomial approximation
c results in an effective smoothing of the data, the degree of which
c depends upon the choice of ipts; the larger its value, the more
c the data are "smoothed". The optimal choice of ipts depends upon
c the sampling rate and the frequency content of the sampled variable.
c The algorithm given below is a generalization of the 7 point
c equation given in Schaum's outline on Numerical Analysis, page 261
c (problem 21.68) and the 5 point equation given on page 237 of the
c same book. For reference, the 5 point equation is:
c
c dy(k)/dt = (1/(10*h))*(-2*y(k-2)-y(k-1)+y(k+1)+2*y(k+2)),
c
c where h is the sampling interval and y is the sampled data.
c
c The seven point formula is:
c
c dy(k)/dt = (1/(28*h))*(-3*y(k-3)-2*y(k-2)-y(k-1)+y(k+1)+2*y(k+2)
c +3*y(k+3))
c
c
c Finally, you may need to change the syntax of the DO loop depending
c upon your implementation of FORTRAN. Good luck!
c
c Robert F. Kirsch 03 Feb 1993, Montreal.

cccccccccccccccccccccccccccccccccccccccccccccccccc cccccccccccccccccccccc

subroutine diff (h,x,dx,n,ipts)
real h,x(n),dx(n),den,c
integer n,ipts,index,i

index=(ipts-1)/2
den=0.
do (i=1,index)
den=den+2.*(float(i)**2)
repeat
den=den*h


do (i=1,index)
dx(i)=0.
dx(n+1-i)=0.
repeat
do (i=index+1,n-index)
dx(i)=0.
do (j=1,index)
c=float(j)
dx(i)=dx(i)+c*x(i+j)-c*x(i-j)
repeat
dx(i)=dx(i)/den
repeat

return
end
-------------------------------------------------------------------------

>From dapena@valeri.hper.indiana.edu Fri Feb 5 00:31:43 1993

In our lab we use quintic spline to do the smoothing and
differentiation, and are very satisfied with the results that it
produces. Have you considered this possibility, or are you committed
to digital filtering?

Jesus Dapena
-------------------------------------------------------------------------

>From awsmith@utcc.utoronto.ca Fri Feb 5 02:09:23 1993

I use a critically-damped, 4th order Butterworth filter in most
of my data smoothing applications and I find that it works
quite well. Once smoothed, the data can be differentiated
using a finite differences algorithim.

The following is the Turbo BASIC source code of the general
data smoothing subroutine. It performs the following:

1. Removes the trend line (reducing the filter transients at
both the start-up and finish of the data window).
2. Smooths the data with a 2nd order Butterworth in normal time
order (this introduces a time lag).
3. Smooths the data with the same filter, but in the reverse
order of time (removing the time lag) and this second pass
results in a 4th order filter.
4. Replaces trend in data.

The minimum start-up and finish data points required are 3 of each,
however, the filter will work better with more points. This, along
with the removal of the trend line, will result in virtually no
spurious data due to filter transients:

'************************************************* ****************************
SUB DigFil (X(2),Y(2),NF,NC,CutOff(2),TIME) 'X(Frames,Markers)
'Y(Frames,Markers)
SHARED pi,Markers$()
LOCAL i,j,temp01,temp02,temp03,temp04,slopex,slopey,l,jj
' Subroutine to smooth XY coordinate data using a critically-damped
' Butterworth-type fourth-order (zero-lag) digital filter.
' NF = number of frames: NC = number of coordinates
' CutOff = cutoff frequency (suggested: 1/5th of sampling frequency at least)
' TIME = time between frames / OR / 1/sampling frequency
' X and Y are 2-D arrays => (NF x NC)

locate 5,4rint "Smoothing coordinate data - "
NF = NF + 40
DIM XF(NF)
DIM YF(NF)
DIM XHOLD(NF)
DIM YHOLD(NF)

FOR I = 1 TO NC

locate 5,33rint Markers$(i);
print using ": X Cutoff = ##.# Y Cutoff = ##.#";CutOff(i,1),CutOff(i,2)

' Define critically-damped Butterworth coefficients => A1, B1 and B2

AA = SQR(2^.25 -1)
WAx = 2 * PI * CutOff(I,1)
WAy = 2 * PI * CutOff(I,2)
OAx = TAN(WAx * TIME / 2)
OAy = TAN(WAy * TIME / 2)
AAx = OAx / AA
AAy = OAy / AA
Bx = AAx^2
By = AAy^2
Ax = AAx*2
Ay = AAy*2
A1x = Bx / (1+ Ax + Bx)
A1y = By / (1+ Ay + By)
B1x = (-2 * A1x) + (2 * A1x / Bx)
B1y = (-2 * A1y) + (2 * A1y / By)
B2x = 1 - (4 * A1x + B1x)
B2y = 1 - (4 * A1y + B1y)

RemoveTrendLine:

FOR J = NF-40 TO 1 STEP -1
XHOLD(J+20) = X(J,I)
YHOLD(J+20) = Y(J,I)
NEXT J

FOR J = 1 TO 20
K = 21 - J
XHOLD(K) = XHOLD(K+1) - (XHOLD(21+J) - XHOLD(20+J))
YHOLD(K) = YHOLD(K+1) - (YHOLD(21+J) - YHOLD(20+J))
XHOLD(NF-20+J) = XHOLD(NF-21+J) - (XHOLD(NF-20-J) - XHOLD(NF-19-J))
YHOLD(NF-20+J) = YHOLD(NF-21+J) - (YHOLD(NF-20-J) - YHOLD(NF-19-J))
NEXT J

TEMP01 = XHOLD(1)
TEMP02 = XHOLD(NF)
TEMP03 = YHOLD(1)
TEMP04 = YHOLD(NF)
SLOPEX = (TEMP02 - TEMP01) / (NF-1)
SLOPEY = (TEMP04 - TEMP03) / (NF-1)

FOR J = 1 TO NF
XHOLD(J) = XHOLD(J) - TEMP01 - SLOPEX * (J-1)
YHOLD(J) = YHOLD(J) - TEMP03 - SLOPEY * (J-1)
NEXT J

L = 1

SecondPass:

XF(1) = XHOLD(1)
YF(1) = YHOLD(1)
XF(2) = XHOLD(2)
YF(2) = YHOLD(2)

FOR J = 3 TO NF ' Butterworth-type, critically-damped filter
XF(J) = A1x * (XHOLD(J) + 2 * XHOLD(J-1) + XHOLD(J-2))_
+ B1x * XF(J-1) + B2x * XF(J-2)

YF(J) = A1y * (YHOLD(J) + 2 * YHOLD(J-1) + YHOLD(J-2))_
+ B1y * YF(J-1) + B2y * YF(J-2)
NEXT J

FOR J = 1 TO NF ' Reverses the array in time
JJ = NF - J + 1
XHOLD(J) = XF(JJ)
YHOLD(J) = YF(JJ)
NEXT J

IF L = 2 THEN GOTO ReplaceTrendLine
L = L + 1
GOTO SecondPass

ReplaceTrendLine:

FOR J = 1 TO NF
XHOLD(J) = XHOLD(J) + TEMP01 + SLOPEX * (J-1)
YHOLD(J) = YHOLD(J) + TEMP03 + SLOPEY * (J-1)
NEXT J

FOR J = 21 TO NF - 20
X(J-20,I) = XHOLD(J)
Y(J-20,I) = YHOLD(J)
NEXT J

NEXT I

NF = NF - 40

END SUB
'************************************************* ************************

* Andrew (Drew) Smith, Ph.D. | INTERNET: awsmith@utcc.utoronto.ca *
* Director of Research | BITNET: awsmith@utorgpu *
* Lyndhurst Spinal Cord Centre | COMPUSERVE: 70620,1234 *
--------------------------------------------------------------------------

>From dapena@valeri.hper.indiana.edu Fri Feb 5 09:23:31 1993

Quintic spline is supposed to be slower than digital filters.
I don't know how much slower, but that is what I hear.


One of the nice things about quintic spline is that the
derivatives of the quintic spline equations allow you to get
instantaneous velocity and acceleration data. I like that a lot.

Jesus Dapena
--------------------------------------------------------------------------

>From BI_SINCLAIR@COCO.CCHS.SU.OZ.AU Fri Feb 5 10:49:31 1993

Formulae for implimenting Butterworth filters are available in Winter (1990). A
table of coefficients for common cut-offs is given as well as formulae for
calculating the coefficients for other cut-offs. I often just implement the
equations in a spreadsheet. They are very simple, however and are easy to code.

Whether the Butterworth filter is most appropriate I am not sure. I am
suspicious of a technique which always smooths towards the mean when you may
not be sure of the direction of the error. I am sure others will be better
qualified than I to answer this.

Peter Sinclair
The University of Sydney
--------------------------------------------------------------------------

>From MEEV213@orange.cc.utexas.edu Mon Feb 8 11:01:44 1993
Reference:nce
Barr, R.E. amd Chan, E.K.Y. (1986) Design and implementation of digital
filters for biomedical signal processing. Journal of Electrophysiology
Tech. 13:73-93.
--------------------------------------------------------------------------

>From sgard@merle.acns.nwu.edu Mon Feb 8 11:02:13 1993

Butterworth filters have been used by many researchers for filtering kinematic
data of human movement. However, the Butterworth filter imparts a nonlinear
phase lag to the signal. This causes the input signal's frequency components
to be shifted in time as a function of frequency, resulting in an output
signal distortion known as 'phase distortion'.

David Winter, in his book "Biomechanics and Motor Control of Human Movement"
(1990), describes a filtering process known as 'bidirectional filtering' in
which nonlinear phase lag is removed from a signal. He passes the signal
through a 2nd order Butterworth filter, reverses the output data sequence in
time, then passes this signal through the filter once again. The final output
sequence has zero phase lag. The order of the effective filter is twice that
of the designed filter since the signal is filtered twice. In this case the
effective filter is a 4th order Butterworth. Note, however, that the cutoff
frequency of the designed filter will be different from the effective filter's
cutoff frequency.

Winter gives 2nd order Butterworth filter coefficients for several ratios of
cutoff frequency/sampling frequency, as well as implementation algorithms and
a more indepth discussion of bidirectional filtering.


Steven Gard
PhD Student
Northwestern University Prosthetics Research Laboratory
sag@merle.acns.nwu.edu
----------------------------------------------------------------------------

>From STMCCAW%ILSTU.BITNET@VMD.CSO.UIUC.EDU Mon Feb 8 11:05:56 1993

For coefficients, see Winter, DA Bioemchanics of Human Movement, 2nd
edition. Wiley Interscience, 1990, isnb 0-471-50908-6, stored in
libraries as QP303.W59 1990. Page 38 and 39 has the algorithm and
coefficients. Graham Wood, another aussie, published a review article
several years ago in Exercise and Sport Science Reviews, comparing the
different filtering techniques, as did Pezzack et al in the Journal of
Biomechanics. Tim Derrick has programmed the digital filter from the
Winter text. Contact Tim at Dept of Exercise Science, U of
Massachusetts, Amherst MA 01003 USA.

Steve McCaw
Dept of HPERD
****************************** End of Posting *****************************