bmailhe@200: function dico = dico_decorr_symetric(dico, mu)
daniele@195:     %DICO_DECORR decorrelate a dictionary
daniele@195:     %   Parameters:
bmailhe@200:     %   dico: the dictionary, either a matrix or a cell array of matrices.
daniele@195:     %   mu: the coherence threshold
daniele@195:     %
daniele@195:     %   Result:
bmailhe@200:     %   dico: if the input dico was a matrix, then a matrix close to the 
bmailhe@200:     %   input one with coherence mu.
bmailhe@200:     %   If the input was a cell array, a cell array of the same size
bmailhe@200:     %   containing matrices such that the coherence between different cells
bmailhe@200:     %   is lower than mu.
daniele@195:     
bmailhe@200:     eps = 1e-3; % 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:     
bmailhe@200:     % if only one dictionary is provided, then decorrelate it
bmailhe@200:     if ~iscell(dico)
bmailhe@200:         % several decorrelation iterations might be needed to reach global
bmailhe@200:         % coherence mu. niter can be adjusted to needs.
bmailhe@200:         niter = 1;
bmailhe@200:         while max(max(abs(dico'*dico -eye(length(dico))))) > mu + eps
bmailhe@200:             % find pairs of high correlation atoms
bmailhe@200:             colors = dico_color(dico, mu);
bmailhe@200:             
bmailhe@200:             % iterate on all pairs
bmailhe@200:             nbColors = max(colors);
bmailhe@200:             for c = 1:nbColors
bmailhe@200:                 index = find(colors==c);
bmailhe@200:                 if numel(index) == 2
bmailhe@200:                     if dico(:,index(1))'*dico(:,index(2)) > 0
bmailhe@200:                         %build the basis vectors
bmailhe@200:                         v1 = dico(:,index(1))+dico(:,index(2));
bmailhe@200:                         v1 = v1/norm(v1);
bmailhe@200:                         v2 = dico(:,index(1))-dico(:,index(2));
bmailhe@200:                         v2 = v2/norm(v2);
bmailhe@200:                         
bmailhe@200:                         dico(:,index(1)) = ctheta*v1+stheta*v2;
bmailhe@200:                         dico(:,index(2)) = ctheta*v1-stheta*v2;
bmailhe@200:                     else
bmailhe@200:                         v1 = dico(:,index(1))-dico(:,index(2));
bmailhe@200:                         v1 = v1/norm(v1);
bmailhe@200:                         v2 = dico(:,index(1))+dico(:,index(2));
bmailhe@200:                         v2 = v2/norm(v2);
bmailhe@200:                         
bmailhe@200:                         dico(:,index(1)) = ctheta*v1+stheta*v2;
bmailhe@200:                         dico(:,index(2)) = -ctheta*v1+stheta*v2;
bmailhe@200:                     end
bmailhe@200:                 end
bmailhe@200:             end
bmailhe@200:             niter = niter+1;
bmailhe@200:         end
bmailhe@200:         %if a cell array of dictionaries is provided, decorrelate among
bmailhe@200:         %different dictionaries only
bmailhe@200:     else
bmailhe@200:         niter = 1;
bmailhe@200:         numDicos = length(dico);
bmailhe@200:         G = cell(numDicos);
bmailhe@200:         maxCorr = 0;
bmailhe@200:         for i = 1:numDicos
bmailhe@200:             for j = i+1:numDicos
bmailhe@200:                 G{i,j} = dico{i}'*dico{j};
bmailhe@200:                 maxCorr = max(maxCorr,max(max(abs(G{i,j}))));
bmailhe@200:             end
bmailhe@200:         end
daniele@195:         
bmailhe@200:         while maxCorr > mu + eps
bmailhe@200:             % find pairs of high correlation atoms
bmailhe@200:             [colors nbColors] = dico_color_separate(dico, mu);
bmailhe@200: 
bmailhe@200:             % iterate on all pairs
bmailhe@200:             for c = 1:nbColors
bmailhe@200:                 for tmpI = 1:numDicos
bmailhe@200:                     index = find(colors{tmpI}==c);
bmailhe@200:                     if ~isempty(index)
bmailhe@200:                         i = tmpI;
bmailhe@200:                         m = index;
bmailhe@200:                         break;
bmailhe@200:                     end
bmailhe@200:                 end
bmailhe@200:                 for tmpJ = i+1:numDicos
bmailhe@200:                     index = find(colors{tmpJ}==c);
bmailhe@200:                     if ~isempty(index)
bmailhe@200:                         j = tmpJ;
bmailhe@200:                         n = index;
bmailhe@200:                         break;
bmailhe@200:                     end
bmailhe@200:                 end
bmailhe@200:                 
bmailhe@200:                 if dico{i}(:,m)'*dico{j}(:,n) > 0
daniele@195:                     %build the basis vectors
bmailhe@200:                     v1 = dico{i}(:,m)+dico{j}(:,n);
daniele@195:                     v1 = v1/norm(v1);
bmailhe@200:                     v2 = dico{i}(:,m)-dico{j}(:,n);
daniele@195:                     v2 = v2/norm(v2);
daniele@195:                     
bmailhe@200:                     dico{i}(:,m) = ctheta*v1+stheta*v2;
bmailhe@200:                     dico{j}(:,n) = ctheta*v1-stheta*v2;
daniele@195:                 else
bmailhe@200:                     v1 = dico{i}(:,m)-dico{j}(:,n);
daniele@195:                     v1 = v1/norm(v1);
bmailhe@200:                     v2 = dico{i}(:,m)+dico{j}(:,n);
daniele@195:                     v2 = v2/norm(v2);
daniele@195:                     
bmailhe@200:                     dico{i}(:,m) = ctheta*v1+stheta*v2;
bmailhe@200:                     dico{j}(:,n) = -ctheta*v1+stheta*v2;
bmailhe@200:                 end
bmailhe@200:             end
bmailhe@200:             niter = niter+1;
bmailhe@200:             
bmailhe@200:             % Remove noegative components and renormalize
bmailhe@200:             for i = 1:length(dico)
bmailhe@200:                 dico{i} = max(dico{i},0);
bmailhe@200:                 for m = 1:size(dico{i},2)
bmailhe@200:                     dico{i}(:,m) = dico{i}(:,m)/norm(dico{i}(:,m));
bmailhe@200:                 end
bmailhe@200:             end
bmailhe@200:             
bmailhe@200:             maxCorr = 0;
bmailhe@200:             for i = 1:numDicos
bmailhe@200:                 for j = i+1:numDicos
bmailhe@200:                     G{i,j} = dico{i}'*dico{j};
bmailhe@200:                     maxCorr = max(maxCorr,max(max(abs(G{i,j}))));
daniele@195:                 end
daniele@195:             end
daniele@195:         end
daniele@195:     end
daniele@195: end