%
% Withers modified from Aster to use in class by including pauses
%
% This script demonstrates deconvolution by spectral division with
% water level regularization.
%
% First reset MATLAB and its random number generators to a known state.
%
clear all
close all
rand('state',0);
randn('state',0);
%
% We'll use 100 Hz sampling, and a ten second time period.
%
r=100;
deltat=1/r;
t=0.0:deltat:10.0;
%
% For our input signal, we'll use m(t)=t*exp(-t)
%
m=t.*exp(-t);
%
% Plot the input signal.
%
figure(1);
clf;
plot(t,m,'k');
xlabel('t');
ylabel('m(t)');
%print -deps fig1.eps
disp(' ')
disp('We use m(t) = t*exp(-1) as an input to our linear system');
disp('press enter to continue');
pause;
%
% For our impulse response we'll use g(t)=exp(-5t)*sin(10t)
%
g=exp(-5*t).*sin(10.*t);
%
% Plot the impulse response.
%
figure(2);
clf;
plot(t,g,'k');
xlabel('t');
ylabel('g(t)');
%print -deps fig2.eps
disp(' ')
disp('We use g(t)=exp(-5t)*sin(10t) as the impulse response of our linear system');
disp(' Note that this impulse response doesnt have unit gain.  Although');
disp(' most of the output comes in a single pulse in the first half');
disp(' second or so, the integral of g(t) is about 0.08.');
disp('press enter to continue');
pause;
%
% Note that this impulse response doesn't have unit gain.  Although
% most of the output comes in a single pulse in the first half
% second or so, the integral of g(t) is about 0.08.
% 
integralg=sum(g)*deltat;
%
% Also note that we've extended both of these signals in time so that we
% won't have to do any extra 0 padding.
%
%
% Do the convolution using the fft and produce the data.  Add a bit
% of noise.
%
G=fft(g);
M=fft(m);
D=G.*M;
d=ifft(D);
%
% Scale the data by deltat to get a correct level.
%
d=d*deltat;
%
% Plot the clean data.
%
figure(3);
clf;
plot(t,d,'k');
xlabel('t');
ylabel('d(t)');
%print -deps fig3.eps
disp(' ')
disp('d(t) = g(t)*m(t) where * is convolution');
disp('press enter to continue');
pause;
%
% Add the noise.
%
dn=d+0.0001*randn(size(d));
%
% Plot the noisy data.
%
figure(4);
clf;
plot(t,dn,'k');
xlabel('t');
ylabel('dn(t)');
%print -deps fig4.eps
disp(' ')
disp('d(t) with some added noise');
disp('press enter to continue');
pause;
%
% Now try to do the deconvolution by spectral division.
%
DN=fft(dn);
Minv=DN./(G*deltat);
minv=ifft(Minv);
%
% Plot the inverse solution.
%
figure(5);
clf;
plot(t,minv,'k');
xlabel('t');
ylabel('minv(t)');
%print -deps fig5.eps
disp(' ')
disp('m(t) recovered from noisy d(t) using spectral division');
disp('press enter to continue');
pause;
%
% Argh- that's too noisy.  Try some Water Level regularization.
%
w=0.1;
Ghat=w*G./abs(G);
for k=1:size(Ghat,2);
  if abs(G(k))>w
     Ghat(k)=G(k);
  elseif abs(G(k)) == 0
     Ghat(k)=0;
  end
end
mhat=ifft(DN./(Ghat*deltat));
%
% Plot the inverse solution.
%
figure(6);
clf;
plot(t,mhat,'k');
xlabel('t');
ylabel('mhat(t)');
title('Water level solution, w=0.1');
%print -deps fig6.eps
disp(' ')
disp('m(t) recovered from noisy d(t) using water level regularization with w=0.1');
disp('press enter to continue');
pause;
%
% Try more regularization.
%
w=1.0;
Ghat=w*G./abs(G);
for k=1:size(Ghat,2);
  if abs(G(k))>w
     Ghat(k)=G(k);
  elseif abs(G(k)) == 0
     Ghat(k)=0;
  end
end
mhat=ifft(DN./(Ghat*deltat));
%
% Plot the inverse solution.
%
figure(7);
clf;
plot(t,mhat,'k');
xlabel('t');
ylabel('mhat(t)');
title('Water level solution, w=1.0');
%print -deps fig7.eps
disp(' ')
disp('m(t) recovered from noisy d(t) using water level regularization with w=1.0');
disp('press enter to continue');
pause;
%
% Try even more regularization.
%
w=10.0;
Ghat=w*G./abs(G);
for k=1:size(Ghat,2);
  if abs(G(k))>w
     Ghat(k)=G(k);
  elseif abs(G(k)) == 0
     Ghat(k)=0;
  end
end
mhat=ifft(DN./(Ghat*deltat));
%
% Plot the inverse solution.
%
figure(8);
clf;
plot(t,mhat,'k');
xlabel('t');
ylabel('mhat(t)');
title('Water level solution, w=10.0');
%print -deps fig8.eps
disp(' ')
disp('m(t) recovered from noisy d(t) using water level regularization with w=10.0');
disp('press enter to continue');
pause;
%
% a little too much regularization
w=100.0;
Ghat=w*G./abs(G);
for k=1:size(Ghat,2);
  if abs(G(k))>w
     Ghat(k)=G(k);
  elseif abs(G(k)) == 0
     Ghat(k)=0;
  end
end
mhat=ifft(DN./(Ghat*deltat));
%
% Plot the inverse solution.
%
figure(9);
clf;
plot(t,mhat,t,m);
xlabel('t');
ylabel('mhat(t)');
title('Original compared with water level solution, w=100.0');
hold on
plot(t,m,'k');
%print -deps fig9.eps
disp(' ')
disp('m(t) recovered from noisy d(t) using water level regularization with w=100.0');
disp('origin m(t) in black')
disp('press enter to continue');
pause;
%
% a ridiculous amount of regularization
w=1000.0;
Ghat=w*G./abs(G);
for k=1:size(Ghat,2);
  if abs(G(k))>w
     Ghat(k)=G(k);
  elseif abs(G(k)) == 0
     Ghat(k)=0;
  end
end
mhat=ifft(DN./(Ghat*deltat));
% 
%
% Plot the inverse solution.
%
figure(10);
clf;
plot(t,mhat,t,m);
xlabel('t');
ylabel('mtik3(t)');
title('Original compared with water level solution, w=1000.0');
hold on
plot(t,m,'k');
%print -deps fig10.eps
disp(' ')
disp('m(t) recovered from noisy d(t) using water level regularization with w=1000.0');
disp('origin m(t) in black')
disp('press enter to continue');
pause;

