filtering of the inverse kinematics results

[Farzaneh Yazdani]

Dear Biomch-L Members,
I am working on gait kinematics using ‘Opensim gait 2392’ model. I filtered the results of inverse kinematics using "butterworth 4th order, cut off 6" then i performed Qubic Spline function on them for the purpose of time normalization.
Would you please take a look at attached file to see the result?
I'm afraid if the data were manipulated too much.

...Any recommendation would be very appreciative.

Farzaneh Yazdani
PhD Candidate, SUMS school of Rehabilitation Sciences

[Allan Carman]

Re: filtering of the inverse kinematics results

Some notes on using a Butterworth filter.
i) A Butterworth filter uses previous smoothed data points (along with the current and previous raw data points) to calculate the current smoothed point. Therefore you need padding data before and after the signal of interest to tune the filter prior to reaching the first data point of interest. If you have pre and post data then this should be left on until after filtering, then the gait cycle can be extracted from filtered data. If you don't have padding data you will need to generate some. There are various options in the literature like linear extrapolation backwards or reflection about the first data point, or if we are looking at a cyclic signal like gait you could repeat the cycle before and after the current data.
ii) Butterworth filters always start from zero (the first smoothed data point). As your data does not start from zero you will introduce oscillations in the signal until the filter adjusts to your data. Therefore extra padding points are needed or alternatively for a low pass filter, remove the offset to the very first padded data point, smooth and then add the offset back on, then remove the padding data.
iii) Butterworth filters introduced a lag in smoothed signal (phase shift). To remove lag, after filtered the data for first time reverse all data points and run the same filter through the data again. This will give you a dual pass or zero lag filter.
iv) Plot residuals, the difference between smoothed and raw data points. Residuals should have a random pattern with near zero correlation, any sign of a pattern emerging in the residuals and you may have over smoothed and removed signal.

Cheers
Allan

[Farzaneh Yazdani]

Re: filtering of the inverse kinematics results

Dear Allan, i am grateful for your guidance.
actually the biggest problem that i have is that the original data does not start from zero, but after filtering the signals start from zero and this leads to unexpected patterns in some trials.
as i am working with MATLAB, can i apply oscillations in signal using this software?

regards

Farzaneh Yazdani
PhD Candidate, SUMS school of Rehabilitation Sciences

[Ton van den Bogert]

Re: filtering of the inverse kinematics results

If you work in Matlab, you should try the "filtfilt" function instead of "filter". The Matlab documentation says "filtfilt minimizes start-up and ending transients by matching initial conditions" (unfortunately they don't say how). Filtfilt does not force your filtered data to be zero initially. But still, it will try to guess the initial conditions from the data, your result will be affected by how good that guess is. The safest strategy is to always have extra data before and after the time period of interest. Then you can clip off the data that was affected by the filter transients.

As an additional benefit, filtfilt does not have the time lag that Allan mentioned, because it filters the signal twice: forward in time and reverse in time.

A good way to see the transient, is to filter a constant signal:

>> x = ones(100,1) + 0.1*rand(100,1);       % one  second of a constant signal, noise added and sampled at 100 Hz
>> [b,a] = butter(2,6/50);                  % a 6 Hz low pass filter
>> y1 = filter(b,a,x);
>> y2 = filtfilt(b,a,x);
>> plot([x y1 y2]);
>> legend('original','filter','filtfilt')

Result:
untitled.jpg
So filtfilt does very well on this test signal. But still be careful, it may not work as well on your signal.

One thing to be aware of is that the signal is now filtered twice, and the cutoff frequency is slightly lowered: 0.802 * 6 Hz instead of 6 Hz. This has been explained in previous Biomch-L postings, and it is also explained in Winter's Biomechanics of Human Movement.

Ton van den Bogert

[Will Robertson]

Re: filtering of the inverse kinematics results

If you work in Matlab, you should try the "filtfilt" function instead of "filter".

Ton beat me to it! I was also going to give the example of biased data to show that filtfilt gives output "like you would expect":

t = linspace(0,10,200);
x = -t(end)/2 + t + 2*rand(size(t));

[b,a] = butter(5,0.5,'low');
y = filtfilt(b,a,x);

figure(1); clf;
plot(t,x,t,y)

Ton also mentioned that the cut-off frequency decreases; that's another way of saying that the attenuation at the cut-off frequency is increased
The following example shows this fairly clearly; with "filter" you get 3dB attenuation (0.707 the amplitude) at the cut-off frequency (plus some lag), and with "filtfilt" you get 6dB attenuation (0.5 the amplitude):

w = 5;
t = linspace(0,10,2000);
dt = t(2)-t(1);
Fs = 1/dt;
Fn = Fs/2;
x = sin(2*pi*w*t);

Ff = 5;
[b,a] = butter(5,Ff/Fn,'low');
y1 = filter(b,a,x);
y2 = filtfilt(b,a,x);

figure(1); clf;
plot(t,x,t,y1,t,y2)



Cheers,
Will

[Tim Wrigley]

Re: filtering of the inverse kinematics results

Nice examples guys !

One could append the following to Will's code to see additionally that by applying the correction to the 2nd order cut-off frequency of 5 Hz for the dual passes (0.802), one gets back to the equivalent of -3dB attenuation (0.707 amplitude) at that original desired cut-off frequency, with no phase shift ... ie the red dotted curve. Note that I have switched to a 2nd order design, for which 0.802 is the correct correction for dual passes.

ORDER = 2;
Ff = 5/0.802; % adjust 2nd order cutoff frequency (up) for dual pass
[b,a] = butter(ORDER,Ff/Fn,'low');
y3 = filtfilt(b,a,x);
hold on
plot(t,y3,'r:')

Tim Wrigley
University of Melbourne

[Will Robertson]

Re: filtering of the inverse kinematics results

...by applying the correction to the 2nd order cut-off frequency of 5 Hz for the dual passes (0.802), one gets back to the equivalent of -3dB attenuation (0.707 amplitude) at that original desired cut-off frequency...

Okay, now I have to ask. Do you collectively think that it's important that the cut-off frequency is defined explicitly in terms of the -3dB point? If you read in a paper: "a 3rd order, bidirectional 5Hz Butterworth filter was applied" -- what would you expect to see? Personally I interpret that as in my original code; i.e., at 5Hz there's a -6dB attenuation. But perhaps I've been wrong all this time!

Should we invent a new filter (e.g., the "Bogert filter") which is equivalent a bidirectional butterworth filter with 1/0.802 the cut-off frequency? This would help cut down some words when trying to explain ourselves in publications

Will

[Ton van den Bogert]

Re: filtering of the inverse kinematics results

I don't think we need a new name. Winter seems to define "cutoff frequency" as the -3dB point, and so does Wikipedia [1].

Winter first introduced zero-lag filtering here [2]:

"a low pass second order Butterworth filter having a cutoff at 6 Hz was designed [...] The initial filtering operation attenuates the high frequency noise. but introduces a phase lag into the output. This lag is cancelled out by again filtering the once-filtered data. but in the reverse direction of time, i.e. the last point first. This also reduces the cutoff frequency to 4.8 Hz."

In his later publications [3], he then specifies the filter like this:

"the coordinates were digitally filtered using a fourth-order zero lag low-pass Butterworth filter cutting off at 5 Hz".

To me that sounds fine, and I interpret this as having designed a 2nd order BW filter with a cutoff of 5/0.802 = 6.23 Hz, and then using it twice. Exactly like Tim Wrigley did it in his earlier posting.

1. https://en.wikipedia.org/wiki/Cutoff_frequency
2. Winter DA, Sidwall HG, Hobson DA.Measurement and reduction of noise in kinematics of locomotion. J Biomech. 1974 Mar;7(2):157-9.
3. Winter DA.Energy generation and absorption at the ankle and knee during fast, natural, and slow cadences. Clin Orthop Relat Res. 1983 May;(175):147-54.

[Tim Wrigley]

Re: filtering of the inverse kinematics results

Winter first introduced zero-lag filtering here [2]:

....

In his later publications [3], he then specifies the filter like this:

"the coordinates were digitally filtered using a fourth-order zero lag low-pass Butterworth filter cutting off at 5 Hz".

That second quoted description (even) from Winter is where I think we do get into confusion. Having been schooled according to Winter I interpret that statement exactly the same way as Ton. Perhaps Winter assumed that everyone would have read his book closely - a not unreasonable expectation. But such potentially ambiguous statements are also common in other papers.

There are several possible interpretations of such a description, and how one might do the same thing:
A. As we agree it should be: a 2nd order design, applied twice in forward and reverse direction, with fc correction applied, effectively achieving a 4th order zero phase response with fc of 5hz.
To be clear ... butter(2,(fc/0.802)/Fnyquist,'low') followed by filtfilt()*
B. Same as A but with fc correction omitted (usually inadvertantly), leading to over-filtering at fc (50% attenuation instead of 70.7%).
C. a 4th order design, applied twice, with fc correction applied, effectively (inadvertantly ?) achieving an 8th order zero phase response.
D. Same as C but with fc correction omitted.

I think we have both in the past suggested here some less ambiguous ways of describing that most commonly desired 'A' variant that Winter popularized: eg "2nd order dual pass Butterworth filter" or "double 2nd order Butterworth filter". Reference to "zero lag" is also helpful, but reference to "4th order" is likely to confuse, without further explanation.

*incidentally, one might have expected that filtfilt(), being a dual pass filter function, would correct fc for the dual pass effect itself. It doesn't. It does however address end-effects, as has been pointed out.

Tim Wrigley
University of Melbourne

[Ton van den Bogert]

Re: filtering of the inverse kinematics results

Like Tim, I actually prefer "double 2nd order Butterworth filter". Winter's wording in his 1983 paper is ambigious (even though I understand it) and possibly confusing if taken literally. His filter was NOT a 4th order Butterworth filter. A 4th order Butterworth filter would have this frequency response:

H(f) = 1 / sqrt(1 + (f/f0)^8)

The frequency response of a double 2nd order Butterworth filter is the square of the frequency response of a single BW filter:

H(f) = 1 / (1 + (f/f0)^4)

And these two are not the same. They are both 4th order, and become the same in the limit of low frequencies (f << f0) and high frequencies (f >> f0) but clearly different in the transition region, even after correcting the f0 by the 0.802 factor.

It would be nice if filtfilt compensated for the dual pass effect, but that would be difficult because it uses the filter coefficients (a,b) generated by the butter function.

Ton

[Will Robertson]

Re: filtering of the inverse kinematics results

Like Tim, I actually prefer "double 2nd order Butterworth filter". Winter's wording in his 1983 paper is ambigious (even though I understand it) and possibly confusing if taken literally. His filter was NOT a 4th order Butterworth filter.

Dear Tim and Ton,

Thank you very much for the continued discussion. I agree with you (now) that the cut-off frequency should be defined by the -3dB point of the resulting filter, and I do imagine there is a fair degree of variation in how the Butterworth filters are applied by different researchers with more or less experience in the area. Tim broke that down very nicely, and I do suspect I've inadvertently been a member of group D in the past (oops). Whether that variation makes too much of a difference would depend on the type of data being analysed — it might be interesting to look into!

To be completely unambiguous, I would imagine writing something verbose like

"a zero-lag 4th-order 4Hz low-pass filter, implemented with a double 2nd-order 4.988Hz Butterworth filter"

... or perhaps "zero-lag 4th-order 4Hz low-pass filter (as per Winter, 19xx)" works.

Cheers,
Will

[Ganesh B]

Re: filtering of the inverse kinematics results

Insightful overall discussion about filtering.
Ton mentioned"The Matlab documentation says "filtfilt minimizes start-up and ending transients by matching initial conditions" (unfortunately they don't say how)."
Matlab does it based on this paper by Gustafsson, F. "Determining the initial states in forward-backward filtering." IEEE® Transactions on Signal Processing. Vol. 44, April 1996, pp. 988–992 which is cited in Matlab documentation.