Find the max width of a peak

Hi

I am trying to find a method to find the max width of peaks from a csv file.

I have used findpeaks to locate the peaks but I now need to know how wide they are.

I understand Matlab uses full width at half amplitude and have been able to use this but I need to know the width at the bottom or treshold of the peak. Basically the lenght of the green lines in the attached peaks.jpeg.

Ideally I need the rough starting and end points (index?) of the of the peaks, so I can reference these points back to my csv file. I need to extract the data from the other column between these points for each peak.

Below is some of the code I have used.

Any help would be greatly appriciated!

Array=csvread('accel test 1 bench.csv');
col1 = Array(:, 1);
col2 = Array(:, 2);
plot(col1, col2)
>> plot(col1)
>> pks = findpeaks(col1)
findpeaks(col1,'MinPeakDistance',5, "MinPeakHeight",3)

[pks,locs,widths,proms] = findpeaks(col1,'MinPeakDistance',5, "MinPeakHeight",3); % I used this to obtain the FWHM

[Jonas]

use the prominence of the peak values. then search for the first time the value peakValue-peakProminence appears behind the peak and the last time peakValue-peakProminence appears before the peak

you can check if this procedure will be successfull by looking into the graph eith the additional name-value pair 'Annotate','extents'

this will work many peaks but not for all

[Nigel Davis]

Hi

Thanks for the quick reply!

I don't fully understand, I'm pretty new to Matlab!

What do you mean by search for the peakValue-peakProminence? And do you mean subtract the prominence from the peak value?

Answer 1

[Mathieu NOE]

hello Nigel

maybe there is some way to make findpeaks compute the prominence in a manner that better works on this example , but I could not "tweak" it as I wanted .

So I opted to go my own way and created myprom , based on looking at n samples left then right from each peak and compute the max of left min and right min .

Finally could get more or less what you wanted : the "width" interval is computed in period1 array.

the root of each peak starts with a red diamond marker and ends with a green one.

period1 =

19.1716

16.2720

11.7574

14.8297

13.6186

15.1826

NB : your data does not mentionned any sampling rate so I created the x axis (time) with increment of 1

Array=csvread('accel test 1 bench.csv');
y = Array(:, 1);
x = (0:numel(y)-1);
[pks,locs,widths,proms] = findpeaks(y,'MinPeakDistance',5, "MinPeakHeight",3); % I used this to obtain the FWHM

figure(1),plot(x, y , x(locs) , pks ,'dk')

hold on 

for ci = 1:numel(pks)
    
    % my own "prom" definition find the max of local minima left and right
    % from this peak and look only n samples left and right 
    n = 20; % samples (left and right from peak)
    left_y = y(x<locs(ci));
    left_y = left_y(end-n+1:end);
    
    right_y = y(x>locs(ci));
    right_y = right_y(1:n);
    
    local_min = max(min(left_y),min(right_y));
    myproms(ci,1) = pks(ci) - local_min;

%     threshold = pks(ci)-proms(ci) + 0.5 ; % your value here
    threshold = pks(ci)-myproms(ci) + 0.5 ; % your value here

    [t0_pos1,s0_pos1,t0_neg1,s0_neg1]= crossing_V7(y,x,threshold,'linear'); % positive (pos) and negative (neg) slope crossing points 

    % find nearest t0_pos1 and t0_neg1 values to selected peak locs
    ind = find(t0_pos1 < locs(ci));
    t0_pos1 = t0_pos1(ind(end));

    ind = find(t0_neg1 > locs(ci));
    t0_neg1 = t0_neg1(ind(1));

    % periods
    period1(ci,1) = (t0_neg1 - t0_pos1); % time delta

    plot(t0_pos1,s0_pos1,'dr',t0_neg1,s0_neg1,'dg','markersize',12);
    

end
hold off 

period1

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [t0_pos,s0_pos,t0_neg,s0_neg] = crossing_V7(S,t,level,imeth)

% [ind,t0,s0,t0close,s0close] = crossing_V6(S,t,level,imeth,slope_sign) % older format

% CROSSING find the crossings of a given level of a signal
%   ind = CROSSING(S) returns an index vector ind, the signal
%   S crosses zero at ind or at between ind and ind+1
%   [ind,t0] = CROSSING(S,t) additionally returns a time
%   vector t0 of the zero crossings of the signal S. The crossing
%   times are linearly interpolated between the given times t
%   [ind,t0] = CROSSING(S,t,level) returns the crossings of the
%   given level instead of the zero crossings
%   ind = CROSSING(S,[],level) as above but without time interpolation
%   [ind,t0] = CROSSING(S,t,level,par) allows additional parameters
%   par = {'none'|'linear'}.
%	With interpolation turned off (par = 'none') this function always
%	returns the value left of the zero (the data point thats nearest
%   to the zero AND smaller than the zero crossing).
%


% check the number of input arguments
error(nargchk(1,4,nargin));

% check the time vector input for consistency
if nargin < 2 | isempty(t)
	% if no time vector is given, use the index vector as time
    t = 1:length(S);
elseif length(t) ~= length(S)
	% if S and t are not of the same length, throw an error
    error('t and S must be of identical length!');    
end

% check the level input
if nargin < 3
	% set standard value 0, if level is not given
    level = 0;
end

% check interpolation method input
if nargin < 4
    imeth = 'linear';
end


% make row vectors
t = t(:)';
S = S(:)';

% always search for zeros. So if we want the crossing of 
% any other threshold value "level", we subtract it from
% the values and search for zeros.
S   = S - level;

% first look for exact zeros
ind0 = find( S == 0 ); 

% then look for zero crossings between data points
S1 = S(1:end-1) .* S(2:end);
ind1 = find( S1 < 0 );

% bring exact zeros and "in-between" zeros together 
ind = sort([ind0 ind1]);

% and pick the associated time values
t0 = t(ind); 
s0 = S(ind);

if ~isempty(ind)

    if strcmp(imeth,'linear')
        % linear interpolation of crossing
        for ii=1:length(t0)
            %if abs(S(ind(ii))) >= eps(S(ind(ii)))    % MATLAB V7 et +
            if abs(S(ind(ii))) >= eps*abs(S(ind(ii)))    % MATLAB V6 et -    EPS * ABS(X)

                % interpolate only when data point is not already zero
                NUM = (t(ind(ii)+1) - t(ind(ii)));
                DEN = (S(ind(ii)+1) - S(ind(ii)));
                slope =  NUM / DEN;
                slope_sign(ii) = sign(slope);
                t0(ii) = t0(ii) - S(ind(ii)) * slope;
                s0(ii) = level;
            end
        end
    end

    % extract the positive slope crossing points 
    ind_pos = find(sign(slope_sign)>0);
    t0_pos = t0(ind_pos);
    s0_pos = s0(ind_pos);

    % extract the negative slope crossing points 
    ind_neg = find(sign(slope_sign)<0);
    t0_neg = t0(ind_neg);
    s0_neg = s0(ind_neg);

else

    % empty output
    ind_pos = [];
    t0_pos = [];
    s0_pos = [];

    % extract the negative slope crossing points 
    ind_neg = [];
    t0_neg = [];
    s0_neg = [];
    
end

end