Creating events from continuous data

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

[John D'Errico]

You did not provide the data. So it may be diffficult for someone to help you in filtering data they cannot see or use.

Answer 1

[Mathieu NOE]

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   

[Mathieu NOE]

hello

problem solved?

[Perri Johnson]

@Mathieu NOE Thanks for the comment. If you are familiar with how to how I could Normalize each step to 101 cells, I could use some help there. Since each step has a different number of data points, I can't effectively average each step together.

[Mathieu NOE]

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