bmailhe@200
|
1 function dico = dico_decorr_symetric(dico, mu)
|
daniele@195
|
2 %DICO_DECORR decorrelate a dictionary
|
daniele@195
|
3 % Parameters:
|
bmailhe@200
|
4 % dico: the dictionary, either a matrix or a cell array of matrices.
|
daniele@195
|
5 % mu: the coherence threshold
|
daniele@195
|
6 %
|
daniele@195
|
7 % Result:
|
bmailhe@200
|
8 % dico: if the input dico was a matrix, then a matrix close to the
|
bmailhe@200
|
9 % input one with coherence mu.
|
bmailhe@200
|
10 % If the input was a cell array, a cell array of the same size
|
bmailhe@200
|
11 % containing matrices such that the coherence between different cells
|
bmailhe@200
|
12 % is lower than mu.
|
daniele@195
|
13
|
bmailhe@200
|
14 eps = 1e-3; % define tolerance for normalisation term alpha
|
daniele@195
|
15
|
daniele@195
|
16 % convert mu to the to the mean direction
|
daniele@195
|
17 theta = acos(mu)/2;
|
daniele@195
|
18 ctheta = cos(theta);
|
daniele@195
|
19 stheta = sin(theta);
|
daniele@195
|
20
|
daniele@195
|
21 % compute atom weights
|
daniele@195
|
22 % if nargin > 2
|
daniele@195
|
23 % rank = sum(amp.*amp, 2);
|
daniele@195
|
24 % else
|
daniele@195
|
25 % rank = randperm(length(dico));
|
daniele@195
|
26 % end
|
daniele@195
|
27
|
bmailhe@200
|
28 % if only one dictionary is provided, then decorrelate it
|
bmailhe@200
|
29 if ~iscell(dico)
|
bmailhe@200
|
30 % several decorrelation iterations might be needed to reach global
|
bmailhe@200
|
31 % coherence mu. niter can be adjusted to needs.
|
bmailhe@200
|
32 niter = 1;
|
bmailhe@200
|
33 while max(max(abs(dico'*dico -eye(length(dico))))) > mu + eps
|
bmailhe@200
|
34 % find pairs of high correlation atoms
|
bmailhe@200
|
35 colors = dico_color(dico, mu);
|
bmailhe@200
|
36
|
bmailhe@200
|
37 % iterate on all pairs
|
bmailhe@200
|
38 nbColors = max(colors);
|
bmailhe@200
|
39 for c = 1:nbColors
|
bmailhe@200
|
40 index = find(colors==c);
|
bmailhe@200
|
41 if numel(index) == 2
|
bmailhe@200
|
42 if dico(:,index(1))'*dico(:,index(2)) > 0
|
bmailhe@200
|
43 %build the basis vectors
|
bmailhe@200
|
44 v1 = dico(:,index(1))+dico(:,index(2));
|
bmailhe@200
|
45 v1 = v1/norm(v1);
|
bmailhe@200
|
46 v2 = dico(:,index(1))-dico(:,index(2));
|
bmailhe@200
|
47 v2 = v2/norm(v2);
|
bmailhe@200
|
48
|
bmailhe@200
|
49 dico(:,index(1)) = ctheta*v1+stheta*v2;
|
bmailhe@200
|
50 dico(:,index(2)) = ctheta*v1-stheta*v2;
|
bmailhe@200
|
51 else
|
bmailhe@200
|
52 v1 = dico(:,index(1))-dico(:,index(2));
|
bmailhe@200
|
53 v1 = v1/norm(v1);
|
bmailhe@200
|
54 v2 = dico(:,index(1))+dico(:,index(2));
|
bmailhe@200
|
55 v2 = v2/norm(v2);
|
bmailhe@200
|
56
|
bmailhe@200
|
57 dico(:,index(1)) = ctheta*v1+stheta*v2;
|
bmailhe@200
|
58 dico(:,index(2)) = -ctheta*v1+stheta*v2;
|
bmailhe@200
|
59 end
|
bmailhe@200
|
60 end
|
bmailhe@200
|
61 end
|
bmailhe@200
|
62 niter = niter+1;
|
bmailhe@200
|
63 end
|
bmailhe@200
|
64 %if a cell array of dictionaries is provided, decorrelate among
|
bmailhe@200
|
65 %different dictionaries only
|
bmailhe@200
|
66 else
|
bmailhe@200
|
67 niter = 1;
|
bmailhe@200
|
68 numDicos = length(dico);
|
bmailhe@200
|
69 G = cell(numDicos);
|
bmailhe@200
|
70 maxCorr = 0;
|
bmailhe@200
|
71 for i = 1:numDicos
|
bmailhe@200
|
72 for j = i+1:numDicos
|
bmailhe@200
|
73 G{i,j} = dico{i}'*dico{j};
|
bmailhe@200
|
74 maxCorr = max(maxCorr,max(max(abs(G{i,j}))));
|
bmailhe@200
|
75 end
|
bmailhe@200
|
76 end
|
daniele@195
|
77
|
bmailhe@200
|
78 while maxCorr > mu + eps
|
bmailhe@200
|
79 % find pairs of high correlation atoms
|
bmailhe@200
|
80 [colors nbColors] = dico_color_separate(dico, mu);
|
bmailhe@200
|
81
|
bmailhe@200
|
82 % iterate on all pairs
|
bmailhe@200
|
83 for c = 1:nbColors
|
bmailhe@200
|
84 for tmpI = 1:numDicos
|
bmailhe@200
|
85 index = find(colors{tmpI}==c);
|
bmailhe@200
|
86 if ~isempty(index)
|
bmailhe@200
|
87 i = tmpI;
|
bmailhe@200
|
88 m = index;
|
bmailhe@200
|
89 break;
|
bmailhe@200
|
90 end
|
bmailhe@200
|
91 end
|
bmailhe@200
|
92 for tmpJ = i+1:numDicos
|
bmailhe@200
|
93 index = find(colors{tmpJ}==c);
|
bmailhe@200
|
94 if ~isempty(index)
|
bmailhe@200
|
95 j = tmpJ;
|
bmailhe@200
|
96 n = index;
|
bmailhe@200
|
97 break;
|
bmailhe@200
|
98 end
|
bmailhe@200
|
99 end
|
bmailhe@200
|
100
|
bmailhe@200
|
101 if dico{i}(:,m)'*dico{j}(:,n) > 0
|
daniele@195
|
102 %build the basis vectors
|
bmailhe@200
|
103 v1 = dico{i}(:,m)+dico{j}(:,n);
|
daniele@195
|
104 v1 = v1/norm(v1);
|
bmailhe@200
|
105 v2 = dico{i}(:,m)-dico{j}(:,n);
|
daniele@195
|
106 v2 = v2/norm(v2);
|
daniele@195
|
107
|
bmailhe@200
|
108 dico{i}(:,m) = ctheta*v1+stheta*v2;
|
bmailhe@200
|
109 dico{j}(:,n) = ctheta*v1-stheta*v2;
|
daniele@195
|
110 else
|
bmailhe@200
|
111 v1 = dico{i}(:,m)-dico{j}(:,n);
|
daniele@195
|
112 v1 = v1/norm(v1);
|
bmailhe@200
|
113 v2 = dico{i}(:,m)+dico{j}(:,n);
|
daniele@195
|
114 v2 = v2/norm(v2);
|
daniele@195
|
115
|
bmailhe@200
|
116 dico{i}(:,m) = ctheta*v1+stheta*v2;
|
bmailhe@200
|
117 dico{j}(:,n) = -ctheta*v1+stheta*v2;
|
bmailhe@200
|
118 end
|
bmailhe@200
|
119 end
|
bmailhe@200
|
120 niter = niter+1;
|
bmailhe@200
|
121
|
bmailhe@200
|
122 % Remove noegative components and renormalize
|
bmailhe@200
|
123 for i = 1:length(dico)
|
bmailhe@200
|
124 dico{i} = max(dico{i},0);
|
bmailhe@200
|
125 for m = 1:size(dico{i},2)
|
bmailhe@200
|
126 dico{i}(:,m) = dico{i}(:,m)/norm(dico{i}(:,m));
|
bmailhe@200
|
127 end
|
bmailhe@200
|
128 end
|
bmailhe@200
|
129
|
bmailhe@200
|
130 maxCorr = 0;
|
bmailhe@200
|
131 for i = 1:numDicos
|
bmailhe@200
|
132 for j = i+1:numDicos
|
bmailhe@200
|
133 G{i,j} = dico{i}'*dico{j};
|
bmailhe@200
|
134 maxCorr = max(maxCorr,max(max(abs(G{i,j}))));
|
daniele@195
|
135 end
|
daniele@195
|
136 end
|
daniele@195
|
137 end
|
daniele@195
|
138 end
|
daniele@195
|
139 end
|