%
%
% Optical Imaging and Spectroscopy
%
% David J. Brady
% Duke University
% www.opticalimaging.org
%
% Figure 10.43
%
% SVD for random and hadamard codes
%
figure(1);set(gcf,'color','white');
colormap gray;
hr=round(rand(511));
[ar br cr]=svd(hr);
%
% truncated SVD inversion matrix
%
rr=300;
hrTSVD=cr(:,1:rr)*pinv(br(1:rr,1:rr))*ar(:,1:rr)';
subplot(3,1,1);semilogy(diag(br),'-k');hold on;
%
% hadamard sampling
%
h=(1-hadamard(512))/2;
hMat=h(2:512,2:512);
[a b c2]=svd(hMat);
subplot(3,1,1);semilogy(diag(b),'-k');
title('(a)Singular value spectrum');
axis([1 511 0.1 100]);
%
%
%
% compare Random sampling and Hadamard noise 
%
% Definition of spectral signature
%
%
lambda=linspace(400,700,511);
spec=exp(-(lambda-450).^2/10^2);
spec=spec+.7*exp(-(lambda-480).^2/5^2);
spec=spec+.5*exp(-(lambda-510).^2/7^2);
spec=spec+.25*exp(-(lambda-550).^2/3^2);
spec=spec+.1*exp(-(lambda-600).^2/5^2);
spec=spec+.8*exp(-(lambda-650).^2/1^2);
spec1=hrTSVD*ar(:,1:rr)*ar(:,1:rr)'*poissrnd(hr*(1000*spec'))/1000;
spec1=nullSpaceSmoothing(null(hrTSVD'),spec1);spec1=spec1';
spec2=nnls(hMat,poissrnd(hMat*(1000*spec')))/1000;spec2=spec2';
lev = 5;
spec1Den = wden(spec1,'minimaxi','s','sln',lev,'sym8');
spec2Den = wden(spec2,'minimaxi','s','sln',lev,'sym8');
subplot(3,1,2);plot(lambda,[spec;spec1+.5;spec2+1],'-k');axis([400,700, 0,2.2]); 
title(['(b) \sigma_{random}^2=' num2str((spec1-spec)*(spec1-spec)',2) '     \sigma_{Hadamard}^2=' num2str((spec2-spec)*(spec2-spec)',2)]);
xlabel('\lambda (nm)');ylabel('Spectral density');
set(gca,'Ytick',[]);
subplot(3,1,3);plot(lambda,[spec1Den+.1;spec2Den+.6],'-k');axis([400,700, 0,1.8]); 
title(['(c) \sigma_{random}^2=' num2str((spec1Den-spec)*(spec1Den-spec)',2) '     \sigma_{Hadamard}^2=' num2str((spec2Den-spec)*(spec2Den-spec)',2)]);
xlabel('\lambda (nm)');ylabel('Spectral density');
set(gca,'Ytick',[]);