-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathfindMultiTresholds.m
More file actions
119 lines (92 loc) · 2.24 KB
/
findMultiTresholds.m
File metadata and controls
119 lines (92 loc) · 2.24 KB
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
118
119
function T = findMultiTresholds(sortedList,th, orderLabel)
tTask = 1;
Q = sort(sortedList);
N = length(Q);
t_exp = th; % experical
G = [];
temp_D2_Q = 0;
temp_D1_Q = 0;
D1_Q = [];
D2_Q = [];
D = [];
for i = 1:N
%-----------the 1st and 2nd order derivative.------------------
if i>1 && i < N
temp_D1_Q = (Q(i+1) - Q(i-1))/2;
end
D1_Q = [D1_Q, temp_D1_Q];
if i>2 && i < N-1 % the second order
temp_D2_Q = abs((Q(i+1)+Q(i-2) - Q(i-1) - Q(i)))/4;
end
D2_Q = [D2_Q, temp_D2_Q];
%-----------------------------------------------------
end
if orderLabel == 1
D = Q'.*D1_Q;
elseif orderLabel == 2
D = Q'.*Q'.* D2_Q;
end
Qm = [];
for i = 1:N
m = mean(Q(1:i));
Qm = [Qm,D(i)./m];
end
figure
plot(Q,'LineWidth',1)
hold on
plot(Qm,'-r','LineWidth',1)
set(get(gca, 'Title'), 'String', 'Q & Qm');
hold on
legend('Q','Qm');
grid on
%-----------------find T------------------
if tTask == 0
indx = find(Qm >= t_exp,1,'first');
if indx == N
T = Q(indx);
else
T = 0.5*(Q(indx)+Q(indx+1)); %(indx+1)
end
if ~isempty(T)
T = T(1);
else
T = Q(end); % if T dosen't exsit, T = max(Q).
end
elseif tTask==1
indx = find(Qm >= t_exp);
fInd = [];
for j = 1:length(indx)-1
jd = indx(j)-1;
jb = indx(j);
if Qm(jd) <= th && Qm(jb) >= th
fInd = [fInd indx(j)];
end
end
for k = 1:length(fInd)
if fInd(k) == N
T(k) = Q(N);
else
T(k) = 0.5*(Q(fInd(k))+Q(fInd(k)+1)); %(indx+1)
end
end
T= [0,T];
else
disp('wrong task!!')
end
%--------------------plot T in Q graph-------
%{
xt = [find(Q >= T, 1 ) find(Q >= T, 1 )]; % actually, no need to +1;
yt = [0 T];
plot(xt,yt,'--m','LineWidth',2)
plot(xt,yt,'.k','MarkerSize',20)
hold on
xt = [0 find(Q >= T, 1 )]; % no need to +1
yt = [T T];
plot(xt,yt,'--m','LineWidth',1.5)
plot(xt,yt,'.k','MarkerSize',20)
xt = [0 find(Q >= T, 1 )]; % no need to +1; th = 0.12
yt = [t_exp t_exp];
plot(xt,yt,'--g','LineWidth',1.5)
plot(xt,yt,'.k','MarkerSize',20)
hold off
%}