%
% -----------------------------------------------------------------
% Estimation of a Sparse Reflectivity Series
% 
% Mauricio D. Sacchi, SAIG, 2012
%
% Please refer this article when using the code: 
%       
% M.D. Sacchi, M.D., 1997, Re-weighting strategies in seismic 
% deconvolution: Geophysical Journal International, 129, 651-656.
%


% -------------------- Start deconvolution ------------------------

% Load data and wavelet (Marine trace, Gulf of Mexico)

 load w.mat;
 load s.mat;

 tol      = 1e-4;
 Max_iter = 50;
 lambda   = 0.4 
 silent   = 0; 

 [r,sp,J] = l1_l2_irls(s, w, lambda, Max_iter, tol, silent);


% -------------------- End deconvolution --------------------------

% The rest is for plotting 

 r = max(s)*r/max(r); 
 taxis = [0:1:length(r)-1]*4/1000;

 figure(1); 

 subplot(211); 

  plot(taxis, s, taxis, r+1 , taxis, sp+2);  
  legend('Seismogram', 'Inverted Reflectivity', 'Predicted Seismogram'); 
  xlabel('Time (s)'); 
  title('IRLS \lambda = 0.4');

 subplot(223); 

  plot(J(:,1)); 
  xlabel('Iteration'); 
  ylabel('Cost function'); 

 subplot(224); 

  plot(taxis(1:length(w)),w/max(w));
  xlabel('Time (s)'); 
  title('Wavelet'); axis tight;