daniele@170: function [dic mus srs] = iterativeprojections(dic,mu,Y,X,params)
daniele@169: % grassmanian attempts to create an n by m matrix with minimal mutual
daniele@169: % coherence using an iterative projection method.
daniele@169: %
daniele@170: % REFERENCE
daniele@169: %
daniele@169: 
daniele@169: %% Parameters and Defaults
daniele@170: if ~nargin, testiterativeprojections; return; end
daniele@169: 
daniele@170: if ~exist('params','var') || isempty(param), params = struct; end
daniele@170: if ~isfield(params,'nIter'), params.nIter = 10; end %number of iterations
daniele@170: if ~isfield(params,'eps'),	 params.eps = 1e-9;	 end %tolerance level
daniele@170: [n m] = size(dic);
daniele@169: 
daniele@170: SNR = @(dic) snr(Y,dic*X);									 %SNR function
daniele@170: MU  = @(dic) max(max(abs((dic'*dic)-diag(diag(dic'*dic))))); %coherence function
daniele@169: 
daniele@170: %% Main algorithm
daniele@170: dic = normc(dic);	%normalise columns
daniele@170: alpha = m/n;		%ratio between number of atoms and ambient dimension
daniele@170: 
daniele@170: mus = zeros(params.nIter,1);		%coherence at each iteration
daniele@170: srs = zeros(params.nIter,1);		%signal to noise ratio at each iteration
daniele@170: iIter = 1;
daniele@170: while iIter<=params.nIter && MU(dic)>mu
daniele@170: 	fprintf(1,'Iteration number %u\n', iIter);
daniele@170: 	% calculate snr and coherence
daniele@170: 	mus(iIter) = MU(dic);
daniele@170: 	srs(iIter) = SNR(dic);
daniele@170: 	
daniele@170: 	% calculate gram matrix
daniele@170: 	G = dic'*dic;
daniele@170: 	
daniele@170: 	% project into the structural constraint set
daniele@170: 	H = zeros(size(G));				%initialise matrix
daniele@170: 	ind1 = find(abs(G(:))<=mu);		%find elements smaller than mu
daniele@170: 	ind2 = find(abs(G(:))>mu);		%find elements bigger than mu
daniele@170: 	H(ind1) = G(ind1);				%copy elements belonging to ind1
daniele@170: 	H(ind2) = mu*sign(G(ind2));		%threshold elements belonging to ind2
daniele@170: 	H(1:m+1:end) = 1;				%set diagonal to one
daniele@170: 	
daniele@170: 	% project into spectral constraint set
daniele@170: 	[~ , S, V] = svd(H);
daniele@170: 	%G = alpha*(V(:,1:n)*V(:,1:n)');
daniele@170: 	G = V(:,1:n)*S(1:n,1:n)*V(:,1:n)';
daniele@170: 	
daniele@170: 	% calculate dictionary
daniele@170: 	[~, S V] = svd(G);
daniele@170: 	dic = sqrt(S(1:n,:))*V';
daniele@170: 	
daniele@170: 	% rotate dictionary
daniele@170: 	options = struct('nIter',100,'step',0.001);
daniele@170: 	[~, ~, W] = rotatematrix(Y,dic*X,'conjgradlie',options);
daniele@170: 	dic = W*dic;
daniele@170: 	
daniele@170: 	iIter = iIter+1;
daniele@169: end
daniele@170: if iIter<params.nIter
daniele@170: 	mus(iIter:end) = mus(iIter);
daniele@170: 	srs(iIter:end) = srs(iIter);
daniele@170: end
daniele@170: % Test function
daniele@170: function testiterativeprojections
daniele@170: clc
daniele@170: %define parameters
daniele@170: n = 256;							%ambient dimension
daniele@170: m = 512;							%number of atoms
daniele@170: N = 1024;							%number of signals
daniele@170: mu_min = sqrt((m-n)/(n*(m-1)));		%minimum coherence
daniele@169: 
daniele@170: %initialise data
daniele@170: X = sprandn(m,N,1);					%matrix of coefficients
daniele@170: phi = normc(randn(n,m));			%dictionary
daniele@170: temp = randn(n);
daniele@170: W = expm(0.5*(temp-temp'));			%rotation matrix
daniele@170: Y = W*phi*X;						%observed signals
daniele@169: 
daniele@170: %optimise dictionary
daniele@170: [~, mus srs] = iterativeprojections(phi,0.2,Y,X);
daniele@169: 
daniele@170: %plot results
daniele@170: nIter = length(mus);
daniele@169: 
daniele@170: figure, 
daniele@170: subplot(2,1,1)
daniele@170: plot(1:nIter,srs,'kd-');
daniele@170: xlabel('nIter');
daniele@170: ylabel('snr (dB)');
daniele@170: grid on
daniele@169: 
daniele@170: subplot(2,1,2), hold on
daniele@170: plot(1:nIter,mus,'ko-');
daniele@170: plot([1 nIter],[mu_min mu_min],'k')
daniele@170: grid on
daniele@170: legend('\mu','\mu_{min}');
daniele@169: