How to cut out specific segment from findpeaks function?

조회 수: 18 (최근 30일)
Matlaber
Matlaber 2019년 2월 7일
댓글: Giacomo Echevers 2022년 1월 18일
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:
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
1.png
Then, using findpeaks function,
findpeaks(PeakSig,x,'Annotate','extents','WidthReference','halfheight')
title('Signal Peak Widths')
2.png
I want to cut the segment of the border into 3 segment.
4.PNG
The method I am thinking is find the X value corresponding the Y value, which is not an effective way to do it.

채택된 답변

Star Strider
Star Strider 2019년 2월 8일
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.
  댓글 수: 16
Matlaber
Matlaber 2019년 2월 14일
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
Star Strider 2019년 2월 14일
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).

댓글을 달려면 로그인하십시오.

추가 답변 (3개)

João
João 2019년 2월 7일
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)) = [];
  댓글 수: 3
João
João 2019년 2월 8일
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
Matlaber 2019년 2월 9일
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?

댓글을 달려면 로그인하십시오.


Image Analyst
Image Analyst 2019년 2월 9일
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.
  댓글 수: 7
Image Analyst
Image Analyst 2019년 2월 14일
Looks like you got it working since you've accepted an answer.
Matlaber
Matlaber 2019년 2월 14일
편집: Matlaber 2019년 2월 14일
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

댓글을 달려면 로그인하십시오.


Giacomo Echevers
Giacomo Echevers 2022년 1월 18일
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.
  댓글 수: 2
Image Analyst
Image Analyst 2022년 1월 18일
Helps me anyway. I never knew findpeaks() had those nice annotation lines available for display.
Giacomo Echevers
Giacomo Echevers 2022년 1월 18일
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.

댓글을 달려면 로그인하십시오.

카테고리

Help CenterFile Exchange에서 Axis Labels에 대해 자세히 알아보기

태그

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by