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