% ------------------------------------------------------------------------- %
  %                                                                           %
  %                                                                           %
  %             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 rotation + 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));