%
% Optical Imaging and Spectroscopy
%
% David J. Brady
% Duke University
% www.opticalimaging.org
%
% Figure 8.16
%
%
% This code is derived from the file l1eq_example.m
%
%(l1 minimization with equality constraints).
%
% The code calls functions in the l1-magic code on line at
% http://www.acm.caltech.edu/l1magic/
%
% Written by: Justin Romberg, Caltech
% Email: jrom@acm.caltech.edu
% Created: October 2005
%

% put key subdirectories in path if not already there
path(path, '.\Optimization');
%path(path, '.\Data');

% load random states for repeatable experiments
% load RandomStates
% rand('state', rand_state);
% randn('state', randn_state);
hold off;
% signal length
N = 765;
% number of spikes in the signal
T = 4;
% number of observations to make
K = 130;
load Xenonspectrum

% random +/- 1 signal
%x = threshold(spectrum,.05)';
x = spectrum';
% q = randperm(N);
% x(q(1:T)) = sign(randn(T,1));

% measurement matrix
disp('Creating measurment matrix...');
A = randn(K,N);
A = orth(A')';
disp('Done.');
	
% observations
y = A*x;

% initial guess = min energy
x0 = A'*y;

% solve the LP
tic
xp = l1eq_pd(x0, A, [], y, 1e-3);
toc

figure(1); set(gcf,'color','white');

subplot(3,1,1);plot(wavelengths, x,'-k'); axis([860 930 -.05 1]); title('(a)'); 
subplot(3,1,2);plot(wavelengths, xp,'-k');axis([860 930 -.05 1]);  title('(b)');
subplot(3,1,3);plot(wavelengths, 10*circshift(xp,-10)+2,'-k');axis([860 930 0 7]);title('(c)'); xlabel('\lambda (nm)');
hold on
subplot(3,1,3);plot(wavelengths, 10*x+1,'-k');  

% spectrum without noise
x = threshold(spectrum,.05)';
%x = spectrum';
% observations
y = A*x;

% initial guess = min energy
x0 = A'*y;

% solve the LP
tic
xp = l1eq_pd(x0, A, [], y, 1e-3);
toc
subplot(3,1,3);plot(wavelengths, 10*circshift(xp,-20)+3,'-k');
%
% noise but more samples
%
x = spectrum';
K=200;
A = randn(K,N);
A = orth(A')';
% observations
y = A*x;

% initial guess = min energy
x0 = A'*y;

% solve the LP
tic
xp = l1eq_pd(x0, A, [], y, 1e-3);
toc
subplot(3,1,3);plot(wavelengths, 10*circshift(xp,-30)+4,'-k');
%
%
%
% noise but fewer samples
%
x = spectrum';
K=100;
A = randn(K,N);
A = orth(A')';
% observations
y = A*x;

% initial guess = min energy
x0 = A'*y;

% solve the LP
tic
xp = l1eq_pd(x0, A, [], y, 1e-3);
toc
subplot(3,1,3);plot(wavelengths, 10*circshift(xp,-40)+5,'-k');
%
% noise but still fewer samples
%
x = spectrum';
K=90;
A = randn(K,N);
A = orth(A')';
% observations
y = A*x;

% initial guess = min energy
x0 = A'*y;

% solve the LP
tic
xp = l1eq_pd(x0, A, [], y, 1e-3);
toc
subplot(3,1,3);plot(wavelengths, 10*circshift(xp,-50)+6,'-k');