%
% This script demonstrates deconvolution by spectral division with
% Tikhonov regularization.
%
% First reset MATLAB and its random number generators to a known state.
%
clear
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
%
% 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
%
% 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 date 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
%
% 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
%
% 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
%
% Argh- that's too noisy.  Try some Tikonov regularization.
%
Mtik1=(conj(G).*DN)./((conj(G).*G+0.01)*deltat);
mtik1=ifft(Mtik1);
%
% Plot the inverse solution.
%
figure(6);
clf;
plot(t,mtik1,'k');
xlabel('t');
ylabel('mtik1(t)');
title('Tikhonov solution, \lambda=0.01');
print -deps fig6.eps
%
% Try more regularization.
%
Mtik2=(conj(G).*DN)./((conj(G).*G+0.1)*deltat);
mtik2=ifft(Mtik2);
%
% Plot the inverse solution.
%
figure(7);
clf;
plot(t,mtik2,'k');
xlabel('t');
ylabel('mtik2(t)');
title('Tikhonov solution, \lambda=0.1');
print -deps fig7.eps
%
% Try even more regularization.
%
Mtik3=(conj(G).*DN)./((conj(G).*G+1.0)*deltat);
mtik3=ifft(Mtik3);
%
% Plot the inverse solution.
%
figure(8);
clf;
plot(t,mtik3,'k');
xlabel('t');
ylabel('mtik3(t)');
title('Tikhonov solution, \lambda=1.0');
print -deps fig8.eps
