Count impulses into sequence of impulses and segments of sequences of impulses from signal

For example I have signal of current:

And I count the all of impulses such same:

numberOfPulses = nnz(diff(outSig_1_curr > 1) > 0);

But I don't know how to count numbers of marked sequences of impulses in signal and count nums of impulses for each sample (markes as black rectangle on png)?

So, the some of signals contain the such abnormal samples (sequences of impulses)

How can I count this abnormal samples and impulses for each abnormal samples?

Maybe I could findpeaks function, but I don't know ho to separate each samples at the other impulses in signal..

Signals .mat files available in GDrive: Google Drive

Thanks in advance for your help!

[Mathieu NOE]

hello

could you please precise

1/ sampling freq

2/ wich portion of the data are you showing ?

[Alex Milns]

1/ sampling freq - 1 000 000 Hz

2/ wich portion of the data are you showing ?

I demonstrate time fragment from signals (all signal available in Google Drive). First screen - sig1 data, Second screen - sig16 data

Answer 1

[Mathieu NOE]

hello again

so this is my attempt... not maybe the shortest and most efficienct code , but seems to do what I wanted

so we can slect blocks of data based on two selection parameters :

• the minimal window length (in samples) where we want peaks : here 1000 samples = 1 ms
• the min distance between blocks of data : here we want separation above 400 samples = 0.4 ms

The code shows here 3 blocks of data which contains this amount of peaks :

number_of_peaks = 18 111 15

hope it helps

code

clc
clearvars

load('sig1.mat') % outSig_1_curr2      32800000x1
samples = numel(outSig_1_curr2);
Fs = 1000000; %  Hz
dt = 1/Fs;
t = (0:samples-1)*dt;

% extract signal between t = 21.9112 and t = 21.94 s
ind1 = round(21.9112*Fs);
ind2 = round(21.94*Fs);

t = t(ind1:ind2);
outSig_1_curr2 = outSig_1_curr2(ind1:ind2);

%% find peaks
% find positive local maxima
[tf, P] = islocalmax(outSig_1_curr2,'MinProminence',200);
t_peak = t(tf);
outSig_1_peak = outSig_1_curr2(tf);


%% localize / extract valid data blocks 
blocks_distance_samples = 100; % select data blocks  that are distant above value = blocks_distance_samples (in samples)
min_contiguous_samples = 1000; % select data blocks  that are contiguous with min amount of samples = min_contiguous_samples

detect_signal = zeros(size(t))+0.5; % initialisation

    % now define start en end point of  segments above threshold
        %%%%%%%%%%
        % This locates the beginning /ending points of data groups
        ind = find(tf); % logical to index conversion
        D = diff([0;ind;0]);
        ind_ends = find(D > blocks_distance_samples) - 1;
        ends = ind(ind_ends);
        begin = ind(ind_ends+1);
        % remove incorrect begin / ends points (like end point before begin  point and vice versa)
        ends = ends(ends>begin(1));
        begin = begin(begin<ends(end));
        m = min(numel(begin),numel(ends));
        ends = ends(1:m);
        begin = begin(1:m);
        
        
        %%%%%%%%%%
        
    length_ind = ends - begin; % compute block length
    ind2= length_ind>min_contiguous_samples; % check if their length is valid (above min_contiguous_samples value)
    begin = begin(ind2); % selected  points 
    ends = ends(ind2); % selected  points 
    
    t_peak2 = t(ind);
    outSig_1_peak2 = outSig_1_curr2(ind);
    % define the begin / ending x, y values of raw data 
    t_peak2_begin = t(begin);
    outSig_1_curr2_begin = outSig_1_curr2(begin);
    t_peak2_ends = t(ends);
    outSig_1_curr2_ends = outSig_1_curr2(ends);

    
for ci = 1:length(begin)
    ind = (t>=t_peak2_begin(ci) & t<=t_peak2_ends(ci));
    detect_signal(ind) = 200;
    % include here your code if you want to extract the data blocks
    % individually
    ind4 = (t_peak>=t_peak2_begin(ci) & t_peak<=t_peak2_ends(ci));
    number_of_peaks(ci) = numel(find(ind4));
    
end

number_of_peaks

% plot
figure(1)
plot(t,outSig_1_curr2,t_peak,outSig_1_peak,'dr',t(begin),outSig_1_curr2(begin),'dg',t(ends),outSig_1_curr2(ends),'dk',t,detect_signal,'k');
legend('signal',' peaks ','begin points','end points','selected signal')

[Alex Milns]

This is very very good applies, but if I want to count every blocks

I didn't choose all signal, or I did?

And if sequence of pulses lools like a screen two this code inspect and detect sequences?

Could you answer of this questions, please!

[Alex Milns]

So, if I try to run script I got error:

Array indices must be positive integers or logical values.

Error in testing_counter (line 37)

ends = ind(ind_ends);

So, please, let me clarify what can be the minimum window size and the size between blocks when using this script?

[Mathieu NOE]

hello

I changed a bit the parameter to select more blocks (less long)

min_contiguous_samples changed from 1000 to 100 :

min_contiguous_samples = 100; % select data blocks  that are contiguous with min amount of samples = min_contiguous_samples

this is the result now :

number_of_peaks = 18 4 4 3 5 111 10 15 7 6 8

also I made a slight change in the code to make it more robust

clc
clearvars

load('sig1.mat') % outSig_1_curr2      32800000x1
samples = numel(outSig_1_curr2);
Fs = 1000000; %  Hz
dt = 1/Fs;
t = (0:samples-1)*dt;

% extract signal between t = 21.9112 and t = 21.94 s
ind1 = round(21.9112*Fs);
ind2 = round(21.94*Fs);

t = t(ind1:ind2);
outSig_1_curr2 = outSig_1_curr2(ind1:ind2);

%% find peaks
% find positive local maxima
[tf, P] = islocalmax(outSig_1_curr2,'MinProminence',200);
t_peak = t(tf);
outSig_1_peak = outSig_1_curr2(tf);


%% localize / extract valid data blocks 
blocks_distance_samples = 100; % select data blocks  that are distant above value = blocks_distance_samples (in samples)
min_contiguous_samples = 100; % select data blocks  that are contiguous with min amount of samples = min_contiguous_samples

detect_signal = zeros(size(t))+0.5; % initialisation

    % now define start en end point of  segments above threshold
        %%%%%%%%%%
        % This locates the beginning /ending points of data groups
        ind = find(tf); % logical to index conversion
        D = diff([0;ind;0]);
        ind_ends = find(D > blocks_distance_samples) - 1;
        ind_ends = ind_ends(ind_ends>0); % remove potential zero value
        ends = ind(ind_ends);
        begin = ind(ind_ends+1);
        % remove incorrect begin / ends points (like end point before begin  point and vice versa)
        ends = ends(ends>begin(1));
        begin = begin(begin<ends(end));
        m = min(numel(begin),numel(ends));
        ends = ends(1:m);
        begin = begin(1:m);
        
        
        %%%%%%%%%%
        
    length_ind = ends - begin; % compute block length
    ind2= length_ind>min_contiguous_samples; % check if their length is valid (above min_contiguous_samples value)
    begin = begin(ind2); % selected  points 
    ends = ends(ind2); % selected  points 
    
    t_peak2 = t(ind);
    outSig_1_peak2 = outSig_1_curr2(ind);
    % define the begin / ending x, y values of raw data 
    t_peak2_begin = t(begin);
    outSig_1_curr2_begin = outSig_1_curr2(begin);
    t_peak2_ends = t(ends);
    outSig_1_curr2_ends = outSig_1_curr2(ends);

    
for ci = 1:length(begin)
    ind = (t>=t_peak2_begin(ci) & t<=t_peak2_ends(ci));
    detect_signal(ind) = 200;
    % include here your code if you want to extract the data blocks
    % individually
    ind4 = (t_peak>=t_peak2_begin(ci) & t_peak<=t_peak2_ends(ci));
    number_of_peaks(ci) = numel(find(ind4));
    
end

number_of_peaks

% plot
figure(1)
plot(t,outSig_1_curr2,t_peak,outSig_1_peak,'dr',t(begin),outSig_1_curr2(begin),'dg',t(ends),outSig_1_curr2(ends),'dk',t,detect_signal,'k');
legend('signal',' peaks ','begin points','end points','selected signal')

[Alex Milns]

Hello!

Now I check this script, so for sinal very large size (for example ind1 = round(0.001*Fs); ind2 = round(36.001*Fs);) i's work slow or not work. I can count the pulses only if I take a long window and a large distance between the pulses .. Probably I should try to increase speed of evaluation inds...

Let me ask the last question: if a piece of the pulse sequence does not start from zero, but, for example, in the middle in amplitude and above 6000 in signal amplitude, is it possible to make a separate calculation of the same, i.e. set the parameter islocalmax() in the range (search for a peak from the average value to the maximum in amplitude, i.e. the threshold value is exactly greater than the middle and greater than one more value, and count only for such a fragment). Those like in the second picture?

[Mathieu NOE]

hello again

maybe it would be good to split the large data in smaller chuncks , because from a code development and debugging perspective , I cannot work over the full length of data , which is impractical to plot and look for eventual bugs and side effects (once the window is filled with millions of peaks)

your data is very long so maybe think of a method to select a portion of interest before applying the sauce to the full monthy.

[Alex Milns]

You are right, this signal is very copmlex for processing, maybe I try to optimize speed, one of the problem: calculate separatly this abnormal segments

and so I ask is it possible to set a parameter 'MinProminence' for counting peaks in the form of a band from the beginning to the end of the amplitudes of such fragments... Have any idea?

Otherwise, the program is simply magnificent, for which many thanks to you!

[Alex Milns]

For example, if I extract such fragments sig1.mat from sig16.mat when I using this counter script I get the same error:

Array indices must be positive integers or logical values.

Error in testing_counter2 (line 31)

ends = ind(ind_ends);

I can assume that this is due to the almost zero difference between the pulses in such a section, but there is still a problem with just such a component.

[Mathieu NOE]

sorry, where do the fragments sig11.mat from sig16.mat come from ?

[Alex Milns]

sorry, Mathieu NOE, that I did not specify: sig 1.mat and sig16.mat - files of data points (signals), but I can't attach to post, so files is very large (over 250 mb)

[Alex Milns]

Now I have a idea create envelope for signals:

and calculate local minima peaks for fragment:

[tf, P] = islocalmin(outSig_1_curr2,'MinProminence',min(envelope));

but again I get the error:

Array indices must be positive integers or logical values.

Error in movavg (line 43)

ends = ind(ind_ends);

[Mathieu NOE]

can you somehow make a short extract of the data and share it ?

[Alex Milns]

yes

signal no 1 (extract signal between t = 18 and t = 20 s) sig1.mat

signal no 16 (extract signal between t = 4.2 and t = 4.8 s) sig16.mat

[Mathieu NOE]

hello again

so this is an improved code tested on new sig1 file

I changed a bit how I work with islocalmax. I prefer this new version where the test is done on a modified signal (the negative portion is replaced by zero , so the "prominence" value acts more like a fixed y threshold)

clc
clearvars

load('sig1.mat') % outSig_1_curr2
samples = numel(outSig_1_curr2);
Fs = 1000000; %  Hz
dt = 1/Fs;
t = (0:samples-1)*dt;

% % extract signal between t = 18 and t = 20 s
% ind1 = round(0.001*Fs);
% ind2 = round(0.1*Fs);
% t = t(ind1:ind2);
% outSig_1_curr2 = outSig_1_curr2(ind1:ind2);

% use tmp instead of outSig_1_curr2 data for islocalmax 
% tmp is the same as outSig_1_curr2 , except that negative values are
% replaced by zeroes
tmp = outSig_1_curr2;
tmp(tmp<0) = 0;

%% find peaks
% find positive local maxima
[tf, P] = islocalmax(tmp,'MinProminence',200); % see remark above
t_peak = t(tf);
outSig_1_peak = outSig_1_curr2(tf);


%% localize / extract valid data blocks 
blocks_distance_samples = 100; % select data blocks  that are distant above value = blocks_distance_samples (in samples)
min_contiguous_samples = 100; % select data blocks  that are contiguous with min amount of samples = min_contiguous_samples

detect_signal = zeros(size(t))+0.5; % initialisation

    % now define start en end point of  segments above threshold
        %%%%%%%%%%
        % This locates the beginning /ending points of data groups
        ind = find(tf); % logical to index conversion
        D = diff([0;ind;0]);
        ind_ends = find(D > blocks_distance_samples) - 1;
        ind_ends = ind_ends(ind_ends>0); % remove potential zero value
        ends = ind(ind_ends);
        begin = ind(ind_ends+1);
        % remove incorrect begin / ends points (like end point before begin  point and vice versa)
        ends = ends(ends>begin(1));
        begin = begin(begin<ends(end));
        m = min(numel(begin),numel(ends));
        ends = ends(1:m);
        begin = begin(1:m);
        
        
        %%%%%%%%%%
        
    length_ind = ends - begin; % compute block length
    ind2= length_ind>min_contiguous_samples; % check if their length is valid (above min_contiguous_samples value)
    begin = begin(ind2); % selected  points 
    ends = ends(ind2); % selected  points 
    
    t_peak2 = t(ind);
    outSig_1_peak2 = outSig_1_curr2(ind);
    % define the begin / ending x, y values of raw data 
    t_peak2_begin = t(begin);
    outSig_1_curr2_begin = outSig_1_curr2(begin);
    t_peak2_ends = t(ends);
    outSig_1_curr2_ends = outSig_1_curr2(ends);

    
for ci = 1:length(begin)
    ind = (t>=t_peak2_begin(ci) & t<=t_peak2_ends(ci));
    detect_signal(ind) = 200;
    % include here your code if you want to extract the data blocks
    % individually
    ind4 = (t_peak>=t_peak2_begin(ci) & t_peak<=t_peak2_ends(ci));
    number_of_peaks(ci) = numel(find(ind4));
    
end

number_of_peaks

% plot
figure(1)
plot(t,outSig_1_curr2,t_peak,outSig_1_peak,'dr',t(begin),outSig_1_curr2(begin),'dg',t(ends),outSig_1_curr2(ends),'dk',t,detect_signal,'k');
legend('signal',' peaks ','begin points','end points','selected signal')

[Alex Milns]

Hi, let me clarify, i.e. if a fragment occurs in the signal, when the sequence of pulses is raised relative to normal fragments, then it automatically takes the maximum value for the normal section as a threshold, and then already counts the number of pulses for abnormal sections? Did I understand you correctly?

It's just that on the example of signal sig16.mat, he cannot count such non-characteristic sections (when the minimum amplitude of the pulse clearly lies to the example of about 1500) ...

[Mathieu NOE]

ok

I didn't say the code above works for sig16

this is a different topic I still need to address separately

[Alex Milns]

Ok, I understand! Thanks a lot!

[Mathieu NOE]

hello Alex

finally, tried a slightly different approach for sig16 signals

I do first a lower envelope that allows me to rougly find the interesting portion of the signal, then i will extract the portion and apply islocalmin to count the peaks

you have a plot for each identified segment and the number of peaks is recorded (and shown in the title)

try it and let me know

load('sig16.mat') % outSig_16_curr2 

samples = numel(outSig_16_curr2);
Fs = 1000000; %  Hz
dt = 1/Fs;
t = (0:samples-1)*dt;
% 
% % extract signal between t = 18 and t = 20 s
% ind1 = 1%round(0.05*Fs);
% ind2 = round(0.15*Fs);
% t = t(ind1:ind2);
% outSig_16_curr2 = outSig_16_curr2(ind1:ind2);

% lower envelope
n_env = 300;
% tmp is the same as outSig_1_curr2 , except that negative values are
% replaced by zeroes
tmp = outSig_16_curr2;
tmp(tmp<0) = 0;

[up,down] = envelope(tmp,n_env,'peak');

%% positive threshold crossing
thresh = 300;  % Specify threshold
% Essentially need to pick out values which exceed threshold with the condition that the previous value  
% needs to be below the threshold
idxl = down>=thresh;
idxl(1) = 0;
idx = find(idxl);
yest = down(idx-1)<thresh; 
pos_trig_ind = idx(yest);  % Final output

%% negative threshold crossing
% Essentially need to pick out values which are below threshold with the condition that the previous value  
% needs to be above the threshold
idx2 = down<=thresh;
idx2(1) = 0;
idx = find(idx2);
yest = down(idx-1)>thresh; 
neg_trig_ind = idx(yest);  % Final output

%% remove incorrect begin / ends points (like end point before begin  point and vice versa)
begin = pos_trig_ind(pos_trig_ind<neg_trig_ind(end));
ends = neg_trig_ind(neg_trig_ind>pos_trig_ind(1));

%% localize / extract valid data blocks 
min_contiguous_samples = 1000; % select data blocks  that are contiguous with min amount of samples = min_contiguous_samples
length_ind = ends - begin; % compute block length
ind2= length_ind>min_contiguous_samples; % check if their length is valid (above min_contiguous_samples value)
begin = begin(ind2); % selected  points 
ends = ends(ind2); % selected  points 

%% plot    
figure(1)
plot(t,outSig_16_curr2,t,down,t(begin),down(begin),'*',t(ends),down(ends),'*','markersize',14);
legend('signal','down envelope',' begin ',' ends ')

prepost = n_env/Fs; % samples/Fs = seconds (add some samples before start / end points defined on envelop signal)

for ci = 1:length(begin)
    ind = (t>=t(begin(ci))-prepost & t<=t(ends(ci))+prepost);
    t_extract = t(ind);
    y_extract = outSig_16_curr2(ind);
    
    %% find peaks
    % find all negative local minima
    [tf, P] = islocalmin(y_extract,'MinProminence',200); % see remark above
    ti = find(tf);
    t_peak = t_extract(ti);
    data1_peak = y_extract(ti);

    % find the negative peaks between amplitude 500 and 3000
    id = find(data1_peak>500 & data1_peak<3000);
    t_peak = t_peak(id);
    data1_peak = data1_peak(id);
    number_of_peaks(ci) = numel(id);
    figure(1+ci)
    plot(t_extract,y_extract,t_peak,data1_peak,'dr');
    legend('signal',' peaks ')
    title([' Peaks found :' num2str(number_of_peaks(ci))]);

    
end

number_of_peaks

[Alex Milns]

Hello, Mathieu NOE!

It's very very cool solution to the counting problem. Thus, the first code can be used to count the working impulses, and this code even does an excellent job for the entire duration of the signal (the program will run for about 3 minutes)! Thank you very much for such great help!

[Mathieu NOE]

Glad I could be of some help !

All the best

Answer 2

[Karim]

If the rectangles were drawn using matlab, then you know the start and stop time.

Thus, if this works for the total number of pulses:

numberOfPulses = nnz(diff(outSig_1_curr > 1) > 0);

then you can add an extra logic operation to only count in a time interval:

numberOfPulses = nnz( [0 diff(outSig_1_curr > 1)] > 0 & t>=t_start & t<=t_stop);

note that I added a 0 to the first logic check, this is since the diff(outSig) will contain 1 element less than the time vector due to the diff operation.

also, note that you can add .mat files directly to your question by using the paperclip symbol

[Alex Milns]

This signal extract from TDMS file, AE and Current sensor, this signal is random, and all I know, more precisely, I assume that the duration of each pulse is 5 microseconds, but the number of pulses in each block, and the number of such blocks, I still do not know how to calculate. Moreover, it can be seen that some fragments (sequences of pulses are located above the midline)

The mat file 252 mb, so I attach link!

This fragmets has diffirent length and arbitrary start and end time!

[Alex Milns]

So I understand, that in my situation I should set 200 threshold value.

numberOfPulses = nnz(diff(outSig_1_curr > 200) > 0);

And I consider, that I use logical indexing for mark sequences of impulses, and count this numbers segments:

idx1 = [0 0 0 1 1 1 0 0 1 1 1 0 0 0 0 1 1 1 1 1 1 0 0 0 0 0 0 0 1 1 1 0 0];
diffs = diff([0 idx1 0]);

% find start index of the blocks
loc = find(diffs == 1);
% calc block lengths by subtracting end - start indices
Ns = find(diffs == -1) - loc;

C = mat2cell([loc' Ns'],ones(size(loc)))

  4×1 cell array

    {[ 4 3]}
    {[ 9 3]}
    {[16 6]}
    {[29 3]}

And I have a idea create two idx for count normal and abnormal segments, and count pulse for each segments, but I don't have idea for separate normal and abnormal fragments..

So this method not working:

 numberOfPulses = nnz( [0 diff(outSig_1_curr > 1)] > 0 & t>=t_start & t<=t_stop);
 

Error using horzcat

Dimensions of arrays being concatenated are not consistent.