![](/uploads/1/2/5/7/125744261/573347697.jpg)
I ended up duplicating one of the legs and attached it again the the world frame and the torso, effectively closing the 'loop'. A contact force between the left foot and the floor. A contact force between the right foot and the floor. But I would recommend you look at this library on the MATLAB Central File. Hacked spotify apk. Linefewermarkers - line with controlled amount of markers and correct legend. 1: offset all markers such that one marker on first max of y(x). Triangle (left).
72 Comments
on 30 Sep 2016
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!
on 30 Sep 2016
The easiest way to do that is to do asecond call tofindpeaks, but with thenegative 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 thesecondfindpeaks call (to get the ‘troughs’) will be thenegative of the actual value returned, so consider that in your plot. The‘locs’ value will not be affected.
[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)
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)
grid
legend('Force Plate Data', 'Right Foot', 'Left Foot', 'Location','E')
ylabel('Force')
on 30 Sep 2016
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!
on 30 Sep 2016
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
lfoot = []; % Initialise Vector
rfoot = [rfoot; d(tloc(k1):tloc(k1+1))];
end
If you want to separate the individual records, the easiest way is to putNaN values between them:
tloc = [1; tloc; length(d)]; % Augment ‘tloc’ Vector
lfoot = []; % Initialise Vector
rfoot = [rfoot; d(tloc(k1):tloc(k1+1)); NaN];
end
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.
on 30 Sep 2016
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!
on 30 Sep 2016
The loop changes to:
t = linspace(0, 1, length(d))'*Ts; % Time Vector
tloc = [1; tloc; length(d)]; % Augment ‘tloc’ Vector
lfoot = []; % Initialise Vector
lfoot_t = []; % Initialise Vector
rfoot = [rfoot; d(tloc(k1):tloc(k1+1)); NaN];
rfoot_t = [rfoot_t; t(tloc(k1):tloc(k1+1)); NaN];
lfoot_t = [lfoot_t; t(tloc(k1+1):tloc(k1+2)); NaN];
Note the transpose(') in thelinspace call to force a column vector.
You can now also plot them individually with respect to time using those vectors:
plot(rfoot_t, rfoot, 'r', 'LineWidth',1)
plot(lfoot_t, lfoot, 'g', 'LineWidth',1)
grid
legend('Right Foot', 'Left Foot', 'Location','E')
ylabel('Force')
If my Answer is solving your problems, pleaseAccept it.
on 30 Sep 2016
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!
on 30 Sep 2016
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.)
[d,s,r] = xlsread('Anthony 1 Column Gait Data.xlsx');
Ts = 1/1200; % Sampling Time Interval
[pks,pkt] = findpeaks( d, t, 'MinPeakDistance',t(500), 'MinPeakHeight',1000);
[trf,trt] = findpeaks(-d, t, 'MinPeakDistance',t(500), 'MinPeakHeight',-100);
figure(1)
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)
grid
legend('Force Plate Data', 'Right Foot', 'Left Foot', 'Location','E')
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
lfoot = []; % Initialise Vector
lfoot_t = []; % Initialise Vector
rfoot = [rfoot; d(tloc(k1):tloc(k1+1)); NaN];
rfoot_t = [rfoot_t; t(tloc(k1):tloc(k1+1)); NaN];
lfoot_t = [lfoot_t; t(tloc(k1+1):tloc(k1+2)); NaN];
plot(rfoot_t, rfoot, 'r', 'LineWidth',1)
plot(lfoot_t, lfoot, 'g', 'LineWidth',1)
grid
legend('Right Foot', 'Left Foot', 'Location','E')
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.
on 30 Sep 2016
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.
on 30 Sep 2016
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 about50 N was to subtract50 from‘d’ to create the zero-crossings. You may need to start the subscripts from1 rather than2 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 thedependent variable, so indexing into the data such as I did with the‘pos’ matrix is necessary. Then plotting it or recording it is straightforward.
zci = @(v) find(v(:).*circshift(v(:), [-1 0]) <= 0); % Returns Approximate Zero-Crossing Indices Of Argument Vector
zx = zci(d - 50);
plot(rfoot_t, rfoot, 'r', 'LineWidth',1)
plot(lfoot_t, lfoot, 'g', 'LineWidth',1)
plot(pos(:,1), pos(:,2), 'bp', 'MarkerSize',9, 'MarkerFaceColor','b')
grid
legend('Right Foot', 'Left Foot', 'Initial Contact (50 N)', 'Location','E')
ylabel('Force')
on 1 Oct 2016
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!!! :)
on 1 Oct 2016
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.
zci = @(v) find(v(:).*circshift(v(:), [-1 0]) <= 0); % Returns Approximate Zero-Crossing Indices Of Argument Vector
zx = zci(d - 50);
init_cont_rf = init_cont(2:2:end,:);
plot(rfoot_t, rfoot, 'r', 'LineWidth',1)
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')
grid
legend('Right Foot', 'Left Foot', 'Right Foot Initial Contact (50 N)', 'Left Foot Initial Contact (50 N)', 'Location','E')
ylabel('Force')
on 1 Oct 2016
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!!
on 1 Oct 2016
I certainlyintend to be here.
This was an interesting and enjoyable challenge. I never thought to use my little‘zci’ function this way.
on 1 Oct 2016
I'm glad you enjoyed helping! Thank you for the nice wishes!
on 1 Oct 2016
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 apaired 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 advisorwill (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) isPubMed. 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.
on 3 Oct 2016
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!
on 3 Oct 2016
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 theend 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)
grid
ylabel('Force Rate (N/sec)')
figure(5)
hold on
hold off
axis([0 t(20000) ylim])
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!
on 3 Oct 2016
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!
on 3 Oct 2016
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 interesting0.2 millisecond periodicity. I look forward to your analysis of it, because I have no idea how to explain it.
[d,s,r] = xlsread('Anthony 1 Column Gait Data.xlsx');
Ts = 1/1200; % Sampling Time Interval
[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)
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)
grid
legend('Force Plate Data', 'Right Foot', 'Left Foot', 'Location','E')
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
lfoot = []; % Initialise Vector
lfoot_t = []; % Initialise Vector
rfoot = [rfoot; d(tloc(k1):tloc(k1+1)); NaN];
rfoot_t = [rfoot_t; t(tloc(k1):tloc(k1+1)); NaN];
lfoot_t = [lfoot_t; t(tloc(k1+1):tloc(k1+2)); NaN];
plot(rfoot_t, rfoot, 'r', 'LineWidth',1)
plot(lfoot_t, lfoot, 'g', 'LineWidth',1)
grid
legend('Right Foot', 'Left Foot', 'Location','E')
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_rf = init_cont(2:2:end,:);
plot(rfoot_t, rfoot, 'r', 'LineWidth',1)
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')
grid
legend('Right Foot', 'Left Foot', 'Right Foot Initial Contact (50 N)', 'Left Foot Initial Contact (50 N)', 'Location','E')
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)
hold on
plot(init_cont_lf(:,1), lf_frate, '-g')
grid
xlabel('Time (sec)')
legend('Right Foot', 'Left Foot', 'Location', 'NE')
figure(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)
grid
xlabel('Time (sec)')
legend('Force Data', 'Right Foot', 'Left Foot', 'Location', 'NE')
on 3 Oct 2016
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!
on 3 Oct 2016
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.
on 3 Oct 2016
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!
on 3 Oct 2016
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, specificallyfinddelay,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. TheR2015a documentation forfft 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.
on 3 Oct 2016
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!
on 3 Oct 2016
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 forapproximate entropy brought up a couple contributions that may work with your data.
on 3 Oct 2016
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!
on 3 Oct 2016
I printed out the original1990PNAS 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 inScience (I’m an AAAS member).
on 3 Oct 2016
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!
on 3 Oct 2016
There are some linked-toPubMed 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.)
on 4 Oct 2016
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!
on 4 Oct 2016
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.
on 4 Oct 2016
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
on 4 Oct 2016
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.
on 4 Oct 2016
Ok let me know if it works out. If not, I will try to figure something out on my end!
on 4 Oct 2016
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).
on 4 Oct 2016
on 4 Oct 2016
on 7 Oct 2016
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!
on 8 Oct 2016
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 usinggradient twice, and usingdel2 once, multiplying it by4. 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 (thesgolayfilt 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 theR2015a documentation forfft (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 thedesignfilt function. I prefer to useIIR filters, but that’s my personal preference. Many timesFIR filters are the only design that will work effectively. MyIIR 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 usesgolayfilt, 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.
on 8 Oct 2016
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?
-Anthony
on 8 Oct 2016
You have aserious 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 thefft to see if the noise is band-limited. It’s the end of my day here (that being2230 LCL (UCT-6)) so I’ll work on this tomorrow.
1. Fourier transform.
2. If the noise is band-limited, design & implement low-pass filter (or other appropriate filter).
3. Experiment with the4*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.
![For loop left right leg matlab function For loop left right leg matlab function](/uploads/1/2/5/7/125744261/733249091.jpg)
on 8 Oct 2016
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!
on 9 Oct 2016
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 infigure(7) will tell you if you designed the filter correctly. This codes for a ButterworthIIR filter, so the Bode plot should have a flat top, with no irregularities in the passband, and a smooth rolloff.
Append this code to theend of the previous code I posted.
t = linspace(0, 1, length(d))'*Ts; % Time Vector
L = size(d,1);
Fs = 1/Ts; % Sampling Frequency
Fv = linspace(0, 1, fix(L/2)+1)*Fn; % Frequency Vector
semilogy(Fv, 1+abs(FTd(Iv))*2)
axis([0 2.5E+6 ylim])
Ws = [3.5E+6]/Fn; % Stopband Frequency
Rs = 40; % Stopband Ripple
[b,a] = butter(n,Wn); % Filter Transfer Function Coefficients
[sos,g] = tf2sos(b,a); % Second-Order Section For Stability
freqz(sos)
figure(8)
plot(t, d, 'LineWidth',0.6)
ylabel('Force (N)')
grid
plot(t, df, '-r', 'LineWidth',0.6)
ylabel('Force (N)')
grid
plot(t, del2(df)*4, '-g', 'LineWidth',0.6)
xlabel('Time (s)')
title('Second Derivative Of Filtered Signal')
The Plot:
on 10 Oct 2016
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!
on 10 Oct 2016
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.
on 10 Oct 2016
There is absolutely no need to apologize. I greatly appreciate the help! Thanks again!
Anthony
on 10 Oct 2016
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.
on 10 Oct 2016
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!
on 10 Oct 2016
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 MATLABmovefiles 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.
AS always, my pleasure!
on 11 Oct 2016
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.
-Anthony
on 11 Oct 2016
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!
on 11 Oct 2016
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!
Anthony
on 11 Oct 2016
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 thatwe are not powerless.
on 12 Oct 2016
That is a very important lesson. Nice to see that these companies can't get away with absolutely everything!!!
on 12 Oct 2016
on 19 Oct 2016
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!
on 20 Oct 2016
Rasonably so. I hope all is going well with you as well.
I don’t see that it makes a difference. My code here:
init_cont = [t(zx(2:2:end)) d(zx(2:2:end))];
init_cont_lf = init_cont(1:2:end,:);
relies on theindices (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));
on 20 Oct 2016
Thank you very much for the fast response! I seem to be getting an error when running this code. This is the error:
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));
-Anthony
on 20 Oct 2016
My pleasure.
The‘r_pks’ array shouldn’t change, so I doubt it’s the problem.
Check thesize of‘init_cont_rf’. It should be a(Nx2) matrix. It has to have‘t’ (time) as thefirst column and initial force‘d’ at the same time as thesecond column.
![For loop left right leg matlab function For loop left right leg matlab function](https://bcsecure01-a.akamaihd.net/6/62009828001/201909/1295/62009828001_6083246113001_6022520889001-vs.jpg?pubId=62009828001&videoId=6022520889001)
The row size (dimension1 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).
on 25 Oct 2016
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!
on 25 Oct 2016
No worries. Iwould like to know what the problem was and how you solved it.
on 28 Oct 2016
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!
on 28 Oct 2016
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.
on 31 Oct 2016
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.
on 31 Oct 2016
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.
on 31 Oct 2016
Absolutely! I completely understand your interest!
on 29 Jan 2018
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%?
Felipe
on 29 Jan 2018
I would first subtract the minimum to create a new record, then divide by the maximum of those data.
on 29 Jan 2018
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!
on 4 Mar 2018
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.
Felipe
on 4 Mar 2018
Sign in to comment.
![](/uploads/1/2/5/7/125744261/573347697.jpg)