How to cut out specific segment from findpeaks function?

Hi All,

I am trying to figure out how to cut out the specific segment of the signal where I draw the border using findpeaks function.

I am using example of findpeaks link below:

Example Link of findpeaks function

x = linspace(0,1,1000);

Pos = [1 2 3 5 7 8]/10;
Hgt = [4 4 2 2 2 3];
Wdt = [3 8 4 3 4 6]/100;

for n = 1:length(Pos)
    Gauss(n,:) =  Hgt(n)*exp(-((x - Pos(n))/Wdt(n)).^2);
end

PeakSig = sum(Gauss);

plot(x,Gauss,'--',x,PeakSig)
grid

Then, using findpeaks function,

findpeaks(PeakSig,x,'Annotate','extents','WidthReference','halfheight')
title('Signal Peak Widths')

I want to cut the segment of the border into 3 segment.

The method I am thinking is find the X value corresponding the Y value, which is not an effective way to do it.

Answer 1

[Star Strider]

You can get the ‘Borders’ x-coordinates using:

(Previous Code)
findpeaks(PeakSig,x,'Annotate','extents','WidthReference','halfheight')
title('Signal Peak Widths')

Ax = gca;
Kids = Ax.Children;
Borders = Kids(1:2);
Line = Borders.XData;

You have to search the ‘Axes.Children’ result manually to find the ‘(Border)’ elements. There appears to be no way to get them automatically, or even to specifically search for a ‘Line (Border)’ or ‘Border’ entry in the ‘Kids’ variable. (There may be a way, but it would take more patience than I have tonight to find it. I’ve tried every cell and structure addressing approach I can think of or can find in the documentation, including using strcmp and related functions, and couldn’t get a good result.) It would be very nice if this was more straightforward.

[Matlaber]

Thanks for your reply.

Do you know what is the name of the documentation in matlab to read more information on using "Axes", "Border"?

I run the code, it give me some "Ax", "Borders", "Kids" information in my Workspace", but I have no idea how to read it.

[Star Strider]

The ‘Kids’ variable is simply ‘Ax.Children’, the ‘Children’ field of the ‘Ax’structure. The ‘Line’ variable has the x-axis coordinates of the ‘Border’ lines.

Do you know what is the name of the documentation in matlab to read more information on using "Axes", "Border"?

I can’t specifically access the ‘Border’ Line objects by name, so I need to simply refer to them by subscript reference. See the documentation on Axes Properties (link) and Access Data in a Structure Array (link).

An easier way to find the ‘borders’ is likely to invert (actually negate) the ‘PeakSig’ array(here) and then find those locations:

[Pks,Locs] = findpeaks(-PeakSig,x);

The value of the peaks returned is now ‘-Pks’, ‘Locs’ is correct as returned.

[Matlaber]

Thanks for your reply.

I tried:

[Pks,Locs] = findpeaks(-PeakSig,x);

result:

>> Pks

Pks =

   -2.9151   -2.9962   -0.0039   -0.0011   -1.8394

>> Locs

Locs =

    0.1451    0.2663    0.4184    0.5866    0.7397

While the original peak:

pks =
    4.8799    4.0040    3.0328    1.9994    2.2135    3.0039
locs =
    0.1031    0.2002    0.2833    0.4995    0.7057    0.7998

I still couldn't get the meaning.

[Star Strider]

The ‘Locs’ corresponding to the result of

[Pks,Locs] = findpeaks(-PeakSig,x);

are the ‘borders’ you want.

[Matlaber]

Thank you for your reply.

Yes, correct. I tried the code below and check with the plotted graph, the "Locs" is the width value.

[Pks,Locs] = findpeaks(-PeakSig,x)

Can I know the reason behind? "-PeakSig" mean the opposite/negative/reflect side of the original "PeakSig"?

Thanks.

[Star Strider]

Taking the negative of ‘PeakSig’ inverts it, so the peaks are valleys amd the valleys are peaks. The findpeaks function then does what it si designed to do, and finds the peaks of the inverted (negated) signal.

[Matlaber]

Thanks. Is inverted means Up side down?

Original signal

Inverted signal

[Star Strider]

Yes, although inverting it and then using the original code means that the peaks are the previous borders, and the borders are the previous peaks.

It’s not really that complicated. It’s just different.

[Matlaber]

Thanks for your reply.

Looks complicated to me. As border is a line beside peaks. Is not directly opposite/reflect on peak.

Do you look at the findpeaks.m and analysis it?

[Star Strider]

My pleasure.

I did not look at the findpeaks code. I have no desire to, since it is obvious that the 'borders’ are simply the ‘peaks’ of the negated (inverted) signal.

Run this, and you will find that the red ‘+’ markers exactly correspond to the ‘borders’, and are plotted at the ends of the ‘border’ lines:

x = linspace(0,1,1000);
Pos = [1 2 3 5 7 8]/10;
Hgt = [4 4 2 2 2 3];
Wdt = [3 8 4 3 4 6]/100;
for n = 1:length(Pos)
    Gauss(n,:) =  Hgt(n)*exp(-((x - Pos(n))/Wdt(n)).^2);
end
PeakSig = sum(Gauss);
plot(x,Gauss,'--',x,PeakSig)
grid


findpeaks(PeakSig,x,'Annotate','extents','WidthReference','halfheight')
title('Signal Peak Widths')
hold on
[pks,locs] = findpeaks(-PeakSig,x)
plot([locs; locs], [zeros(size(pks)); -pks], '+r', 'MarkerSize',10)
hold off

QUOD ERAT DEMONSTRANDUM

[Matlaber]

Thanks a lot. I am much appreciated.

g

[Star Strider]

As always, my pleasure.

[Matlaber]

Thanks.

I am not sure I am correctly checking this.

It seemed like missing the first border.

load sunspot.dat
year = sunspot(:,1);
avSpots = sunspot(:,2);

findpeaks(avSpots,year,'Annotate','extents','WidthReference','halfheight')
title('Signal Peak Widths')
hold on

[pks,locs] = findpeaks(-avSpots,year)
plot([locs; locs], [zeros(size(pks)); -pks], '+r', 'MarkerSize',10)
hold off

[Star Strider]

‘It seemed like missing the first border.’

The ends of a signal are not ‘peaks’ (or ‘valleys’) because there is nothing on the other sides of them to define them as such. This is the expected behaviour of findpeaks. Note that the other end of the signal is also not identified as either a ‘peak’ or a ‘valley’.

[Matlaber]

Thanks! I am appreciate it.

Let say I want to plot the first peak (range from 1700 to 1711), it turn out to be a full range of "pks" and "locs" value.

I am not sure can I plot only specific range and then find the pks in that specific range?

load sunspot.dat
year = sunspot(:,1);
avSpots = sunspot(:,2);

findpeaks(avSpots,year,'Annotate','extents','WidthReference','halfheight')
title('Signal Peak Widths')
axis([1700 1711 0 100])
hold on

[pks,locs] = findpeaks(-avSpots,year)
plot([locs; locs], [zeros(size(pks)); -pks], '+r', 'MarkerSize',10)
hold off

Results:

pks =

         0
  -11.0000
   -5.0000
   -5.0000
  -47.7000
   -9.6000
  -11.4000
   -7.0000
  -10.2000
   -4.1000
  -43.1000
         0
   -1.8000
   -8.5000
  -10.7000
   -4.3000
  -44.0000
   -7.3000
  -11.3000
   -3.4000
   -6.3000
  -26.2000
   -2.7000
  -53.8000
   -1.4000
   -5.8000
   -5.7000
   -9.6000
   -4.4000
  -10.2000
  -66.6000
  -12.6000
  -13.4000


locs =

        1711
        1723
        1733
        1744
        1751
        1755
        1766
        1775
        1784
        1798
        1803
        1810
        1823
        1833
        1843
        1856
        1863
        1867
        1876
        1878
        1889
        1897
        1901
        1906
        1913
        1923
        1933
        1944
        1954
        1964
        1971
        1976
        1986

Maybe I should change the code in ploting in a specific range, but i run out of idea.

[Star Strider]

‘I am not sure can I plot only specific range and then find the pks in that specific range?’

You can. Starting with:

[pks,locs] = findpeaks(-avSpots,year)

(or something similar) then calculate the indices as as:

idx = find((year >= yr1) & (year <= yr2));

where ‘yr1’ is the beginning year and ‘yr2’ is the end year, then analyse your data as:

spots_sub = avSpots(idx);
year_sub = year(idx);

That should work (I did not test it).

Answer 2

[João]

Try this,

[pks,locs] = findpeaks(PeakSig,x,'Annotate','extents','WidthReference','halfheight')

The variable locs have the indexes of the peaks found. Use them like this to remove a specific peak (your case the third):

PeakSig(locs(3)) = [];

[Matlaber]

Thanks for your reply.

Can you elebrate more on the following code meaning:

PeakSig(locs(3)) = [];

Is that mean PeakSig find the points of "locs" at 3rd point?

If I want to find 1st point, it will be:

PeakSig(locs(1)) = [];

However, when I run it, it give error:

>> PeakSig(locs(3)) = []
Subscript indices must either be real positive integers or logicals.
 

[João]

I didn't notice your were using x in findpeaks. You need to find the index that corresponds to locs. Do it like this then:

[pks,locs] = findpeaks(PeakSig,x,'Annotate','extents','WidthReference','halfheight')
peakToFind = 1; % the peak you want to delete (in this case the first one)
% where you get the postion in x corresponding to locs
idx_x = find(x==locs(peakToFind));
PeakSig(idx_x) = [];

[Matlaber]

Thanks for your reply.

I managed to find the peak point, which is below:

pks =

    4.8799    4.0040    3.0328    1.9994    2.2135    3.0039


locs =

    0.1031    0.2002    0.2833    0.4995    0.7057    0.7998

Thus,

peakToFind = 1; % the peak you want to delete (in this case the first one)
% where you get the postion in x corresponding to locs
idx_x = find(x==locs(peakToFind));
PeakSig(idx_x) 

Give the answer:

ans =

    4.8765

which is the first peak answer.

However, how can I plot of first segment of the peak?

Answer 3

[Image Analyst]

The borders appear to be available in findpeaks() at line 558 (in R2018b version):

% get the x-coordinates of the half-height width borders of each peak
[wxPk,iLBh,iRBh] = getPeakWidth(yFinite,x,iPk,bPk,iLB,iRB,refW);

See it with edit, then search for border

>> edit findpeaks.m

Make a copy of findpeaks.m and call it findpeaksandborders() in some utilities folder on your path with all your other m-files. BE SURE NOT TO ALTER THE ORIGINAL ONE!!!

Then have findpeaksandborders() return the borders (iLB and iRB) in the list of output arguments.

[Matlaber]

Thanks for your reply.

Yes, I also thinking of the boders appear in built in function of "findpeaks".

I managed to find the findpeaks.m and the line 558:

% get the x-coordinates of the half-height width borders of each peak
[wxPk,iLBh,iRBh] = getPeakWidth(yFinite,x,iPk,bPk,iLB,iRB,refW);

And make a copy of findpeaks.m.

So I have "findpeaksandborders.m" and "findpeaks.m" in my "C:\Program Files\MATLAB\R2017b\toolbox\signal\signal" folder.

However, I do not understand what should I change in "findpeaksandborders.m" file to get have findpeaksandborders() return the borders (iLB and iRB) in the list of output arguments.

[Image Analyst]

Well the quick way is to just declare it global and send it out. Down at line 558:

% get the x-coordinates of the half-height width borders of each peak
global iLB;
global iRB;
[wxPk,iLBh,iRBh] = getPeakWidth(yFinite,x,iPk,bPk,iLB,iRB,refW);

Then at line 1 of findpeaksandborders.m:

function [Ypk,Xpk,Wpk,Ppk, iLB, iRB] = findpeaksandborders(Yin,varargin)
global iLB;
global iRB;

Then when youi call findpeaksandborders() just make sure you accept all 6 returned variables.

[Matlaber]

Thanks for your reply.

I have 2 more short questions:

1. After I edit my findpeaksandborders.m file by adding in the global value (following your steps), I can now use/call my findpeaksandborders() function.

[Ypk,Xpk,Wpk,Ppk, iLB, iRB] = findpeaksandborders(Yin,varargin)
iLB;
iRB;

2. I am still figuring out why I am unable to edit (Access is denied). Maybe I am using university matlab licence PC.

Thank you.

[Image Analyst]

You're not allowed to change anything under Program Files. That's a Windows rule.

You'll have to save a copy of the m-file where you have the rest of your m-files, and edit it there.

[Matlaber]

Thanks. I tried edit the "findpeaksandborders.m" in my matlab directory (not in matlab file under Program File) following your suggestion.

[pks,locs,Wpk,Ppk, iLB, iRB] = findpeaksandborders(PeakSig,x)
iLB;
iRB;

It said:

Error using findpeaksandborders
Too many output arguments.

Error in Find_All_Peak (line 315)
[pks,locs,Wpk,Ppk, iLB, iRB] = findpeaksandborders(PeakSig,x)

In your words "make sure you accept all 6 returned variables". Any idea which part should I check?

Thanks.

[Image Analyst]

Looks like you got it working since you've accepted an answer.

[Matlaber]

I used the above option. I still unable to figure out the error:

Error using findpeaksandborders
Too many output arguments.

Error in Find_All_Peak (line 315)
[pks,locs,Wpk,Ppk, iLB, iRB] = findpeaksandborders(PeakSig,x)

However, it cannot find the first border.

I think it is due to invert of peak.

load sunspot.dat

year = sunspot(:,1);
avSpots = sunspot(:,2);

findpeaks(avSpots,year,'Annotate','extents','WidthReference','halfheight')
title('Signal Peak Widths')
hold on
[pks,locs] = findpeaks(-avSpots,year)
plot([locs; locs], [zeros(size(pks)); -pks], '+r', 'MarkerSize',10)
hold off

Answer 4

[Giacomo Echevers]

Well, I kind of needed exactly this feature so I ventured and fought with the findpeaks code, luckily I found what I needed. I attached the modified version of the findpeaks function (now called Findpeaks.m).

Note: It seems that because MATLAB recognized that this is a custom function, the plotting features are not working, but that's no issue, because you can simply use the orginial function and that'll do, this mod is just to extract the x position for the borders.

I'm going to include here a small example using the previous presented case.

x = linspace(0,1,1000);
Pos = [1 2 3 5 7 8]/10;
Hgt = [4 4 2 2 2 3];
Wdt = [3 8 4 3 4 6]/100;
for n = 1:length(Pos)
   Gauss(n,:) =  Hgt(n)*exp(-((x - Pos(n))/Wdt(n)).^2);
end
PeakSig = sum(Gauss);
findpeaks(PeakSig,x,'Annotate','extents','WidthReference','halfheight');

Now, let's use the modified version of the findpeaks function to get the borders x values. Please note the use of the pair parameters 'BorderExtraction' and logical 1. Basically, logical 1 to get the x values for borders in the loc output variable position and logical 0 to tell it to behave exactly the same as the original findpeaks function.

[~,X] = Findpeaks(PeakSig,x,'Annotate','extents','WidthReference','halfheight','BorderExtraction',1)

X = 6×2

0 0.1451 0.1451 0.2663 0.2663 0.4184 0.4184 0.5866 0.5866 0.7397 0.7397 1.0000

And now let's use this to get the second peak just for fun.

plot(PeakSig(x>=X(2) & x<=X(3)))

Hope this helps

A94 out.

[Image Analyst]

Helps me anyway. I never knew findpeaks() had those nice annotation lines available for display.

[Giacomo Echevers]

Glad to hear that! Just noticed those borders, literally yesterday and got quite triggered that it didn't have a native way to extract those x coordinates.