Skip to content

Commit

Permalink
Added initial code
Browse files Browse the repository at this point in the history
  • Loading branch information
Janis Wojtusch committed May 4, 2015
1 parent 883637c commit c32ffbe
Show file tree
Hide file tree
Showing 47 changed files with 8,095 additions and 28 deletions.
2 changes: 2 additions & 0 deletions .gitignore
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,2 @@
# Temporary files
*.*~
887 changes: 887 additions & 0 deletions AntropometricParameterGeneration.m
View file

Large diffs are not rendered by default.

563 changes: 563 additions & 0 deletions CenterOfPressureEstimation.m
View file

Large diffs are not rendered by default.

Binary file added EventDetection.fig
View file
Binary file not shown.
1,060 changes: 1,060 additions & 0 deletions EventDetection.m
View file

Large diffs are not rendered by default.

203 changes: 203 additions & 0 deletions EventVisualization.m
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,203 @@
% ------------------------------------------------------
% This script visualizes ground reaction force and events data.
% ------------------------------------------------------
% Technische Universität Darmstadt
% Department of Computer Science
% Simulation, Systems Optimization and Robotics Group
% Janis Wojtusch (wojtusch@sim.tu-darmstadt.de), 2015
% Licensed under BSD 3-Clause License
% ------------------------------------------------------

% Clean up workspace
clc;
clear all;
close all;

% Add functions to search path
addpath('Scripts');

% Set parameters
savePath = [getPath, '../Grafiken/Events/'];
datasets = {
'1.1', ...
'1.2', ...
'1.3', ...
'2.1', ...
'2.2', ...
'2.3', ...
'3', ...
'4', ...
'5.1', ...
'5.2', ...
'6', ...
'7', ...
'8', ...
'9.1', ...
'9.2', ...
'9.3' ...
};
subjects = {
'A',...
'B'...
};

for subjectIndex = 1:length(subjects)
for datasetIndex = 1:length(datasets)

% Set parameters
dataset = datasets{datasetIndex};
subject = subjects{subjectIndex};

% Load data file
file = getFile(subject, dataset);
if file
variables = load(file);
if isfield(variables, 'force') && isfield(variables, 'events')
force = variables.force;
events = variables.events;
name = regexp(file, '[^/]*(?=\.[^.]+($|\?))', 'match');
name = name{1};
else
fprintf('WARNING: No matching data found!\n');
continue;
end
else
fprintf('WARNING: No matching data file found!\n');
continue;
end

% Plot force and events data
visualization = figure('Name', 'Events', 'NumberTitle', 'off', 'Color', 'white', 'Position', [0, 0, 1400, 600]);
time = 0:(1 / force.frameRate):((force.frames - 1) / force.frameRate);
subplot(3, 1, 1);
if isfield(force, 'groundReactionForceX_L')
plot(time, force.groundReactionForceX, 'Color', [0.7, 0.7, 0.7]);
hold on;
plot(time, force.groundReactionForceX_L, 'r-');
plot(time, force.groundReactionForceX_R, 'b-');
title('Ground reaction force X (gray) / X_L (red) / X_R (blue)');
else
plot(time, force.groundReactionForceX, 'k-');
title('Ground reaction force X');
end
xlabel('Time in s');
ylabel('Force in N');
grid on;
maximumValue = max(force.groundReactionForceX);
minimumValue = min(force.groundReactionForceX);
if ~strcmp(dataset, '6') && ~strcmp(dataset, '7')
for eventIndex = 1:length(events.eventStart_L)
if events.groundReactionForceCorrection_L(eventIndex)
eventColor = [1, 0.7, 0];
else
eventColor = 'r';
end
patch([events.eventStart_L(eventIndex), events.eventStart_L(eventIndex), events.eventEnd_L(eventIndex), events.eventEnd_L(eventIndex)], [minimumValue, maximumValue, maximumValue, minimumValue], eventColor, 'EdgeColor', 'none', 'FaceAlpha', 0.3);
end
for eventIndex = 1:length(events.eventStart_R)
if events.groundReactionForceCorrection_R(eventIndex)
eventColor = [0, 0.7, 1];
else
eventColor = 'b';
end
patch([events.eventStart_R(eventIndex), events.eventStart_R(eventIndex), events.eventEnd_R(eventIndex), events.eventEnd_R(eventIndex)], [minimumValue, maximumValue, maximumValue, minimumValue], eventColor, 'EdgeColor', 'none', 'FaceAlpha', 0.3);
end
else
for eventIndex = 1:length(events.eventStart)
if events.groundReactionForceCorrection(eventIndex)
eventColor = [1, 0.7, 0];
else
eventColor = 'r';
end
patch([events.eventStart(eventIndex), events.eventStart(eventIndex), events.eventEnd(eventIndex), events.eventEnd(eventIndex)], [minimumValue, maximumValue, maximumValue, minimumValue], eventColor, 'EdgeColor', 'none', 'FaceAlpha', 0.3);
end
end
subplot(3, 1, 2);
plot(time, force.groundReactionForceY_L, 'r-');
hold on;
plot(time, force.groundReactionForceY_R, 'b-');
title('Ground reaction force Y_L (red) / Y_R (blue)');
xlabel('Time in s');
ylabel('Force in N');
grid on;
maximumValue = max(max(force.groundReactionForceY_L), max(force.groundReactionForceY_R));
minimumValue = min(min(force.groundReactionForceY_L), min(force.groundReactionForceY_R));
if ~strcmp(dataset, '6') && ~strcmp(dataset, '7')
for eventIndex = 1:length(events.eventStart_L)
if events.groundReactionForceCorrection_L(eventIndex)
eventColor = [1, 0.7, 0];
else
eventColor = 'r';
end
patch([events.eventStart_L(eventIndex), events.eventStart_L(eventIndex), events.eventEnd_L(eventIndex), events.eventEnd_L(eventIndex)], [minimumValue, maximumValue, maximumValue, minimumValue], eventColor, 'EdgeColor', 'none', 'FaceAlpha', 0.3);
end
for eventIndex = 1:length(events.eventStart_R)
if events.groundReactionForceCorrection_R(eventIndex)
eventColor = [0, 0.7, 1];
else
eventColor = 'b';
end
patch([events.eventStart_R(eventIndex), events.eventStart_R(eventIndex), events.eventEnd_R(eventIndex), events.eventEnd_R(eventIndex)], [minimumValue, maximumValue, maximumValue, minimumValue], eventColor, 'EdgeColor', 'none', 'FaceAlpha', 0.3);
end
else
for eventIndex = 1:length(events.eventStart)
if events.groundReactionForceCorrection(eventIndex)
eventColor = [1, 0.7, 0];
else
eventColor = 'r';
end
patch([events.eventStart(eventIndex), events.eventStart(eventIndex), events.eventEnd(eventIndex), events.eventEnd(eventIndex)], [minimumValue, maximumValue, maximumValue, minimumValue], eventColor, 'EdgeColor', 'none', 'FaceAlpha', 0.3);
end
end
subplot(3, 1, 3);
if isfield(force, 'groundReactionForceZ_L')
plot(time, force.groundReactionForceZ, 'Color', [0.7, 0.7, 0.7]);
hold on;
plot(time, force.groundReactionForceZ_L, 'r-');
plot(time, force.groundReactionForceZ_R, 'b-');
title('Ground reaction force Z (gray) / Z_L (red) / Z_R (blue)');
else
plot(time, force.groundReactionForceZ, 'k-');
title('Ground reaction force Z');
end
xlabel('Time in s');
ylabel('Force in N');
grid on;
maximumValue = max(force.groundReactionForceZ);
minimumValue = min(force.groundReactionForceZ);
if ~strcmp(dataset, '6') && ~strcmp(dataset, '7')
for eventIndex = 1:length(events.eventStart_L)
if events.groundReactionForceCorrection_L(eventIndex)
eventColor = [1, 0.7, 0];
else
eventColor = 'r';
end
patch([events.eventStart_L(eventIndex), events.eventStart_L(eventIndex), events.eventEnd_L(eventIndex), events.eventEnd_L(eventIndex)], [minimumValue, maximumValue, maximumValue, minimumValue], eventColor, 'EdgeColor', 'none', 'FaceAlpha', 0.3);
end
for eventIndex = 1:length(events.eventStart_R)
if events.groundReactionForceCorrection_R(eventIndex)
eventColor = [0, 0.7, 1];
else
eventColor = 'b';
end
patch([events.eventStart_R(eventIndex), events.eventStart_R(eventIndex), events.eventEnd_R(eventIndex), events.eventEnd_R(eventIndex)], [minimumValue, maximumValue, maximumValue, minimumValue], eventColor, 'EdgeColor', 'none', 'FaceAlpha', 0.3);
end
else
for eventIndex = 1:length(events.eventStart)
if events.groundReactionForceCorrection(eventIndex)
eventColor = [1, 0.7, 0];
else
eventColor = 'r';
end
patch([events.eventStart(eventIndex), events.eventStart(eventIndex), events.eventEnd(eventIndex), events.eventEnd(eventIndex)], [minimumValue, maximumValue, maximumValue, minimumValue], eventColor, 'EdgeColor', 'none', 'FaceAlpha', 0.3);
end
end
suplabel([dataset, ' ', name, ' - HuMoD Database']);

% Save figure
saveTightFigure(visualization, [savePath, subject, filesep, dataset, '.png']);
fprintf('STATUS: Figure for dataset %s %s was saved.\n', subject, dataset);
close all;

end
end
182 changes: 182 additions & 0 deletions ForceFilter.m
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,182 @@
% ------------------------------------------------------
% This script processes the raw force data and transforms the coordinate
% system according to [Wu2002], [Wu2005]. Also motion and force data is
% synchronized by compensating the time delay between motion and force
% measurements.
% ------------------------------------------------------
% Technische Universität Darmstadt
% Department of Computer Science
% Simulation, Systems Optimization and Robotics Group
% Janis Wojtusch (wojtusch@sim.tu-darmstadt.de), 2015
% Licensed under BSD 3-Clause License
% ------------------------------------------------------

% Clean up workspace
clc;
clear all;
close all;

% Set parameters
intercept = 0.017767380272882; % Intercept of the linear regression for force measurement delay
gradient = -5.713883360543501e-05; % Gradient of the linear regression for force measurement delay
filterHalfOrder = 3; % Half order of the zero-phase high-pass and low-pass filters
filterCutOff = 50; % Cut-off frequency of the zero-phase low-pass filter
datasets = {
'1.1', ...
'1.2', ...
'1.3', ...
'2.1', ...
'2.2', ...
'2.3', ...
'3', ...
'4', ...
'5.1', ...
'5.2', ...
'6', ...
'7', ...
'8', ...
'9.1', ...
'9.2', ...
'9.3' ...
};
subjects = {
'A',...
'B'...
};
signals = {
'Fx1', ...
'Fx2', ...
'Fy1', ...
'Fy2', ...
'FzL1', ...
'FzL2', ...
'FzL3', ...
'FzL4', ...
'FzR1', ...
'FzR2', ...
'FzR3', ...
'FzR4' ...
};

% Add functions to search path
addpath('Scripts');

for subjectIndex = 1:length(subjects)
for datasetIndex = 1:length(datasets)

% Set parameters
dataset = datasets{datasetIndex};
subject = subjects{subjectIndex};
fprintf('STATUS: Processing dataset %s %s.\n', subject, dataset);

% Load data file
saveFile = getFile(subject, dataset);
if saveFile
[path, file, extention] = fileparts(saveFile);
dataFile = [path, filesep, 'Raw data', filesep, 'Force ', file, extention];
variables = load(dataFile);
frameRate = variables.fs;
frames = length(variables.vBden);
else
fprintf('WARNING: No matching data file found!\n');
continue;
end

% Apply zero-phase low-pass filter
filterPasses = 2;
filterCorrectedCutOff = filterCutOff / ((2^(1 / filterPasses) - 1)^(1 / 4));
[filterB, filterA] = butter(filterHalfOrder, (filterCorrectedCutOff / (0.5 * frameRate)));
for signalIndex = 1:length(signals)
variables.(signals{signalIndex}) = filtfilt(filterB, filterA, variables.(signals{signalIndex}));
end

% Synchronize force and motion measurements
originalSamplePoints = 0:(1 / frameRate):((frames - 1) / frameRate);
compensatedSamplePoints = originalSamplePoints - (intercept + gradient * originalSamplePoints);
for signalIndex = 1:length(signals)
variables.(signals{signalIndex}) = interp1(compensatedSamplePoints, variables.(signals{signalIndex}), originalSamplePoints, 'pchip', 0);
end

% Compute forces with transforming the coordinate system
Fx = variables.Fy1 + variables.Fy2;
Fx1 = variables.Fy1;
Fx2 = variables.Fy2;
Fy_L = variables.FzL1 + variables.FzL2 + variables.FzL3 + variables.FzL4;
Fy_L1 = variables.FzL1;
Fy_L2 = variables.FzL2;
Fy_L3 = variables.FzL3;
Fy_L4 = variables.FzL4;
Fy_R = variables.FzR1 + variables.FzR2 + variables.FzR3 + variables.FzR4;
Fy_R1 = variables.FzR1;
Fy_R2 = variables.FzR2;
Fy_R3 = variables.FzR3;
Fy_R4 = variables.FzR4;
Fz = variables.Fx1 + variables.Fx2;
Fz1 = variables.Fx1;
Fz2 = variables.Fx2;

% Compensate force sensor offset and drift by linear regression
compensationFx = [ones(16000, 1), [1:8000, (frames - 7999):frames]'] \ [Fx(1:8000), Fx((frames - 7999):frames)]';
compensationFx1 = [ones(16000, 1), [1:8000, (frames - 7999):frames]'] \ [Fx1(1:8000), Fx1((frames - 7999):frames)]';
compensationFx2 = [ones(16000, 1), [1:8000, (frames - 7999):frames]'] \ [Fx2(1:8000), Fx2((frames - 7999):frames)]';
compensationFy_L = [ones(16000, 1), [1:8000, (frames - 7999):frames]'] \ [Fy_L(1:8000), Fy_L((frames - 7999):frames)]';
compensationFy_L1 = [ones(16000, 1), [1:8000, (frames - 7999):frames]'] \ [Fy_L1(1:8000), Fy_L1((frames - 7999):frames)]';
compensationFy_L2 = [ones(16000, 1), [1:8000, (frames - 7999):frames]'] \ [Fy_L2(1:8000), Fy_L2((frames - 7999):frames)]';
compensationFy_L3 = [ones(16000, 1), [1:8000, (frames - 7999):frames]'] \ [Fy_L3(1:8000), Fy_L3((frames - 7999):frames)]';
compensationFy_L4 = [ones(16000, 1), [1:8000, (frames - 7999):frames]'] \ [Fy_L4(1:8000), Fy_L4((frames - 7999):frames)]';
compensationFy_R = [ones(16000, 1), [1:8000, (frames - 7999):frames]'] \ [Fy_R(1:8000), Fy_R((frames - 7999):frames)]';
compensationFy_R1 = [ones(16000, 1), [1:8000, (frames - 7999):frames]'] \ [Fy_R1(1:8000), Fy_R1((frames - 7999):frames)]';
compensationFy_R2 = [ones(16000, 1), [1:8000, (frames - 7999):frames]'] \ [Fy_R2(1:8000), Fy_R2((frames - 7999):frames)]';
compensationFy_R3 = [ones(16000, 1), [1:8000, (frames - 7999):frames]'] \ [Fy_R3(1:8000), Fy_R3((frames - 7999):frames)]';
compensationFy_R4 = [ones(16000, 1), [1:8000, (frames - 7999):frames]'] \ [Fy_R4(1:8000), Fy_R4((frames - 7999):frames)]';
compensationFz = [ones(16000, 1), [1:8000, (frames - 7999):frames]'] \ [Fz(1:8000), Fz((frames - 7999):frames)]';
compensationFz1 = [ones(16000, 1), [1:8000, (frames - 7999):frames]'] \ [Fz1(1:8000), Fz1((frames - 7999):frames)]';
compensationFz2 = [ones(16000, 1), [1:8000, (frames - 7999):frames]'] \ [Fz2(1:8000), Fz2((frames - 7999):frames)]';
for dataIndex = 1:frames
Fx(dataIndex) = Fx(dataIndex) - (compensationFx(2) * dataIndex + compensationFx(1));
Fx1(dataIndex) = Fx1(dataIndex) - (compensationFx1(2) * dataIndex + compensationFx1(1));
Fx2(dataIndex) = Fx2(dataIndex) - (compensationFx2(2) * dataIndex + compensationFx2(1));
Fy_L(dataIndex) = Fy_L(dataIndex) - (compensationFy_L(2) * dataIndex + compensationFy_L(1));
Fy_L1(dataIndex) = Fy_L1(dataIndex) - (compensationFy_L1(2) * dataIndex + compensationFy_L1(1));
Fy_L2(dataIndex) = Fy_L2(dataIndex) - (compensationFy_L2(2) * dataIndex + compensationFy_L2(1));
Fy_L3(dataIndex) = Fy_L3(dataIndex) - (compensationFy_L3(2) * dataIndex + compensationFy_L3(1));
Fy_L4(dataIndex) = Fy_L4(dataIndex) - (compensationFy_L4(2) * dataIndex + compensationFy_L4(1));
Fy_R(dataIndex) = Fy_R(dataIndex) - (compensationFy_R(2) * dataIndex + compensationFy_R(1));
Fy_R1(dataIndex) = Fy_R1(dataIndex) - (compensationFy_R1(2) * dataIndex + compensationFy_R1(1));
Fy_R2(dataIndex) = Fy_R2(dataIndex) - (compensationFy_R2(2) * dataIndex + compensationFy_R2(1));
Fy_R3(dataIndex) = Fy_R3(dataIndex) - (compensationFy_R3(2) * dataIndex + compensationFy_R3(1));
Fy_R4(dataIndex) = Fy_R4(dataIndex) - (compensationFy_R4(2) * dataIndex + compensationFy_R4(1));
Fz(dataIndex) = Fz(dataIndex) - (compensationFz(2) * dataIndex + compensationFz(1));
Fz1(dataIndex) = Fz1(dataIndex) - (compensationFz1(2) * dataIndex + compensationFz1(1));
Fz2(dataIndex) = Fz2(dataIndex) - (compensationFz2(2) * dataIndex + compensationFz2(1));
end
clear compensationFx compensationFx1 compensationFx2 compensationFy_L compensationFy_L1 compensationFy_L2 compensationFy_L3 compensationFy_L4 compensationFy_R compensationFy_R1 compensationFy_R2 compensationFy_R3 compensationFy_R4 compensationFz compensationFz1 compensationFz2 dataIndex

% Save processed data
fprintf('STATUS: Saving dataset %s %s.\n', subject, dataset);
force.frameRate = frameRate;
force.frames = frames;
force.treadmillVelocity = variables.vBden;
force.groundReactionForceX = Fx;
force.forceSensorX1 = Fx1;
force.forceSensorX2 = Fx2;
force.groundReactionForceY_L = Fy_L;
force.forceSensorY_L1 = Fy_L1;
force.forceSensorY_L2 = Fy_L2;
force.forceSensorY_L3 = Fy_L3;
force.forceSensorY_L4 = Fy_L4;
force.groundReactionForceY_R = Fy_R;
force.forceSensorY_R1 = Fy_R1;
force.forceSensorY_R2 = Fy_R2;
force.forceSensorY_R3 = Fy_R3;
force.forceSensorY_R4 = Fy_R4;
force.groundReactionForceZ = Fz;
force.forceSensorZ1 = Fz1;
force.forceSensorZ2 = Fz2;
if exist(saveFile, 'file') == 2
save(saveFile, 'force', '-append');
else
save(saveFile, 'force');
end
end
end
Loading

0 comments on commit c32ffbe

Please sign in to comment.