Creating events from continuous data

2 ビュー (過去 30 日間)
Perri Johnson
Perri Johnson 2021 年 12 月 9 日
コメント済み: Mathieu NOE 2021 年 12 月 17 日
Hi,
I have a data set that contains force over time duing walking trials. Trying to figure out how to create an event or trigger that indicates heel strike (for instance when the force value reaches 25. Will do the same for the force value falling below 25 to indicate toe off). Have tried using the find function but I need a function that will help me seperate the x number of steps during the stance phase so I can later average each step against each other. I've attached the script I have so far.
Appreciate any advise that anyone could offer for how I can go about this. I have also had trouble filtering the data if anyone is familiar with that command as well?
Thank you in advance.
%This code is designed to import and filter loadsol insole force data
clear all;
close all;
clc;
%Establish the file path for data folder
DataFolder='E:\BME Program\UTHSC Lab work\DAO project\Data processing\';
%Establishing file path for subject folders
subj_folder = ([DataFolder 'Data\']);
%Establishing file path for output folder
OutputFolder = ([DataFolder 'Output\']);
%Define the subject data you're interested in - displays oddly, need to fix
prompt = {'Enter subject number(ex. S1)'};
dlgtitle = 'Subject Number';
definput = {'S1'};
dims = [1 2];
subj = inputdlg(prompt,dlgtitle,dims,definput);
subj = char(subj);
%define subject weight and convert to newtons
%prompt = {'Enter subject weight(lbs)'};
%dlgtitle = 'Subject Weight(lbs)';
%definput = {' '};
%dims = [1 2];
%lbs_weight = inputdlg(prompt,dlgtitle,dims,definput);
%weight = lbs_weight(1,1);
%N_weight = char(weight);
cond_list = {'C1'};
for i_cond = 1:length(cond_list)
cond = char(cond_list(i_cond));
%Identify the files in the subject folder
InFolder = [subj_folder subj '\'];
FileList = dir([InFolder '*.txt']);
for i_file = 1:length(FileList)
clearvars raw_data alldata data raw_time raw_Rforce
%identify filename & create input path/filename
filename = char(FileList(i_file).name);
input_file = [subj_folder subj '\' filename];
%read raw data to array
raw_data = readmatrix(input_file);
alldata = raw_data(5:end,2:5);
%resample data to 240Hz to match force treadmill GRF sampling rate
data = resample(alldata,24,10);
raw_time = data(:,1);
raw_Rforce = data(:,4);
%threshold_data = find(raw_Rforce>25)
%plot raw data
plot(raw_time, raw_Rforce);
title([subj 'Loadsol Force-Time Graph']);
xlabel('Time(s)');
ylabel('Force(N)');
legend('Right insole force data');
%yline (25, 'Color', 'r', 'LineWidth',2);
hold on
%Lowpass filter the force data using 25Hz similar to V3D code
%time = lowpass(raw_time, 25, 240);
%Rforce = lowpass(raw_Rforce, 25, 240);
outputfile = [OutputFolder subj '\' 'LS' subj 'C' char(string(i_cond)) '_Processed.xlsx'];
header1 = {'subject weight(lbs)'};
header2 = {'Time(sec)' 'Right Insole Force(N)'};
%output array
outputArray = [raw_time raw_Rforce];
%write data to excel
xlswrite(outputfile,header1,cond,'A1');
xlswrite(outputfile,lbs_weight,cond,'B1');
xlswrite(outputfile,header2,cond,'A3');
xlswrite(outputfile,outputArray,cond,'A4');
end %file list
end %condition list
  1 件のコメント
John D'Errico
John D'Errico 2021 年 12 月 9 日
You did not provide the data. So it may be diffficult for someone to help you in filtering data they cannot see or use.

サインインしてコメントする。

回答 (1 件)

Mathieu NOE
Mathieu NOE 2021 年 12 月 13 日
hello
I took only the central portion of your code and applied my recipe - see below
had to do some preliminary data processing before doing the "threshold crossing" stuff
(remove NaN's , get a time vector with constant increment , smooth the force signal...)
code
clc; clearvars; close all
%read raw data to array
raw_data = readmatrix('S1C1.txt');
alldata = raw_data(5:end,2:5);
%resample data to 240Hz to match force treadmill GRF sampling rate
data = resample(alldata,24,10);
raw_time = data(:,1);
raw_Rforce = data(:,4);
% remove NaN
ind_nan = isnan(raw_time);
raw_time(ind_nan) = [];
raw_Rforce(ind_nan) = [];
% resample (interp1) the data to get a truly uniformly increasing
% time vector
samples = length(raw_time);
dt = (max(raw_time) - min(raw_time))/(samples-1);
%time = linspace(min(raw_time),max(raw_time),samples);
time = (0:samples-1)'*dt;
Rforce = interp1(raw_time,raw_Rforce,time);
% add a bit of smoothing because signal exibit sometimes some zig
% zags
Rforce = smoothdata(Rforce,'gaussian',10);
%% my code below
threshold = 50; % your value here
[t0_pos,s0_pos,t0_neg,s0_neg]= crossing_V7(Rforce,time,threshold,'linear'); % positive (pos) and negative (neg) slope crossing points
% ind => time index (samples)
% t0 => corresponding time (x) values
% s0 => corresponding function (y) values , obviously they must be equal to "threshold"
figure(4)
plot(time,Rforce,time,threshold*ones(size(time)),'k--',t0_pos,s0_pos,'dr',t0_neg,s0_neg,'dg','linewidth',2,'markersize',12);grid on
legend('signal','threshold','positive slope crossing points','negative slope crossing points');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [t0_pos,s0_pos,t0_neg,s0_neg] = crossing_V7(S,t,level,imeth)
% [ind,t0,s0,t0close,s0close] = crossing_V6(S,t,level,imeth,slope_sign) % older format
% CROSSING find the crossings of a given level of a signal
% ind = CROSSING(S) returns an index vector ind, the signal
% S crosses zero at ind or at between ind and ind+1
% [ind,t0] = CROSSING(S,t) additionally returns a time
% vector t0 of the zero crossings of the signal S. The crossing
% times are linearly interpolated between the given times t
% [ind,t0] = CROSSING(S,t,level) returns the crossings of the
% given level instead of the zero crossings
% ind = CROSSING(S,[],level) as above but without time interpolation
% [ind,t0] = CROSSING(S,t,level,par) allows additional parameters
% par = {'none'|'linear'}.
% With interpolation turned off (par = 'none') this function always
% returns the value left of the zero (the data point thats nearest
% to the zero AND smaller than the zero crossing).
%
% check the number of input arguments
error(nargchk(1,4,nargin));
% check the time vector input for consistency
if nargin < 2 | isempty(t)
% if no time vector is given, use the index vector as time
t = 1:length(S);
elseif length(t) ~= length(S)
% if S and t are not of the same length, throw an error
error('t and S must be of identical length!');
end
% check the level input
if nargin < 3
% set standard value 0, if level is not given
level = 0;
end
% check interpolation method input
if nargin < 4
imeth = 'linear';
end
% make row vectors
t = t(:)';
S = S(:)';
% always search for zeros. So if we want the crossing of
% any other threshold value "level", we subtract it from
% the values and search for zeros.
S = S - level;
% first look for exact zeros
ind0 = find( S == 0 );
% then look for zero crossings between data points
S1 = S(1:end-1) .* S(2:end);
ind1 = find( S1 < 0 );
% bring exact zeros and "in-between" zeros together
ind = sort([ind0 ind1]);
% and pick the associated time values
t0 = t(ind);
s0 = S(ind);
if ~isempty(ind)
if strcmp(imeth,'linear')
% linear interpolation of crossing
for ii=1:length(t0)
%if abs(S(ind(ii))) >= eps(S(ind(ii))) % MATLAB V7 et +
if abs(S(ind(ii))) >= eps*abs(S(ind(ii))) % MATLAB V6 et - EPS * ABS(X)
% interpolate only when data point is not already zero
NUM = (t(ind(ii)+1) - t(ind(ii)));
DEN = (S(ind(ii)+1) - S(ind(ii)));
slope = NUM / DEN;
slope_sign(ii) = sign(slope);
t0(ii) = t0(ii) - S(ind(ii)) * slope;
s0(ii) = level;
end
end
end
% extract the positive slope crossing points
ind_pos = find(sign(slope_sign)>0);
t0_pos = t0(ind_pos);
s0_pos = s0(ind_pos);
% extract the negative slope crossing points
ind_neg = find(sign(slope_sign)<0);
t0_neg = t0(ind_neg);
s0_neg = s0(ind_neg);
else
% empty output
ind_pos = [];
t0_pos = [];
s0_pos = [];
% extract the negative slope crossing points
ind_neg = [];
t0_neg = [];
s0_neg = [];
end
end
  3 件のコメント
Mathieu NOE
Mathieu NOE 2021 年 12 月 17 日
hello
see below - this is the avrage curve of all events (dotted black line)
clc; clearvars; close all
%read raw data to array
raw_data = readmatrix('S1C1.txt');
alldata = raw_data(5:end,2:5);
%resample data to 240Hz to match force treadmill GRF sampling rate
data = resample(alldata,24,10);
raw_time = data(:,1);
raw_Rforce = data(:,4);
% remove NaN
ind_nan = isnan(raw_time);
raw_time(ind_nan) = [];
raw_Rforce(ind_nan) = [];
% resample (interp1) the data to get a truly uniformly increasing
% time vector
samples = length(raw_time);
dt = (max(raw_time) - min(raw_time))/(samples-1);
%time = linspace(min(raw_time),max(raw_time),samples);
time = (0:samples-1)'*dt;
Rforce = interp1(raw_time,raw_Rforce,time);
% add a bit of smoothing because signal exibit sometimes some zig
% zags
Rforce = smoothdata(Rforce,'gaussian',10);
%% my code below
threshold = 25; % your value here
[t0_pos,s0_pos,t0_neg,s0_neg]= crossing_V7(Rforce,time,threshold,'linear'); % positive (pos) and negative (neg) slope crossing points
% ind => time index (samples)
% t0 => corresponding time (x) values
% s0 => corresponding function (y) values , obviously they must be equal to "threshold"
figure(4)
plot(time,Rforce,time,threshold*ones(size(time)),'k--',t0_pos,s0_pos,'dr',t0_neg,s0_neg,'dg','linewidth',2,'markersize',12);grid on
legend('signal','threshold','positive slope crossing points','negative slope crossing points');
%% now creating a loop for averaging all events
% tip : we are going to resample with the same time length all patterns
% define by the start / stop indexes computed above
m = min(length(t0_pos),length(t0_neg));
nb_points = floor(length(Rforce)/m); % nb of points between start / stop indexes (average)
for ci = 1:m
% extract data
new_time = linspace(t0_pos(ci),t0_neg(ci),nb_points);
data_extrc(ci,:) = interp1(time,Rforce,new_time);
new_time = new_time-new_time(1); % so new_time starts always at zero;
end
% and we can now make the averaged curve
Rforce_all_avg = mean(data_extrc,1);
figure(5)
plot(new_time,data_extrc,'linewidth',2);hold on
plot(new_time,Rforce_all_avg,'k--','linewidth',4); hold off
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [t0_pos,s0_pos,t0_neg,s0_neg] = crossing_V7(S,t,level,imeth)
% [ind,t0,s0,t0close,s0close] = crossing_V6(S,t,level,imeth,slope_sign) % older format
% CROSSING find the crossings of a given level of a signal
% ind = CROSSING(S) returns an index vector ind, the signal
% S crosses zero at ind or at between ind and ind+1
% [ind,t0] = CROSSING(S,t) additionally returns a time
% vector t0 of the zero crossings of the signal S. The crossing
% times are linearly interpolated between the given times t
% [ind,t0] = CROSSING(S,t,level) returns the crossings of the
% given level instead of the zero crossings
% ind = CROSSING(S,[],level) as above but without time interpolation
% [ind,t0] = CROSSING(S,t,level,par) allows additional parameters
% par = {'none'|'linear'}.
% With interpolation turned off (par = 'none') this function always
% returns the value left of the zero (the data point thats nearest
% to the zero AND smaller than the zero crossing).
%
% check the number of input arguments
error(nargchk(1,4,nargin));
% check the time vector input for consistency
if nargin < 2 | isempty(t)
% if no time vector is given, use the index vector as time
t = 1:length(S);
elseif length(t) ~= length(S)
% if S and t are not of the same length, throw an error
error('t and S must be of identical length!');
end
% check the level input
if nargin < 3
% set standard value 0, if level is not given
level = 0;
end
% check interpolation method input
if nargin < 4
imeth = 'linear';
end
% make row vectors
t = t(:)';
S = S(:)';
% always search for zeros. So if we want the crossing of
% any other threshold value "level", we subtract it from
% the values and search for zeros.
S = S - level;
% first look for exact zeros
ind0 = find( S == 0 );
% then look for zero crossings between data points
S1 = S(1:end-1) .* S(2:end);
ind1 = find( S1 < 0 );
% bring exact zeros and "in-between" zeros together
ind = sort([ind0 ind1]);
% and pick the associated time values
t0 = t(ind);
s0 = S(ind);
if ~isempty(ind)
if strcmp(imeth,'linear')
% linear interpolation of crossing
for ii=1:length(t0)
%if abs(S(ind(ii))) >= eps(S(ind(ii))) % MATLAB V7 et +
if abs(S(ind(ii))) >= eps*abs(S(ind(ii))) % MATLAB V6 et - EPS * ABS(X)
% interpolate only when data point is not already zero
NUM = (t(ind(ii)+1) - t(ind(ii)));
DEN = (S(ind(ii)+1) - S(ind(ii)));
slope = NUM / DEN;
slope_sign(ii) = sign(slope);
t0(ii) = t0(ii) - S(ind(ii)) * slope;
s0(ii) = level;
end
end
end
% extract the positive slope crossing points
ind_pos = find(sign(slope_sign)>0);
t0_pos = t0(ind_pos);
s0_pos = s0(ind_pos);
% extract the negative slope crossing points
ind_neg = find(sign(slope_sign)<0);
t0_neg = t0(ind_neg);
s0_neg = s0(ind_neg);
else
% empty output
ind_pos = [];
t0_pos = [];
s0_pos = [];
% extract the negative slope crossing points
ind_neg = [];
t0_neg = [];
s0_neg = [];
end
end

サインインしてコメントする。

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by