Separating left and right leg data given ground reaction force gait data?

Hi,

I have an Excel file with several columns of ground reaction force data where each column is a different test subject. My question is, how can I separate each column into left leg data and right leg data. Where stride 1 is assumed to be the right leg, therefore all odd numbered strides are right leg strides and all even numbered strides are left leg strides. I have attached the first column of my data to give you an idea of what I am looking at. Any guidance would be great. Thank you so much!

[jennifer aniston]

Hi Anthony,

Is the data that you have provided an excel sheet for, COP data collected using a force plate? I have a similar set of raw data and i am a bit lost on how to proceed with it. I will really appreciate your guidance. thanks, Jen

[Star Strider]

@Jen —

If you want to open an new Question, I’ll help you as I can when we get Anthony’s data analysis and coding sorted. Anthony’s Question has all the information needed to Answer it (including the data), so consider a similar presentation with your Question.

[jennifer aniston]

Thanks Star Strider,

I am new on this group so finding my way around. If you find the time please take a look at my question. Link: https://in.mathworks.com/matlabcentral/answers/305313-advice-on-processing-raw-gait-data

Thanks!

[Star Strider]

My pleasure!

Answer 1

[Star Strider]

You need to use the Signal Processing Toolbox findpeaks function. It retuens the values of the peaks and their locations. I only had the indices to work with, but you can easily use the indices to convert to time if you need to.

The code:

[d,s,r] = xlsread('Anthony 1 Column Gait Data.xlsx');

[pks,locs] = findpeaks(d, 'MinPeakDistance',500, 'MinPeakHeight',1000);

figure(1)
plot(d)
hold on
plot(locs(1:2:end), pks(1:2:end), '^r', 'MarkerFaceColor','r')
plot(locs(2:2:end), pks(2:2:end), '^g', 'MarkerFaceColor','g')
hold off
grid
axis([0  5000    ylim])
legend('Force Plate Data', 'Right Foot', 'Left Foot', 'Location','E')
xlabel('Index')
ylabel('Force')

I plotted a portion of the data to show the details here. The code works for the entire data set.

The plot:

[Anthony]

Hi Star Strider,

Thank you very much for the response! I think you are on right on track here. This is already very informative and useful for me, but I was wondering if there was a way to filter out the data for each entire stride. I copied your picture and marked lines to explain what I mean.I basically want to separate all of the data that is coming from the right leg from all the data that is coming from the left leg. Thank you so much again!

[Star Strider]

My pleasure!

The easiest way to do that is to do a second call to findpeaks, but with the negative of your signal. Those peaks will be the ‘troughs’, and with essentially the same code (you will only need to modify the 'MinPeakHeight' value), you should be able to identify the low points separating the strides. Remember that the ‘pks’ value for the second findpeaks call (to get the ‘troughs’) will be the negative of the actual value returned, so consider that in your plot. The ‘locs’ value will not be affected.

The code:

[d,s,r] = xlsread('Anthony 1 Column Gait Data.xlsx');

[pks,locs] = findpeaks( d, 'MinPeakDistance',500, 'MinPeakHeight',1000);
[trf,tloc] = findpeaks(-d, 'MinPeakDistance',500, 'MinPeakHeight',-100);

figure(1)
plot(d)
hold on
plot(locs(1:2:end), pks(1:2:end), '^r', 'MarkerFaceColor','r')
plot(locs(2:2:end), pks(2:2:end), '^g', 'MarkerFaceColor','g')
plot([tloc tloc]', repmat(ylim,size(tloc,1),1)', '-k', 'LineWidth',1.5)
hold off
grid
axis([0  5000    ylim])
legend('Force Plate Data', 'Right Foot', 'Left Foot', 'Location','E')
xlabel('Index')
ylabel('Force')

The plot:

[Anthony]

Hi Star Strider,

Thank you again! I'm glad you understood my question! For one last question, now that we have the left and right leg data separated, how can I put that into an Excel file in separate columns? Thank you so much again. You have no idea how helpful this is!

[Star Strider]

My pleasure.

Among other things (one being a retired Board Certified Internal Medicine Physician), I’m also a M.S. biomedical engineer, and have done some work (long ago) as part of a gait analysis project.

This is how I would separate the left and right foot force records. There are no gaps in the records:

tloc = [1; tloc; length(d)];                                % Augment ‘tloc’ Vector
rfoot = [];                                                 % Initialise Vector
lfoot = [];                                                 % Initialise Vector
for k1 = 1:2:length(tloc)-2
    rfoot = [rfoot; d(tloc(k1):tloc(k1+1))];
    lfoot = [lfoot; d(tloc(k1+1):tloc(k1+2))];
end

If you want to separate the individual records, the easiest way is to put NaN values between them:

tloc = [1; tloc; length(d)];                                % Augment ‘tloc’ Vector
rfoot = [];                                                 % Initialise Vector
lfoot = [];                                                 % Initialise Vector
for k1 = 1:2:length(tloc)-2
    rfoot = [rfoot; d(tloc(k1):tloc(k1+1)); NaN];
    lfoot = [lfoot; d(tloc(k1+1):tloc(k1+2)); NaN];
end

Then, just test for the NaN values.

They are now individual vectors, so you can write them to your Excel file as individual columns. (Note that the columns are not equal length.)

If you want the corresponding times for each, you can do that similarly, in the same loop. The easiest way is to create two new vector assignments, one for each, including the initialization assignment before the loop. If your time vector is ‘t’, substitute ‘t’ for ‘d’, and add a ‘t’ to the beginning or end of the variable names here.

[Anthony]

Thank you again very much. I have been able to write the 2 data sets into two separate columns. Regarding the time vector, if my sampling rate is 1200 Hz (therefore my time vector is 0:1/1200:end), how would I match these time points with the data I have from each foot? I'm sorry for all of the questions but you have been incredibly helpful! I think you alluded to that in your previous answer but I am not exactly sure how to incorporate that into the code. Thank you so much again!

[Star Strider]

My pleasure.

The loop changes to:

Ts = 1/1200;                                                % Sampling Time Interval
t = linspace(0, 1, length(d))'*Ts;                          % Time Vector
tloc = [1; tloc; length(d)];                                % Augment ‘tloc’ Vector
rfoot = [];                                                 % Initialise Vector
lfoot = [];                                                 % Initialise Vector
rfoot_t = [];                                               % Initialise Vector
lfoot_t = [];                                               % Initialise Vector
for k1 = 1:2:length(tloc)-2
    rfoot = [rfoot; d(tloc(k1):tloc(k1+1)); NaN];
    lfoot = [lfoot; d(tloc(k1+1):tloc(k1+2)); NaN];
    rfoot_t = [rfoot_t; t(tloc(k1):tloc(k1+1)); NaN];
    lfoot_t = [lfoot_t; t(tloc(k1+1):tloc(k1+2)); NaN];
end

Note the transpose (') in the linspace call to force a column vector.

You can now also plot them individually with respect to time using those vectors:

figure(2)
plot(rfoot_t, rfoot, 'r', 'LineWidth',1)  
hold on
plot(lfoot_t, lfoot, 'g', 'LineWidth',1)
hold off
grid
axis([0  t(5000)    ylim])
legend('Right Foot', 'Left Foot', 'Location','E')
xlabel('Time (s)')
ylabel('Force')

If my Answer is solving your problems, please Accept it.

[Anthony]

I am so sorry about that. I meant to accept it previously. Thank you again for the help! I have one more question if you can help I would be eternally grateful. Would I be able to then find the time of each peak for each stride and then export them to Excel in addition to the separated data? My end goal is to have the time of each peak for each stride (with the corresponding foot) and the time of each initial contact. If you could guide me in the right direction for that, it would be incredible! Thank you so much again. If there is anything else I could possibly do to help, let me know please!

[Star Strider]

No worries, and thanks for the additional vote! Yours is a fun question for me because of my interest in the subject.

This is a complete rewrite of the code, using time rather than indices in all computations. (You still have the index code if you need it. In this version, the independent variable is in units of time, here seconds, throughout the code.)

Specifically, ‘pkt’ (peaks times) and ‘trt’ (troughs times) are now in units of seconds. (Note that the variables in the plots have changed.) You are of course free to change the variable names to something more appropriate in the context of your research. (I like short variable names because the probability of my mistyping them is less.)

The Code:

[d,s,r] = xlsread('Anthony 1 Column Gait Data.xlsx');

Ts = 1/1200;                                                % Sampling Time Interval
t = linspace(0, 1, length(d))'*Ts;                          % Time Vector

[pks,pkt] = findpeaks( d, t, 'MinPeakDistance',t(500), 'MinPeakHeight',1000);
[trf,trt] = findpeaks(-d, t, 'MinPeakDistance',t(500), 'MinPeakHeight',-100);

figure(1)
plot(t, d)
hold on
plot(pkt(1:2:end), pks(1:2:end), '^r', 'MarkerFaceColor','r')
plot(pkt(2:2:end), pks(2:2:end), '^g', 'MarkerFaceColor','g')
plot([trt trt]', repmat(ylim,size(trt,1),1)', '-k', 'LineWidth',1.5)
hold off
grid
axis([0  t(5000)    ylim])
legend('Force Plate Data', 'Right Foot', 'Left Foot', 'Location','E')
xlabel('Time (s)')
ylabel('Force')

% Ts = 1/1200;                                                % Sampling Time Interval
% t = linspace(0, 1, length(d))'*Ts;                          % Time Vector
t2i = round(trt*length(d)/Ts);                              % Convert ‘t’ To Indices
tloc = [1; t2i; length(d)];                                 % Augment ‘tloc’ Vector
rfoot = [];                                                 % Initialise Vector
lfoot = [];                                                 % Initialise Vector
rfoot_t = [];                                               % Initialise Vector
lfoot_t = [];                                               % Initialise Vector
for k1 = 1:2:length(trt)-2
    rfoot = [rfoot; d(tloc(k1):tloc(k1+1)); NaN];
    lfoot = [lfoot; d(tloc(k1+1):tloc(k1+2)); NaN];
    rfoot_t = [rfoot_t; t(tloc(k1):tloc(k1+1)); NaN];
    lfoot_t = [lfoot_t; t(tloc(k1+1):tloc(k1+2)); NaN];
end

figure(2)
plot(rfoot_t, rfoot, 'r', 'LineWidth',1)  
hold on
plot(lfoot_t, lfoot, 'g', 'LineWidth',1)
hold off
grid
axis([0  t(5000)    ylim])
legend('Right Foot', 'Left Foot', 'Location','E')
xlabel('Time (s)')
ylabel('Force')

The plots are the same, so I won’t reproduce them here, except for the independent variable now being in seconds, in both plots.

[Anthony]

Thank you again so much! It's amazing that you are able to understand exactly what I am asking. I apologize if I my questions seem redundant. I do not have much experience with Matlab. I was wondering if it would be possible to find the times the correspond with the initial contact of each stride? For instance if I let 50 N of force be the threshold for which I assumed that was initial contact, how would I isolate each of these points with their corresponding time. I have attached a picture to show what I mean to hopefully clarify.

[Star Strider]

As always, my pleasure!

That turned out to be easier than I thought it would be, once I remembered ‘zci’, a little utility routine I wrote a while back. It calculates zero-crossings, so all I needed to do to detect values of about 50 N was to subtract 50 from ‘d’ to create the zero-crossings. You may need to start the subscripts from 1 rather than 2 in other of your records, but that should be easy to determine when you plot your data. My ‘zci’ function returns the (z)ero (c)rossing (i)ndices of the dependent variable, so indexing into the data such as I did with the ‘pos’ matrix is necessary. Then plotting it or recording it is straightforward.

The Code:

zci = @(v) find(v(:).*circshift(v(:), [-1 0]) <= 0);                    % Returns Approximate Zero-Crossing Indices Of Argument Vector

zx = zci(d - 50);
pos = [t(zx(2:2:end))  d(zx(2:2:end))];

figure(3)
plot(rfoot_t, rfoot, 'r', 'LineWidth',1)  
hold on
plot(lfoot_t, lfoot, 'g', 'LineWidth',1)
plot(pos(:,1), pos(:,2), 'bp', 'MarkerSize',9, 'MarkerFaceColor','b')
hold off
grid
axis([0  t(5000)    ylim])
legend('Right Foot', 'Left Foot', 'Initial Contact (50 N)', 'Location','E')
xlabel('Time (s)')
ylabel('Force')

The Plot:

[Anthony]

Thank you very much again! Does that code separate the right initial contact values from the left initial contact values? I'm sorry, that code is a little over my head. I can't seem to pick out where it separates the odd-numbered initial contacts from the even-numbered initial contacts. I already have my program writing the peak values of left and right strides and the entire data from each foot as well. This would never have been possible without you! Thank you again for everything! I think that will be my last question so I won't bother you anymore!!! :)

[Star Strider]

My pleasure!

This code separates them into right and left initial contact times and forces. Again, you may need to tweak the initial indices by ±1 for other data.

The Code:

zci = @(v) find(v(:).*circshift(v(:), [-1 0]) <= 0);                    % Returns Approximate Zero-Crossing Indices Of Argument Vector

zx = zci(d - 50);
init_cont = [t(zx(2:2:end))  d(zx(2:2:end))];
init_cont_rf = init_cont(2:2:end,:);
init_cont_lf = init_cont(1:2:end,:);

figure(3)
plot(rfoot_t, rfoot, 'r', 'LineWidth',1)  
hold on
plot(lfoot_t, lfoot, 'g', 'LineWidth',1)
plot(init_cont_rf(:,1), init_cont_rf(:,2), 'bp', 'MarkerSize',9, 'MarkerFaceColor','y')
plot(init_cont_lf(:,1), init_cont_lf(:,2), 'mp', 'MarkerSize',9, 'MarkerFaceColor','b')
hold off
grid
axis([0  t(5000)    ylim])
legend('Right Foot', 'Left Foot', 'Right Foot Initial Contact (50 N)', 'Left Foot Initial Contact (50 N)', 'Location','E')
xlabel('Time (s)')
ylabel('Force')

The Plot:

[Anthony]

Awesome! Thank you so much again! I really can't thank you enough for all of this help! You have been so much help! I will keep you updated on how everything goes. Hopefully I won't have anymore questions, but I hope you are there to help if I do! Thanks again!!

[Star Strider]

As always, my pleasure!

I certainly intend to be here.

This was an interesting and enjoyable challenge. I never thought to use my little ‘zci’ function this way.

I wish you well in your research!

[Anthony]

I'm glad you enjoyed helping! Thank you for the nice wishes!

[Star Strider]

I like medical and biomedical engineering problems! I really do have fun with them!

I noticed that the left leg was generating less force that the right leg. The left leg may be slightly shorter than the right leg. The swing times may also be slightly different if that hypothesis is correct, since the ‘pendulum period’ would be different. (You can also assess this with EMGs from the hip flexors and extensors. The shorter leg would have greater EMG activity to keep it synchronised with the longer leg.)

Measure leg lengths. Leg lengths are rarely equal, and if you can account for the difference in force by a difference in leg lengths, that could be an important covariate in your eventual statistical analysis. By the same token, if you can equalise the length (and force) with an appropriately-designed orthotic in the shoe of the shorter leg, you might have an additional paper out of your current study. (This would be a paired study, so consider the appropriate statistical analysis.) An orthopaedic surgeon or physical therapist can assist you with these measurements if they’re not part of your background.

If you really want to get precise, get a lower-segment (lumbar, lumbosacral, pelvic and distal) MRI of each of your participants. (CT scans generate too much ionising radiation.) You can not only measure precise leg lengths, but also muscle volumes, and search for specific pathology (arthritis, atrophy, peripheral neuropathy, ad infinitum). Avoid gadolinium contrast if possible. It may be centrally neurotoxic.

Just a thought to make your project more interesting. Your committee or faculty advisor will (if they are good) ask about the force discrepancy (and hypothesise a leg-length discrepancy as have I), so be prepared to explain it with data. If you can offer a paired study with intervention, you get serious points for having thought about it, investigated it, and have an explanation and possibly and intervention for it. Never overlook details! They may amount to nothing, but they may be the next Nobel Prize. (Hyperbole here, but you get the idea.)

Your other friend in all this (if you have not met it yet) is PubMed. You will learn to rely on each other.

I look forward to following your research. You’re off to a good start with the way your stated your Question and presented your data. That made it very easy for me.

[Anthony]

I cannot thank you enough for your help and words of encouragement! I think you bring up a tremendous amount of great idea that can be used in my research!

I love your enthusiasm and willingness to help. I would like to give you some background on what I may be trying to accomplish. For one, I think I will be including the calculation of slope between the initial contact point and peak point (I may need some help on that one!)

But in terms of outside of the Matlab program, I will be focusing a lot on fatigue and its corresponding effects on neuromuscular control and stability. Looking at gait data before and after a fatigue protocol will give a great look in to the effects fatigue can have on an individual's biomechanics. Particularly, their lower extremity control. I have looked a lot into PubMed! It has been a very good friend indeed thus far.

I love many of your suggestions! I definitely will try to implement them into what I am doing. You seem to have a brilliant mind and your words have been inspiring. The length of leg comparison brings up a very intriguing possibility. I think I will have to give a good look into that. I'm actually fascinated by that possibility now that I am thinking about it. I can't wait to get this program totally finished (all I need is the slope calculation now!), and apply it to my studies. You have given me a whole new level of excitement! Thank you again!

-Anthony

[Star Strider]

As always, my pleasure!

The slope between the initial contact and the peak force (or load), is relatively straightforward. Divide the difference between the peak force and the initial contact force by the time difference between them. I’m not certain this is what you want. If it is, add it to the end of the last code I posted:

f_rate = (pks(1:size(init_cont,1)) - init_cont(:,2)) ./ (init_cont(:,1) - pkt(1:size(init_cont,1)));    % Force Rate (N/sec)

figure(4)
plot(init_cont(:,1), f_rate)
grid
xlabel('Time (sec)')
ylabel('Force Rate (N/sec)')

figure(5)
plot(t, d*5E+5)
hold on
plot(init_cont(:,1), f_rate)
hold off
grid
axis([0  t(20000)    ylim])
xlabel('Time (sec)')
ylabel('Force Rate (N/sec)')

I called it ‘f_rate’ for ‘force rate’, but that’s for my convenience. It produces an interesting plot that reveals information I didn’t realise was in the data. (It uses the initial contact time, ‘init_cont(:,1)’, as the time base for the plot. Both vectors have to be equal length, so you may have to tweak the code for different records.

Fatigue has many components, the depletion of ATP and oxygen (myoglobin oxygen saturation) and accumulation of lactate in the muscle, and to a lesser extent, depletion of neurotransmitters, particularly acetylcholine, or neurotransmitter insensitivity (consider the effects of succinylcholine). Some aspects of fatigue may be accessible with surface electromyography. It’s been a while since I studied muscle fatigue, so you will want to see what the literature has on that.

You happened to hit one of my research interest ‘eigenfrequencies’, so I characteristically resonated.

I very much appreciate your compliments! I am always happy to help, within the scope of my expertise. Again, my pleasure!

[Anthony]

Thank you so much again! That is exactly what I wanted! Is it possible to sort those force rates for each stride into their correlating strides and feet? Meaning can I match each force rate with each corresponding stride and therefore separate them by which foot is producing that stride? I may have worded that a little unclear. I apologize for that.

Regarding my research, I am indeed looking into a component (or multiple components) that constitutes fatigue.I am leaning towards certain EMG techniques and looking at VO2 max change. Previous studies have shown these two to be reliable in fatigue studies. I will once again keep you updated! I am really looking forward to getting my research going! Thank you again!

[Star Strider]

As always, my pleasure!

I added another two variables ‘r_pks’ and ‘l_pks’, tweaked the code to accommodate them, and created different force rates for the right and left foot as ‘rf_frate’ and ‘lf_frate’ respectively. I’m at a loss to explain the off- scale-high final value for the right foot force rate. It may just be a data anomaly, since the rest are consistent.

I find the results of the force rate analysis fascinating. I had no idea that the force rates would display such interesting 0.2 millisecond periodicity. I look forward to your analysis of it, because I have no idea how to explain it.

The Code:

[d,s,r] = xlsread('Anthony 1 Column Gait Data.xlsx');

Ts = 1/1200;                                                % Sampling Time Interval
t = linspace(0, 1, length(d))'*Ts;                          % Time Vector

[pks,pkt] = findpeaks( d, t, 'MinPeakDistance',t(500), 'MinPeakHeight',1000);
[trf,trt] = findpeaks(-d, t, 'MinPeakDistance',t(500), 'MinPeakHeight',-100);

r_pks = [pkt(1:2:end), pks(1:2:end)];                       % Right Foot Times & Peaks
l_pks = [pkt(2:2:end), pks(2:2:end)];                       % Left Foot Times & Peaks

figure(1)
plot(t, d)
hold on
plot(r_pks(:,2), r_pks(:,2), '^r', 'MarkerFaceColor','r')
plot(l_pks(:,2), l_pks(:,2), '^g', 'MarkerFaceColor','g')
plot([trt trt]', repmat(ylim,size(trt,1),1)', '-k', 'LineWidth',1.5)
hold off
grid
axis([0  t(5000)    ylim])
legend('Force Plate Data', 'Right Foot', 'Left Foot', 'Location','E')
xlabel('Time (s)')
ylabel('Force')

% Ts = 1/1200;                                                % Sampling Time Interval
% t = linspace(0, 1, length(d))'*Ts;                          % Time Vector
t2i = round(trt*length(d)/Ts);                              % Convert ‘t’ To Indices
tloc = [1; t2i; length(d)];                                 % Augment ‘tloc’ Vector
rfoot = [];                                                 % Initialise Vector
lfoot = [];                                                 % Initialise Vector
rfoot_t = [];                                               % Initialise Vector
lfoot_t = [];                                               % Initialise Vector
for k1 = 1:2:length(trt)-2
    rfoot = [rfoot; d(tloc(k1):tloc(k1+1)); NaN];
    lfoot = [lfoot; d(tloc(k1+1):tloc(k1+2)); NaN];
    rfoot_t = [rfoot_t; t(tloc(k1):tloc(k1+1)); NaN];
    lfoot_t = [lfoot_t; t(tloc(k1+1):tloc(k1+2)); NaN];
end

figure(2)
plot(rfoot_t, rfoot, 'r', 'LineWidth',1)  
hold on
plot(lfoot_t, lfoot, 'g', 'LineWidth',1)
hold off
grid
axis([0  t(5000)    ylim])
legend('Right Foot', 'Left Foot', 'Location','E')
xlabel('Time (s)')
ylabel('Force')

zci = @(v) find(v(:).*circshift(v(:), [-1 0]) <= 0);        % Returns Approximate Zero-Crossing Indices Of Argument Vector

zx = zci(d - 50);
init_cont = [t(zx(2:2:end))  d(zx(2:2:end))];
init_cont_rf = init_cont(2:2:end,:);
init_cont_lf = init_cont(1:2:end,:);

figure(3)
plot(rfoot_t, rfoot, 'r', 'LineWidth',1)  
hold on
plot(lfoot_t, lfoot, 'g', 'LineWidth',1)
plot(init_cont_rf(:,1), init_cont_rf(:,2), 'bp', 'MarkerSize',9, 'MarkerFaceColor','y')
plot(init_cont_lf(:,1), init_cont_lf(:,2), 'mp', 'MarkerSize',9, 'MarkerFaceColor','b')
hold off
grid
axis([0  t(5000)    ylim])
legend('Right Foot', 'Left Foot', 'Right Foot Initial Contact (50 N)', 'Left Foot Initial Contact (50 N)', 'Location','E')
xlabel('Time (s)')
ylabel('Force')

f_rate = @(fi,fp,ti,tp) (fp - fi) ./ (tp - ti);             % Force Rate Function (N/sec)
rf_frate = f_rate(init_cont_rf(:,2), r_pks(2:end,2), init_cont_rf(:,1), r_pks(2:end,1));
lf_frate = f_rate(init_cont_lf(:,2), l_pks(:,2), init_cont_lf(:,1), l_pks(:,1));

figure(4)
plot(t, d*4E+5)
hold on
plot(init_cont_rf(:,1), rf_frate, '-r')
plot(init_cont_lf(:,1), lf_frate, '-g')
hold off
grid
axis([0  t(20000)    ylim])
xlabel('Time (sec)')
ylabel('Force Rate (N/sec)')
legend('Right Foot', 'Left Foot', 'Location', 'NE')

figure(5)
plot(t, d*4E+5)
hold on
plot(init_cont_rf(:,1), rf_frate, '-r', 'LineWidth',1.2)
plot(init_cont_lf(:,1), lf_frate, '-g', 'LineWidth',1.2)
hold off
grid
% axis([0  t(20000)    ylim])
xlabel('Time (sec)')
ylabel('Force Rate (N/sec)')
legend('Force Data', 'Right Foot', 'Left Foot', 'Location', 'NE')

The Plot:

[Anthony]

That is perfect! Once again, I am very glad you understood my question. Thank you again. I find the force rates very interesting as well. I will have to apply this code to other data I have. Im curious to see how force rates are affected under certain conditions. The periodic nature of the force rates seems to be a very interesting find. I'm very excited about moving forward with this code! Thank you again very much!

[Star Strider]

My pleasure.

You may be making a significant original contribution with the force rate data. I would consult with a statistician at your university to devise appropriate statistical approaches to analysing all your data before you proceed much further. I’m sure you’ve already developed a statistical design before you began, but the force rate analysis and perhaps some other results may require you to devise new analysis methods for them.

I’m interested in following your research, and I would appreciate a PDF of your paper when you publish it, as well.

[Anthony]

Thank you again for your interest and help! I am looking into different ways of interpreting the force rate data. I will be sure to keep you in the loop every step of the way! Thank you very much!

[Star Strider]

I’d definitely like to know where that strange peroidicity comes from. There are some functions that will help analyse the left-right phase differences, specifically finddelay, alignsignals, xcorr, xcov, and if you want to do system identification, invfreqz and the System Identification Toolbox. Most biomedical signals are sufficiently nonlinear as to preclude the use of linear techniques, but your gait data may be well-suited to it.

I’m not certain how much information a Fourier transform would give you. It’s certainly worth doing the analysis to see if there’s anything else hidden in your data. The R2015a documentation for fft is the easiest to understand, so I recommend it, and particularly the code between the first two (top two) plot figures.

As always, my pleasure! I’m very much interested in what you find.

[Anthony]

Thank you again for your suggestions and help! Regarding the use of linear techniques, a lot of research has been published recently showing that in order to measure stability (especially in biological systems such as human movement), nonlinear techniques should be used. This is due to the idea of an "optimal range of variability" which promotes the idea that human locomotion is actually most stable when it has some degree of variability. There are arguments that when there is too small of a degree of variability, the movement becomes too robust and robotic and is not able to adapt to perturbation and other external forces. Some linear techniques have shown to be insufficient in measuring the oscillations in human movement in terms of their stability so I think we will be focusing on nonlinear techniques such as Approximate Entropy. I think that will help us tremendously in analyzing gait data.

It has been a pleasure going back and forth with you! I'm really enjoying this conversation and feel like I am getting a TON out of this. Thank you so much again!

[Star Strider]

As always, my pleasure. Your literature review is interesting. When I last worked with gait, the ‘inverted pendulum’ model ruled, and it could be modeled using linear techniques.

The approximate entropy approach is intriguing. I remember studying it in connection with pituitary hormone secretory dynamics (part of my checkered past is having trained as an endocrinologist). I’ve not reviewed it in a while and I’ve not worked with it, so I’ll do some reading.

A brief File Exchange search for approximate entropy brought up a couple contributions that may work with your data.

[Anthony]

If you are really interested in the use of ApEn in measuring dynamic stability, I can point you in the right direction of a few articles. I really appreciate the File Exchange search. That is another great contribution from you! Thank you again for the help!

[Star Strider]

My pleasure!

I printed out the original 1990 PNAS Pincus paper to read later. I would appreciate any references you can provide, since ‘ApEn’ and other ways to characterise nonlinear and specifically chaotic processes are always useful. It would help if they’re free, or in Science (I’m an AAAS member).

[Anthony]

No problem at all! I came across this article as one of the free articles available on PubMed: https://www.ncbi.nlm.nih.gov/pubmed/22052133 There are a few more I wish I could link you to but I can't seem to find free versions of them anywhere. I will work on finding more! Thanks again!

[Star Strider]

Thank you!

There are some linked-to PubMed Central articles on that page that I’ll look through as well. (My interest was in looking at gait disturbances in diabetic myopathy and peripheral neuropathy.)

[Anthony]

No problem at all! That is a very interesting research topic! I actually have read a paper that may interest you. It is a great paper written by Jonathan B. Dingwell which was published in 2000. It is entitled "Nonlinear time series analysis of normal and pathological human ". I cannot find a free source for the paper but I here is the link in case you want to take a look at the abstract: http://scitation.aip.org/content/aip/journal/chaos/10/4/10.1063/1.1324008

In this paper they use Lyapunov exponents instead of ApEn, however it is still a great read for this topic of discussion! Let me know what you think!

[Star Strider]

I need to review nonlinear dynamics and chaos. I have the books, and I’ve worked with Lyapunov exponents. I’ve just not worked with them recently.

I’ll see if I can sneak into the library of the nearby branch campus of the state university and download the PDF of the paper. If I can’t email it to me (this has worked elsewhere), I may be able to download it to a SD card.

I didn’t consider chaotic dynamics back when I was doing my research, since we were interested in detection and classification of the extent of impairment rather than the dynamics. That was a clinical study, and we had to keep it simple.

[Anthony]

I wish I could be of more help. I would love to give you a copy of the paper if I was able to. Using nonlinear dynamic models to look at human motion seems to be a relatively new idea so I am very excited to get going on my research using these ideas.

Let me know if you are able to get a copy of that paper or anything else you might think is interesting! I wish I could be more help! I apologize

[Star Strider]

I understand. I need to visit the university library. I’m not formally associated with the university, so I can’t do this from home.

[Anthony]

Ok let me know if it works out. If not, I will try to figure something out on my end!

[Star Strider]

That works. This is a busier week than usual for me, so it may have to wait until next week. It’s definitely an area of my interest, and nonlinear dynamics is something I need to review (a long time since my one course in it).

[Anthony]

Ok just keep me updated! Thanks again!

[Star Strider]

I will. My pleasure.

[Anthony]

Hi Star Strider,

I hope you are still checking this! I have another question I think you may be able to help with. I would like to know how to take the double derivative of the initial data column to produce an acceleration data column. I would then like to find the peaks of each stride in that acceleration column. I believe each stride will have a maximum acceleration at initial contact which will give a more clear picture of when initial contact occurs. Let me know if you can help! Thank you so much again!

-Anthony

[Star Strider]

Hi Anthony,

As always, my pleasure!

Questions I Answer with the most recent activity pop to the top of my profile Answers page. I check it every morning. I don’t check it during the day unless I’ve been away for a while. I scroll through Answers occasionally otherwise. I look at all new activity and catch most of them.

There are two functions that will give you the second derivative. The first is using gradient twice, and using del2 once, multiplying it by 4. Either works.

Your data have a leading ‘spike’ (that I interpret as a heel strike), followed by what seems to be a transition to the stance phase, and then toe-off. The ‘spike’ could give you serious problems with either function.

Note that derivatives amplify noise, and although your data seem reasonably noise free, you’ll soon find out where the noise is when you take the derivatives! There are ways to deal with that, from designing a discrete filter if the noise is band-limited, or using a Savitzky-Golay filter (the sgolayfilt function) if it’s broadband. These can be either straightforward to implement, or a real challenge, depending on the signal.

To design a discrete filter, you first have to take the Fourier transform to determine the frequency content. Use the R2015a documentation for fft (link) as your guide, particularly the code between the first (top) two plot figures. If you have band-limited noise, the easiest way to design a filter is to use the designfilt function. I prefer to use IIR filters, but that’s my personal preference. Many times FIR filters are the only design that will work effectively. My IIR filter design procedure is here: How to design a lowpass filter for ocean wave data in Matlab? It’s probably worth a look even if you don’t use it. Your signals don’t have any baseline drift, so you can probably use a low-pass filter rather than a bandpass filter. These are easier to design.

If your data turn out to have significant band-limited noise, give the discrete filter design a go. I have a fair amount of experience with discrete filter design and implementation, so I’m here if you have questions or problems.

If you need to use sgolayfilt, the implementation is much more heuristic (and I can provide less effective help). You may need to experiment with that until you get the result you want.

[Anthony]

Hi Star Strider,

Very nice to know! Thank you again for the response!

I went ahead and used the 'gradient' function and I think I see what you mean with the noise being amplified. I have attached the picture of the plot I achieved by plotting the second derivative vs. time. Are the peaks in that plot true peaks? Meaning, if all I want is the peaks and their corresponding locations, do I necessarily need to design a filter? I am mainly asking because I have no experience designing filters. For example, can I just extract the second derivative (acceleration data) and use the 'findpeaks' function like you did in the beginning of the code? Is there a simple way to do this? And how would I go about doing this in the most simple way?

Thank you again for the help!

-Anthony

[Star Strider]

Hi Anthony,

You have a serious noise problem with the second derivatives. They are ‘true peaks’ only in the most permissive sense. They are, in fact, likely garbage.

The first approach is to do the fft to see if the noise is band-limited. It’s the end of my day here (that being 2230 LCL (UCT-6)) so I’ll work on this tomorrow.

My approach will be:

1. Fourier transform.

2. If the noise is band-limited, design & implement low-pass filter (or other appropriate filter).

3. Experiment with the 4*del2 function on the filtered data.

I’ll do the first two and post the code, You can reproduce the third to see if it meets your needs. (I can’t determine that.)

Note that what you want to do may not be possible in practice. You can only do so much with real-world signals. Reality (sometimes) sucks.

My pleasure!

[Anthony]

Hi Star Strider,

Ok thank you very much! I look forward to seeing your approach. Even if my goal is not entirely possible, I look forward to seeing what you can do. I appreciate the help as always! Look forward to hearing from you!

-Anthony

[Star Strider]

Hi Anthony,

This is about as good as it gets. Further lowering the cutoff frequency begins to slur the details of the original heel-strike ‘spike’, that then compromises the derivative function and resolution, so I left it as is.

If you experiment with the passband and stopband frequencies, the plot in figure(7) will tell you if you designed the filter correctly. This codes for a Butterworth IIR filter, so the Bode plot should have a flat top, with no irregularities in the passband, and a smooth rolloff.

Append this code to the end of the previous code I posted.

The Code:

Ts = 1/1200;                                                % Sampling Time Interval
t = linspace(0, 1, length(d))'*Ts;                          % Time Vector

L = size(d,1);
Ts = mean(diff(t));                                         % Sampling Interval 
Fs = 1/Ts;                                                  % Sampling Frequency 
Fn = Fs/2;                                                  % Nyquist Frequency 

FTd = fft(d)/L;                                             % Fourier Transform
Fv = linspace(0, 1, fix(L/2)+1)*Fn;                         % Frequency Vector
Iv = 1:length(Fv);                                          % Index Vector

figure(6)
semilogy(Fv, 1+abs(FTd(Iv))*2)
grid
axis([0  2.5E+6    ylim])
xlabel('Frequency (Hz)')

Wp = [2.5E+6]/Fn;                                           % Passband Frequency
Ws = [3.5E+6]/Fn;                                           % Stopband Frequency
Rp = 20;                                                    % Passband Ripple
Rs = 40;                                                    % Stopband Ripple
[n,Wn] = buttord(Wp,Ws,Rp,Rs);                              % Filter Order
[b,a] = butter(n,Wn);                                       % Filter Transfer Function Coefficients
[sos,g] = tf2sos(b,a);                                      % Second-Order Section For Stability
df = filtfilt(sos,g,d);                                     % Phase-Neutral Filtering

figure(7)
freqz(sos)

figure(8)
subplot(3,1,1)
plot(t, d, 'LineWidth',0.6)
axis([0  5E-5    -100  2000])
ylabel('Force (N)')
title('Original Signal')
grid
subplot(3,1,2)
plot(t, df, '-r', 'LineWidth',0.6)
axis([0  5E-5    -100  2000])
ylabel('Force (N)')
title('Lowpass-Filtered Signal')
grid
subplot(3,1,3)
plot(t, del2(df)*4, '-g', 'LineWidth',0.6)
axis([0  5E-5    ylim])
xlabel('Time (s)')
ylabel('\Delta Force (d^2N/dt^2)')
title('Second Derivative Of Filtered Signal')
grid

The Plot:

[Anthony]

Hi Star Strider,

This is fantastic! This gives such a clear picture as to where heel-strike occurs. Thank you so much again! I appreciate the help so much! I'm pretty sure this gives me exactly what I want but I will let you know if I need anything else. Thank you so much again!

-Anthony

[Star Strider]

Hi Anthony,

As always, my pleasure.

I apologise for the delay, but this has been a busy week. I was catching up, and since filter design and implementation can take a while, I wanted to be sure I had the uninterrupted time to devote to it.

[Anthony]

Hi Star Strider,

There is absolutely no need to apologize. I greatly appreciate the help! Thanks again!

Best,

Anthony

[Star Strider]

My pleasure.

New problem here is that my primary Win 8.1 MATLAB computer is down and it’ll be about three weeks before I’ll have it back (power management chip failed again). Neither Firefox nor IE will allow me to download and save .mat files on this Win 10 machine (so I don’t have access to yours), another reason to hate Win 10 (as though there aren’t enough already). They want to save them as .txt files. I’ve done everything I can to troubleshoot this, including checking the file associations, trying to get both browsers to accept .mat files specifically, and everything else, all to no avail. I asked TMW for help (since it affects MATLAB Answers for me), but it’s too early to have heard back.

[Anthony]

Hi Star Strider,

I am very sorry you are having computer troubles. I actually hate Windows 10 as well and have kept my personal computer on 8.1. I wish you luck with getting everything fixed! I am assuming several people are having similar problems so hopefully you will have an answer soon! Thank you again!

-Anthony

[Star Strider]

Hi Anthony,

Thank you. It turns out that your .xlsx file isn’t affected, only .mat files. I was able to download your file to my Win 10 machine without problems. I can convert inappropriately-saved .txt files to .mat files with the MATLAB movefiles function (something to remember if you have to change file extensions), so while inconvenient, it works.

I didn’t have a choice with my utility computer — it came with Win 10, as did my upgraded MATLAB-&-Gaming machine that I’ve not yet migrated to.

So we’re good for your data.

AS always, my pleasure!

[Anthony]

Hi Star Strider,

Nice to hear it won't have any effect on my data. I still feel for you regarding Windows 10. We actually had a Windows 10 computer in our lab crash today. Even the troubleshooting in Windows 10 is awful. So, I feel your pain.

Thank you again!

-Anthony

[Star Strider]

Hi Anthony,

I completely agree. Win 7 was a terrific OS, and I still have one seven-year-old laptop running it. (That one came with Vista.) Everything since was a ‘downgrade’ IMO. Someone won a rather significant award in a lawsuit against Micro$oft for ‘upgrading’ her computer to Win 10 without her permission. I believe the suit also made M$ roll it back to the previous OS. Justice prevails!

My pleasure!

[Anthony]

Hi Star Strider,

Wow I didn't know that! I'm not surprised someone took legal action though. I would have been incredibly frustrated if Windows 10 was forced upon me!

Best,

Anthony

[Star Strider]

Hi Anthony,

As the successful plaintiff in at least one action (others pending) with a terrific legal team, justice can definitely prevail. I’d like to know more about the Micro$oft Win 10 litigation (not that I’d understand the technical details of it, since I’m not a lawyer), but to learn how to prevent such ‘involuntary servitude’ in the future. The important lesson to learn from this is that we are not powerless.

[Anthony]

Hi Star Strider,

That is a very important lesson. Nice to see that these companies can't get away with absolutely everything!!!

-Anthony

[Star Strider]

Quite definitely!

[Anthony]

Hi Star Strider,

Hope all is going well! I had another question regarding the code. Now that I have identified the initial contact using the peaks of the second derivative, I would like to once again calculate the rate change of the force between the initial contact time and the peak force time. Is this done differently now because one data point is from the second derivative and one is from the normal data. Let me know if you need any other information! Thank you again!

Anthony

[Star Strider]

Hi Anthony,

Rasonably so. I hope all is going well with you as well.

I had to go back and remember what I was doing.

I don’t see that it makes a difference. My code here:

zx = zci(d - 50);
init_cont = [t(zx(2:2:end))  d(zx(2:2:end))];
init_cont_rf = init_cont(2:2:end,:);
init_cont_lf = init_cont(1:2:end,:);

relies on the indices (here ‘zx’) to create ‘init_cont’, so if you have your data references as indices, there should be no problem. Just substitute them into ‘zx’.

If you’ve defined the ‘init_cont’ matrix differently, it should not make any difference so long as the time column and initial contact force column are appropriate to your design, with alternating left and right foot times and contact forces in each column.

So long as ‘init_cont’ has the appropriate row size to be compatible with the force rate calculations:

f_rate = @(fi,fp,ti,tp) (fp - fi) ./ (tp - ti);             % Force Rate Function (N/sec)
rf_frate = f_rate(init_cont_rf(:,2), r_pks(2:end,2), init_cont_rf(:,1), r_pks(2:end,1));
lf_frate = f_rate(init_cont_lf(:,2), l_pks(:,2), init_cont_lf(:,1), l_pks(:,1));

there should be no problem.

[Anthony]

Hi Star Strider,

Thank you very much for the fast response! I seem to be getting an error when running this code. This is the error:

Index exceeds matrix dimensions.

Error in practice (line 81) rf_frate = f_rate(init_cont_rf(:,2), r_pks(2:end,2), init_cont_rf(:,1), r_pks(2:end,1));

Is there any simple fix to this? Thank you again!

-Anthony

[Star Strider]

Hi Anthony,

My pleasure.

The ‘r_pks’ array shouldn’t change, so I doubt it’s the problem.

Check the size of ‘init_cont_rf’. It should be a (Nx2) matrix. It has to have ‘t’ (time) as the first column and initial force ‘d’ at the same time as the second column.

The row size (dimension 1 or ‘N’ here) of ‘r_pks’ and ‘init_cont_rf’ must be the same (although if they didn’t that would throw a ‘matrix dimensions must agree’ error, not the error you’re seeing).

[Anthony]

Hi Star Strider,

Sorry for the delayed response. I had a few exams last week that took up most of my time. I think I figured out what I was doing wrong. Thank you again so much for all of your help!

-Anthony

[Star Strider]

My pleasure!

No worries. I would like to know what the problem was and how you solved it.

[Anthony]

Hi Star Strider,

For some reason I didn't have the time data as my first column of the init_cont_rf matrix. It was basically what you said to check. I mustve mixed something up somewhere in the code. But once I double checked the size of init_cont_rf like you told me to, it was clear what my problem was! Thanks again!

-Anthony

[Star Strider]

My pleasure!

I was curious because I wanted to be certain that I didn’t overlook something, and to be sure my code didn’t glitch.

I’m still interested in your project, so I’ll be here as long as necessary. It seems to me to be publishable, and if so, I’d appreciate a PDF of the published paper.

[Anthony]

No problem at all! It might be a while before I get to that point but I will keep you in the loop as I progress! Thank you again for all of the help.

-Anthony

[Star Strider]

As always, my pleasure.

I’m curious about your research, especially since I’ve become so involved in it. The paper will let me read about all the details.

[Anthony]

Absolutely! I completely understand your interest!

[Felipe Mattioni Maturana]

Hello guys! That is amazing, I am currently running into a similar problem and this helped me quite a lot. I have one more question, though. How would you guys normalize the force and emg data from 0 to 100%?

Thank you in advance,

Felipe

[Star Strider]

I would first subtract the minimum to create a new record, then divide by the maximum of those data.

[Christina Kaltenbach]

That is what I thought. I was wondering if Matlab had a specific function for that. One more thing: can I use your code above to automatically analyze the gait cycles (left foot + right foot) in a Trial of approximately 10 steps? I mean, I have these data from different speeds on the treadmill, and each Trial contains around 10 steps for each foot. I wanted to automatically retrieve each gait cycle (right+left, using the 50-N threshold as above), and normalize them (0-100%). Thank you in advance!

[Felipe Mattioni Maturana]

Dear Star Strider,

I have been working on my force data using the codes you posted here. Everything worked perfectly so far. The only problem I have now is that I am analyzing data from an old Treadmill and it is quite noisy. The 50-N threshold is not working on this data set because between each strike there are usually two bumps on the data and it is often above 50 N. Would you have any suggestion on how I can overcome this issue? I am posting a figure and an example of my data (please, note that commas instead of dots might appear as decimal separators as I am using a German computer) to help you understand what is going on. Any help would be very much appreciated.

Best,

Felipe

[Star Strider]

Post this as a new Question.

Answer 2

[KSSV]

[data,txt,raw] = xlsread('1 Column Gait Data.xlsx');
% Right / odd numbered 
right = data(1:2:end) ;
% left / even numbered 
left = data(2:2:end) ;