-
Notifications
You must be signed in to change notification settings - Fork 0
/
findpeaks.m
117 lines (114 loc) · 2.99 KB
/
findpeaks.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
function [acme] = findpeaks(y, thresh, mode, strtitle)
%trim vortex core data
% first find the peak in height direction
% loc is in wind frame, first valley starting from zmax to zmin
% zmin is an update of zmin (last valley), should decrease monotonically
acme = findacme(y);
thresh = min(thresh, autothresh(acme));
while 1
[flag, acme] = trimacme(acme, thresh);
if flag==0
break;
end
end
if nargin>=3 && mode>0
figure;
plot(y);
hold on
plot(acme(:,1), acme(:,2),'or');
hold off;
title(strtitle);
end
end
function acme = findacme(h)
nsize = length(h);
flag = 0;
acme(1,1)=1;
acme(1,2)=h(1);
acount = 2;
for ii=2:1:nsize-1
if h(ii) <= h(ii-1) && h(ii)<=h(ii+1) && (h(ii)<h(ii-1) || h(ii)<h(ii+1));
if flag==-1;
acme(acount-1, 1) = round((acme(acount-1, 1) + ii)/2);
else
acme(acount, 1) = ii;
acme(acount, 2) = h(ii);
acount = acount + 1;
end
flag = -1.;
end
if h(ii) >= h(ii-1) && h(ii)>=h(ii+1) && (h(ii)>h(ii-1) || h(ii)>h(ii+1));
if flag==1;
acme(acount-1, 1) = round((acme(acount-1, 1) + ii)/2);
else
acme(acount, 1) = ii;
acme(acount, 2) = h(ii);
acount = acount + 1;
end
flag = 1.;
end
end
if nsize>1;
acme(acount, 1) = nsize;
acme(acount, 2) = h(nsize);
end
end
function thresh = autothresh(acme)
len = length(acme(:,2));
thresh = 0.8*( max( abs( acme(2:len,2)-acme(1:len-1,2) ) ) );
end
function [flag, acme2] = trimacme(acme, thresh)
nsize = length(acme);
flag = 0;
for ii=1:1:nsize-1
if abs(acme(ii+1,2) - acme(ii,2)) >= thresh;
sign = 1;
if acme(ii+1,2) < acme(ii, 2);
sign = -1;
end
if ii<=nsize-3 && abs(acme(ii+1,2)-acme(ii+2,2))<thresh;
if sign*acme(ii+3,2)>sign*acme(ii+1,2);
acme(ii+1,1) = -1;
acme(ii+2,1 ) =-1;
else
acme(ii+2,1) = -1;
acme(ii+3,1 ) =-1;
end
flag = 1;
break;
end
if ii >=3 && abs(acme(ii,2)-acme(ii-1,2))<thresh;
if sign*acme(ii-2,2)<sign*acme(ii,2);
acme(ii,1) = -1;
acme(ii-1,1 ) =-1;
else
acme(ii-1,1) = -1;
acme(ii-2,1 ) =-1;
end
flag = 1;
break;
end
if ii==nsize-2 && abs(acme(ii+1,2)-acme(ii+2,2))<thresh;
acme(ii+2,1 ) =-1;
flag = 1;
break;
end
if ii==2 && abs(acme(ii,2)-acme(ii-1,2)) < thresh;
acme(ii-1, 1) = -1;
flag = 1;
break
end
end
end
if flag==1
count = 1;
for ii=1:1:nsize
if acme(ii,1)>0;
acme2(count, :) = acme(ii, :);
count = count + 1;
end
end
else
acme2 = acme;
end
end