%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Given: x(t)
% Find inverse of x(t) by spectral decon
%  dt     -input- Temporal sampling interval
%  eps0   -input- Damping parameter
%  fr     -input- Dominant frequency of Ricker
%  np     -input- # of samples for input wavelet
%  x(t)   -input- Ricker wavelet
%  IX(f)  -output- Spectral decon filter
%  ix(t)  -output- Result after decon filter applied to x(t)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all
dt=.002;fr=17;np=1024;T=round(1/20/dt);
npt=np*dt;t=(-npt/2):dt:npt/2;;
x=(1-t .*t * fr^2 *pi^2  ) .*exp(- t.^2 * pi^2 * fr^2 ) ; % x(t) is the input wavelet
nt=length(x);x=x(round(nt/2)-T:nt-1);nt=length(x);
subplot(221);
plot([0:nt-1]*dt,x);title('Input Wavelet');
xlabel('Time (s)');axis([0 nt*dt min(x) max(x) ]);

nt2=round(nt/2);df=1/((nt-1)*dt);
subplot(222);
X=fft(x);                                                 % X(f) is Spectrum of x(t)
plot([0:nt2-1]*df,abs(X(1:nt2)));
xlabel('Frequency (Hz)');
title('Input Wavelet');

eps0=.05;
eps1=max(abs(X.^.2))*eps0;
IX=conj(X)./(conj(X).*X+eps1);                            % Inverse Filter for X(f)
subplot(224);plot([0:nt2-1]*df,abs(IX(1:nt2)));
xlabel('Frequency (Hz)');
title('Deconvolved Spectrum');

subplot(223);
ix=ifft(IX.*X);                                           % Deconvolved Wavelet = Inverse Filter * X(f)
nt=length(ix);nt2=round(nt/2);
ix=[flipud(ix(nt2:nt)) ix(1:nt2-1)];
plot([-nt2:1:nt2-2]*dt,real(ix));
axis([-nt2*dt (nt2-2)*dt min(real(ix)) max(real(ix)) ]);
xlabel('Time (s)');
title('Deconvolved Wavelet');

