Can someone help me with codes or function combination to evaluate the Average force from tdms file with four channels? Pls, see codes I applied "manually" to attempt it.

#1 codes used to extract the tdms file

data1 = tdmsread("20221128131019RUN1force.tdms");

data1 = tdmsread("20221128131019RUN1force.tdms", ...

ChannelGroupName = "Force", ...

ChannelNames = ["Feed force" "Main force" "Passive force" "AERMS"], ...

RowTimes = "ArrTime");

#2 codes used to separate each forces from the channels eg for "Main force" - Cutting force

>> plot(data{1, 1}.(2))

>> plot(data1{1, 1}.(2))

>> main_force1 = data1{1, 1}.(2);

# codes used to compute the average Forces after picking the points Ranges for Fc1, Fc2 and Fc3 as shown the Main force plot figure at the bottom of codes;

>> sep_m1b=main_force1(950000:1592000);

>> sep_Fc1=main_force1(950000:1592000);

>> mean(sep_Fc1)

ans =

-75.6126 % the negative values shows the direction of the force due to position of dynamometer

>> sep_Fc1=main_force1(1741000:1971000);

>> sep_Fc2=main_force1(1741000:1971000);

>> mean(sep_Fc2)

ans =

-63.5206

>> sep_Fc3=main_force1(2051000:2221000);

>> mean(sep_Fc3)

ans =

-64.0661 75.6,

#Question : Can someone assist or suggest codes/functions that will COMBINE codes #1, #2 and #3, such that it Reads the 20221128131019RUN1force.tdms file, Extract any/each of the forces (Feed, Main, Passive or AERMS), Evaluate the Average values from the range and give the values 75.6N, 63.5N and 64.1N respectively skipping the ROUGH JUMPS and NON-cutting signal zones. The link to the tdms file, I have provided in the link. https://drive.google.com/file/d/1as21ZQWwY5RsKSNISNI8m761oG84i89I/view?usp=sharing

[Mathieu NOE]

hello

not everyone (like me) has the DAQ Toolbox, so it's better to share the data as mat file

the FEX submission did not help me either as you have raw TDMS file format which is unsupported by this fex submission

TDMS Reader - File Exchange - MATLAB Central (mathworks.com)

by sharing the data in matlab format more people on this forum may answer you

[Ugifada]

Thank you @Mathieu NOE for you kind advice. I have attached the mat file now. Thanks

[Mathieu NOE]

hello

no , I don't see any mat file here

you attached a m file that contains the beginning or your code (that we know already from your post above)

[Ugifada]

Oh my error. However, I tried to upload the mat file but is very large (39MB) even as a zip or .rar file exceeding the max of 5MB. So I shared with the link from google drive. Can this be helpful.? https://drive.google.com/file/d/1_ooEOQGzYIj_rO88AdiqWKGeV9t77XsQ/view?usp=sharing

Answer 1

[Mathieu NOE]

Ok so this is what I can offer you today

1/ I noticed your sampling rate could be much lower (at least factor 1000) without any loss of information !! so this could reduce your file size and further processing time by large amount

2/ the rest is quite simple : we look for the positive and negative transitions and use that time stamps to extract the "flat" sections of data , and take the mean value of them . a bit of smoothing is applied before derivation

last , I prefer peakseek over the built in findpeaks but it's everyone taste so pick either one or the other....

results :

force_mean = -75.6448 -63.6153 -64.4441

plot

code

load('ExpRUN1.mat');
ForceData = table2array(ForceData{1,1});  % ["Feed force" "Main force" "Passive force" "AERMS"], ...

mainforce = ForceData(:,2);
samples = numel(mainforce);
t = (1:samples);

%% task 1 : decimate (because we don't need that much data !!)
decim = 1000;
mainforce2 = decimate(mainforce,decim);
samples2 = numel(mainforce2);
t2 = (decim/2:decim:samples);

% figure(1),
% plot(t,mainforce,t2,mainforce2);

%% task 2 : find start / stop indexes for flat signal areas
clear ForceData mainforce samples t  % we don't need the hude initial data anymore so let's clear some memory

% smooth the signal 
mainforce2s = smoothdata(mainforce2,'gaussian',50);
% derivation for finding transients
d = gradient(mainforce2s);


% positive peaks
minpeakdist = 100;
minpeakh = max(d)/3;
[locspos, pkspos]=peakseek(d,minpeakdist,minpeakh);
t2pos = t2(locspos);

% negative peaks
minpeakdist = 100;
minpeakh = -min(d)/3;
[locsneg, pksneg]=peakseek(-d,minpeakdist,minpeakh);
t2neg = t2(locsneg);
pksneg = -pksneg;

figure(2),
subplot(2,1,1),plot(t2,mainforce2,t2,mainforce2s);
hold on 

for k = 1:numel(pksneg)
    ind_start = locsneg(k);
    ind_stop = locspos(k);
    length=(ind_stop-ind_start);
    samples2remove = round(0.15*length); % remove this amount of samples at the beginning and end of this data buffer
    ind_start = ind_start+samples2remove;
    ind_stop = ind_stop-samples2remove;   
    %duration(k)=(ind_stop-ind_start);
    
    time = t2(ind_start:ind_stop);
    data = mainforce2s(ind_start:ind_stop);

    subplot(2,1,1),plot(time,data,'*-');
    time_mean(k) = mean(time);
    force_mean(k) = mean(data);

end
subplot(2,1,1),plot(time_mean,force_mean,'dr','Markersize',15);
legend('signal decimated','signal decimated and smoothed','data segment 1','data segment 2','data segment 3','force mean')

subplot(2,1,2),plot(t2,d,t2pos,pkspos,'dg',t2neg,pksneg,'dk','Markersize',10);
legend('signal derivative','positive peaks','negative peaks')

hold off 

force_mean

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [locs, pks]=peakseek(x,minpeakdist,minpeakh)
% x is a vector input (generally a timecourse)
% minpeakdist is the minimum desired distance between peaks (optional, defaults to 1)
% minpeakh is the minimum height of a peak (optional)
%
% (c) 2010
% Peter O'Connor
% peter<dot>ed<dot>oconnor .AT. gmail<dot>com
if size(x,2)==1, x=x'; end
% Find all maxima and ties
locs=find(x(2:end-1)>=x(1:end-2) & x(2:end-1)>=x(3:end))+1;
if nargin<2, minpeakdist=1; end % If no minpeakdist specified, default to 1.
if nargin>2 % If there's a minpeakheight
    locs(x(locs)<=minpeakh)=[];
end
if minpeakdist>1
    while 1
        del=diff(locs)<minpeakdist;
        if ~any(del), break; end
        pks=x(locs);
        [garb, mins]=min([pks(del) ; pks([false del])]); %#ok<ASGLU>
        deln=find(del);
        deln=[deln(mins==1) deln(mins==2)+1];
        locs(deln)=[];
    end
end
if nargout>1
    pks=x(locs);
end
end

[Ugifada]

Thank you so much @Mathieu NOE, sincerely appreciate your expert and kind assistance.

[Mathieu NOE]

My pleasure !

have a good day

[Ugifada]

Good afternoon @Mathieu NOE, I tried the codes force the main force codes and plot (forceEvaluation1.m, Plot_mainforce.jpg) and it is perfectly correct ! However, I tried evaluating the feed force (forceEvaluation2.m, Plot_feedforce.jpg) and passive force (forceEvaluation3.m, Plot_passiveforce.jpg), unfortunately the force_mean returned = NaN NaN NaN in both cases. Please, check my codes to know what I am doing wrong. Thank you

FEED FORCE

force_mean = NaN NaN NaN ???

load ExpRUN1.mat

ForceData = table2array(ForceData{1,1}); % ["Feed force" "Main force" "Passive force" "AERMS"], ...

feedforce = ForceData(:,1);

samples = numel(feedforce);

t = (1:samples);

%% task 1 : decimate (because we don't need that much data !!)

decim = 1000;

feedforce2 = decimate(feedforce,decim);

samples2 = numel(feedforce2);

t2 = (decim/2:decim:samples);

% figure(1),

% plot(t,feedforce,t2,feedforce2);

%% task 2 : find start / stop indexes for flat signal areas

clear ForceData feedforce samples t % we don't need the hude initial data anymore so let's clear some memory

% smooth the signal

feedforce2s = smoothdata(feedforce2,'gaussian',50);

% derivation for finding transients

d = gradient(feedforce2s);

% positive peaks

minpeakdist = 100;

minpeakh = max(d)/3;

[locspos, pkspos]=peakseek(d,minpeakdist,minpeakh);

t2pos = t2(locspos);

% negative peaks

minpeakdist = 100;

minpeakh = -min(d)/3;

[locsneg, pksneg]=peakseek(-d,minpeakdist,minpeakh);

t2neg = t2(locsneg);

pksneg = -pksneg;

figure(2),

subplot(2,1,1),plot(t2,feedforce2,t2,feedforce2s);

hold on

for k = 1:numel(pksneg)

ind_start = locsneg(k);

ind_stop = locspos(k);

length=(ind_stop-ind_start);

samples2remove = round(0.15*length); % remove this amount of samples at the beginning and end of this data buffer

ind_start = ind_start+samples2remove;

ind_stop = ind_stop-samples2remove;

%duration(k)=(ind_stop-ind_start);

time = t2(ind_start:ind_stop);

data = feedforce2s(ind_start:ind_stop);

subplot(2,1,1),plot(time,data,'*-');

time_mean(k) = mean(time);

force_mean(k) = mean(data);

end

subplot(2,1,1),plot(time_mean,force_mean,'dr','Markersize',15);

legend('signal decimated','signal decimated and smoothed','data segment 1','data segment 2','data segment 3','force mean')

subplot(2,1,2),plot(t2,d,t2pos,pkspos,'dg',t2neg,pksneg,'dk','Markersize',10);

legend('signal derivative','positive peaks','negative peaks')

hold off

force_mean

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [locs, pks]=peakseek(x,minpeakdist,minpeakh)

% x is a vector input (generally a timecourse)

% minpeakdist is the minimum desired distance between peaks (optional, defaults to 1)

% minpeakh is the minimum height of a peak (optional)

%

% (c) 2010

% Peter O'Connor

% peter<dot>ed<dot>oconnor .AT. gmail<dot>com

if size(x,2)==1, x=x'; end

% Find all maxima and ties

locs=find(x(2:end-1)>=x(1:end-2) & x(2:end-1)>=x(3:end))+1;

if nargin<2, minpeakdist=1; end % If no minpeakdist specified, default to 1.

if nargin>2 % If there's a minpeakheight

locs(x(locs)<=minpeakh)=[];

end

if minpeakdist>1

while 1

del=diff(locs)<minpeakdist;

if ~any(del), break; end

pks=x(locs);

[garb, mins]=min([pks(del) ; pks([false del])]); %#ok<ASGLU>

deln=find(del);

deln=[deln(mins==1) deln(mins==2)+1];

locs(deln)=[];

end

end

if nargout>1

pks=x(locs);

end

end

[Mathieu NOE]

hello

this is an improved code that works for all 4 cases

there is no modification to do , except change which channel you want to proceed

force = ForceData(:,1); % choose channel 1 / 2 / 3 / 4 corresponding to ["Feed force" "Main force" "Passive force" "AERMS"]

new code :

load('ExpRUN1.mat');
ForceData = table2array(ForceData{1,1});  % ["Feed force" "Main force" "Passive force" "AERMS"], ...

force = ForceData(:,1); % choose channel 1 / 2 / 3 / 4 corresponding to ["Feed force" "Main force" "Passive force" "AERMS"]
samples = numel(force);
t = (1:samples);

%% task 1 : decimate (because we don't need that much data !!)
decim = 1000;
force2 = decimate(force,decim);
samples2 = numel(force2);
t2 = (decim/2:decim:samples);

% figure(1),
% plot(t,mainforce,t2,mainforce2);

%% task 2 : find start / stop indexes for flat signal areas
clear ForceData mainforce samples t  % we don't need the hude initial data anymore so let's clear some memory

% smooth the signal 
force2s = smoothdata(force2,'gaussian',50);
% derivation for finding transients
d = gradient(force2s);

% positive peaks
minpeakdist = 100;
minpeakh = max(d)/10;
[locspos, pkspos]=peakseek(d,minpeakdist,minpeakh);
t2pos = t2(locspos);

% negative peaks
minpeakdist = 100;
minpeakh = -min(d)/10;
[locsneg, pksneg]=peakseek(-d,minpeakdist,minpeakh);
t2neg = t2(locsneg);
pksneg = -pksneg;

figure(2),
% subplot(2,1,1),
plot(t2,force2,t2,force2s);
hold on 

% init some variables
duration_min = 100; % minimal duration in samples (short signal are not considered valid)
time_mean = [];
force_mean = [];

for k = 1:numel(pksneg)
    % check if positive or negative peak comes first
    if locsneg(1) < locspos(1)
        ind_start = locsneg(k);
        ind_stop = locspos(k);
    else
        ind_start = locspos(k);
        ind_stop = locsneg(k);
    end
    
    length=(ind_stop-ind_start);
    samples2remove = round(0.15*length); % remove this amount of samples at the beginning and end of this data buffer
    ind_start = ind_start+samples2remove;
    ind_stop = ind_stop-samples2remove;   
    duration(k)=(ind_stop-ind_start);
    
    if duration(k)>= duration_min % only longer buffer are valid
    
    time = t2(ind_start:ind_stop);
    data = force2s(ind_start:ind_stop);

%     subplot(2,1,1),
    plot(time,data,'*-');
    time_mean = [time_mean mean(time)];
    force_mean = [force_mean mean(data)];
    
    end

end
% subplot(2,1,1),
plot(time_mean,force_mean,'dr','Markersize',15);
legend('signal decimated','signal decimated and smoothed','data segment 1','data segment 2','data segment 3','force mean')

% subplot(2,1,2),plot(t2,d,t2pos,pkspos,'dg',t2neg,pksneg,'dk','Markersize',10);
% legend('signal derivative','positive peaks','negative peaks')
% 
hold off 

force_mean

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [locs, pks]=peakseek(x,minpeakdist,minpeakh)
% x is a vector input (generally a timecourse)
% minpeakdist is the minimum desired distance between peaks (optional, defaults to 1)
% minpeakh is the minimum height of a peak (optional)
%
% (c) 2010
% Peter O'Connor
% peter<dot>ed<dot>oconnor .AT. gmail<dot>com
if size(x,2)==1, x=x'; end
% Find all maxima and ties
locs=find(x(2:end-1)>=x(1:end-2) & x(2:end-1)>=x(3:end))+1;
if nargin<2, minpeakdist=1; end % If no minpeakdist specified, default to 1.
if nargin>2 % If there's a minpeakheight
    locs(x(locs)<=minpeakh)=[];
end
if minpeakdist>1
    while 1
        del=diff(locs)<minpeakdist;
        if ~any(del), break; end
        pks=x(locs);
        [garb, mins]=min([pks(del) ; pks([false del])]); %#ok<ASGLU>
        deln=find(del);
        deln=[deln(mins==1) deln(mins==2)+1];
        locs(deln)=[];
    end
end
if nargout>1
    pks=x(locs);
end
end

[Ugifada]

@Mathieu NOE . . . This works PERFECTLY 100%, Sir! Much appreciated, especially your patience and prompt replies!

[Mathieu NOE]

my pleasure !!

[Ugifada]

@Mathieu NOE Also, in some cases (experimental RUN) it returns this error message and display no results for the Mean_force (just for a particular channel selection)

>> Mean_force

Index exceeds the number of array elements. Index must not exceed 2.

Error in Mean_force (line 44)

ind_start = locspos(k);

What should I do in such cases ?

[Mathieu NOE]

Can you share the data file with this issue ?

[Ugifada]

@Mathieu NOE ExpRUN2, gave the error message in channel 4 (AERMS), others were perfect. https://drive.google.com/file/d/1cHMuETIKcjSkN6WPo1KG6odspF9AfDu6/view?usp=sharing and ExpRUN3, gave the error message in channel 3 (Passive Force), others were perfect. - https://drive.google.com/file/d/1whWVsmLWURWDV-u1cQaJyyS9ID4B_VIf/view?usp=sharing

[Mathieu NOE]

hello again

ok, so I cahnged a bit the logic to make the code more robust

basically, what we are now doing : look for segments that has a fairly moderate slope (there is a logical test with the derivative of the signal) and also that have a value above a minimal threshold (to not select almost zero force data)

it seems to work pretty well on all 3 files and 4 channels

the only parameter you may have to tweak it the minimal length of the segments, so if you cycles are shorter or longer you may have to change that value to have always the right number of valid segments

hope it helps !

clc
clearvars
% close all

% select file and channel : 
% load('ExpRUN1.mat'); % ok all 4 channels
% duration_min = 10; % minimal duration in samples (shorter signal are not considered valid)

% load('ExpRUN2.mat');  % ok all 4 channels
% duration_min = 20; % minimal duration in samples (shorter signal are not considered valid)

load('ExpRUN3.mat'); % ok all 4 channels
duration_min = 20; % minimal duration in samples (shorter signal are not considered valid)

channel = 4;% choose channel 1 / 2 / 3 / 4 corresponding to ["Feed force" "Main force" "Passive force" "AERMS"]

%% main code 
ForceData = table2array(ForceData{1,1});  % ["Feed force" "Main force" "Passive force" "AERMS"], ...
force = ForceData(:,channel); 
samples = numel(force);
t = (1:samples);

%% task 1 : decimate (because we don't need that much data !!)
decim = 5000;% init some variables

force2 = decimate(force,decim);
samples2 = numel(force2);
t2 = (round(decim/2):decim:samples);

%% task 2 : find start / stop indexes for flat signal areas
clear ForceData mainforce samples t  % we don't need the hude initial data anymore so let's clear some memory

% smooth the signal 
force2s = smoothdata(force2,'gaussian',15);
% derivation for finding transients
d = gradient(force2s);

% find segments of certain legth and certain slope (d)
ms = max(abs(d));  % max absolute slope
ind_slope = abs(d)< ms/10; % check if factor 10 is enough or too high

% find segment of non zero force 
mf = max(abs(force2s));  % max  force
ind_force = abs(force2s)> mf/10; % check if factor 10 is enough or too high

% combine both conditions 
ind = ind_slope&ind_force;

figure(1),
% subplot(2,1,1),
plot(t2,force2,t2,force2s);
hold on

% subplot(2,1,2),plot(t2,d,t2,double(ind_slope),t2,double(ind_force),t2,double(ind));
% legend('signal derivative','ind slope','ind force','ind')
% subplot(2,1,2),plot(t2,d,t2,double(ind));
% legend('signal derivative','ind')

[begin,ends] = find_start_end_group(ind);
length=(ends-begin);
id = (length>duration_min);
begin = begin(id);
ends = ends(id);


for k = 1:numel(begin)
    % check if positive or negative peak comes first
    ind_start = begin(k);
    ind_stop = ends(k);
    
    samples2remove = round(0.1*length(k)); % remove this amount of samples at the beginning and end of this data buffer
    ind_start = ind_start+samples2remove;
    ind_stop = ind_stop-samples2remove;   
    duration(k)=(ind_stop-ind_start);
    
    time = t2(ind_start:ind_stop);
    data = force2s(ind_start:ind_stop);

%     subplot(2,1,1),
    plot(time,data,'*-');
    time_mean(k) = mean(time);
    force_mean(k) = mean(data);
%     leg_str{k} = ['Segment' num2str(n)];
    


end
% % subplot(2,1,1),
plot(time_mean,force_mean,'dr','Markersize',15);
legend('signal decimated','signal decimated and smoothed')
title(['Force mean = ' num2str(force_mean,3), ' N']);

% subplot(2,1,2),plot(t2,d,t2,double(ind));
% legend('signal derivative','ind')
% 
hold off 

force_mean

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [begin,ends] = find_start_end_group(ind)
% This locates the beginning /ending points of data groups
% Important : ind must be a LOGICAL array

    D = diff([0;ind(:);0]);
    begin = find(D == 1);
    ends = find(D == -1) - 1;
end

[Mathieu NOE]

hello again

so does this new code work for you ?

[Ugifada]

Dear @Mathieu NOE, please pardon me, I have been on exams at school and other end of semester tasks. However, I have about 27 RUNS of the experiment and I have tested on RUN1 to RUN4 before the tight engagments. I discovered as you said that I need to vary these particular lines duration time, and segment ....

duration_min = (10-30); % minimal duration in samples (shorter signal are not considered valid)

ind_slope = abs(d)< ms/(10-200); % check if factor 10 is enough or too high

% find segment of non zero force

ind_force = abs(force2s)> mf/(10-200); % check if factor 10 is enough or too high

% combine both conditions

For instance in RUN4 ( channel 3 was giving mean force = 116 127 NaN ) but after adjusting I got value 116, although when manually one is 204. Similarly(channel 2 gave -34.7770 -38.7909 NaN) after adjusting a got value 116 but when manually done -80N.

Therefore, I would say the last codes you provided seems to be perfect but ONE NEED TO CAREFULLY ADJUST TIME AND SEGMENTS, and I will complete the remaining RUN5-27 and give a feedback once am done with it.

Thank you so much, and I appreciate the follow up! You are a good teacher.

[Mathieu NOE]

hello again

hope you have succeeded your exams !

I was just wondering if this new approach was better , more robust than the previous one

yes for sure, you may have to tweak some parameters according to your data

as soon as the main code structure is fine, that's a good point and i am glad I could help you on that

good luck for the future !