clear;

echo on
%Convolution principles demo featuring smooth functions.

%Adapted from Example 8.2 of Aster, Borchers, and Thurber, "Parameter Estimation and Inverse Problems" (Elsevier, 2004).

%Learning Goals

% Show basic syntax for those new to Matlab 

% Demonstrate Frequency-domain deconvolution

% Demonstrate Time-domain devonvolution

% Demonstrate deconvolution instability in the presence of noise

% Illustrate fundamental smoothness vs. data fit tradeoffs with both freqeuncy-and time-domain regularization.

%set noise level
noise=0.05;

%time series length parameter
N=211;

%sampling interval (s)
dt=0.5;

echo off

%frequency spacing of Fourier transform terms (Hz)
dF=1/dt;

%time vector for plots
t=linspace(-5,100,N);

%frequency vector for plots
F=linspace(0,dF/2,N-1);

%Convolution kernel option 1
%synthetic instrument response is a critically-damped pulse (e.g,, a critically-damped seismometer)
%tau is the time constant
tau=10;

%generate system response here (actual time series are N-1 samples long)
for i=1:N-1
if (t(i)<0)
g(i) = 0;
else
g(i) = t(i)*exp(-t(i)/tau);
end
end

gmax=max(g);

%final response is a unit-height column vector
g=g'/gmax;

%synthetic signal to be estimated via deconvolution consists of two pulses of standard deviation sig
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);

%Kernel option 2: uncomment to produce a two-pulse convolution kernel
sig=4;
g = (exp(-((t(1:N-1)-8).^2/(sig^2*2)))+0.5*exp(-((t(1:N-1)-25).^2/(sig^2*2))))';
g=g/max(g);

%Kernel option 3: uncomment to produce a three-pulse convolution kernel
sig=1;
g = (exp(-((t(1:N-1)-8).^2/(sig^2*2)))+0.5*exp(-((t(1:N-1)-20).^2/(sig^2*2)))+0.25*exp(-((t(1:N-1)-30).^2/(sig^2*2))))';
g=g/max(g);

echo on
%perform a time domain (serial) convolution to get the data (observed) time series vector
echo off

%note that the length of d will be equal to 2*(N-1)-1
d=conv(g,mtrue);

%generate Gaussian noise of amplitude noise
nn=noise*randn(length(d),1);

%generate the noisy data vector
dn=d+nn;

echo on
%first, we will plot the true deconvolution solution
%
echo off
figure(1)
clf
plot(t(1:N-1),mtrue,'k')
set(gca,'FontSize',18,'LineWidth',1.0);
xlabel('Time (s)')
axis tight
title('True Model')
pause

echo on
%Here is a plot of the convolution kernal that smooths the model to produce
%the data vector
echo off
figure(2)
clf
plot(t(1:N-1),g,'k')
set(gca,'FontSize',18,'LineWidth',1.0);
xlabel('Time (s)')
axis tight
title('System Impulse Response')
pause

echo on
%and here is the result of convolving the previous two plotted functions to
%produce a noiseless data vector 
echo off
figure(3)
clf
plot(t(1:N-1),d(1:N-1),'k')
set(gca,'FontSize',18,'LineWidth',1.0);
xlabel('Time (s)')
axis tight
title('Noise-free Data')
pause

echo on
%Now we generate spectra for the data and impulse response using the fft
%pad all time series to 2*(N-1) points to avoid wrap-around effects from circular convolution 
%note that in a large problem we would want to use highly composite (ideally powers of 2) length time series
echo off

%FFT of the padded noise-free data
dspec=fft([d; 0]);
%FFT of the padded noisy data
dnspec=fft([dn; 0]);
gspec=fft([g; zeros(length(g),1)]);
mperf=real(ifft(dspec./gspec));
mn=real(ifft(dnspec./gspec));

echo on
%The recovered model for the noiseless data vector using spectral
%division.  The result should be almost exact thanks to Matlab's
%double-precision numerical accuracy.
%
echo off
figure(4)
clf
plot(t(1:N-1),mperf(1:N-1),'k')
set(gca,'FontSize',18,'LineWidth',1.0);
xlabel('Time (s)')
axis tight
title('Noise-Free Data Frequency-Domain Deconvolution')
pause

echo on
%We add a bit of noise (as specified in the noise level set at the
%beginning of the program)...
%
echo off
figure(5)
clf
plot(t(1:N-1),dn(1:N-1),'k')
set(gca,'FontSize',18,'LineWidth',1.0);
xlabel('Time (s)')
axis tight
title(['Noisy Data (',num2str(100*noise),'%)'])
pause

%This is the deconvolution result using the noisy data.  If the data noise
%is large enough, the result will bear no resemblence to the true model
%because of the ill-posed nature of the deconvolution problem.
%Essentially, frequencies that were highly attenuated in the convolution
%will be greatly amplified in the deconvolution.  If unmodeled signal
%('noise' in this context) contains appreciable energy at these
%frequencies, the influence of noise will dominate the deconvolution
%result.
%
echo off
figure(6)
clf
plot(t(1:N-1),mn(1:N-1),'k')
set(gca,'FontSize',18,'LineWidth',1.0);
xlabel('Time (s)')
axis tight
title('Noisy Data Frequency Domain Deconvolution')
pause

%generate spectral parameters for noise investigation
gs=abs(gspec(1:length(gspec)/2));
ds=abs(dspec(1:length(dspec)/2));
dns=abs(dnspec(1:length(dspec)/2));

echo on
%Examining the convolving kernel (impulse response), and the noise free and
%noisy data spectra shows the nature of the deconvolution instability more
%clearly.  The noisy data has lots of energy in bands that are greatly
%attenuated by the convolution with the impulse response (at high
%frequencies in this example).  Note the log scale for spectral amplitude.
%
echo off
figure(7)
clf
loglog(F,gs,'k');
xlabel('f (Hz)')
ylabel('Spectral Amplitude')
ylim([1e-3 1e2]);
xlim([1/N 1]);
hold on
loglog(F,ds,'k-.')
loglog(F,dns,'k-')
hold off
legend('Impulse Response','Noise-Free Data','Noisy data')
title('Frequency-Domain Comparison')
pause

echo on
%when we divide the data in the frequency domain by the impulse response
%spectrum, we get amplification of the high-frequency part of the data
%spectrum.  In the case of the noisy data, this component will dominate the
%result in this example if the noise level is even a few percent of the signal.
%
echo off
figure(8)
clf
loglog(F,ds./gs,'k')
hold on
loglog(F,dns./gs,'k-')
xlabel('f (Hz)')
ylabel('Spectral Amplitude')
axis tight
title('Deconvolution Spectra')
legend('Noise-Free Data','Noisy data','Location','SouthWest')
pause

echo on
%A straightforward way to improve the deconvolution result in this simple example would be to low-pass
%filter the result.  A more adaptive way to stabilize a spectral divisio deconvolution (that will work for
%peaks in the deconvolution spectrum as well as high-frequency
%amplification effects such as observed here), is to "regularize" the
%deconvolution by filling in the especially low amplitude features of the convolution
%spectrum before spectral division.  This will prevent high amplifications
%of the data spectrum in the deconvolution at those frequencies.  There
%will be a tradeoff, as in all regularization methodologies, between
%stabilizing the result and degrading the accuracy of the deconvolution by
%over-regularization.  In the absurd limit as we "flood" the entire spectrum with
%a very high water level, the deconvolution will equal the data vector divided
%by the amplitude of the water level!
%We will examine this tradeoff by regularizing with a range of water
%levels and looking at the tradeoff curve between model norm and the residual norm
%(how well the deconvolution result fits the data when plugged back into the forward
%problem.
%
%Here is an animation of the evolving solution as the water level is
%increased...
%
echo off

wcount=0;
%select an exponentially-spaced range of water levels 
for wlevexp=-1:.1:2
  wcount=wcount+1;
  wlev=10^wlevexp;
%
%apply water level damping
%by constructing successive water level patched spectra, gwspec (wlev increasing)
 for i=1:length(gspec)
    if abs(gspec(i)) < wlev
      gwspec(i,1)=wlev*(gspec(i)./abs(gspec(i)));
    else
      gwspec(i,1)=gspec(i);
    end
  end

%recovered signal for this particular water level
  mw=real(ifft(dnspec./gwspec));
  wval(wcount)=wlev;
  mnorm(wcount)=norm(mw);
  resid(wcount)=norm((dnspec./gwspec).*gspec-dnspec);
  mstore(:,wcount)=mw;

  figure(9)
  clf
  plot(t(1:N-1),mw(1:N-1),'k')
set(gca,'FontSize',18,'LineWidth',1.0);
  xlabel('Time (s)')
legend(['Water level: ',num2str(wlev)]);
%plot successive water level deconvolution results
  drawnow;
title('Frequency Domain Deconvoution Result')
  pause(0.2);
end
pause

echo on
%and here is the tradeoff curve for the water level regularized
%deconvolution.  Generally, a "best" solution is selected based on the
%corner of the tradeoff curve.
%
echo off

figure(10)
clf
plot(resid,mnorm,'ko')
set(gca,'FontSize',18,'LineWidth',1.0);
xlabel('Misfit 2-Norm');
ylabel('Model 2-Norm');
for i=1:10:wcount
  H=text(resid(i),mnorm(i)+.2,['  w=',num2str(wval(i))]);
  set(H,'Fontsize',18);
end
title('Model Norm vs. Misfit Tradeoff Curve, Water Level Regularization')
pause

echo on
%here we display all of the solutions in another manner
%
echo off
figure(11)
clf
for i=1:31
hold on
plot(t(1:N-1),mstore(1:N-1,i)+10*log10(wval(i)),'k')
set(gca,'FontSize',18,'LineWidth',1.0);
plot(t(1:N-1),mtrue+10*log10(wval(i)),'-.r')
set(gca,'FontSize',18,'LineWidth',1.0);
end
hold off
xlabel('Time (s)')
ylabel('10 log_{10} (w)')
axis tight
legend('Frequency-Domain Deconvolution','True Model')
pause

echo on
%finally, we show the "preferred" solution, using the water level value
%that lies at the approximate corner of the model norm vs. residual curve
%
echo off

figure(12)
clf
%preferred model from tradeoff corner
mpref=16;
plot(t(1:N-1),mstore(1:N-1,mpref),'k')
set(gca,'FontSize',18,'LineWidth',1.0);
hold on
plot(t(1:N-1),mtrue(1:N-1),'-.r')
set(gca,'FontSize',18,'LineWidth',1.0);
hold off
xlabel('Time (s)')
axis tight
legend(['Frequency Domain Deconvolution for w=',num2str(wval(mpref))],'True Model');
pause

echo on
%Next, we will do this deconvolution in the time domain.
%
%The G matrix for the forward problem, G*m=d, has columns that are simply time-shifted Green's
%functions
%
echo off

%preallocate G
M=2*(N-1);
G=zeros(M,N);
gpad=[g;zeros(size(g))];
for i=1:N
G(:,i)=circshift(gpad,[i-1,0]);
end

echo on
%now find the least-squares solution using time domin deconvolution (and
%keep track of the time to execute)
tic;
%least-squares deconvolution solution for noise-free data (using the Matlab
%slash function
mtd=G\[d;0];
runtime1=toc;
disp(['Time domain deconvolution run time is: ',num2str(runtime1),' s'])
tic
%reevaluate the spectral deconvolution for runtime comparison
mfd=real(ifft(dspec./gspec));
runtime2=toc;
echo off

disp(['Frequency domain deconvolution run time is: ',num2str(runtime2),' s (',num2str((1-runtime2/runtime1)*100),'% faster)'])
tcpu0=cputime;
figure(13)
plot(t,mtd,'k')
set(gca,'FontSize',18,'LineWidth',1.0);
axis tight
xlabel('Time (s)')
title('Noise Free Data Time-Domain Deconvolution')

%least-squares deconvolution solution for noisy data
mtdn=G\[dn;0];

echo on
%here's the time-domain deconvolution for the noisy data with no
%regularization.  No surprise, it's the same outrageous solution that we
%obtained in the frequency domain (and it took a lot more crunching to get
%it).
%
echo off
pause

figure(14)
plot(t,mtdn,'k')
set(gca,'FontSize',18,'LineWidth',1.0);
axis tight
xlabel('Time (s)')
title('Noisy Data Time-Domain Deconvoluion');
pause

echo on
%A desirable feature of time-domain deconvolution is that we can regularize
%the deconvolution in various ways that are easily interpred in the time
%domain.  Here we will demonstrate a deconvolution solution for noisy data
%that minimizes a discrete approximation to the integral of the second derivative (second-order
%Tikhonov regularization).  The basic methodology is to augment the system
%matrix (forward problem matrix) with additional constraints that force the
%integral of the second derivative to be small).  The balance between
%making this "roughness" small and fitting the data is set by a
%regularization parameter that is analogous to the water level parameter used earlier.
%
%Here is an animation of the evolving second-order Tikhonov solution as the regularization parameter, alpha, is
%increased...
%
%
echo off

%preallocate and generate the second-order roughening matrix
L=zeros(N-2,N);
lvec=zeros(1,N);
lvec(1:3)=[1 -2 1];
for i=1:N-2
L(i,:)=circshift(lvec,[0,i-1]);
end

alphacount=0;
%select an exponentially-spaced range of tradeoff parameters
for alphaexp=-1:.1:2
  alpha=10^alphaexp;
  alphacount=alphacount+1;
%
malpha=[G;alpha*L]\[dn;0;zeros(N-2,1)];
  alphaval(alphacount)=alpha;
  mnormtd(alphacount)=norm(malpha);
  dpred=G*malpha;
  residtd(alphacount)=norm(dpred(1:M-1)-dn);
  mstoretd(:,alphacount)=malpha;

figure(15)
clf
plot(t,mstoretd(:,alphacount),'k')
set(gca,'FontSize',18,'LineWidth',1.0);
axis tight
xlabel('Time (s)')
title('Time Domain Deconvolution Result');
legend(['2^{nd} Order Regularization Parameter: ',num2str(alpha)]);
pause(0.2)
end
pause

echo on
%Here is the tradeoff curve for the time-domain deconvolution, this time
%with the parameter alpha
%
echo off

figure(16)
clf
plot(residtd,mnormtd,'ko')
set(gca,'FontSize',18,'LineWidth',1.0);
xlabel('Misfit 2-Norm');
ylabel('Model 2-Norm');
for i=1:10:alphacount
  H=text(residtd(i),mnormtd(i)+.05,['\alpha=',num2str(alphaval(i))]);
  set(H,'Fontsize',18);
end
title('Model Norm vs. Misfit Tradeoff Curve (Time Domain Deconvolution, 2^{nd} Order Regularization)')
pause

echo on
%A display of time-domain regularized deconvolution results.  
%
echo off

figure(17)
clf
for i=1:31
hold on
plot(t(1:N-1),mstoretd(1:N-1,i)+10*log10(alphaval(i)),'k')
set(gca,'FontSize',18,'LineWidth',1.0);
plot(t(1:N-1),mtrue+10*log10(alphaval(i)),'-.r')
set(gca,'FontSize',18,'LineWidth',1.0);
set(gca,'FontSize',18,'LineWidth',1.0);
end
hold off
xlabel('Time (s)')
ylabel('10 log_{10} (\alpha)')
axis tight
legend('Time Domain Deconvolution','True Model')
pause

echo on
%And a plot of the preferred model from the tradeoff corner in the previous
%figure.  Note that there are higher-order water level methods as well that
%can approximate higher-order Tikonov solutions to combine the speed of the
%frequency domain deconvolution with the desireable features of Tikonov
%regularization (Aster et al., ex. 8.4).
%
%In this particular case, the Tikonov result especially good because the true model
%actualy consists of two fairly broad (relatively small second derivative)
%pulses and the unmodeled noise was white broadband noise.
%
%Moral: A priori constraints (or prejudices!) about the model are inherent in the
%choice of a specific regularization methodology in a decolvolution or other ill-posed
%inverse problem!
%
echo off

figure(18)
clf
mpref=18;
plot(t(1:N-1),mstoretd(1:N-1,mpref),'k')
set(gca,'FontSize',18,'LineWidth',1.0);
hold on
plot(t(1:N-1),mtrue(1:N-1),'-.r')
set(gca,'FontSize',18,'LineWidth',1.0);
hold off
xlabel('Time (s)')
axis tight
legend(['Time Domain Deconvolution for \alpha=',num2str(alphaval(mpref))],'True Model');
