clear all

n=1024;	% Signal length.

% Save figures: 1/0 (default 0).
savfig = 0;

% Signal.
% Four cosines have close frequencies, not in the dictionary
fineness = 4;
t = ((1:n)' - .5) / n;freq = pi * ((1:(4*n))' - 1) / (fineness*n);
freq1 = pi * (126.55-1) / (fineness*n);
freq2 = pi * (127.55-1) / (fineness*n);
freq3 = pi * (128.55-1) / (fineness*n);
freq4 = pi * (129.55-1) / (fineness*n);
const = (2/n) ^ .5;
x1 = const * cos(pi * ((126.55  - 1) / fineness) * t);
x2 = const * cos(pi * ((127.55  - 1) / fineness) * t);
x3 = const * cos(pi * ((128.55  - 1) / fineness) * t);
x4 = const * cos(pi * ((129.55  - 1) / fineness) * t);
zerosn = zeros(n,1);
EnergyRatio = 1E1;
y0dct = EnergyRatio*(x1 + x2 + x3 + x4);
y0dir = SparseVector(n, 3, 'UNIFORM', true);
y0dir(find(y0dir)) = normsignal(y0dir(find(y0dir)),1,2);
pos = find(y0dir);
y0 = y0dct + y0dir;

% Representation dictionary: overcomplete DCT+DIRAC.
dict  = MakeList('DCT','DIRAC');
pars1 = MakeList(fineness,0);
pars2 = MakeList(0,0);
pars3 = MakeList(0,0);

% Coefficient vector length.
p = SizeOfDict1(n,dict,pars1,pars2,pars3);
tightFrame = p/n; % The dictrionary is a tight frame with constant p/n.

% Observed noisy data.
SNR = 15;
sigma = std(y0)*10^(-SNR/20);
y = y0 + sigma*randn(size(y0));

% Solve Lasso in l1-penalized form with FISTA. 
% lambda chosen such that ||r|| <= epsilon.
options.niter = 300;
options.method = 'fista';
options.verb  = 0;
lambda = 0.9;
K  = @(x)FastS1(x,n,dict,pars1,pars2,pars3);
KS = @(x)FastA1(x,dict,pars1,pars2,pars3);
ProxF = @(x,mu)SoftThresh(x,mu*lambda);
GradG = @(x)KS(K(x) - y);
options.report = @(x)norm(y-K(x))^2/2+lambda*norm(x,1);
tic;[xhatlassofista,Efista] = fb(zeros(p,1), ProxF, GradG, 1/tightFrame, options);timelassofista=toc

% Solve BP with DR.
options.gamma = 0.1;    % Stepsize parameter for Douglas-Rachford iteration.
options.tau   = 1;    % No relaxation.
options.niter = 300;
options.verb  = 0;
ProxF = @(x,ga)SoftThresh(x,ga);
ProxG = @(x,ga)x + KS(y0 - K(x))/tightFrame;
tic;[xhatBPDR,EBP] = dr(zeros(p,1),ProxF,ProxG,options);timeBPDR=toc

% Solve BPDN with DR.
options.gamma = 0.1;    % Stepsize parameter for Douglas-Rachford iteration.
options.tau   = 1;    % No relaxation.
options.niter = 300;
options.verb  = 0;
epsilon = sqrt(n)*sigma*sqrt(1+2*sqrt(2/n));	      
ProxF = @(x,ga)SoftThresh(x,ga);
ProxG = prox_F_tightframe(@(x,ga)compute_l2ball_projection(x,epsilon),K,KS,tightFrame,y);
tic;[xhatBPDNDR,EBPDN] = dr(zeros(p,1),ProxF,ProxG,options);timeBPDNDR=toc

% Solve Dantzig-Selector with PD.
options.mu    = .99*n/p;      % Stepsize on the dual.
options.nu    = .99*n/p;      % Stepsize on the primal.
options.tau   = 1;	      % No relaxation.
options.niter = 300;
options.verb  = 0;
delta = lambda;       % Same parameter as Lasso.      
ProxG = @(x,mu)SoftThresh(x,mu);
z = KS(y);
ProxF = @(x,nu)z+compute_linfball_projection(x-z,delta);
ProxFS= prox_conjugate(ProxF);
K  = @(x)FastA1(FastS1(x,n,dict,pars1,pars2,pars3),dict,pars1,pars2,pars3);
KS = K;
GradH = @(x)0;
options.report = @(x,u)norm(x,1)+sum(z.*u)+delta*norm(u,1);
tic;[xhatDSPD,EDS] = primaldual(zeros(p,1), K,  KS, ProxFS, ProxG, GradH, options);timeDSPD=toc


% Compute recovered signals and respective morphological components.
NumberOfDicts = LengthList(dict);
ind = 0;
for i = 1:NumberOfDicts,
	NAME = NthList(dict, i);
	PAR1 = NthList(pars1, i);
	PAR2 = NthList(pars2, i);
	PAR3 = NthList(pars3, i);
	DimDomain = SizeOfDict1(n, NAME, PAR1, PAR2, PAR3);
	C = xhatlassofista((ind+1):(ind+DimDomain));
	yhatlassofista(:,i) = eval(['Fast' NAME, 'Synthesis(C, n, PAR1, PAR2, PAR3)']);
	C = xhatBPDR((ind+1):(ind+DimDomain));
	yhatBPDR(:,i) = eval(['Fast' NAME, 'Synthesis(C, n, PAR1, PAR2, PAR3)']);
	C = xhatBPDNDR((ind+1):(ind+DimDomain));
	yhatBPDNDR(:,i) = eval(['Fast' NAME, 'Synthesis(C, n, PAR1, PAR2, PAR3)']);
	C = xhatDSPD((ind+1):(ind+DimDomain));
	yhatDSPD(:,i) = eval(['Fast' NAME, 'Synthesis(C, n, PAR1, PAR2, PAR3)']);
	ind = ind + DimDomain;
end


fprintf('%s\n','*'*ones(1,90));
fprintf('%40sSummary\n',' ');
fprintf('%s\n','*'*ones(1,90));
fprintf('%-10s%-15s\t%-15s\t%-15s\t%-15s\t%-15s\n',' ','Original','LASSO-FISTA','BP-DR','BPDN-DR','DS-PD');
fprintf('%-10s%-15g\t%-15g\t%-15g\t%-15g\t%-15g\n','||x||_0:',4+length(pos),  ...
							      length(find(abs(xhatlassofista))), ...
							      length(find(abs(xhatBPDR))), ...
							      length(find(abs(xhatBPDNDR))), ...
							      length(find(abs(xhatDSPD))));

fprintf('%-10s%-15g\t%-15g\t%-15g\t%-15g\t%-15g\n','||x||_1:',4*EnergyRatio+norm(y0dir,1), ...
							      norm(xhatlassofista,1), 	  ...
							      norm(xhatBPDR,1), 	  ...
							      norm(xhatBPDNDR,1), 	  ...
							      norm(xhatDSPD,1));

fprintf('%-25s\t%-15g\t%-15g\t%-15g\t%-15g\n','CPU (s):', timelassofista, timeBPDR, timeBPDNDR, timeDSPD);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(1);clf
% Plot recovered signals.
subplot(5,1,1); 
plot([y0 y]);axis tight
legend('Noiseless',sprintf('Noisy SNR_{in}=%g dB',SNR));legend boxoff

subplot(5,3,4); 
plot([sum(yhatlassofista,2) y0]);axis tight;
title(sprintf('LASSO-FISTA PSNR=%g dB',psnr(y0,sum(yhatlassofista,2))));

subplot(5,3,7); 
plot([sum(yhatBPDR,2) y0]);axis tight;
title(sprintf('BP-DR PSNR=%g dB',psnr(y0,sum(yhatBPDR,2))));

subplot(5,3,10); 
plot([sum(yhatBPDNDR,2) y0]);axis tight;
title(sprintf('BPDN-DR PSNR=%g dB',psnr(y0,sum(yhatBPDNDR,2))));

subplot(5,3,13); 
plot([sum(yhatDSPD,2) y0]);axis tight;
title(sprintf('DS-PD PSNR=%g dB',psnr(y0,sum(yhatDSPD,2))));

% Plot DCT recovered component.
subplot(5,3,5); 
plot([yhatlassofista(:,1) y0dct]);axis([1 n min(y0dct) max(y0dct)]);
title(sprintf('LASSO-FISTA PSNR=%g dB',psnr(y0dct(:),yhatlassofista(:,1))));

subplot(5,3,8); 
plot([yhatBPDR(:,1) y0dct]);axis([1 n min(y0dct) max(y0dct)]);
title(sprintf('BP-DR PSNR=%g dB',psnr(y0dct(:),yhatBPDR(:,1))));

subplot(5,3,11); 
plot([yhatBPDNDR(:,1) y0dct]);axis([1 n min(y0dct) max(y0dct)]);
title(sprintf('BPDN-DR PSNR=%g dB',psnr(y0dct(:),yhatBPDNDR(:,1))));

subplot(5,3,14); 
plot([yhatDSPD(:,1) y0dct]);axis([1 n min(y0dct) max(y0dct)]);
title(sprintf('DS-PD PSNR=%g dB',psnr(y0dct(:),yhatDSPD(:,1))));

% Plot Dirac recovered component.
subplot(5,3,6); 
plot(1:n,yhatlassofista(:,2),'-b',pos,y0dir(pos),'+r');axis([1 n min(y0dir) max(y0dir)]);
title(sprintf('LASSO-FISTA PSNR=%g dB',psnr(y0dir(:),yhatlassofista(:,2))));

subplot(5,3,9); 
plot(1:n,yhatBPDR(:,2),'-b',pos,y0dir(pos),'+r');axis([1 n min(y0dir) max(y0dir)]);
title(sprintf('BP-DR PSNR=%g dB',psnr(y0dir(:),yhatBPDR(:,2))));

subplot(5,3,12); 
plot(1:n,yhatBPDNDR(:,2),'-b',pos,y0dir(pos),'+r');axis([1 n min(y0dir) max(y0dir)]);
title(sprintf('BPDN-DR PSNR=%g dB',psnr(y0dir(:),yhatBPDNDR(:,2))));

subplot(5,3,15); 
plot(1:n,yhatDSPD(:,2),'-b',pos,y0dir(pos),'+r');axis([1 n min(y0dir) max(y0dir)]);
title(sprintf('DS-PD PSNR=%g dB',psnr(y0dir(:),yhatDSPD(:,2))));


if savfig,
  saveas(gcf,'1D/Datasets/testsSparseApproxSinesDiracsSynthesis1D_sig.fig','fig');
  saveas(gcf,'1D/Datasets/testsSparseApproxSinesDiracsSynthesis1D_sig.eps','epsc');
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(2);clf
% Plot recovered coeffs.
subplot(4,3,1); 
stem(xhatlassofista,'.');axis([1 p min(xhatlassofista) max(xhatlassofista)]);
title(sprintf('Coeffs of LASSO-FISTA'));

subplot(4,3,4); 
stem(xhatBPDR,'.');axis([1 p min(xhatBPDR) max(xhatBPDR)]);
title(sprintf('Coeffs of BP-DR'));

subplot(4,3,7); 
stem(xhatBPDNDR,'.');axis([1 p min(xhatBPDNDR) max(xhatBPDNDR)]);
title(sprintf('Coeffs of BPDN-DR'));

subplot(4,3,10); 
stem(xhatDSPD,'.');axis([1 p min(xhatDSPD) max(xhatDSPD)]);
title(sprintf('Coeffs of DS-PD'));

% Plot DCT recovered coeffs.
subplot(4,3,2); 
stem(freq,log(abs(xhatlassofista(1:n*fineness))+1),'.');
axis1 = axis;
axis([freq(120) freq(140) axis1(3) axis1(4)]);
X1 = [freq1 freq2 freq3 freq4]; X1 = [X1;X1];
Y1 = [axis1(3)*ones(1,4);axis1(4)*ones(1,4)];
hold on
plot(X1, Y1, ':b');
hold off
title(sprintf('DCT coeffs log scale LASSO-FISTA'));

subplot(4,3,5); 
stem(freq,log(abs(xhatBPDR(1:n*fineness))+1),'.');
axis1 = axis;
axis([freq(120) freq(140) axis1(3) axis1(4)]);
X1 = [freq1 freq2 freq3 freq4]; X1 = [X1;X1];
Y1 = [axis1(3)*ones(1,4);axis1(4)*ones(1,4)];
hold on
plot(X1, Y1, ':b');
hold off
title(sprintf('DCT coeffs log scale BP-DR'));

subplot(4,3,8); 
stem(freq,log(abs(xhatBPDNDR(1:n*fineness))+1),'.');
axis1 = axis;
axis([freq(120) freq(140) axis1(3) axis1(4)]);
X1 = [freq1 freq2 freq3 freq4]; X1 = [X1;X1];
Y1 = [axis1(3)*ones(1,4);axis1(4)*ones(1,4)];
hold on
plot(X1, Y1, ':b');
hold off
title(sprintf('DCT coeffs log scale BPDN-DR'));

subplot(4,3,11); 
stem(freq,log(abs(xhatDSPD(1:n*fineness))+1),'.');
axis1 = axis;
axis([freq(120) freq(140) axis1(3) axis1(4)]);
X1 = [freq1 freq2 freq3 freq4]; X1 = [X1;X1];
Y1 = [axis1(3)*ones(1,4);axis1(4)*ones(1,4)];
hold on
plot(X1, Y1, ':b');
hold off
title(sprintf('DCT coeffs log scale DS-PD'));

% Plot Dirac recovered coeffs.
xmin = min(pos)-30;xmax = max(pos)+30;
subplot(4,3,3); 
stem(log(abs(xhatlassofista(end-n+1:end))+1),'.');
hold on;
plot(pos,log(abs(y0dir(pos))+1),'+r');axis tight
hold off
axis1 = axis;
axis([xmin xmax axis1(3) axis1(4)]);
title(sprintf('Dirac coeffs log scale LASSO-FISTA'));

subplot(4,3,6); 
stem(log(abs(xhatBPDR(end-n+1:end))+1),'.');
hold on;
plot(pos,log(abs(y0dir(pos))+1),'+r');axis tight
hold off
axis1 = axis;
axis([xmin xmax axis1(3) axis1(4)]);
title(sprintf('Dirac coeffs log scale BP-DR'));

subplot(4,3,9); 
stem(log(abs(xhatBPDNDR(end-n+1:end))+1),'.');
hold on;
plot(pos,log(abs(y0dir(pos))+1),'+r');axis tight
hold off
axis1 = axis;
axis([xmin xmax axis1(3) axis1(4)]);
title(sprintf('Dirac coeffs log scale BPDN-DR'));

subplot(4,3,12); 
stem(log(abs(xhatDSPD(end-n+1:end))+1),'.');
hold on;
plot(pos,log(abs(y0dir(pos))+1),'+r');axis tight
hold off
axis1 = axis;
axis([xmin xmax axis1(3) axis1(4)]);
title(sprintf('Dirac coeffs log scale DS-PD'));

if savfig,
  saveas(gcf,'1D/Datasets/testsSparseApproxSinesDiracsSynthesis1D_coeffs.fig','fig');
  saveas(gcf,'1D/Datasets/testsSparseApproxSinesDiracsSynthesis1D_coeffs.eps','epsc');
end

tilefigs