I have a set of data that I need to separate into unequal intervals in order to calculate the RMS value of these data (The data points are continuous over time)

I have a set of torque data and I need to calculate the rms (root-mean-square) torque. The problem is that I have almost 16000 data points with negative and positive values over time, These data must be seperated into intervals first and the length (unequal intervals) of each interval must be found.

The formula I am following is taken from here and can be written as:

What possible code can I use to perform a similar operation but on the full range of data as in the example picture below?

The torque data and time values saved as an Excel file and are attached below.

I know that MATLAB is able to find the rms value for a set of independent data as I read here but my data is continuous over time.

Note: some torque values do not repeat over a period of time and show up once in between two intervals, these can be omitted as I am intrested in the overall average of the intervals.

Thank you in advance.

[dpb]

This certainly isn't clear to me what the rest of the data are and what the result would be for those from this description, sorry.

Attach a section of the real data and show what (and how) the expected result would be from it. Finding a grouping variable wouldn't be any issue, see discretize or findgroups or any of several alternate ways to do that; the Q? is what then, in general? It's not clear to me from what was written above why the answer is not just a constant; the result doesn't seem to depend on the grouping, or I'm simply not recognizing how.

[Mathieu NOE]

can you share your excel file ?

[ABDALLA AL KHALEDI]

So what I am trying to do is seperate the torque data into intervals; however, these intervals are not equally spaced. The reason being that I want to calculate the rms-torque (root-mean-square) value on a continuous time basis according to the formula found here. I know that MATLAB is able to find this value for a set of independent data as I read here but my data are continuous over time. In the picture above I wrote a simplified version of what I wanted to calculate. Note: some torque values do not repeat over a period of time and show up once in between two intervals, these can be omitted as I am intrested in the overall average.

Answer 1

[dpb]

For the solution for all torque levels without collapsing any,

tq=readmatrix('Torque Data.xlsx');
d=[true; diff(tq(:,2)) ~= 0];         % TRUE if values change (down to machine eps)
TQ=tq(d,2);                           % torque values at change points
ix=find([d(:); true]);                % location of changes in logical vector
N = diff(ix);                         % difference in locations --> number
% so calculate the torque rms by group weighted by time
dt=mean(diff(tq(:,1)));               % remove roundoff in dt calculationto get one dt overall
T=dt*(size(tq,1)-1);                  % total T
tq_rms=sqrt(sum(TQ.^2.*(N-1)*dt)/T)

tq_rms = 0.7850

@Jan coded <an FEX submission as a mex function RunLength> that is quite a bit faster if series lengths were to get to be really, really long...and is the root source of the M-code outline above although the general technique has wide applicability in finding crossing, etc., etc, etc., ...

ADDENDUM: The above takes care of the one-sample elements automatically in that the time interval is zero for a count of one...you don't even have to remove them first.

[Mathieu NOE]

Hello @dpb

i just wondered why you have a different result compared to mine (the second is by definition a true rms value whereas the first one is obtained by splitting data in bins)

Torque_rms = 0.885089037584030 (rms "bins")

Torque_rms_alt = 0.885089037584036 (true rms)

the accuracy of my first result depends of course of the number of bins (levels in my code)

[Mathieu NOE]

Could it be that ?

tq_rms=sqrt(sum(TQ.^2.*(N-1)*dt)/T)

replaced by

tq_rms=sqrt(sum(TQ.^2.*N*dt)/T)

gives

tq_rms = 0.885116049594405

still there is a slight difference (?)

[dpb]

The OP's calculation is time-weighted, not an unweighted global estimate.

[ABDALLA AL KHALEDI]

I was also wondering about this myself, the two codes give different results.

Also I tried the following code which gave close reults

%Torque Post processor No. 3

clc;

%Import the Excel File

file = uigetfile;

data = readmatrix(file);

t = data(:,1);

Torque = data(:,2);

%Total cycle time in [s]

t_tot = 81.92;

%The unique torque values in the data set

Torque_int = sort(unique(Torque, "stable"));

[Rep, Torque_int] = groupcounts(Torque);

T_rms = sqrt(sum(Torque_int.^2 .*(Rep-1).*0.005)/t_tot)

However, I faced another small problem, that is the readmatrix command will give an error unless the m-file and the excel file have the same path, how can I resolve this?

Error using readmatrix (line 148)

Unable to find or open 'Torque Data.xlsx'. Check the path and filename or file permissions.

Error in Torque_Post_Processor_No3 (line 6)

data = readmatrix(file);

I want to store the excel files I have in a seperate folder than the m-files.

Thank you for your help.

[dpb]

Use fully-qualified filenames to address the data paths, don't rely on everything being in your local working directory (and don't just keep adding stuff to the MATLABPATH). fullfile is useful for building qualified file names...

[dpb]

"...the two codes give different results."

Of course, they all calculate slightly different things/make different assumptions about what time-weighting to use and the number of binning levels (the global estimate doesn't bin/time-weight at all).

The Answer above computes precisely what the formula says with your initial additional caveat of not counting the one-element samples. Although the differences are small enough it would seem unlikely that it would really make any practical difference.

[Mathieu NOE]

just to be clear

my code gives both time weighted and global rms outputs , and both values matches as soon as the number of bins / levels are reasonnably high (> 5), which is what we expect

for the single isolated values , it's simply stated that they could be omited , not that we are asked to remove them

but maybe here I am pin pointing and my goal is certainly not to make long debates about details that don't matter from a practical point of view

[dpb]

OP went on to state "...these can be omitted as I am int[e]rested in the overall average of the intervals." (emphasis added--dpb) from which it can be inferred the single sample spikes should not be counted. <vbg>

[dpb]

"...about that don't matter from a practical point of view"

For the sample dataset, all the above return essentially the same result. The time-weighted value could be somewhat more different depending on the particular transient behavior, however...

[ABDALLA AL KHALEDI]

Thank you all. Yes the sample spikes are to be neglected.

Answer 2

[Joe Vinciguerra]

Here is one way to do it based on your example:

% import the data from your file using readmatrix(), but using your example:
t = [0.2 0.4 0.6 0.8 1.0]'
data = [3 3 2 2 2]'

% find the unique torque values in data set to create an index for use below
[C, ia, ic] = unique(data, "stable")

% find the max and min times at which the unique values occure
t_max = accumarray(ic, t, [] ,@max)
t_min = accumarray(ic, t, [], @min)

% date the delta
t_delta = t_max - t_min

% perform you calcultion on these arrays
P = sum(C.^2 .*  t_delta)

[ABDALLA AL KHALEDI]

The problem with this code is that the same torque value appears at several intervals, and when this code is applied just the first and last time value related to the indexed torque are considered resulting in an error in the t_delta values. The same unique torque index can be observed over several intervals.

For example,

The instantaneous time = [ 0.2, 0.4, 0.6, 0.8, 1.0, 1.2]

The Torque intervals:

from 0.2 to 0.4 ---> T= 2

from 0.4 to 0.6 ---> T= -5

at t = 0.6 ---> T= -6

from 0.6 to 0.8 ---> T= 2

from 0.8 to 1 ---> T= 3

from 1.0 to 1.2 ---> T= 3

So you see torque intervals are repeating for various intervals, plus some torque values appear as a single data point at a certain time but do not continue to show up as time passes like (-6) above, these I can ommit if possible but I didn't know how.

Please check the excel file I have attached with the question to fully understand the data.

I wrote the following code, how can I apply what I mentioned above?

also what is ia and ic? I didn't understand that.

%Data Analizer - RMS Torque Calculator

clc;

%Import the Torque data from your Excel File

%Ex: 'Torque_Data.xlsx' with the quotation marks

prompt = 'What is The File Name ? ';

fileName= input(prompt);

%Example 'C4:C10' with the quotation marks

prompt = 'What is The Range of Data ? ';

xlRange = input(prompt);

%Instantaneous Torque Data

Data = xlsread(fileName,xlRange);

%Instantaneous time in [s]

t = (0.005:0.005:81.92)';

%Total cycle time in [s]

t_tot = 81.92;

%The unique torque values in the data set

[T_int] = sort(unique(Data, "stable"));

[C, ia, ic] = unique(Data, "stable")

% find the max and min times at which the unique values occure

t_max = accumarray(ic, t, [], @max)

t_min = accumarray(ic, t, [], @min)

% date the delta (error)

t_delta = t_max - t_min

% perform you calcultion on these arrays

T_rms = (sum(C.^2.*t_delta))./t_tot

Answer 3

[Mathieu NOE]

hello

maybe this ? I know I reinvented the histogram ( :)) , in few lines...

% Torque post processing 

data = readmatrix("Torque Data.xlsx");
t = data(:,1);
torque = data(:,2);
dt = mean(diff(t));


% create levels (edges, like for a histogram)
Tmin = floor(min(torque));
Tmax = ceil(max(torque));
levels = 10;
Trange = linspace(Tmin,Tmax,levels+1);

for ci = 1:levels
    id = (torque>=Trange(ci) & torque<Trange(ci+1));
    torq_extr = torque(id);
    samples = numel(find(id));
    duration(ci) = samples*dt;
    % T² computation
    torq2(ci) = mean(torq_extr.^2); % T1², T2², etc
end

% finaly :
Torque_rms = sqrt(sum(torq2.*duration)/sum(duration));

 % direct rms computation
Torque_rms_alt = sqrt(mean(torque.^2));

Answer 4

[Star Strider]

It is not obvious to me exactly what you want to calculate. This isolates the regions between the spikes (that I refer to as ‘segments’), determines the times at the beginning and end of each segment, and calculates the RMS value of the signal in each region. I am not certain what you want to do beyond that.

One approach —

T1 = readtable('Torque Data.xlsx', 'VariableNamingRule','preserve')

T1 = 16384×2 table
    Time [s]    Torque [N.m]
    ________    ____________

     0.005        -0.61183  
      0.01         -0.7138  
     0.015         -0.7138  
      0.02         -0.7138  
     0.025         -0.7138  
      0.03         -0.7138  
     0.035         -0.7138  
      0.04         -0.7138  
     0.045         -0.7138  
      0.05         -0.7138  
     0.055         -0.7138  
      0.06         -0.7138  
     0.065         -0.7138  
      0.07         -0.7138  
     0.075         -0.7138  
      0.08         -0.7138  

VN = T1.Properties.VariableNames;
[pks,locs] = findpeaks(T1{:,2}, 'MinPeakProminence',1.5);               %  Quality Check

idx = [1; find(diff(sign(T1{:,2}-0.5))); size(T1,1)];                   % Crossing Indices
blen = 2*ceil(numel(idx)/2);                                            % Calculate Column Length Of 'idxmtx'
idxmtx = zeros(blen,1)+idx(end);                                        % Preallocate
idxmtx(1:numel(idx)) = idx;                                             % Assign Values To Vector
idxmtx = reshape(idxmtx, 2, [])                                         % Convert To (2xN) Matrix Of Segment Beginning & End Indices

idxmtx = 2×33

1 774 1798 2276 2772 3286 3780 4307 4787 5280 5762 5809 6288 6769 7264 7792 8272 8766 9295 9773 10304 10783 11278 11806 12267 12780 13310 13787 14282 14778 742 1752 2243 2754 3251 3745 4259 4752 5248 5745 5777 6240 6734 7246 7743 8238 8750 9260 9740 10255 10748 11244 11756 12250 12763 13258 13740 14248 14745 15241

% idxmtx(:,end)                                                           % Check (Optional)
for k = 1:size(idxmtx,2)
    idxrng = idxmtx(1,k) : idxmtx(2,k);                                 % Vector Index Range
    [b1,b2] = bounds(T1{idxrng,2});                                     % Range Of Segment Values (Optional
    rangmtx(:,k) = [b1; b2; b2-b1];                                     % Store As Matrix (Optional)
    timev{k} = T1{idxrng,1};                                            % Segment Time Vector
    times(:,k) = T1{idxmtx(:,k),1};                                     % Segment Time Limits
    rmsv(k) = rms(T1{idxrng,2});                                        % Segment RMS Value
end

times                                                                   % Segment Beginning & End Times

times = 2×33

0.0050 3.8700 8.9900 11.3800 13.8600 16.4300 18.9000 21.5350 23.9350 26.4000 28.8100 29.0450 31.4400 33.8450 36.3200 38.9600 41.3600 43.8300 46.4750 48.8650 51.5200 53.9150 56.3900 59.0300 61.3350 63.9000 66.5500 68.9350 71.4100 73.8900 3.7100 8.7600 11.2150 13.7700 16.2550 18.7250 21.2950 23.7600 26.2400 28.7250 28.8850 31.2000 33.6700 36.2300 38.7150 41.1900 43.7500 46.3000 48.7000 51.2750 53.7400 56.2200 58.7800 61.2500 63.8150 66.2900 68.7000 71.2400 73.7250 76.2050

rangmtx                                                                 % Segment Minimum, Maximum & Difference

rangmtx = 3×33

-1.5296 -1.5296 -1.5296 -1.1217 -1.0197 -1.0197 -1.4276 -1.0197 -1.2237 -1.0197 0.3059 -1.0197 -1.0197 -1.2237 -1.4276 -1.0197 -1.0197 -1.5296 -1.2237 -1.0197 -1.0197 -1.1217 -1.3256 -1.1217 -1.2237 -0.9177 -0.9177 -1.5296 -0.8158 -0.9177 0 0.7138 0.6118 0.7138 0.8158 0.6118 0.5099 0.8158 1.3256 0.6118 0.6118 0.9177 0.8158 0.7138 0.7138 0.6118 0.6118 0.7138 0.6118 1.0197 0.9177 0.5099 0.7138 0.7138 0.8158 0.5099 0.6118 1.0197 0.6118 0.5099 1.5296 2.2434 2.1414 1.8355 1.8355 1.6316 1.9375 1.8355 2.5493 1.6316 0.3059 1.9375 1.8355 1.9375 2.1414 1.6316 1.6316 2.2434 1.8355 2.0394 1.9375 1.6316 2.0394 1.8355 2.0394 1.4276 1.5296 2.5493 1.4276 1.4276

rmsv                                                                    % Segment RMS Values

rmsv = 1×33

0.7204 0.6390 0.6854 0.6242 0.5899 0.5981 0.6026 0.6295 0.6766 0.6445 0.3467 0.6364 0.6269 0.5885 0.6655 0.6218 0.6183 0.6538 0.6334 0.6280 0.6363 0.6604 0.6280 0.6270 0.6467 0.5686 0.6178 0.6929 0.5897 0.5984

figure
hp{1} = plot(T1{:,1}, T1{:,2}, 'DisplayName','Data');
hold on
% plot(T1{locs,1}, pks, 'r+', 'DisplayName','Peaks')
hp{2} = plot(times, [1;1]*rmsv,'-m', 'DisplayName','Segment RMS Values');
hold off
grid
xlabel(VN{1})
ylabel(VN{2})
legend([hp{1},hp{2}(1)], 'Location','best')

I believe this gest close to what you asked for, however I do not understand how the times weight the RMS values. Incorporating that would be straightforward.

.