From 7646546de7b29d9e304fdba413fc6b2bd5b2ce5b Mon Sep 17 00:00:00 2001 From: tvdbogert Date: Mon, 30 Mar 2015 10:35:16 -0400 Subject: [PATCH 1/3] myfiltfilt no longer calls rtfilter.m for each sample This should be faster, hopefully it makes a significant difference. I have not tested this code (sorry) but it's no more than putting the code from rtfilter.m in the for loop in myfiltfilt.m. It would be good to try this version in Python too. --- .../octave/2d_inverse_dynamics/myfiltfilt.m | 38 +++++++++++++------ 1 file changed, 27 insertions(+), 11 deletions(-) diff --git a/gaitanalysis/octave/2d_inverse_dynamics/myfiltfilt.m b/gaitanalysis/octave/2d_inverse_dynamics/myfiltfilt.m index 01b1765..ba0dc78 100644 --- a/gaitanalysis/octave/2d_inverse_dynamics/myfiltfilt.m +++ b/gaitanalysis/octave/2d_inverse_dynamics/myfiltfilt.m @@ -9,28 +9,44 @@ % f0 (scalar) Corner frequency % % Outputs -% y (Nsamples x 1) Filtered signal -% yd (Nsamples x 1) First derivative -% ydd (Nsamples x 1) Second derivative +% y (Nsamples x Nchannels) Filtered signal +% yd (Nsamples x Nchannels) First derivative +% ydd (Nsamples x Nchannels) Second derivative - C = 0.802; % correction factor for dual pass 2nd order filter (Winter book) - y = rtfilter_batch(-flipud(t), flipud(x), f0/C); % filter backwards in time - [y, yd, ydd] = rtfilter_batch(t, flipud(y), f0/C); % filter forward in time + C = 0.802; % correction factor for dual pass 2nd order filter (Winter book) + y = filter_batch(-flipud(t), flipud(x), f0/C); % filter backwards in time + [y, yd, ydd] = filter_batch(t, flipud(y), f0/C); % filter forward in time end %=================================================================================== -function [y, yd, ydd] = rtfilter_batch(t,x,f0) +function [y, yd, ydd] = filter_batch(t,x,f0) % filters a time series and also returns derivatives % uses real-time second order Butterworth filter (rtfilter.m) + + % some constants we will need + a = (2*pi*f0)^2; + b = sqrt(2)*(2*pi*f0); + % allocate memory for the results n = size(x,1); y = zeros(size(x)); yd = zeros(size(x)); ydd = zeros(size(x)); - for i = 1:n - [y(i,:), yd(i,:), ydd(i,:)] = rtfilter(t(i),x(i,:),f0); + % Integrate the filter state equation using the midpoint Euler method with step h + % initial conditions are y=0 and yd=0 + for i = 2:n + h = t(i)-t(i-1); % time step + denom = 4 + 2*h*b + h^2*a; + A = (4 + 2*h*b - h^2*a)/denom; + B = 4*h/denom; + C = -4*h*a/denom; + D = (4 - 2*h*b - h^2*a)/denom; + E = 2*h^2*a/denom; + F = 4*h*a/denom; + y(i) = A*y(i-1) + B*yd(i-1) + E*(x(i)+x(i-1))/2; + yd(i) = C*y(i-1) + D*yd(i-1) + F*(x(i)+x(i-1))/2; + ydd = (yd(i)-yd(i-1))/h; end - end -%=================================================================================== \ No newline at end of file +%=================================================================================== From acd711d3ad6882afc4cb6282b4cb900a3c27603a Mon Sep 17 00:00:00 2001 From: tvdbogert Date: Mon, 30 Mar 2015 10:37:14 -0400 Subject: [PATCH 2/3] Deleted rtfilter.m Is no longer needed. --- .../octave/2d_inverse_dynamics/rtfilter.m | 52 ------------------- 1 file changed, 52 deletions(-) delete mode 100644 gaitanalysis/octave/2d_inverse_dynamics/rtfilter.m diff --git a/gaitanalysis/octave/2d_inverse_dynamics/rtfilter.m b/gaitanalysis/octave/2d_inverse_dynamics/rtfilter.m deleted file mode 100644 index 67964e1..0000000 --- a/gaitanalysis/octave/2d_inverse_dynamics/rtfilter.m +++ /dev/null @@ -1,52 +0,0 @@ -function [y, yd, ydd] = rtfilter(t, x, f0) - -% Real-time low pass Butterworth filter, processes one data sample at a time. -% -% Inputs: -% t time stamp of input x (s) -% x input data at time t (may be a scalar, array, or matrix) -% f0 cutoff frequency (Hz) -% -% Outputs: -% y filtered x -% yd filtered derivative of x -% ydd filtered second derivative of x - - % Internal state of the filter must be preserved between function calls - persistent state - - % If this is the first call to this function, or we have gone back in time, reset the filter - if isempty(state) || (t <= state.t) - y = x; - yd = zeros(size(x)); - ydd = zeros(size(x)); - % Otherwise, solve the state equation by Euler integration - else - % Calculate time between current and previous sample - h = t - state.t; - - % Compute coefficients of the state equation - a = (2*pi*f0)^2; - b = sqrt(2)*(2*pi*f0); - - % Integrate the filter state equation using the midpoint Euler method with step h - denom = 4 + 2*h*b + h^2*a; - A = (4 + 2*h*b - h^2*a)/denom; - B = 4*h/denom; - C = -4*h*a/denom; - D = (4 - 2*h*b - h^2*a)/denom; - E = 2*h^2*a/denom; - F = 4*h*a/denom; - y = A*state.y + B*state.yd + E*(x+state.x)/2; - yd = C*state.y + D*state.yd + F*(x+state.x)/2; - ydd = (yd-state.yd)/h; - - end - - % Store the filter state - state.t = t; - state.x = x; - state.y = y; - state.yd = yd; - -end From 124567c16cc0f9451563e9792f07ca37173b5b6a Mon Sep 17 00:00:00 2001 From: "Jason K. Moore" Date: Mon, 30 Mar 2015 11:59:39 -0700 Subject: [PATCH 3/3] Added ydd index and removed unecessary semicolon. --- gaitanalysis/octave/2d_inverse_dynamics/myfiltfilt.m | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/gaitanalysis/octave/2d_inverse_dynamics/myfiltfilt.m b/gaitanalysis/octave/2d_inverse_dynamics/myfiltfilt.m index ba0dc78..aee08b3 100644 --- a/gaitanalysis/octave/2d_inverse_dynamics/myfiltfilt.m +++ b/gaitanalysis/octave/2d_inverse_dynamics/myfiltfilt.m @@ -1,4 +1,4 @@ -function [y, yd, ydd] = myfiltfilt(t, x, f0); +function [y, yd, ydd] = myfiltfilt(t, x, f0) % performs low pass filtering and differentiation % method is the same as used in HBM but bidirectional to eliminate lag @@ -46,7 +46,7 @@ F = 4*h*a/denom; y(i) = A*y(i-1) + B*yd(i-1) + E*(x(i)+x(i-1))/2; yd(i) = C*y(i-1) + D*yd(i-1) + F*(x(i)+x(i-1))/2; - ydd = (yd(i)-yd(i-1))/h; + ydd(i) = (yd(i)-yd(i-1))/h; end end %===================================================================================