daniele@167: function dico = dico_decorr_symetric(dico, mu, amp)
daniele@167:     %DICO_DECORR decorrelate a dictionary
daniele@167:     %   Parameters:
daniele@167:     %   dico: the dictionary
daniele@167:     %   mu: the coherence threshold
daniele@167:     %   amp: the amplitude coefficients, only used to decide which atom to
daniele@167:     %   project
daniele@167:     %
daniele@167:     %   Result:
daniele@167:     %   dico: a dictionary close to the input one with coherence mu.
daniele@167:     
daniele@167:     eps = 1e-6; % define tolerance for normalisation term alpha
daniele@167:     
daniele@167:     % convert mu to the to the mean direction
daniele@167:     theta = acos(mu)/2;
daniele@167:     ctheta = cos(theta);
daniele@167:     stheta = sin(theta);
daniele@167:     
daniele@167:     % compute atom weights
daniele@167:     %     if nargin > 2
daniele@167:     %         rank = sum(amp.*amp, 2);
daniele@167:     %     else
daniele@167:     %         rank = randperm(length(dico));
daniele@167:     %     end
daniele@167:     
daniele@167:     % several decorrelation iterations might be needed to reach global
daniele@167:     % coherence mu. niter can be adjusted to needs.
daniele@167:     niter = 1;
daniele@169:     while max(max(abs(dico'*dico -eye(length(dico))))) > mu + 0.01
daniele@167:         % find pairs of high correlation atoms
daniele@167:         colors = dico_color(dico, mu);
daniele@167:         
daniele@167:         % iterate on all pairs
daniele@167:         nbColors = max(colors);
daniele@167:         for c = 1:nbColors
daniele@167:             index = find(colors==c);
daniele@167:             if numel(index) == 2
daniele@167:                 if dico(:,index(1))'*dico(:,index(2)) > 0               
daniele@167:                     %build the basis vectors
daniele@167:                     v1 = dico(:,index(1))+dico(:,index(2));
daniele@167:                     v1 = v1/norm(v1);
daniele@167:                     v2 = dico(:,index(1))-dico(:,index(2));
daniele@167:                     v2 = v2/norm(v2);
daniele@167:                     
daniele@167:                     dico(:,index(1)) = ctheta*v1+stheta*v2;
daniele@167:                     dico(:,index(2)) = ctheta*v1-stheta*v2;
daniele@167:                 else
daniele@167:                     v1 = dico(:,index(1))-dico(:,index(2));
daniele@167:                     v1 = v1/norm(v1);
daniele@167:                     v2 = dico(:,index(1))+dico(:,index(2));
daniele@167:                     v2 = v2/norm(v2);
daniele@167:                     
daniele@167:                     dico(:,index(1)) = ctheta*v1+stheta*v2;
daniele@167:                     dico(:,index(2)) = -ctheta*v1+stheta*v2;
daniele@167:                 end
daniele@167:             end
daniele@167:         end
daniele@167:         niter = niter+1;
daniele@167:     end
daniele@167: end
daniele@167: