-
Notifications
You must be signed in to change notification settings - Fork 4
/
exridge.m
70 lines (60 loc) · 1.93 KB
/
exridge.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
function [c,e] = exridge(Tx,lambda,beta,jump)
% exridge : extracts the ridge curve by maximising some energy.
% The result is an approximation computd by a greedy algorithm.
% The algorithm uses several random initializations and then
% forward/backward search.
%
% INPUTS:
% Tx : synchrosqueezed transform
% lambda : the paremeter associated with first derivative penalization
% beta : the parameter associated with the second derivative penalization
% OUTPUTS:
% c : vector containing the indexes of the ridge. length :size(Tx,2).
% e : energy of the returned ridge, scalar.
%Et = log(abs(Tx)+eps^0.25);
Et = abs(Tx).^2;
[na,N] = size(Tx);
% Parameters of the algorithm.
da = jump;
ng = min(60,floor(N/8)); % Number of random initializations.
ccur = zeros(1,N);
c = ccur;
e = -Inf;
for k = floor(linspace(N/(ng+1),N-N/(ng+1),ng))
[ecur,idx] = max(Et(:,k));
ccur(k) = idx;
Iq = max(1,idx-da):min(na,idx+da);
[ecur,idx] = max(Et(Iq,k-1));
ccur(k-1) = idx;
idx = ccur(k);
% forward step
for b=k+1:N
etmp = -Inf;
for a=max(1,idx-da):min(na,idx+da)
if Et(a,b)-lambda*(N/(2*na))^2*(a-ccur(b-1))^2-beta*(N/(2*na))^2*(a-2*ccur(b-1)+ccur(b-2))^2 > etmp
etmp = Et(a,b)-lambda*(N/(2*na))^2*(a-ccur(b-1))^2-beta*(N/(2*na))^2*(a-2*ccur(b-1)+ccur(b-2))^2;
idx = a;
end
end
ccur(b)=idx;
ecur = ecur + etmp;
end
% backward step
idx = ccur(k);
for b=k-1:-1:1
etmp = -Inf;
for a=max(1,idx-da):min(na,idx+da)
if Et(a,b)-lambda*(N/(2*na))^2*(a-ccur(b+1))^2-beta*(N/(2*na))^2*(a-2*ccur(b+1)+ccur(b+2))^2 > etmp
etmp = Et(a,b)-lambda*(N/(2*na))^2*(a-ccur(b+1))^2-beta*(N/(2*na))^2*(a-2*ccur(b+1)+ccur(b+2))^2;
idx = a;
end
end
ccur(b)=idx;
ecur = ecur + etmp;
end
if ecur> e
e = ecur;
c = ccur;
end
end
end