clear;
noise=0.05;
N=211;
M=211;
t=linspace(-5,100,N);
%instrument response is a critically-damped pulse
sigi=10;
for i=1:N-1
  if (t(i)<0)
    g(i) = 0;
  else
    g(i) = t(i)*exp(-t(i)/sigi);
  end
end
gmax=max(g);
g=g/gmax;

for i=2:M
  for j=1:N-1
    tp=t(j)-t(i);
    if (tp > 0)
      G(i-1,j)=0;
    else
      G(i-1,j)=-tp*exp(tp/sigi);
    end
  end
end
deltat=t(2)-t(1);
G=G/gmax*deltat;

[u,s,v]=svd(G);

figure(1)
function void=bookfonts
set(gca,'FontSize',18,'LineWidth',1.0);
colormap(gray)
imagesc(G)
H=colorbar;
set(H,'FontSize',18);
xlabel('j')
ylabel('i')
%print -deps2 gimage.eps

figure(2)
function void=bookfonts
set(gca,'FontSize',18,'LineWidth',1.0);
semilogy(diag(s),'ko')
axis tight
xlabel('i')
ylabel('s_i')
%print -deps2 svalues.eps

%true signal is two pulses of sig deviation
sig=2;
mtrue = (exp(-((t(1:N-1)-8).^2/(sig^2*2)))+0.5*exp(-((t(1:N-1)-25).^2/(sig^2*2))))';
mtrue=mtrue/max(mtrue);
d=G*mtrue;
nn=noise*randn(M-1,1);
dn=G*mtrue+nn;

nkeep=N-1;
up=u(:,1:nkeep);
vp=v(:,1:nkeep);
sp=s(1:nkeep,1:nkeep);
mn=vp*inv(sp)*up'*dn;
mperf=vp*inv(sp)*up'*d;

figure(3)
function void=bookfonts
set(gca,'FontSize',18,'LineWidth',1.0);
plot(t(1:N-1),mtrue,'k')
axis tight
xlabel('Time (s)')
ylabel('Acceleration (m/s)');
%print -deps2 truemodel.eps

figure(4)
function void=bookfonts
set(gca,'FontSize',18,'LineWidth',1.0);
plot(t(1:N-1),g,'k')
axis tight
xlabel('Time (s)')
ylabel('V')
%print -deps2 impresp.eps

figure(5)
function void=bookfonts
set(gca,'FontSize',18,'LineWidth',1.0);
plot(t(1:N-1),d,'k')
xlabel('Time (s)')
ylabel('V')
axis tight
%print -deps2 dsynfig.eps

figure(6)
function void=bookfonts
set(gca,'FontSize',18,'LineWidth',1.0);
plot(t(1:N-1),mperf,'k')
xlabel('Time (s)')
ylabel('Acceleration (m/s)');
%print -deps2 perfectsol.eps

figure(7)
function void=bookfonts
set(gca,'FontSize',18,'LineWidth',1.0);
plot(t(1:N-1),dn,'k')
xlabel('Time (s)')
ylabel('V');
axis tight
%print -deps2 dnoisefig.eps

figure(8)
function void=bookfonts
set(gca,'FontSize',18,'LineWidth',1.0);
plot(t(1:N-1),mn,'k')
xlabel('Time (s)')
axis tight
ylabel('Acceleration (m/s)');
%print -deps2 m2noise.eps

%svd truncation to 26 singluar values
nkeep=26;
up=u(:,1:nkeep);
vp=v(:,1:nkeep);
sp=s(1:nkeep,1:nkeep);
m2=vp*inv(sp)*up'*dn;

figure(9)
function void=bookfonts
set(gca,'FontSize',18,'LineWidth',1.0);
plot(t(1:N-1),m2,'k')
xlabel('Time (s)')
ylabel('Acceleration (m/s)');
axis tight
%print -deps2 srej.eps

figure(10)
function void=bookfonts
set(gca,'FontSize',18,'LineWidth',1.0);
Rm=vp*vp';
colormap(gray)
CLIM=[min(min(Rm)),max(max(Rm))];
imagesc(Rm,CLIM)
xlabel('j')
ylabel('i')
H=colorbar;
set(H,'FontSize',18);
%print -deps2 mres.eps

figure(11)
function void=bookfonts
set(gca,'FontSize',18,'LineWidth',1.0);
mresrow=Rm(80,:);
plot(t(1:N-1),mresrow,'k')
axis tight
xlabel('Time (s)')
ylabel('Element Value')
%print -deps2 mresrow.eps