# HG changeset patch
# User Daniele Barchiesi <daniele.barchiesi@eecs.qmul.ac.uk>
# Date 1314795143 -3600
# Node ID 88578ec2f94aab5af0d74de61aebcf5838225a3c
# Parent  e3035d45d0145a03239a9c8efb15ad3f633312c3
Updated grassmannian function and minor debugs

diff -r e3035d45d014 -r 88578ec2f94a DL/two-step DL/SMALL_two_step_DL.m
--- a/DL/two-step DL/SMALL_two_step_DL.m	Wed Aug 31 10:53:10 2011 +0100
+++ b/DL/two-step DL/SMALL_two_step_DL.m	Wed Aug 31 13:52:23 2011 +0100
@@ -119,7 +119,7 @@
                 dico = dico_decorr(dico, mu, solver.solution);
             case 'tropp'
                 [n m] = size(dico);
-                dico = grassmanian(n,m,[],[],[],dico,true);
+                dico = grassmannian(n,m,[],[],[],dico,true);
             otherwise
         end
     
diff -r e3035d45d014 -r 88578ec2f94a util/classes/dictionaryMatrices/grassmannian.m
--- a/util/classes/dictionaryMatrices/grassmannian.m	Wed Aug 31 10:53:10 2011 +0100
+++ b/util/classes/dictionaryMatrices/grassmannian.m	Wed Aug 31 13:52:23 2011 +0100
@@ -12,33 +12,41 @@
 if ~exist('initA','var') || isempty(initA), initA = randn(n,m); end %initial matrix
 if ~exist('dd2','var')   || isempty(dd2),   dd2 = 0.95; end %shrinking factor
 if ~exist('dd1','var')   || isempty(dd1),   dd1 = 0.9; end %percentage of coherences to be shrinked
-if ~exist('nIter','var') || isempty(nIter), nIter = 5; end %number of iterations
+if ~exist('nIter','var') || isempty(nIter), nIter = 10; end %number of iterations
 
-%% Compute svd and gramian
+%% Main algo
 A = normc(initA); %normalise columns
-[Uinit Sigma] = svd(A); %calculate svd of the matrix
-G = A'*A; %gramian matrix
+G = A'*A; %gram matrix
 
-muMin = sqrt((m-n)/(n*(m-1)));              %Lower bound on mutual coherence
+muMin = sqrt((m-n)/(n*(m-1)));              %Lower bound on mutual coherence (equiangular tight frame)
 res = zeros(nIter,1);
-for iIter = 1:nIter
-   gg  = sort(abs(G(:))); %sort inner products from less to ost correlated
-   pos = find(abs(G(:))>gg(round(dd1*(m^2-m))) & abs(G(:)-1)>1e-6); 
-   G(pos) = G(pos)*dd2;
-   [U S V] = svd(G);
-   S(n+1:end,1+n:end) = 0;
-   G = U*S*V';
-   G = diag(1./abs(sqrt(diag(G))))*G*diag(1./abs(sqrt(diag(G))));
-   gg = sort(abs(G(:)));
-   pos = find(abs(G(:))>gg(round(dd1*(m^2-m))) & abs(G(:)-1)>1e-6);
-   res(iIter) = max(abs(G(pos)));
-   if verb
-    fprintf(1,'%6i  %12.8f   %12.8f  %12.8f \n',...
-            [iIter,muMin,mean(abs(G(pos))),max(abs(G(pos)))]);
-   end
+if verb
+	fprintf(1,'Iter		mu_min		mu \n');
 end
 
-[~, Sigma_gram V_gram] = svd(G); %calculate svd decomposition of gramian
-Sigma_new = sqrt(Sigma_gram(1:n,:)).*sign(Sigma);
-A = Uinit*Sigma_new*V_gram';
-A = normc(A);
+% optimise gram matrix
+for iIter = 1:nIter
+	gg  = sort(abs(G(:))); %sort inner products from less to most correlated
+	pos = find(abs(G(:))>=gg(round(dd1*(m^2-m))) & abs(G(:)-1)>1e-6); %find large elements of gram matrix
+	G(pos) = G(pos)*dd2;	%shrink large elements of gram matrix
+	[U S V] = svd(G);	%compute new SVD of gram matrix
+	S(n+1:end,1+n:end) = 0; %set small eigenvalues to zero (this ensures rank(G)<=d)
+	G = U*S*V';			%update gram matrix
+	G = diag(1./abs(sqrt(diag(G))))*G*diag(1./abs(sqrt(diag(G)))); %normalise gram matrix diagonal
+	if verb
+		Geye = G - eye(size(G));
+		fprintf(1,'%6i    %12.8f  %12.8f  \n',iIter,muMin,max(abs(Geye(:))));
+	end
+end
+
+% [~, Sigma_gram V_gram] = svd(G); %calculate svd decomposition of gramian
+% Sigma_new = sqrt(Sigma_gram(1:n,:)).*sign(Sigma); %calculate singular values of dictionary
+% A = Uinit*Sigma_new*V_gram';	%update dictionary
+% A = normc(A);					%normalise dictionary
+
+[U S] = svd(G);				%calculate svd decomposition of gramian
+A = sqrt(S(1:n,1:n))*U(:,1:n)';		%calculate valid frame, s.t. A'*A=G
+
+% %% Debug visualization function
+% function plotcart2d(A)
+% compass(A(1,:),A(2,:));
diff -r e3035d45d014 -r 88578ec2f94a util/snr.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/util/snr.m	Wed Aug 31 13:52:23 2011 +0100
@@ -0,0 +1,6 @@
+function x = snr(s,r)
+% SNR calculates the signal-to-noise ratio between the signal s and its
+% representation r defined as:
+% SNR = 20*log10(||s||_2/||s-r||_2)
+err = s-r;
+x = mag2db(norm(s(:))/norm(err(:)));
\ No newline at end of file