 % ------------------------------------------------------------------------- %
 %                                                                           %
 %                                                                           %
 %             Script for Maximum Kurtosis Phase Deconvolution               %
 %        using a data set from the Western Canadian Sedimentary Basin       %
 %                                                                           %
 %                 M.D.Sacchi - for Geoph 426 - UofA                         %
 %                                                                           %
 % From SeismicLab:                                                          %
 %                                                                           %
 %           readsegy.m                                                      %
 %           spiking.m                                                       %
 %           kurtosis_of_traces.m                                            %
 %           phase_correction.m                                              %
 %           delay.m                                                         %
 %           bp_filter.m                                                     %
 %           clip.m                                                          %
 %           seismic.m                                                       %
 %                                                                           %
 % ------------------------------------------------------------------------- %



 clear all; 
 close all;

% ---------------------------------
% Read a small stack from the WCB 
% ---------------------------------

 [d,header] = readsegy('wcb_small_stack.su');

 cdp = [header.cdp];
 dt = header(1).dt/1000/1000;


% Extract a window

 d = taper(d(250:900,:),10,10);

 [nt,nx] = size(d);

% ---------------------------------
% Apply spiking deconvolution
% ---------------------------------

 [f,d1] = spiking(d,25,1.);


% ---------------------------------------------------------------
% Assume that the  output of spiking deconvolution has the
% wrong phase. Try to fix it via a phase rotation correction.
%
% Measure the Kurtosis for tentative phase rotations an pick
% the one that maximizes the Kurtosis.
% ---------------------------------------------------------------


 phase_rot = linspace(0,180,20);

 K = [];

   for j = 1:20;
     temp = phase_correction(d1,phase_rot(j));
     K(j) = kurtosis_of_traces(temp);
   end;

 [Kmax,Jmax] = max(K);

  
% ---------------------------------------------------
% Apply the phase rotation that maximizes the Kurtosis
% ---------------------------------------------------

     d2 = phase_correction(d1, phase_rot(Jmax));

% ------------------------------
% Apply band-pass filtering
% Why?
% ------------------------------


     f1 = 2;
     f2 = 4;
     f3 = 55;
     f4 = 70;

     d2 = bp_filter(d2,dt,f1,f2,f3,f4);

% ---------------------------------------
% Delay output by comparison with input
% ---------------------------------------

     d2 = delay(d,d2,40);



% ----------
% Figures 
% ----------

% Display Kurtosis vs phase rotation

 figure(1);
 plot(phase_rot, K); xlabel('Phase rotation in degrees'); ylabel('Kurtosis');


% Display power spectrum before and after process

 figure(2);
 [P_in,freq] = smooth_spectrum(d,dt,10,'li');
 [P_out,freq] = smooth_spectrum(d2,dt,10,'li');

  plot(freq,P_in,freq,P_out); legend('Orignal', 'After Decon');
  xlabel('Frequency [Hz]');
  ylabel('Power Spectral Density');

 
% Display data, data after decon, data after decon + phase rotatuon + filtering

 figure(3);

 taxis = [0:1:nt-1]*dt;

  subplot(131); imagesc(cdp,taxis,clip(d,50,50));  title('Orignal stack');
  subplot(132); imagesc(cdp,taxis,clip(d1,50,50)); title('Spiking decon');
  subplot(133); imagesc(cdp,taxis,clip(d2,50,50)); title('Spiking Decon + Phase Correction');

  colormap(seismic(3));


