Selecting which prominence to use with findpeaks

Hello

I am using findpeaks to find the prominence of my peaks in the attached imaged. However the prominence for peak 4 is instead the entire height of the function. How can I change this to match the other peaks (1,2,3,5) and why is this behaviour happening?

Thanks

HP

Code:

prominance = max(xDose(:,2))/3;

[pks,locs,widths,proms] = findpeaks(xDose(:,2),xDose(:,1),'MinPeakProminence', prominance);

findpeaks(xDose(:,2),xDose(:,1), 'MinPeakProminence', prominance,'Annotate','extents')

text(locs+1,pks,num2str((1:numel(pks))'))

Answer 1

[Mathieu NOE]

hello again @Henry

I am 100% sure to have fully understood all the maths you want to do on your data but I have some suggestions for you below

FYI I often prefer peakseek to use peakseek vs findpeaks . It's also in attachment if you are lazzy like me and want save some time ... PeakSeek - File Exchange - MATLAB Central

so here's my proposal , you can maybe refine it according to your own needs

the peak width is computed based on the lower subplot - distances between the red and green diamonds

peak_width = 5.6650 5.3281 5.9644 5.5536 6.6400

I did not compute the vertical distance between the baseline and the peaks but looking at the top subplot you can easily compute the y delta between the red and black diamonds

x   = linspace(0, 100, 1001);
c   = 30:10:70;
y1  = exp(- 1/2*(0.23*(x - c(1))).^2);
y2  = exp(- 1/2*(0.23*(x - c(2))).^2);
y3  = exp(- 1/2*(0.23*(x - c(3))).^2);
y4  = exp(- 1/2*(0.23*(x - c(4))).^2);
y5  = exp(- 1/2*(0.23*(x - c(5))).^2);
p1  = 0.8;      % peak 1 (leftmost)
p2  = 1.1;      % peak 2
p3  = 1.3;      % peak 3
p4  = 1.2;      % peak 4
p5  = 1.0;      % peak 5
y   = p1*y1.^2 + p2*y2.^2 + p3*y3.^2 + p4*y4.^2 + p5*y5.^2;

%% baseline correction 
Base = env_secant(x, y, 150, 'bottom'); 
Corrected_y = y - Base;

%% find positive and negative peaks
minpeakdist = 50;
[locsP, pksP]=peakseek(y,minpeakdist,0.3);
[locsN, pksN]=peakseek(-y,minpeakdist,-1);
pksN = -pksN;
% corresponding baseline points at x locs =  locsP
Base_locsP = interp1(x,Base,x(locsP));

%% measure peak width at given height
[locsPc, pksPc]=peakseek(Corrected_y,minpeakdist,0.5*max(Corrected_y));
threshold = 0.5*min(pksPc);
[ZxP,ZxN] = find_zc(x,Corrected_y,threshold);
peak_width = (ZxN) - (ZxP);

%% plots
subplot(2,1,1),plot(x,y)
hold on
plot(x(locsP),pksP,'dr','markersize',8)
plot(x(locsN),pksN,'dg','markersize',8)
plot(x(locsP),Base_locsP,'dk','markersize',8)
plot(x,Base,'--')
for k =1:numel(pksP)
    text(x(locsP(k)),pksP(k)*1.1,sprintf('  Peak W =%g',peak_width(k)));
end
hold off

%baseline corrected plot
subplot(2,1,2),plot(x,Corrected_y)
hold on
plot(ZxP,threshold*ones(size(ZxP)),'dr',ZxN,threshold*ones(size(ZxN)),'dg','markersize',8);
hold off

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [ZxP,ZxN] = find_zc(x,y,threshold)

% put data in rows
x = x(:);
y = y(:);

% positive slope "zero" crossing detection, using linear interpolation
y = y - threshold;
zci = @(data) find(diff(sign(data))>0);  %define function: returns indices of +ZCs
ix=zci(y);                      %find indices of + zero crossings of x
ZeroX = @(x0,y0,x1,y1) x0 - (y0.*(x0 - x1))./(y0 - y1); % Interpolated x value for Zero-Crossing 
ZxP = ZeroX(x(ix),y(ix),x(ix+1),y(ix+1));

% negative slope "zero" crossing detection, using linear interpolation
zci = @(data) find(diff(sign(data))<0);  %define function: returns indices of +ZCs
ix=zci(y);                      %find indices of + zero crossings of x
ZeroX = @(x0,y0,x1,y1) x0 - (y0.*(x0 - x1))./(y0 - y1); % Interpolated x value for Zero-Crossing 
ZxN = ZeroX(x(ix),y(ix),x(ix+1),y(ix+1));

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [env] = env_secant(x_data, y_data, view, side)
% Function call: env_secant(x_data, y_data, view, side)
% Calculates the top envelope of data <y_data> over <x_data>.
% Method used: 'secant-method'
% env_secant() observates the max. slope of about <view> points,
% and joints them to the resulting envelope.
% An interpolation over original x-values is done finally.
% <side> ('top' or 'bottom') defines which side to evolve.
% Author: Andreas Martin, Volkswagen AG, Germany

    if nargin == 0
        test( false );
        test( true );
        return
    end

    if nargin < 4
        error( '%s needs at least 4 input arguments!', mfilename );
    end

    assert( isnumeric(view) && isscalar(view) && view > 1, ...
           'Parameter <view> must be a value greater than 1!' );
    assert( isvector(x_data) && isnumeric(x_data) && all( isfinite( x_data ) ), ...
           'Parameter <x_data> has to be of vector type, holding finite numeric values!' );
    assert( isvector(y_data) && isnumeric(y_data) && all( isfinite( y_data ) ), ...
           'Parameter <y_data> has to be of vector type, holding finite numeric values!' );
    assert( isequal( size(x_data), size(y_data) ), ...
           'Parameters <x_data> and <y_data> must have same size and dimension!' );
    assert( ischar(side), ...
            'Parameter <side> must be ''top'' or ''bottom''!' );

    switch lower( side )
        case 'top'
            side = 1;
        case 'bottom'
            side = -1;
        otherwise
            error( 'Parameter <side> must be ''top'' or ''bottom''!' );
    end

    sz = size( x_data );

    x_data = x_data(:);
    x_diff = diff( x_data );
    x_diff = [min(x_diff), max(x_diff)];
    assert( x_diff(1) > 0, '<x_data> must be monotonic increasing!' );

    y_data = y_data(:);
    data_len = length( y_data );
    x_new = [];
    y_new = [];

    if diff( x_diff ) < eps( max(x_data) ) + eps( min(x_data) )
        % x_data is equidistant
        search_fcn = @( y_data, ii, i ) ...
                     max( ( y_data(ii) - y_data(i) ) ./ (ii-i)' * side );
    else
        % x_data is not equidistant
        search_fcn = @( y_data, ii, i ) ...
                     max( ( y_data(ii) - y_data(i) ) ./ ( x_data(ii) - x_data(i) ) * side );
    end


    i = 1;
    while i < data_len;
        ii = i+1:min( i + view, data_len );
        [ m, idx ] = search_fcn( y_data, ii, i );

        % New max. slope: store new "observation point"
        i = i + idx;
        x_new(end+1) = x_data(i);
        y_new(end+1) = y_data(i);
    end;

    env = interp1( x_new, y_new, x_data, 'linear', 'extrap' );
    env = reshape( env, sz );




function test( flagMonotonic )
    npts = 100000;
    y_data = cumsum( randn( npts, 1 ) ) .* cos( (1:npts)/50 )' + 100 * cos( (1:npts)/6000 )';
    if flagMonotonic
        x_data = (1:npts)' + 10;
    else
        x_diff = rand( size( y_data ) );
        x_data = cumsum( x_diff );
    end

    view = ceil( npts * 0.01 ); % 1 Percent of total length
    env_up = env_secant( x_data, y_data, view, 'top' );
    env_lo = env_secant( x_data, y_data, view, 'bottom' );

    figure
    plot( x_data, y_data, '-', 'Color', 0.8 * ones(3,1) );
    hold all
    h(1) = plot( x_data, env_up, 'b-', 'DisplayName', 'top' );
    h(2) = plot( x_data, env_lo, 'g-', 'DisplayName', 'bottom' );

    grid
    legend( 'show', h )
end

end

Answer 2

[Sam Chak]

Hi @Henry

Since you did not provide the data for the 5-peak graph, I created one that resembles the Five-Finger Mountain. According to Chinese mythology and the novel Journey to the West, the Monkey King was imprisoned there after a series of rebellions against the heavens. The documentation of findpeaks() explains how the prominence of a peak is measured (see Peak #6).

x   = linspace(0, 100, 1001)';
c   = 30:10:70;

y1  = exp(- 1/2*(0.23*(x - c(1))).^2);
y2  = exp(- 1/2*(0.23*(x - c(2))).^2);
y3  = exp(- 1/2*(0.23*(x - c(3))).^2);
y4  = exp(- 1/2*(0.23*(x - c(4))).^2);
y5  = exp(- 1/2*(0.23*(x - c(5))).^2);

p1  = 0.8;      % peak 1 (leftmost)
p2  = 1.1;      % peak 2
p3  = 1.3;      % peak 3
p4  = 1.2;      % peak 4
p5  = 1.0;      % peak 5
y   = p1*y1.^2 + p2*y2.^2 + p3*y3.^2 + p4*y4.^2 + p5*y5.^2;

figure
plot(x, y)
xlabel('x'), ylabel('y')
title('Five-Finger Mountain')
grid on

figure
prominance = max(y)/5;
[pks, locs, widths, proms] = findpeaks(y, x, 'MinPeakProminence', prominance);
findpeaks(y, x, 'MinPeakProminence', prominance, 'Annotate', 'extents')
text(locs+1, pks, num2str((1:numel(pks))'))
xlabel('x'), ylabel('y')

[Henry]

Thanks for your answer. Unfortunatley I cannot include the data since, even compressed, it exceeds the upload file size limit.

I can see in your graph that there are some prominences that are measure from the minima to the peaks (1,2,4,5) and one where the prominance is from the peak to the 0 value (peak 3). Is there a way to choose how these are measure?

For example, on your graph could you choose that all the peaks promencnaces are measured from the peak to the zero value, or all are measure from the peak to the minima?

Thanks

[Mathieu NOE]

@Henry

I have the feeling you would have wished to have a "MaxProminence" parameter with findpeaks ?

but that does not exist or you have to use the other available parameters or use a differen function - depends on what your fianl goal is - which is what btw ?

[Sam Chak]

Hi @Henry

I aimed to explain the observation and how prominence is measured according to the documentation and the diagram. Here, you can manually estimate the prominence at the 3rd peak.

x   = linspace(0, 100, 1001)';
c   = 30:10:70;

y1  = exp(- 1/2*(0.23*(x - c(1))).^2);
y2  = exp(- 1/2*(0.23*(x - c(2))).^2);
y3  = exp(- 1/2*(0.23*(x - c(3))).^2);
y4  = exp(- 1/2*(0.23*(x - c(4))).^2);
y5  = exp(- 1/2*(0.23*(x - c(5))).^2);

p1  = 0.8;      % peak 1 (leftmost)
p2  = 1.1;      % peak 2
p3  = 1.3;      % peak 3
p4  = 1.2;      % peak 4
p5  = 1.0;      % peak 5
y   = p1*y1.^2 + p2*y2.^2 + p3*y3.^2 + p4*y4.^2 + p5*y5.^2;

figure
prominance = max(y)/5;
[pks, locs, widths, proms] = findpeaks(y, x, 'MinPeakProminence', prominance);
proms(3) = pks(3) - (pks(4) - proms(4))

proms = 5×1

0.3096 0.4745 0.6462 0.5462 0.4238

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

hold on
plot(x, y), grid on

% plot Peaks
d   = 0.03;
plot(locs(1), pks(1)+d, "v", 'Color', "#0072BD", 'MarkerFaceColor', "blue")
plot(locs(2), pks(2)+d, "v", 'Color', "#0072BD", 'MarkerFaceColor', "blue")
plot(locs(3), pks(3)+d, "v", 'Color', "#0072BD", 'MarkerFaceColor', "blue")
plot(locs(4), pks(4)+d, "v", 'Color', "#0072BD", 'MarkerFaceColor', "blue")
plot(locs(5), pks(5)+d, "v", 'Color', "#0072BD", 'MarkerFaceColor', "blue")

% plot Prominences
p1  = line([locs(1) locs(1)], [(pks(1)-proms(1)) pks(1)], 'Color', 'red'); 
p2  = line([locs(2) locs(2)], [(pks(2)-proms(2)) pks(2)], 'Color', 'red'); 
p3  = line([locs(3) locs(3)], [(pks(3)-proms(3)) pks(3)], 'Color', 'red'); 
p4  = line([locs(4) locs(4)], [(pks(4)-proms(4)) pks(4)], 'Color', 'red'); 
p5  = line([locs(5) locs(5)], [(pks(5)-proms(5)) pks(5)], 'Color', 'red'); 

text(locs+1, pks, num2str((1:numel(pks))'))
xlabel('x'), ylabel('y'), title('Five-Finger Mountain')
hold off

[Henry]

@Mathieu NOE

Yes a MaxProminance would be very helpful. My final aim is to automatically extract the peak-to-valley dose ratio (average peak height divided by averge valley height) as well as the width of the peaks at a given height, and peak to peak distance. Its likeley that these paramaters could vary alot as I try and optimise them.

I will be running many simulations to do this which is why I want to automate it, i understand @Sam Chak suggestion of manually calculating the prominance, and this would work if there was a way to auomate it.

Im open to work arounds such as finding maxima and minima without using findpeaks function if needed. I may be able to find out the number of peaks I would expect to see if thats helpful?

[Image Analyst]

@Sam Chak I think @Henry is wanting a verbal explanation of prominence. For example, explain using the colored peaks, like why does peak 7 stop at point g, while peak 6 does not stop at point f? Is it because peak 6 is taller in absolute (and relative) height than peak 7? Is the prominence equal to the area of the colored peaks+bases?

[Sam Chak]

Thanks @Image Analyst. To be honest, I believe that the procedure for measuring the prominence of the local highest peak is flawed. I did not notice this issue until @Henry posed the question.

The prominence is measured against a reference level, which is defined as the highest minimum between two minimum levels determined by evaluating the left and right intervals from the peak until a higher peak is encountered or until the end of the signal.

Let us take Peak 6 as an example and start by evaluating the left interval. Peak 2 is higher than Peak 6 in the left interval. The minimum point in the left interval must lie between Peak 6 and the crossing due to Peak 2, as indicated at trough 'd'. In the right interval, there is no peak higher than Peak 6; thus, the minimum point on the right side is at trough 'h'. The reference level is determined by identifying the highest minimum between trough 'd' and trough 'h'. Since trough 'd' is higher than trough 'h', it is selected as the reference level. The prominence of Peak 6 is measured against the reference level at 'd', which is nearly at the base of the signal.

Should there be any If–Else conditions to identify the reference level between trough 'e' and trough 'f'?