% ------------------------------------------------------------------------- %
%                                                                           %
%            Elimination of Radon artifacts using the KL transform          %
%                                                                           %
%                  M.D.Sacchi - for Geoph 426 - UofA                        %
%                                                                           %
% I have used a large trade-off parameter in the Parabolic Radon            %
% transform. This is to guarantee stability. The problem is                 %
% that a large trade-off will produce artifacts in the estimated            %
% primaries (residual energy from the multiples is left                     %
% in the panel of primaries). The artifacts are attenuated                  %
% with the KL transform.                                                    %
%                                                                           %
% From SeismicLab: pradon_demultiple.m                                      %
%                  readsegy.m                                                     %
%                  kl.m                                                     %
%                  pimage.m                                                 %
%                  clip.m                                                   %
%                  seismic.m                                                %
%                                                                           %
% ------------------------------------------------------------------------- %



  clear;
  close all;


% Read data, a nmo-ed corrected marine gather
% from the Gulf of Mexico (./SeismicLab_data)

  [D,H] = readsegy('gom_cdp_nmo.su');

% We will use a subset of the data

  D = D(800:1200,:);

% Extract areas that were muted

  Mutes = ones(size(D)); 
  I = find(D==0); 
  Mutes(I) = 0;

% Extract offet and sampling interval

  h = [H.offset];
  dt = H(1).dt/1000/1000;

% Define parameters for Radon Demultiple

  qmin = -0.3;
  qmax = .8;
  nq = 120;
  flow = 2.;
  fhigh = 90;
  mu = 1;              % Trade-off parameter
  q_cut = 0.05;


  [P,M,tau,q] = pradon_demultiple(D,dt,h,qmin,qmax,nq,flow,fhigh,mu,q_cut);

% Apply data mutes back to the estimated  primaries

  P = P.*Mutes;

% KL transform filter using a subspace of 15 eigenvectors

  P_kl = kl(P,15);
  P_kl = P_kl.*Mutes;
 
% Plot data, Radon panel, Primaries and Primaries after
% fitering with the kl transform

 Big =  [D,10*M,P,P_kl,P-P_kl];
 [nt,nx]=size(Big);

 pimage([1:1:nx],tau,clip(Big,50,50));

 colormap(seismic(2));
