peak average of consecutive values in a matrix array

Hello, I would like to calculate the average of the peaks (positive and negative values) of consecutive values from a predifined number of elements from a Nx1 matrix, where N is the number of rows. For example, let's say I have a matrix with the following form (24x1):

A[1;2;3;4;5;6;7;8;9;10;11;12;13;14;15;16;17;18;19;20;21;22;23;24]

I want to obtain the peak and the average for every 5 elements,

B=[mean(peak(1,2,3,4,5)); mean(peak((2,3,4,5,6)); mean(peak((3,4,5,6,7));...........mean(peak((20,21,22,23,24))]

I am chosing every 5 elements to simplify the example. But, this number number of elements can be higher than 100 and the points vary dramatically from positive to negative as they come from acceleration earthquake signals from N higher than 245000.

I would appreciate the help.

[Dyuman Joshi]

Following up from your previous questions, are movmean and movmax not working for you?

[Jorge Luis Paredes Estacio]

Hi, yes it worked. But peaks estimations involve as many peak values (positive or negative) not just one when using movmax. Thank you.

[Dyuman Joshi]

"But peaks estimations involve as many peak values (positive or negative) not just one when using movmax."

Can you provide an example for this?

[Jorge Luis Paredes Estacio]

Of course. I have attached an example file where there are three columns. Each column represents the component X, Y and Z of the acceleration earthquake signal in cm/sec2. I could say based on the previous description I would like to estimate the average of the peaks (positive and negative) for every 51 elements from the total signal that contains 9600 points for each component of the signal X, Y and Z. Thank you for your support.

[Dyuman Joshi]

Suppose this is the data in hand -

x = randi([-20 20],1,15)

x = 1×15

-5 -16 14 -18 19 -3 18 1 -10 -6 -12 16 5 8 0

Assume the moving window is 4, what should be the output in thise case? and what is the logic behind it?

[Jorge Luis Paredes Estacio]

Thank you for your help.

Answer 1

[Mathieu NOE]

hello

maybe this ?

I am using peakseek , this fex submission is available here

PeakSeek - File Exchange - MATLAB Central (mathworks.com)

, but you can use the standard (slower) findpeaks function offered by TMW

the code provided below will buffer your data with length = 51 samples , then split the signal for each cahnnel into a positive and negative signal and find the positive and negative peaks , then take the mean of them (separately not mixing pos and neg peaks !)

example for X channel :

data = readmatrix('Signal.txt'); % acceleration earthquake signal in cm/sec2 ,X, Y and Z
[samples, channels] = size(data);

dt = 1; % sample rate (s)
t = dt*(0:samples-1);  %  time vector 


%% home made solution (you choose the amount of overlap)
buffer_size = 51; % how many samples 
overlap = 0; % overlap expressed in samples

%%%% main loop %%%%
[new_time,data_out_pos,data_out_neg] = my_peakmean(t,data,buffer_size,overlap);

figure(1),
plot(t,data(:,1),new_time,data_out_pos(:,1),'*-g',new_time,data_out_neg(:,1),'*-r');
title('X');
legend('raw data','pos peaks mean','neg peaks mean');
xlabel('Time(s)');
ylabel('cm/sec2');

figure(2),
plot(t,data(:,2),new_time,data_out_pos(:,2),'*-g',new_time,data_out_neg(:,2),'*-r');
title('Y');
legend('raw data','pos peaks mean','neg peaks mean');
xlabel('Time(s)');
ylabel('cm/sec2');

figure(3),
plot(t,data(:,3),new_time,data_out_pos(:,3),'*-g',new_time,data_out_neg(:,3),'*-r');
title('Z');
legend('raw data','pos peaks mean','neg peaks mean');
xlabel('Time(s)');
ylabel('cm/sec2');

%%%%%%%%%% my functions %%%%%%%%%%%%%%
function [new_time,data_out_pos,data_out_neg] = my_peakmean(t,data_in,buffer_size,overlap)

% NB : buffer size and overlap are integer numbers (samples)
% data (in , out) are 1D arrays (vectors)

shift = buffer_size-overlap;    % nb of samples between 2 contiguous buffers  
[samples,channels] = size(data_in);
nb_of_loops = fix((samples-buffer_size)/shift +1);

    for k=1:nb_of_loops
        start_index = 1+(k-1)*shift;
        stop_index = min(start_index+ buffer_size-1,samples);
        x_index(k) = round((start_index+stop_index)/2);  
        % find peaks 
        for n = 1:channels
            tmp = data_in(start_index:stop_index,n);
            tmp_pos = tmp;
            tmp_pos(tmp_pos<0) = 0;
            tmp_neg = tmp;
            tmp_neg(tmp_neg>0) = 0;
            
            [locsp, pksp]= peakseek(tmp_pos,1,max(tmp_pos)/10); % pos peaks above 10% of max pos amplitude
            [locsn, pksn]= peakseek(-tmp_neg,1,max(-tmp_neg)/10); % neg peaks above 10% of max neg amplitude
            data_out_pos(k,n) = mean(pksp); % mean of pos peaks 
            data_out_neg(k,n) = mean(-pksn); % mean of neg peaks 
       
        end

        
    end

new_time = t(x_index); % time values are computed at the center of the buffer

end

[Jorge Luis Paredes Estacio]

This is a very clever answer straight to the issue. I really appreciate your help.

I was checking the code, it works perfectly fine. However, the estimation of peak mean (defined as 51 as a buffer for this example) should go after each consecutive points. In other words, the averaging should start in time 26, then 27, 28, 29 ...until 9575. This means that the averaging of peaks (positive and negative) are estimated from the following range of points:

for time 26 the averaging of peaks are estimated from points 1 to 51

for time 27 the averaging of peaks are estimated from points 2 to 52

for time 28 the averaging of peaks are estimated from points 3 to 53

for time 29 the averaging of peaks are estimated from points 4 to 54

. .

. .

for time 9575 the averaging of peaks are estimated from points 9550 to 9600.

This can be handle with the overlap value set to 0?

I would appreciate your help

Best

[Mathieu NOE]

hello again

thanks for the nice comments

to get the required behaviour simply change the overlap value , you need now :

overlap = buffer_size-1; % overlap expressed in samples

[Mathieu NOE]

just one point : in my code , the time vector is zero based

so the first buffer (points 1 to 51) , assuming a sampling rate = 1 , correspond to t = 0 to 50 , so my buffer is centered at t = 25 (and not 26 as in your comment above)

this explain why we have a time difference of 1 between you and me.

my new time vectors is 25,26,27,... 9574

[Jorge Luis Paredes Estacio]

I really appreaciate your help. it is sorted now. Thank you very much.

[Mathieu NOE]

as always, my pleasure !