-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy path181003_v5_KNN.m
1992 lines (1679 loc) · 88.9 KB
/
181003_v5_KNN.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
% code for OpenKin datalogger's and moticon insole's data processing--
%% Code lines which may need some modifications for different tests
% if the same code is being used for different walking/running tests
% then the following lines may need some modifications--
% Line no: 30 ->> name or (path+name) of the datalogger file
% Line no: 43 ->> index of datalogger required test data.
% Line no: 66-67 ->> INS frame defination
% Line no: 108-109 ->> recommended
% Line no: 217-233 ->> as per device's orientation on human back
% Line no: 247 ->> name or (path+name) of the moticon's data file
% Line no: 262-287 ->> left and right insoles sync points
% Line no: 590-605 ->> to see if the step segmentation using insoles is correct or not
% NOTE: If any chnages in line no's (add/delete) are to be made before line
% no. 605 then the above list should be modified.
%% Load datalogger's data
% This section parses the complete dataset collected by INS datalogger and
% the foot IMU data (foot IMU is present only in few tests).
% The INS data is logged at the rate of 400Hz. Each row in the data file
% contains the data of all sensors at one time instance.
% Therefore, the sensor data in two consecutive rows are separated by the
% time of 0.0025 sec.
walkdata = cell_array_parser_vn200_mod('2018-10-13-07-58-14.txt', 1);
%The purpose of foot IMU data was to achieve better sync with moticon's
%foot insole data but it doesn't help much so the foot imu data is ignored
%in later sections.
%% INS-data to matrices
% 'ind' is range, of data indices of INS datalogger, in which the walking/
% running test is perfrmed. If the data file contains the logged data from
% ealier time then the differnece between those two times can be very large
% so it better to exclude the older data from processing
ind = 1:size(walkdata,2); % range of indeces of test performed.
ts = arrayfun(@(c) c.time.ros_stamp, walkdata(ind)); % ROS timestamp
% imu data seq. 0->65535->0->65535
imuseq = arrayfun(@(c) c.time.imu_seq, walkdata(ind));
times = arrayfun(@(c) c.time.gpstime, walkdata(ind)); %GPS ticks
plot(times);
%ID of row, 'ins' if row contains the data from datalogger or '/base_imu'
%if it is from foot imu.
ids = arrayfun(@(c) c.imu.id, walkdata(ind));
insind = find(contains(ids(1,:),{'ins'})); % indices of INS data
%indices of foot IMU data (in this test foot imu was not present)
mtiind = find(contains(ids(1,:),{'/base_imu'}));
% Acceleration data from data logger at human back and foot IMU.
acc = cell2mat( arrayfun(@(c) c.imu.acc(:), walkdata(ind), 'UniformOutput', 0) );
acc_ins=acc(:,insind); % Acceleration data from data logger at human back.
% This test was performed using new INS orientation in device.
%Frame definations.
% INS frame axis in new device (when device is on human's back), y:
% forward, z upward, x: towards right hand.
% Anatomical frame: for human movement, X in forward, Y in human's right
% hand's direction and Z in vertically upwards direction.
% Rotation quaternion (in XYZ order)
quat = cell2mat( arrayfun(@(c) c.imu.quat(:), walkdata(ind), 'UniformOutput', 0) );
for k=1:length(quat)
%quaternion to euler angles (by using matlab's aerospace toolbox)
[yaw(k), pitch(k), roll(k)] = quat2angle( [quat(4,k) quat(1,k) quat(2,k) quat(3,k)]);
% quaternion to euler angles (without using matlab's toolbox)
% [yaw(k), pitch(k), roll(k)] = quat2euler([quat1(k) quat2(k) quat3(k) quat4(k)]);
end
% Rotation angles obtained from INS datalogger (around sensor frame axis)
roll_ins = roll(insind);
pitch_ins = pitch(insind);
yaw_ins = yaw(insind);
% unwrap yaw angle to avoid jumps at pi and -pi.
yaw_unwrap=unwrap(yaw_ins);
% extract yaw angle oscillations
yaw_osc= highpass(yaw_unwrap, 0.005,'StopbandAttenuation',30,'Steepness',0.7);
v_d= cell2mat(arrayfun(@(c) c.gps.v_d(:), walkdata(ind), 'UniformOutput', 0));% vertical velocity in downwards direction
v_e= cell2mat(arrayfun(@(c) c.gps.v_e(:), walkdata(ind), 'UniformOutput', 0));
v_n= cell2mat(arrayfun(@(c) c.gps.v_n(:), walkdata(ind), 'UniformOutput', 0));
V_d = -v_d(insind);%vertical velocity in upwards direction from INS/GPS (sensor fusion)
v_e = v_e(insind);%velocity in geographical east direction from INS/GPS (sensor fusion)
v_n = v_n(insind);%velocity in geographical north direction from INSGPS (sensor fusion)
ang_vel= cell2mat(arrayfun(@(c) c.imu.ang_vel(:), walkdata(ind), 'UniformOutput', 0));
ang_vel_ins = ang_vel(:,insind);%angular velocity from IMU in INS datalogger
gpsvel_n= cell2mat(arrayfun(@(c) c.gps.raw_v_n(:), walkdata(ind), 'UniformOutput', 0));
gpsvel_e= cell2mat(arrayfun(@(c) c.gps.raw_v_e(:), walkdata(ind), 'UniformOutput', 0));
gpsvel_n= gpsvel_n(:,insind);%velocity in geographical north direction by only GPS
gpsvel_e= gpsvel_e(:,insind);%velocity in geographical east direction by only GPS
%After this section is executed, it is better to delete variable 'walkdata'
%from matlab workspace for fast processing in sections below.
%% Plot yaw angle?
plot(rad2deg(yaw_ins))
hold on
plot(rad2deg(yaw_unwrap))
plot(rad2deg(yaw_osc))
%% To Generate 'custom_eul2rotm.m' file, if file already available then it is not needed to execute this section
% rotation matrix for transforming acceleration vectors in
% three-dimensional space from frame-A to frame-B.
% let's, first consider orientation X,Y,Z as axis of frame-B and
% x,y,z as of frame-A.
syms psi theta phi R_i1_B R_i2_i1 R_A_i2 R_A_B R_B_A
% psi (yaw) around Z axis, theta (pitch) around Y axis, and phi (roll)
% around X axis
% R_i1_b is rotation matrix from frame-B to first intermediate frame.
% frame-B to first intermediate frame (right hand rotation around Z axis)
% X',Y',Z
R_i1_B = [cos(psi) sin(psi) 0 ;...
-sin(psi) cos(psi) 0 ;...
0 0 1];
% first intermediate frame to second intermediate frame (right hand
% rotation around Y axis)
%i.e. X'',Y',Z'
R_i2_i1 = [cos(theta) 0 -sin(theta);...
0 1 0 ;...
sin(theta) 0 cos(theta)];
% second intermediate frame to a-frame (right hand rotation around x axis)
% X'',Y'',Z'' or x,y,z
R_A_i2 = [1 0 0 ;...
0 cos(phi) sin(phi) ;...
0 -sin(phi) cos(phi)];
R_A_B = R_A_i2 * R_i2_i1 * R_i1_B; %frame-B to frame-A
% inverse of Rotation matrix is same as tranpose matrix
R_B_A = R_A_B.'; % frame-A to frame-B
%save symbolic function to the matlab function in .m file
matlabFunction(R_B_A,'File','custom_eul2rotm','Optimize',false,'Vars',[psi, theta, phi]);
% so, this function is good to use with yaw, pitch, roll to convert
% accelerations from INS to geographical frame meaning that INS frame
% is frame-A and geographical frame is frame-B.
%% Transform accelerations from INS (sensor) frame to geographical frame.
acc_geo=zeros(3,length(insind));
for k=1:length(insind)
%accelerations in geographical frame (n,e,downwards)
%The order of the input angles should be the yaw, pitch and roll.
acc_geo(:,k) = custom_eul2rotm(yaw(k), pitch_ins(k), roll_ins(k)) * acc_ins(:,k);
%k
end
%% Transform accelerations from INS (sensor) frame to geographical frame (using robotics toolbox).
% %alternate, by using inbuilt robotics toolbox.
% acc_geo1=zeros(3,length(insind));
% for k=1:length(insind)
% acc_geo1(:,k) = eul2rotm([yaw(k), pitch_ins(k), roll_ins(k)],'ZYX')*acc_ins(:,k);
% %k
% end
%% Ground track and velocity calculations
GrndTrack= zeros(1,length(insind));
v_lon= zeros(1,length(insind));
v_lat= zeros(1,length(insind));
V_gps_lon= zeros(1,length(insind));
for k=1:length(insind)
% angle between north and forward direction. *not very accurate when
% person is standing or slow walking (<1.5 m/s).
GrndTrack(k)=rad2deg(atan2(v_e(k),v_n(k)));
% INS/GPS (sensor fusion) horizontal speed (forward direction of movement)
v_lon(k) = sqrt(v_n(k)^2 + v_e(k)^2);
% Lateral velocity (in direction of right hand) constraint V_y = 0 (straight ine running?)
v_lat(k) = v_e(k)*cosd(GrndTrack(k)) - v_n(k)*sind(GrndTrack(k));
% GPS horizontal speed (only GPS), forward direction of movement
V_gps_lon(k) = sqrt(gpsvel_n(k)^2 + gpsvel_e(k)^2);
% k
end
%% acceleration calculations in anatomical frame
acc_anatomical = zeros(3,length(insind));
for k=1:length(insind)
% acc. in x direction of anatomical frame (i.e. forward or londitudinal)
acc_anatomical(1,k) = acc_geo(1,k)*cosd(GrndTrack(k)) + acc_geo(2,k)*sind(GrndTrack(k));
% acceleration in lateral direction.
acc_anatomical(2,k) = acc_geo(2,k)*cosd(GrndTrack(k)) - acc_geo(1,k)*sind(GrndTrack(k));
%pure acceleration (without gravity bias) in vertically upward dir.
acc_anatomical(3,k) = -acc_geo(3,k)-9.81;
end
%% plot anatomical accelerations?
plot(acc_anatomical(3,:)) %upward direction
hold on
plot(acc_anatomical(1,:)) % forward or londitudinal direction
%% Euler with reference to Anatomical frame
% Since human body is also oscillating as INS datalogger, when walking/
% running. (assuming device is tightly attached to the body).
% % So, body-frame (anatomical frame), x axis in forward and y axis in
% human right hand's direction and z in upward direction. This frame
% oscillates when person walks/runs around X,Y,Z by means of roll_anatomical,
% roll_anatomical, yaw_anatomical.
%These equation depends on the orientation of the sensor frame with the
%anatomical (body) frame
roll_anatomical= pitch_ins;
pitch_anatomical= roll_ins;
yaw_anatomical= yaw_ins;
% angular velocities are yet to be tranformed/ calculated in geo and anatomical frame.
%% Atomic/ TIA (or UTC) time conversion
%Currently,
% UTC Epoch - GPS Epoch= (315964800 - 18) % 18, leap second correction
% Also, GPS & UTC Epoch (seconds) are in order of 10 digits, until next 2 centuries.
% variable 'times'/10e8 -> GPS time in sec (or GPS epoch)
% times2 = arrayfun(@(x) x/10e8 + 315964800 - 18, times);% gps ticks to UTC time conversion.
times2 = arrayfun(@(x) x/10e8 + 315964800 + 19, times);% gps ticks to atomic time (TAI) conversion.
times2ins=times2(insind);
%% load moticon insoles data
raw = importdata('dharmendra_181003.txt', '\t', 2);
raw = raw.data;
% fill missing values with approximation
raw = fillmissing(raw, 'spline', 1,'EndValues','nearest');
%% plot to sync time between moticon and datalogger
figure;
plot(raw(:,1), raw(:,2)+raw(:,3)+raw(:,13)+raw(:,14));
hold on;
plot(raw(:,1), raw(:,2+19)+raw(:,3+19)+raw(:,13+19)+raw(:,14+19));
figure;
plot(times2ins(1,:), acc_ins(3,:));
%% time sync with two points for each foot at the time of jump landing
% sync points for each foot
%lsync = [left_peak_foot_pressure_time coresponding_peak_vertical_acc_TIA_time; ...]
%These points needs to be changed if GPS/UTC time is used for sync.
%TIA is better since it is not changed with leap sec. addition as UTC.
lsync = [107.58 37+1.538553314535636e+09; ...
533.00 37+1.538553739865062e+09; ...
663.09 37+1.538553869958137e+09; ...
951.44 37+1.538554158354394e+09];
rsync = [107.54 37+1.538553314535636e+09; ...
532.98 37+1.538553739865062e+09; ...
663.06 37+1.538553869958137e+09; ...
951.44 37+1.538554158354394e+09];
ltimes = zeros(length(raw),1);
rtimes = zeros(length(raw),1);
% calculate timestamps for each moticon sample with piecewise linear
% approximation
if (length(lsync) > 1)
i = 2;
for k = 1:length(raw)
rtime = raw(k,1);
while (rtime > lsync(i) && length(lsync) > i)
i = i+1;
end
ltimes(k) = ((rtime - lsync(i-1, 1)) * (lsync(i, 2) - lsync(i-1, 2)) / (lsync(i, 1) - lsync(i-1, 1))) + lsync(i-1, 2);
end
end
if (length(rsync) > 1)
i = 2;
for k = 1:length(raw)
rtime = raw(k,1);
while (rtime > rsync(i) && length(rsync) > i)
i = i+1;
end
rtimes(k) = ((rtime - rsync(i-1, 1)) * (rsync(i, 2) - rsync(i-1, 2)) / (rsync(i, 1) - rsync(i-1, 1))) + rsync(i-1, 2);
end
end
%% various plots for moticon steps (-> uncomment)
% %% plot steps detected by Moticon software
%
% % figure;
% % plot(times2(1,:), -acc(3,:))
% % hold on;
% %
% % for i = 1:length(leftsteps)
% % stime = ((leftsteps(i,2) - lp1rawtime) * (lp2realtime - lp1realtime) / (lp2rawtime - lp1rawtime)) + lp1realtime;
% % plot([stime stime], [-20, 40], 'r');
% % end
% %
% % for i = 1:length(leftsteps)
% % stime = ((leftsteps(i,3) - lp1rawtime) * (lp2realtime - lp1realtime) / (lp2rawtime - lp1rawtime)) + lp1realtime;
% % plot([stime stime], [-20, 40], 'b');
% % end
% %
% % for i = 1:length(rightsteps)
% % stime = ((rightsteps(i,2) - rp1rawtime) * (rp2realtime - rp1realtime) / (rp2rawtime - rp1rawtime)) + rp1realtime;
% % plot([stime stime], [-20, 40], 'm');
% % end
% %
% % for i = 1:length(rightsteps)
% % stime = ((rightsteps(i,3) - rp1rawtime) * (rp2realtime - rp1realtime) / (rp2rawtime - rp1rawtime)) + rp1realtime;
% % plot([stime stime], [-20, 40], 'c');
% % end
%
% %% plot total forces
%
% % figure;
% % plot(times2(1,:), -acc(3,:), 'b')
% % hold on;
% %
% % % left total col. 18
% % plot(ltimes, raw(:,18)/100, 'r')
% %
% % % right total col. 37
% % plot(rtimes, raw(:,37)/100, 'g')
%
% %% add moticon steps to figure
% % for i = 1:length(leftsteps)
% % stime = ((leftsteps(i,2) - lp1rawtime) * (lp2realtime - lp1realtime) / (lp2rawtime - lp1rawtime)) + lp1realtime;
% % plot([stime stime], [-20, 40], 'r');
% % end
% %
% % for i = 1:length(leftsteps)
% % stime = ((leftsteps(i,3) - lp1rawtime) * (lp2realtime - lp1realtime) / (lp2rawtime - lp1rawtime)) + lp1realtime;
% % plot([stime stime], [-20, 40], 'b');
% % end
% %
% % for i = 1:length(rightsteps)
% % stime = ((rightsteps(i,2) - rp1rawtime) * (rp2realtime - rp1realtime) / (rp2rawtime - rp1rawtime)) + rp1realtime;
% % plot([stime stime], [-20, 40], 'm');
% % end
% %
% % for i = 1:length(rightsteps)
% % stime = ((rightsteps(i,3) - rp1rawtime) * (rp2realtime - rp1realtime) / (rp2rawtime - rp1rawtime)) + rp1realtime;
% % plot([stime stime], [-20, 40], 'c');
% % end
%
% %% plot individual sensor values with moticon steps
%
% % figure;
% % plot(times2(1,:), -acc(3,:), 'b')
% % hold on;
% %
% % % left toes
% % plot(ltimes, raw(:,2), 'r')
% % plot(ltimes, raw(:,3), 'r--')
% %
% % % left heels
% % plot(ltimes, raw(:,13), 'g')
% % plot(ltimes, raw(:,14), 'g--')
% %
% % for i = 1:length(leftsteps)
% % stime = ((leftsteps(i,2) - lp1rawtime) * (lp2realtime - lp1realtime) / (lp2rawtime - lp1rawtime)) + lp1realtime;
% % plot([stime stime], [-20, 40], 'r');
% % end
% %
% % for i = 1:length(leftsteps)
% % stime = ((leftsteps(i,3) - lp1rawtime) * (lp2realtime - lp1realtime) / (lp2rawtime - lp1rawtime)) + lp1realtime;
% % plot([stime stime], [-20, 40], 'b');
% % end
% %
% %
% % % right toes
% % plot(rtimes, raw(:,2+19)-8, 'r')
% % plot(rtimes, raw(:,3+19)-8, 'r--')
% %
% % % right heels
% % plot(rtimes, raw(:,13+19)-8, 'g')
% % plot(rtimes, raw(:,14+19)-8, 'g--')
% %
% % for i = 1:length(rightsteps)
% % stime = ((rightsteps(i,2) - rp1rawtime) * (rp2realtime - rp1realtime) / (rp2rawtime - rp1rawtime)) + rp1realtime;
% % plot([stime stime], [-20, 40], 'm');
% % end
% %
% % for i = 1:length(rightsteps)
% % stime = ((rightsteps(i,3) - rp1rawtime) * (rp2realtime - rp1realtime) / (rp2rawtime - rp1rawtime)) + rp1realtime;
% % plot([stime stime], [-20, 40], 'c');
% % end
%% Detect ground contact time from raw pressures
% to find out indices of all heelstrikes and toeoffs by
% left and right foot in moticon's data
% heel sum pressures
lheel = raw(:,13)+raw(:,14);
rheel = raw(:,13+19)+raw(:,14+19);
% % toe sum pressures
% ltoe = raw(:,2)+raw(:,3);
% rtoe = raw(:,2+19)+raw(:,3+19);
ltoe = raw(:,18);
rtoe = raw(:,37);
% plot(ltoe)
% hold on
% plot(rtoe)
% step state, 0 = unknown, 1 = heel stoke, 2 = toe off
lstate = 0;rstate = 0;
% heelstrikes
lheelcount = 0;lheelstrikes = 0;
rheelcount = 0;rheelstrikes = 0;
% toe offs
ltoecount = 0;ltoeoffs = 0;
rtoecount = 0;rtoeoffs = 0;
% Iterate over raw data
for i = 2:length(raw)-3
% Detect heelstrike left
if (lstate ~= 1 && (lheel(i+1) >= 1 && lheel(i+1)-lheel(i) > 1 && lheel(i+2)-lheel(i+1) > 2.5 || lheel(i+1)-lheel(i) > 3))
% if (lstate ~= 1 && (lheel(i-1) >= 30 && lheel(i+1)-lheel(i) >= 0 && ...
% lheel(i+2) - lheel(i+1) >= 0 && lheel(i+3) - lheel(i+2) >= 0 && lheel(i+3) - lheel(i) >= 100))
lstate = 1;
lheelcount = lheelcount +1;
lheelstrikes(lheelcount) = i;
continue;
end
% Detect toe off left
if (lstate == 1 && (ltoe(i-1)> 50 && ltoe(i-1) - ltoe(i) >= 0 && ltoe(i) - ltoe(i+1)>= 0 && ltoe(i+1) <= 50))
lstate = 2;
ltoecount = ltoecount +1;
ltoeoffs(ltoecount) = i+1;
continue;
end
end
% Iterate over raw data
for i = 2:length(raw)-3
% Detect heelstrike right
if (rstate ~= 1 && (rheel(i+1) >= 1 && rheel(i+1)-rheel(i) > 1 && rheel(i+2)-rheel(i+1) > 2.5 || rheel(i+1)-rheel(i) > 3))
% if (rstate ~= 1 && (rheel(i-1) >= 30 && rheel(i+1) - rheel(i) >= 0 && ...
% rheel(i+2) - rheel(i+1) >= 0 && rheel(i+3) - rheel(i+2) >= 0 && rheel(i+3) - rheel(i) >= 100))
rstate = 1;
rheelcount = rheelcount +1;
rheelstrikes(rheelcount) = i;
continue;
end
% Detect toe off right
% if (rstate == 1 && rheel(i+1) < 3.75 && rtoe(i+1) <= 4 && rtoe(i+2) <= 6 && rtoe(i) < rtoe(i-1)-0.3 && rtoe(i+1) < rtoe(i-1)-0.3 && abs(rtoe(i+1)-rtoe(i+2)) < 3 && (rtoe(i)-rtoe(i+1)) > (rtoe(i+1)-rtoe(i+2))*1.7 && ((rtoe(i)-rtoe(i+1)) >= 1 || (rtoe(i-1)-rtoe(i)) >= 2))
if (rstate == 1 && (rtoe(i-1)> 50 && rtoe(i-1) - rtoe(i) >= 0 && rtoe(i) - rtoe(i+1)>= 0 && rtoe(i+1) <= 50 && rtoe(i+3) <= 40))
rstate = 2;
rtoecount = rtoecount +1;
rtoeoffs(rtoecount) = i+1;
continue;
end
end
%% Find ground contact times of left foot steps and right foot steps
lsteps = []; % left foot heelstrike time, left foot toe-off time and left gct
t = 1; % gct-> ground contact time
for i = 1:length(lheelstrikes)
while(t <= length(ltoeoffs))
%raw(ltoeoffs(t),1)
%raw(lheelstrikes(i),1)
gct = raw(ltoeoffs(t),1) - raw(lheelstrikes(i),1);
if (gct > 0 && gct < 2)
stime = ltimes(ltoeoffs(t));
lsteps = [lsteps; [stime-gct, stime, gct]];
t = t+1;
break;
elseif (gct >= 2)
break;
elseif (gct <= 0)
t = t+1;
end
%pause
end
end
% right foot
rsteps = [];% right foot heelstrike time, right foot toe-off time and right gct
t = 1;
for i = 1:length(rheelstrikes)
while(t <= length(rtoeoffs))
%raw(rtoeoffs(t),1)
%raw(rheelstrikes(i),1)
gct = raw(rtoeoffs(t),1) - raw(rheelstrikes(i),1);
if (gct > 0 && gct < 2)
stime = rtimes(rtoeoffs(t));
rsteps = [rsteps; [stime-gct, stime, gct]];
t = t+1;
break;
elseif (gct >= 2)
break;
elseif (gct <= 0)
t = t+1;
end
%pause
end
end
%% plot ground contact times
% to see if the steps segmentation is done perfectly or not. It can be seen
% by the plots of contact time of left foot and right foot. If it not
% correct for some steps have higher contact time (peaks) than
% nearby steps, which can be seen in plots.
%(sync jumps have higher contact time ~1s (peaks), so it should not be
% considered as step segmentation error)
% if it is not correct then the algo should be tweaked in lines (456, 465, 479, 488)
% as per the slope of force curve by intuition (by increasing those differece no's).
% figure;
% plot(lsteps(:, 1), lsteps(:, 3));
% hold on;
% plot(rsteps(:, 1), rsteps(:, 3));
figure;
plot(lsteps(:, 3));
hold on;
plot(rsteps(:, 3));
%% code to lable the data
lforce=raw(:,18)'; %left foot vertical force
rforce=raw(:,37)'; %right foot vertical force
% find index of foot heel strike and toe_off in INS data (TIA time).
% lsteps contains heel strike, toe-off and their time difference
% for left foot and same for rsteps.
for i=1:length(lsteps)
[tempLsVal,tempLsInd]= min(abs(times2ins - lsteps(i,1)));
[tempLeVal,tempLeInd]= min(abs(times2ins - lsteps(i,2)));
%technically 1/2*Fs should work instead of 1 assuming INS data is
%present at 1/Fs time points
lstepinds(i).no= i;
if (tempLsVal < 1 && tempLeVal < 1)%1s
lstepinds(i).start= tempLsInd;
lstepinds(i).end= tempLeInd;
end
%if length(lsteps) is less than length(lstepinds) then datalogger was
% turned OFF before moticon (given, if insoles were turned ON after the
% datalogger) **
end
for i=1:length(rsteps)
[tempRsVal,tempRsInd]= min(abs(times2ins - rsteps(i,1)));
[tempReVal,tempReInd]= min(abs(times2ins - rsteps(i,2)));
%technically 1/2*Fs should work instead of 1 assuming INS data is
%present at 1/Fs time points
rstepinds(i).no= i;
if ( tempRsVal < 1 && tempReVal < 1)
rstepinds(i).start = tempRsInd;
rstepinds(i).end = tempReInd;
end
end
for i=1:length(lstepinds)
if isempty(lstepinds(i).start)
lstepinds(i).start = 0;
end
if isempty(lstepinds(i).end)
lstepinds(i).end = 0;
end
end
for i=1:length(rstepinds)
if isempty(rstepinds(i).start)
rstepinds(i).start = 0;
end
if isempty(rstepinds(i).end)
rstepinds(i).end = 0;
end
end
% initialise lables for both foot
lstepLable=zeros(1,length(insind));% if 1 then left foot is on the ground
rstepLable=zeros(1,length(insind));% if 1 then right foot is on the ground
% put lable=1 if the timepoint is during foot is contact with the ground.
for i=1:length(lstepinds)
if(lstepinds(i).start ~=0 && lstepinds(i).end ~=0)
lstepLable(lstepinds(i).start:lstepinds(i).end)=ones(1,lstepinds(i).end-lstepinds(i).start+1);
end
end
for i=1:length(rstepinds)
if(rstepinds(i).start ~=0 && rstepinds(i).end ~=0)
rstepLable(rstepinds(i).start:rstepinds(i).end)=ones(1,rstepinds(i).end-rstepinds(i).start+1);
end
end
doubleSupportlable=lstepLable & rstepLable; % 1, if both foot on the ground
flightLable=not(lstepLable) & not(rstepLable); %1, if both feet in the air
% force interpolation (50Hz/100Hz --> 400Hz) to get foot force on instance
% of TIA time points of INS datalogger
lforce_intrp=interp1(ltimes,lforce,times2ins);
rforce_intrp=interp1(rtimes,rforce,times2ins);
plot(ltimes,lforce, 'o')
hold on
plot(times2ins,lforce_intrp,'*')
%% To select the range of velocity for further processing
%plot the figure and put two cursors for suitable data window and export
%points to worspace.
figure
plot(times2ins, v_lon)
hold on
plot(times2ins, acc_anatomical(3,:))
pause(0.00001);
frame_h = get(handle(gcf),'JavaFrame');
set(frame_h,'Maximized',1);
%% To extract various index for selected data range using curser points
timepoints=size(cursor_info,2)
%%index for selected run time windows
for i = 1:timepoints/2
dataset_inds(i).rangeseq=i;
dataset_inds(i).windowStartInds=cursor_info(timepoints-2*(i-1)-0).DataIndex;
dataset_inds(i).windowEndInds=cursor_info(timepoints-2*(i-1)-1).DataIndex;
dataset_inds(i).rangeWindow_timepoints=dataset_inds(i).windowEndInds-dataset_inds(i).windowStartInds;
dataset_inds(i).windowStartTimes=cursor_info(timepoints-2*(i-1)-0).Position(1);
dataset_inds(i).windowEndTimes=cursor_info(timepoints-2*(i-1)-1).Position(1);
[~,dataset_inds(i).ltimesStartInds] = min(abs(ltimes- dataset_inds(i).windowStartTimes));
[~,dataset_inds(i).ltimesEndInds] = min(abs(ltimes- dataset_inds(i).windowEndTimes));
[~,dataset_inds(i).rtimesStartInds] = min(abs(rtimes- dataset_inds(i).windowStartTimes));
[~,dataset_inds(i).rtimesEndInds] = min(abs(rtimes- dataset_inds(i).windowEndTimes));
[~,dataset_inds(i).rstepStartInds] = min(abs(rsteps(:,1)- dataset_inds(i).windowStartTimes));
[~,dataset_inds(i).rstepEndInds] = min(abs(rsteps(:,1)- dataset_inds(i).windowEndTimes));
[~,dataset_inds(i).lstepStartInds] = min(abs(lsteps(:,1)- dataset_inds(i).windowStartTimes));
[~,dataset_inds(i).lstepEndInds] = min(abs(lsteps(:,1)-dataset_inds(i).windowEndTimes));
end
%% Run Once
input_data={};
target_data={};
GRFL_strides={};
GRFR_strides={};
dwc=1; % data window count
%% Velocity Segmentation and time shifts for choosen range
% % % % %range using curser points
% windowStart = dataset_inds(1).windowStartInds ;%
% windowEnd = dataset_inds(1).windowEndInds;
% % walking and then jogging
%windowStart = 135305;
%windowEnd = 203301;
% %walking + fast running
%windowStart = 218000;
%windowEnd = 231000;
% % fast running + then some walking
%windowStart = 244000;
%windowEnd = 263644;
% % % % walking
windowStart = 272070;
windowEnd = 376847;
%% dataset for LSTM
A_data = [lstepLable.' lstepLable.' rstepLable.' acc_ins.' ang_vel_ins.' ...
quat.' v_lon.' V_d.' lforce_intrp.' rforce_intrp.'];
A_data1 = A_data(135305:203301, :);
A_data2 = A_data(218000:231000, :);
A_data3 = A_data(244000:263644, :);
A_data4 = A_data(272070:376847, :);
All_data = [A_data1; A_data2; A_data3; A_data4];
csvwrite('181003_testdata.csv',All_data)
%%
% %%8steps
% windowStart =270000;
% windowEnd = 272800;
% % %%12 steps
% windowStart =245000;
% windowEnd = 247000;
GrndTrack_unwrap= unwrap(GrndTrack);
GrndTrack_osc = highpass(GrndTrack_unwrap, 0.005,'StopbandAttenuation',30,'Steepness',0.7);
Fs = 400;
V_d = highpass(V_d, 0.005,'StopbandAttenuation',30,'Steepness',0.7);
vertVel = V_d(1, windowStart:windowEnd);
vertVel_d= vertVel;
% figure
% plot(vertVel_d,':*')
% hold on
vertical_oscillations=[];
vertical_oscillations=[vertical_oscillations, cumtrapz(V_d./Fs)]; % vertical distance for complete window
vertical_oscillations = highpass(vertical_oscillations, 0.005,'StopbandAttenuation',30,'Steepness',0.7);
vertical_oscillations=vertical_oscillations*100; % in cm
vertical_oscillations_d= vertical_oscillations(windowStart:windowEnd);
maxV = max(vertVel_d);
thresh = 0.05 * maxV;
a = false;b = false;
front = 0;back = 0;
while (true)
if (abs(vertVel_d(1, 1)) > thresh || vertVel_d(1, 1) < vertVel_d(1, 2))%
vertVel_d = vertVel_d(1, 3:end);
front = front + 2;
a = false;
else
if (GrndTrack_osc(windowStart + front + Fs/8)<0)
a = true;
else
vertVel_d = vertVel_d(1, 50:end);
front = front + 49;
end
end
if (abs(vertVel_d(1, end)) > thresh || vertVel_d(1, end-1) < vertVel_d(1, end))%
vertVel_d = vertVel_d(1, 1:end - 2);
%vertVel_d = detrend(vertVel_d);
back = back + 2;
a = false;
b = false;
else
b = true;
end
if (a && b)
break;
end
end
% plot(vertVel_d,':*')
velsteppoints = [];
if (vertVel_d(1, 1) < 0)%
intersections = 1;
velsteppoints = 1;
else
intersections = 0;
end
count = 0;
points = [];
for i = 1:length(vertVel_d) - 1
if (vertVel_d(1, i) < 0 && vertVel_d(1, i + 1) > 0)
intersections = intersections + 1;
points = [points, i];
end
if (vertVel_d(1, i) > 0 && vertVel_d(1, i + 1) < 0 && vertical_oscillations_d(front + i) > 0.5)
intersections = intersections + 1;
points = [points, i];
velsteppoints = [velsteppoints, i];%
count = count + 1;%
end
end
velstepCount = length(velsteppoints) - 1;
disp(['Steps missed during segmentation: ', num2str(2*count-intersections-1)]);
e = datenum('01-jan-1970 00:00:00');
datestr(e);
Itime = datestr(e + (times2ins(windowStart + front)-37)/86400,'mmmm dd, yyyy HH:MM:SS.FFF AM');% UTC time for first point for first step
Itime
% start and end index of a step (when calculated according to the vertical velocity)
velstepinds = struct([]);
for i = 1:velstepCount
velstepinds(i).start=windowStart + front + velsteppoints(i);
velstepinds(i).end=windowStart + front + velsteppoints(i+1)-1;
end
v_lon_window = v_lon(1, windowStart:windowEnd);
velstepDuration = []; speed=[]; vel_lon_diff=[];vel_lon_steps=[];
velstepLength = []; cadence = []; vertical_oscillations_window=[];
vertical_oscillations_segmented=[]; vertical_oscillations_amp=[];
vertOscillation_vel_err=[]; % error over a step
vertical_oscillations_window=[vertical_oscillations_window, cumtrapz(vertVel./Fs)]; % vertical distance for complete selected window (with front and back)
vertical_oscillations_window=vertical_oscillations_window.*100;
vertical_oscillations_segmented = vertical_oscillations_window(front:end-back-1); % without front and back
v_dist_shift= mean(vertical_oscillations_segmented); %shift for the mid point of oscillations
vertical_oscillations_segmented= vertical_oscillations_segmented-v_dist_shift; % shifted vertical displacement without front and back
for i = 1:velstepCount
velstepDuration = [velstepDuration, (front + velsteppoints(i + 1) - (front + velsteppoints(i)))/Fs];
cadence = [cadence, 1/velstepDuration(i)];
stepVelsLon = v_lon_window(front + velsteppoints(i):(front + velsteppoints(i + 1)-1));
% stepVelsLon = v_lon(1,(velstepinds(i).start):(velstepinds(i).end));
speed(i)=mean(stepVelsLon);
vel_lon_diff(i)=max(stepVelsLon)-min(stepVelsLon);
velstepLength = [velstepLength, trapz(stepVelsLon./Fs)];% longitudinal displacement in a step
stepDispDown = vertical_oscillations_window((front + velsteppoints(i)):(front + velsteppoints(i+1)-1));
stepDispDown_p_p= max(stepDispDown) - min(stepDispDown);
vertical_oscillations_amp = [vertical_oscillations_amp,stepDispDown_p_p]; % peak to peak
vertOscillation_vel_err = [vertOscillation_vel_err,mean(vertical_oscillations_segmented(velsteppoints(i):(velsteppoints(i + 1)-1)))];
end
for i = 1:velstepCount
contactTimeL(i)=sum(lstepLable(velstepinds(i).start:velstepinds(i).end)==1)/Fs;
contactTimeR(i)=sum(rstepLable(velstepinds(i).start:velstepinds(i).end)==1)/Fs;
flightTime(i)=sum(flightLable(velstepinds(i).start:velstepinds(i).end)==1)/Fs;
doubleSupportTime(i)=sum(doubleSupportlable(velstepinds(i).start:velstepinds(i).end)==1)/Fs;
end
% Stride parameters for velocity segmenattion, left foot stride and right foot stride and their comparision
%start and end points for each stride in window (when calculated according to the vertical velocity)
% one stride is consist of two steps
vel_stride_points=[];
for j=1:ceil(length(velsteppoints)/2)
vel_stride_points=[vel_stride_points velsteppoints(2*j-1)];
end
velStrideCount = length(vel_stride_points) - 1;
% start and end index of the strides (calculated using vertical velocity)
velStrideinds = struct([]);
for j = 1:velStrideCount
velStrideinds(j).start=windowStart + front + vel_stride_points(j);
velStrideinds(j).end=windowStart + front + vel_stride_points(j+1)-1;
end
% To calculate stride parameters in a velocity stride
velStrideDuration = []; velStrideLength = []; velStridecadence=[];
velStrideSpeed=[]; fVel_pp=[];fVel_std=[];
velstrideDispDown = []; velstrideDispDown_p_p = [];
vertOscillation_velstride_err= [];
velStride_vertDisp_pp = []; velStride_vertDisp_std = [];
V_acc_pp=[]; F_acc_pp=[]; vVel_pp=[]; V_acc_std=[]; F_acc_std=[]; vVel_std=[];
yaw_pp=[]; roll_pp=[]; pitch_pp=[]; yaw_std=[]; roll_std=[]; pitch_std=[];
wx_pp=[]; wy_pp=[]; wz_pp=[]; wx_std=[]; wy_std=[]; wz_std=[];
GrndTrack_pp=[]; GrndTrack_std=[];
GCTL_stride=[];GCTR_stride=[];%GCT in stride
FT_stride=[];%flight time
DST_stride=[];% double suppport time
TD_L=[];TO_L=[];footTimeL=[];TD_R=[];TO_R=[];footTimeR=[];
GRFR_stride={}; GRFL_stride={}; GRFR_area=[]; GRFL_area=[]; GRFR_max=[]; GRFL_max=[];
for i = 1:velStrideCount
i
velStrideDuration(i) = (windowStart + front + vel_stride_points(i + 1) - (windowStart + front + vel_stride_points(i)))/Fs;
velStridecadence(i) = 1/velStrideDuration(i);
velStrideSpeed(i) = mean(v_lon(windowStart + front + vel_stride_points(i) : windowStart + front + vel_stride_points(i + 1)-1));
fVel_pp(i) = peak2peak(v_lon(windowStart + front + vel_stride_points(i) : windowStart + front + vel_stride_points(i + 1)-1));
fVel_std(i) = std(v_lon(windowStart + front + vel_stride_points(i) : windowStart + front + vel_stride_points(i + 1)-1));
velStrideLength(i) = trapz(v_lon(windowStart + front + vel_stride_points(i) : windowStart + front + vel_stride_points(i + 1)-1)./Fs);% stride legth in direction of movement (typically 2* velsteplength)
velStride_vertDisp_pp(i) = peak2peak(vertical_oscillations_window((front + vel_stride_points(i)):(front + vel_stride_points(i+1)-1))); % peak to peak
velStride_vertDisp_std(i) = std(vertical_oscillations_window((front + vel_stride_points(i)):(front + vel_stride_points(i+1)-1))); %std
vertOscillation_velstride_err(i) = mean(vertical_oscillations_segmented(vel_stride_points(i):(vel_stride_points(i + 1)-1)));
V_acc_pp(i) = peak2peak(acc_anatomical(3,(windowStart + front + vel_stride_points(i)): (windowStart + front + vel_stride_points(i+1)-1)));
F_acc_pp(i) = peak2peak(acc_anatomical(1,(windowStart + front + vel_stride_points(i)): (windowStart + front + vel_stride_points(i+1)-1)));
vVel_pp(i) = peak2peak(V_d(windowStart + front + vel_stride_points(i) : windowStart + front + vel_stride_points(i+1)-1));
GrndTrack_pp(i) = peak2peak(GrndTrack_unwrap(windowStart + front + vel_stride_points(i) : windowStart + front + vel_stride_points(i+1)-1));
% yaw_pp(i) = peak2peak(yaw_unwrap(windowStart + front + vel_stride_points(i) : windowStart + front + vel_stride_points(i+1)-1));
% roll_pp(i) = peak2peak(roll(windowStart + front + vel_stride_points(i) : windowStart + front + vel_stride_points(i+1)-1));
% pitch_pp(i) = peak2peak(pitch(windowStart + front + vel_stride_points(i) : windowStart + front + vel_stride_points(i+1)-1));
% wx_pp(i) = peak2peak(ang_vel_ins(1, windowStart + front + vel_stride_points(i) : windowStart + front + vel_stride_points(i+1)-1));
% wy_pp(i) = peak2peak(ang_vel_ins(2, windowStart + front + vel_stride_points(i) : windowStart + front + vel_stride_points(i+1)-1));
% wz_pp(i) = peak2peak(ang_vel_ins(3, windowStart + front + vel_stride_points(i) : windowStart + front + vel_stride_points(i+1)-1));
yaw_pp(i) = max(yaw_unwrap(windowStart + front + vel_stride_points(i) : windowStart + front + vel_stride_points(i+1)-1))-...
min(yaw_unwrap(windowStart + front + vel_stride_points(i) : windowStart + front + vel_stride_points(i+1)-1));
roll_pp(i) = max(roll(windowStart + front + vel_stride_points(i) : windowStart + front + vel_stride_points(i+1)-1))-...
min(roll(windowStart + front + vel_stride_points(i) : windowStart + front + vel_stride_points(i+1)-1));
pitch_pp(i) = max(pitch(windowStart + front + vel_stride_points(i) : windowStart + front + vel_stride_points(i+1)-1))-...
min(pitch(windowStart + front + vel_stride_points(i) : windowStart + front + vel_stride_points(i+1)-1));
wx_pp(i) = max(ang_vel_ins(1, windowStart + front + vel_stride_points(i) : windowStart + front + vel_stride_points(i+1)-1))-...
min(ang_vel_ins(1,windowStart + front + vel_stride_points(i) : windowStart + front + vel_stride_points(i+1)-1));
wy_pp(i) = max(ang_vel_ins(2, windowStart + front + vel_stride_points(i) : windowStart + front + vel_stride_points(i+1)-1))-...
min(ang_vel_ins(2,windowStart + front + vel_stride_points(i) : windowStart + front + vel_stride_points(i+1)-1));
wz_pp(i) = max(ang_vel_ins(3, windowStart + front + vel_stride_points(i) : windowStart + front + vel_stride_points(i+1)-1))-...
min(ang_vel_ins(3, windowStart + front + vel_stride_points(i) : windowStart + front + vel_stride_points(i+1)-1));
V_acc_std(i) = std(acc_anatomical(3,(windowStart + front + vel_stride_points(i)): (windowStart + front + vel_stride_points(i+1)-1)));
F_acc_std(i) = std(acc_anatomical(1,(windowStart + front + vel_stride_points(i)): (windowStart + front + vel_stride_points(i+1)-1)));
vVel_std(i) = std(V_d(windowStart + front + vel_stride_points(i) : windowStart + front + vel_stride_points(i+1)-1));
GrndTrack_std(i) = std(GrndTrack_unwrap(windowStart + front + vel_stride_points(i) : windowStart + front + vel_stride_points(i+1)-1));
yaw_std(i) = std(yaw_unwrap(windowStart + front + vel_stride_points(i) : windowStart + front + vel_stride_points(i+1)-1));
roll_std(i) = std(roll(windowStart + front + vel_stride_points(i) : windowStart + front + vel_stride_points(i+1)-1));
pitch_std(i) = std(pitch(windowStart + front + vel_stride_points(i) : windowStart + front + vel_stride_points(i+1)-1));
wx_std(i) = std(ang_vel_ins(1, windowStart + front + vel_stride_points(i) : windowStart + front + vel_stride_points(i+1)-1));
wy_std(i) = std(ang_vel_ins(2, windowStart + front + vel_stride_points(i) : windowStart + front + vel_stride_points(i+1)-1));
wz_std(i) = std(ang_vel_ins(3, windowStart + front + vel_stride_points(i) : windowStart + front + vel_stride_points(i+1)-1));
end
for i = 1:(velStrideCount)
i
GCTL_stride(i)=sum(lstepLable(velStrideinds(i).start:velStrideinds(i).end)==1)/Fs;
GCTR_stride(i)=sum(rstepLable(velStrideinds(i).start:velStrideinds(i).end)==1)/Fs;
FT_stride(i)=sum(flightLable(velStrideinds(i).start:velStrideinds(i).end)==1)/Fs;
DST_stride(i)=sum(doubleSupportlable(velStrideinds(i).start:velStrideinds(i).end)==1)/Fs;
footTimeL = find((lstepLable(velStrideinds(i).start : velStrideinds(i).end)==1));
TD_L(i) = min(footTimeL)/Fs;
TO_L(i) = max(footTimeL)/Fs;%
%%%%%%%%%%%%%%%%%%%%******************************
footTimeR = find((rstepLable(velStrideinds(i).start:velStrideinds(i).end)==0));%%%%%%%%
TO_R(i)= min(footTimeR)/Fs;
TD_R(i)= max(footTimeR)/Fs;%
if(TO_R(i)*Fs==1)
footTimeR = find((rstepLable(velStrideinds(i).start:velStrideinds(i).end)==1));%%%%%%%%
TD_R(i)= min(footTimeR)/Fs;
TO_R(i)= max(footTimeR)/Fs;%
end
GRFL_stride{i} = lforce_intrp((windowStart + front + vel_stride_points(i)) : (windowStart + front + vel_stride_points(i+1)-1));
GRFR_stride{i} = rforce_intrp((windowStart + front + vel_stride_points(i)) : (windowStart + front + vel_stride_points(i+1)-1));
GRFL_area= [GRFL_area, trapz(GRFL_stride{i} ./Fs)];%
GRFR_area= [GRFR_area, trapz(GRFR_stride{i} ./Fs)];%
GRFL_max = [GRFL_max, max(GRFL_stride{i})];%
GRFR_max = [GRFR_max, max(GRFR_stride{i})];%
end
% combine input / target data cells
input_data{dwc} =[
velStrideSpeed; ...
velStrideDuration;...
velStrideLength;...
F_acc_pp; ...
V_acc_pp; ...
fVel_pp; ...
vVel_pp;...
velStride_vertDisp_pp;...
GrndTrack_pp;...
yaw_pp;...
roll_pp;...
pitch_pp;...
wx_pp;...
wy_pp;...
wz_pp;...
F_acc_std; ...
V_acc_std; ...
fVel_std; ...
vVel_std;...
velStride_vertDisp_std;...
GrndTrack_std;...
yaw_std;...
roll_std;...
pitch_std;...
wx_std;...
wy_std;...
wz_std
];
target_data{1,dwc}=[ ...
TD_L; ... % touch down leftfoot
TO_L; ... % toeoff leftfoot
GCTL_stride; ... % GCTLF
GRFL_max; ... % GRFLF max
GRFL_area; ... % Area under GRFLF curve
TD_R; ... % touch down rightfoot
TO_R; ... % toeoff leftfoot
GCTR_stride; ... % GCTRF
GRFR_max; ... % GRFRF max
GRFR_area; ... % Area under GRFRF curve
FT_stride; ... % flight time