Splitting a vector into separate vectors using thresholds

Hello,

So, I have data from a shoe insole that collects ground reaction force during running. Basically, i get one continous vector that is composed of multiple steps. I am trying to separate this vector into multiple separate vectors based off of when foot contact occurs (when force >= 20) and when toe-off occurs (when force <= 10). For the sake of understanding, if I have a vector x: x = [0 10 20 30 40 50 40 30 20 10 0 10 20 30 40 50 40 30 20 10 0], and I want to create a loop that says "when force is greater than or equal to 20, create a vector and end that vector when force is less than or equal to 10. Then when force is greater than or equal to 20 again, create a new vector and end that vector when the froe is less than or equal to 10. And so on and so forth..".

So far, I have this:

x = [0 10 20 30 40 50 40 30 20 10 0 10 20 30 40 50 40 30 20 10 0]

for i = 1:length(x)

y = find(x>=20)

stance = x(y)

end

But I am having trouble trying to figure out how to separate this into multiple vectors according to the above criteria. Can anyone help out with this?

Answer 1

[Stephen23]

Simpler and more robust:

>> x = [9,15,9,23,15,9,15,7,99,0,12]
x =
     9    15     9    23    15     9    15     7    99     0    12
>> idb = find([x(1)>=20,diff(x>=20)>0]);
>> ide = [find(x<=10)-1,numel(x)];
>> ide = arrayfun(@(b)ide(find(ide>=b,1,'first')),idb);
>> out = arrayfun(@(b,e)x(b:e), idb, ide, 'UniformOutput',false);
>> out{:}
ans =
    23    15
ans =
    99

It correctly detects groups at the start and end of the data vector:

>> x = [20,24,20];
>> idb = find([x(1)>=20,diff(x>=20)>0]);
>> ide = [find(x<=10)-1,numel(x)];
>> ide = arrayfun(@(b)ide(find(ide>=b,1,'first')),idb);
>> out = arrayfun(@(b,e)x(b:e), idb, ide, 'UniformOutput',false);
>> out{:}
ans =
    20    24    20

And with your original data vector:

>> x = [0 10 20 30 40 50 40 30 20 10 0 10 20 30 40 50 40 30 20 10 0];
>> idb = find([x(1)>=20,diff(x>=20)>0]);
>> ide = [find(x<=10)-1,numel(x)];
>> ide = arrayfun(@(b)ide(find(ide>=b,1,'first')),idb);
>> out = arrayfun(@(b,e)x(b:e), idb, ide, 'UniformOutput',false);
>> out{:}
ans =
    20    30    40    50    40    30    20
ans =
    20    30    40    50    40    30    20

[Joseph Gonzales]

Thank you both Stephen. This seems much simpler (not to put Jon's answer down). Again Stephen, I apologize if I misinterpreted anything due to my understanding of matlab at this point. However, I chose this as the accepted answer because it is the most simple and easiest for me to understand. I appreciate it.

Thanks,

Joe

Answer 2

[Jon]

While I was working on this I see you already got an answer from Stephen.

My approach is similar although not as vectorized.

Also I think my definition of the epochs might be a little different (I'd have to analyze Stephen's code a little deeper to see). Not sure if this is how you want it, but I basically say the you are interested in the data from when you go over 20 until it falls back to 10, but not the points between 10 and back up to 20.

Also, I don't assume that the data starts with a foot down.

Anyhow in case this is helpful. Here's another way:

upThresh = 20;
lowThresh = 10;

% make some data just to try out idea
fs = [1:30,29:-2:0];
f = repmat(fs,1,8)
plot(f)

% find start of epochs where foot is down
iDown = f >= upThresh; % goes from zero to one (true to false) at foot down
iStart = find(diff(iDown)==1)

% find end of epoch where foot comes up
iUp = f <= lowThresh
iEnd = find(diff(iUp)== 1)

% align the start and ends properly
if iEnd(1) < iStart(1)
    iEnd = iEnd(2:end);
end

% make start and ends same length 
numEpochs = min(length(iStart),length(iEnd));
iStart = iStart(1:numEpochs);
iEnd = iEnd(1:numEpochs);

% make vectors for each individual epoch
% save in structure since the number of data points may be different for
% each epoch
numEpochs = length(iStart);
epoch(numEpochs).data = 0; % preallocate

for k = 1:numEpochs-1
    % just keep the part that is above the lower threshold
    epoch(k).data = f(iStart(k):iEnd(k))
end

[Joseph Gonzales]

Hi Jon,

Thanks for this. This is a very elegant solution that I don't fully understand due to my ill-understanding of Matlab at this point in my career. It seems to work though! I'm going to go through it in more detail along with documentation so I can attempt to understand it. Thanks again.

[Stephen23]

"This is a very elegant solution..."

Except for the fact that the code is buggy and does not do what the question requested.

1 - Apparently both 19 and 0 are >=20, because that is what all of the leading values are:

>> epoch.data
ans =
   19   20   21   22   23   24   25   26   27   28   29   30   29   27   25   23   21   19   17   15   13   11
ans =
   19   20   21   22   23   24   25   26   27   28   29   30   29   27   25   23   21   19   17   15   13   11
ans =
   19   20   21   22   23   24   25   26   27   28   29   30   29   27   25   23   21   19   17   15   13   11
ans =
   19   20   21   22   23   24   25   26   27   28   29   30   29   27   25   23   21   19   17   15   13   11
ans =
   19   20   21   22   23   24   25   26   27   28   29   30   29   27   25   23   21   19   17   15   13   11
ans =
   19   20   21   22   23   24   25   26   27   28   29   30   29   27   25   23   21   19   17   15   13   11
ans =
   19   20   21   22   23   24   25   26   27   28   29   30   29   27   25   23   21   19   17   15   13   11
ans = 0

2 - It does not return the last group. This is easy to show when there is only one group, e.g.:

>> f = [5:25,24:-1:5];
f =
    5    6    7    8    9   10   11   12   13   14   15   16   17   18   19   20   21   22   23   24   25   24   23   22   21   20   19   18   17   16   15   14   13   12   11   10    9    8    7    6    5

incorrectly returns:

>> epoch.data
ans = 0

3 - It does not detect groups right at the start, instead it throws an error: E.g.:

>> f = [20:2:25,24:-2:9]
f =
    20    22    24    24    22    20    18    16    14    12    10

throws this error:

Attempted to access iStart(1); index out of bounds because numel(iStart)=0.
Error in temp0 (line 20)
if iEnd(1) < iStart(1)

As Jon notes "I don't assume that the data starts with a foot down" which is not quite the whole story because the code can't reliably detect a group at the start anyway!

[Jon]

As Stephen points out I had a few "fence post" errors (index off by 1) where I wasn't careful about how the diff operation shifts things. I think the code below catches those. Otherwise, as Stephen indicates there may be edge cases that I don't handle. My goal was just to give you an idea of how you could work with logical indices and the diff operation to find the segments you want, not to write a piece of bullet proof code.

Of course if there is a nice simple way to use the same idea and also have it cover all cases, as Stephen suggests his code does, than that's even better. Stephen stores his results in a cell array, rather than a structure, which is a good approach, definitely simpler. You could easily modify the code below to use a cell array instead.

In the example below I use a periodic force that maybe looks a little more like what I would think your data would look like (sudden jumps when the foot hits and lifts off). You can also vary where the waveform starts, to see if there are any problems if you start at the wrong point. I tried a few and didn't find any problems, but I didn't thoroughly test it or analyze it.

upThresh = 20;
lowThresh = 10;



% make some data just to try out idea
iBegin = 1; % allow data to begin at variable point in overall cycle
fs = [linspace(25,10,30) 5*ones(1,10)];
f = repmat(fs,1,8)
f = f(iBegin:end);


% find start of epochs where foot is down
iDown = f >= upThresh; % goes from zero to one (true to false) at foot down
iStart = find(diff(iDown)==1)+1

% find end of epoch where foot comes up
iUp = f <= lowThresh
iEnd = find(diff(iUp)== 1);

% align the start and ends properly
if iEnd(1) < iStart(1)
    iEnd = iEnd(2:end);
end

% make start and ends same length 
numEpochs = min(length(iStart),length(iEnd));
iStart = iStart(1:numEpochs);
iEnd = iEnd(1:numEpochs);

% make vectors for each individual epoch
% save in structure since the number of data points may be different for
% each epoch
numEpochs = length(iStart);
% preallocate, after clearing any existing values
% if you don't clear it and already have a bigger structure
% you get left with the remmnants
if exist('epoch','var')
    clear('epoch')
end
epoch(numEpochs).data = 0;

for k = 1:numEpochs-1
    % just keep the part that is above the lower threshold

epoch(k).data = f(iStart(k):iEnd(k))
end

% plot the results
figure
plot(f)
figure
for k = 1:numEpochs
    plot(epoch(k).data)
    hold on
end
hold off

[Joseph Gonzales]

Jon,

Thank you for the clear comments. I am still trying to understand some of the aspects of the above code, but it is starting to make more sense. I apprecaite the post.

Thanks,

Joe