clear all
flops(0)
tic % remember the starting time
echo on;
% This demo illustrates some basic principles of the newly developed
% multiple removal methods described in a recent paper by Fan et al. (2006).
% This approach has potential for major improvements in our ability to
% estimate the impulse response of the Earth, particularly for P 
% waves reflected from the free surface and returned as an equivalent
% to the plane-wave reflection response.  That is, it's greatest 
% potential is to go beyond P to S conversions and also estimate 
% P wave reflection sections from teleseismic data.
%
% This demo is intended as a self-taught lesson.  The idea is to 
% run through it a few times to get the ideas.  For further study, 
% take the m-file home and study it.  We hope some of you can pick 
% up from here and continue this line of work to improve this 
% technology.  
%
% Author:  Chengliang Fan
% Adapted by G Pavlis as a demo for this workshop.
%
% Hit return to start the demo
%
pause;

h=[-100,-40,-4,0];
vp=[8.0,8.0,6.0,4.0];
vs=[4.7,4.7,3.5,2.3];
stairs(vp,h);
axis([0,10,-100,0]);
hold
stairs(vs,h,'r');
xlabel('Velocity km/s)')
ylabel('Depth (km)')
% 
%  First, let's look at the simple 1D earth model used for 
% this simulation.  This is a simple crustal model with a 4 km
% thick sedimentary section.  The simulations in this tutorial 
% use and incident P wave at vertical incidence.  In this case,
% the S model is largely irrelevant (students, why?).  This simple
% model has 2 reflection primaries (bottom of sediments and base
% of the crust) and a typical sequence of free-surface and internal
% multiples.  In this simulation we will see transmission and 
% reflection geometry data derived from this model and how the
% inverse scattering series can be used to remove free surface
% multiples from both the reflection and transmission geometry
% impulse response function.
%
pause;
echo off;

%%   Load transmission responses at different ray parameter
trans2n = zeros(3277,5); 
trans2n_0_temp = load('2nc000ZEX.sac'); % p: 0 s/km
trans2n(:,1)  = trans2n_0_temp(:,3);
% trans2n_1_temp = load('D:\syn_t_r\ps_trans\trans2n\2nc003ZEX.sac'); % p: 0.03 s/km
% trans2n(:,2)  = trans2n_1_temp(:,3);
% trans2n_2_temp = load('D:\syn_t_r\ps_trans\trans2n\2nc006ZEX.sac'); 
% trans2n(:,3)  = trans2n_2_temp(:,3);
% trans2n_3_temp = load('D:\syn_t_r\ps_trans\trans2n\2nc009ZEX.sac');
% trans2n(:,4)  = trans2n_3_temp(:,3);
% trans2n_4_temp = load('D:\syn_t_r\ps_trans\trans2n\2nc012ZEX.sac');
% trans2n(:,5)  = trans2n_4_temp(:,3);

%%  extracted from real data  %%
cmpc = load('CMPC.BHZ.sac'); % total length of signal: 180 seconds
 rd(:,4) = cmpc(:,3);
 rd4 = rd(:,4)/2200;
 rds4 = rd4(280:935); % 656 points

% Resample the source signature and make the sampling points equal 3277
rds4_resample_temp = resample(rds4,5,1);

% Move the signal down to around y = 0; the mean of the signal is 0.4965
rds4_resample = rds4_resample_temp(1:3277) - 0.5;

figure
plot(rds4_resample'); axis tight
title('source signature')
echo on
% This is a plot of source wavelet we will use for generating a synthetic
% data set.  It is a real P wave from the Bolivar experiment in Venezuela.
pause;

% define source signature
ss = rds4_resample;

% Select transmission response generated at different ray parameter
for j = 1:1
% Define the initial transmission response in time domain
T_n =  trans2n(:,j);

% Generate the time series for plotting
echo off;
for i = 1:3277
    x(i) = i*0.005*5;
end
echo on;
% generating reflection response from input (transmission) response
echo off;

% Initialize and call function ricker and gaussian for smoothing responses
n = 3277;
dt = 0.15;
nu = 0.3;
gau = gaussian(n,dt,nu);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%  There are three ways to generate the reflection response from transmission response %%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%% Calculate the real part of reflection response in frequency domain
% We have the following relation in frequency domain in theory:
% 2Re{R(w)}=1 - (T_norm(w))'*T_norm(w)
% Fourier transform T_norm into frequency domain
T_n_w = fft(T_n);
Re_R_w_2 = 1 - (abs(T_n_w)).^2;

%%% First: Reconstruct the reflection response by Hilbert transform
R_w = hilbert(0.5*Re_R_w_2);
R_wt = ifft(R_w);
R_t(1) = R_wt(1);
for i = 2:3277
        R_t(i) = R_wt(3279-i);
end
R_t(1:3) = 0;
R_t(3274:3277) = 0;

%%% Second: Inverse Fourier transform Re_R_w_2 into time domain
% This will reconstruct the reflection response in time domain
Re_R_2 = ifft(Re_R_w_2);
Re_R_2(1:4) = 0;
Re_R_2(1600:3277) = 0;

%%% Third: Reconstruct the reflection response by auto-corr in time domain
aucorr = xcorr(T_n);
for i=1:6553
   ReR2(3274:3280) = 1 - aucorr(3274:3280);
   ReR2(i) = -aucorr(i);
end
R(1:4) = 0;
R(5:3277) = ReR2(3281:6553);

figure
subplot(2,2,1)
plot(x,85.5*4.41*T_n); axis tight
title('Transmission response (TR)')
xlabel('travel time(second)')
subplot(2,2,2)
plot(x,85.5*4.41*real(R_t)); axis tight
title('RR constructed by Hilbert Transform')
xlabel('travel time(second)')
subplot(2,2,3)
plot(x,85.5*4.41*real(Re_R_2)); axis tight
title('RR constructed by t domain real part')
xlabel('travel time(second)')
subplot(2,2,4)
plot(x,85.5*4.41*R); axis tight
title('RR constructed by auto-correlation')
xlabel('travel time(second)')
echo on;
%  These plots illustrate how the backscattered (reflection) response
%  can be constructed from the input (Transmission response) data for 
%  radially symmetric Earth models.  The upper left figure shows the
%  transmission data used to construct the reflection sequences in
%  the other 3 frames of the figure.  The three different reflection
%  responses are indistinguishable.  They are the result of the application
%  of three different algorithms that can be used to construct the
%  reflection response from  the transmission data: 
%   1. (upper right) A method using a Hilbert transform
%   2. (lower left) A method using the real part of the complex FFT
%   3. (lower right) A method using the autocorrelation function.
%
%  These methods are founded on relatively standard signal processing
%  algorithms exploiting well known properties of the Fourier Transform. 
%  At your leisure you should review how these are computed.  Note Wapenaar has 
%  derived the (much more complicated) forms for the more realistic
%  situation with fully 3D structure.  How well we can do with radially
%  symmetric models, however, is currently an open question.  
pause;

% Scale the TR and constructed RR properly for free surface multiple removal (FSMR)
Rs = 40.5*9.8*R;
T_n = 9.8*T_n;

%%%%%%% FSMR for reflection response %%%%%%%

y = Rs;  
yy = 0;
figure
% Entering movie loop to show how series progressively removes
% multiples from the data.  Basic algorithm is taking the data
% to successive powers in the frequency domain.  See Fan et al.(2005)
% for details.  
%
%  Turning off echo to improve performance.
%  Watch the sequence as it progressively removes the multiples from
%  the reflection response.  This will show the application of 
%  the first 10 terms in the series.
%
echo off;
 y_c = 0;
% initialize and calculate the first term in the inverse scattering series
 for i = 1:1
      y_c = conv(y,y_c);
       y_cc =y_c(1:3277);
       yy = yy + y_cc;
       y_n = y+yy;
     tri_ffg = conv(y_n,gau);
     pq = tri_ffg(1639:4915);
     h = plot(x,pq,'k');
     xlabel('travel time(second)')
     title('reflection response free surface back-scattered wave removal')
      set(h,'XData',x, 'YData', pq);
       axis([0.025 82 -0.05 0.1])
    Ref(i) = getframe;
    pause(1)
end
 y_c = y;
for i = 2:10 % totally 10 terms in the inverse scattering series are used
      y_c = conv(y,y_c);
       y_cc =y_c(1:3277);
       yy = yy + y_cc;
       y_n = y+yy;
     tri_ffg = conv(y_n,gau);
     pq = tri_ffg(1639:4915);
     h = plot(x, pq, 'k');
     xlabel('travel time(second)')
     title('reflection response free surface back-scattered wave removal')
      set(h,'XData',x, 'YData',pq);
      axis([0.025 82 -0.05 0.1])
    Ref(i) = getframe;
    pause(1)
end

% write out the movie file
% movie2avi(Ref,'D:\syn_t_r\results\mr4\rrfsmr_trans2n000.avi' ,'fps',1,'quality', 100)

echo on;
pause;
echo off;

% constructed reflection response; smooth to the same pulse width as the RR after FSMR
Rs_conv_temp = conv(Rs, gau);
Rs_conv = Rs_conv_temp(1639:4915);

% Calculate the energy of RR before and after FSMR 
ERs_conv = sum((abs(fft(Rs_conv))).^2) % Total amplitudes of constructed reflection response before FSMR
Epq = sum((abs(fft(pq))).^2) % Total amplitudes of constructed reflection response after FSMR

% Reflection data with free surface multiples
rdwm_temp = conv(ss, Rs_conv);
rdwm = rdwm_temp(193:3469);

% Reflection data without free surface multiples
rdwom_temp = conv(ss, pq);
rdwom = rdwom_temp(193:3469);

% Calculate energy of constructed reflection data before and after FSMR 
Erdwm = sum((abs(fft(rdwm))).^2) % Total amplitudes of constructed reflection data before FSMR
Erdwom = sum((abs(fft(rdwom))).^2) % Total amplitudes of constructed reflection data after FSMR

figure
subplot(2,1,1)
plot(x, Rs_conv, x, pq,'r'); axis tight
title('reflection response before(blue) and after(red) FSMR')
subplot(2,1,2)
plot(x, rdwm, x, rdwom,'r'); axis tight
title('reflection data before(blue) and after(red) FSMR')
echo on;
%
%  Now we compare the results before and after free surface multiple
%  removal.  The top figure shows the impulse response.  The bottom
%  compares how real seismic data might compare with and without 
%  free surface multiples removed.  The point of the later is that
%  the wavelet is essential to make all this work.  The difference
%  between data with and without multiples is very subtle in 
%  the simulated data, but obvious for the impulse response.
%  Note that in this model the two primaries are a reflection 
%  from base of sediments and a reflection from the crust-mantle
%  boundary.  The residual small bump at 27 s is an internal multiple
%  that is not removed by this algorithm.
%
pause;
%
% Now we do the same thing for the transmission response.
% The form of the series is difference, but the algorithm is
% similar and can be illustrated as a movie as before.
echo off;

%%% FSMR for transmission response %%%
 figure
% initialize and calculate the first term in the inverse scattering series
 for i = 1:1
     tri_yt = conv(T_n,gau); 
     pqt = tri_yt(1639:4915);
     h = plot(x,pqt, 'k');
     xlabel('travel time(second)')
     title('transmission response free surface back-scattered wave removal')
      set(h,'XData',x, 'YData', pqt);
       axis([0.025 82 -0.05 0.15])
    Trans(i) = getframe;
    pause(1)
end
 y_ct = y;
 yyt = 0;
for i = 2:9 % 9 is the total number of terms in the inverse scattering series used here
      y_ct = conv(y,y_ct);
       y_cct =y_ct(1:3277);
       yyt = yyt + y_cct;
       y_nt = y+yyt;
      y_tc = conv(T_n, y_nt);
     y_t = T_n + y_tc(1:3277)';       
     tri_yt = conv(y_t,gau);
     pqt = tri_yt(1639:4915);
     h = plot(x, pqt, 'k');
     xlabel('travel time(second)')
     title('transmission response free surface back-scattered wave removal')
      set(h,'XData',x, 'YData', pqt);
       axis([0.025 82 -0.05 0.15])
    Trans(i) = getframe;
    pause(1)
end


% write out the movie file to certain directory
% movie2avi(Trans,'D:\syn_t_r\results\trans2n\trfsmr_trans2n000.avi' ,'fps',1,'quality', 100)

echo on;
% Movie finished.  Notice that the transmission response is
% now a very close approximation to a single delta function 
% for this model.  This is because this is a simple, 3 layer model
% with a P wave at normal incidence.  The single residual pulse
% is an internal multiple between base of sediment and the crust-mantle
% boundary.
%
% Next we'll show the results of the multiple removal overlaid with original.
pause;

% Smooth the transmission response to the same as the smoothing used in the
% movie. To make them the same width of the pulse.
T_n_conv_temp = conv(T_n, gau);
T_n_conv = T_n_conv_temp(1639:4915);

% Calculate the energy of TR before and after FSMR
ET_n_conv = sum((abs(fft(T_n_conv))).^2); % Energy of the transmission response before FSMR
Epqt = sum((abs(fft(pqt))).^2); % Energy of transmission response after FSMR

% Transmission data with free surface multiples
tdwm_temp = conv(ss, T_n_conv);
tdwm = tdwm_temp(193:3469);

% Transmission data without free surface multiples
tdwom_temp = conv(ss, pqt);
tdwom = tdwom_temp(193:3469);

% Calculate the energy of the synthetic transmission data before and after FSMR
Etdwm = sum((abs(fft(tdwm))).^2); % Energy of the transmission data before FSMR
Etdwom = sum((abs(fft(tdwom))).^2); % Energy of transmission data after FSMR

figure
subplot(2,1,1)
plot(x, T_n_conv, x, pqt,'r'); axis tight
title('transmission response before(blue) and after(red) FSMR')
subplot(2,1,2)
plot(x, tdwm, x, tdwom,'r'); axis tight
title('transmission data before(blue) and after(red) FSMR')
echo on;
% This figure is like the one shown earlier for the reflection response.
% We see the impact of removing the free surface multiples is clear on 
% the impulse response, but subtle with the simulated data.
% This again indicates the importance of the wavelet in multiple removal
% methods that use the inverse scattering series.  That is, if we applied
% these same algorithms to blue in the lower figure, we would not produce
% the red curve at all.  It would not work at all.  These methods assume
% implicitly that they are given a Green's function type response function.

pause;
echo off;
%%%%% Use Kolmogorov method to construct min-phased TR from original TR %%%%%

% first to transform the TR with free surface multiples into min-phased TR
pqt_wc = fft(T_n_conv);
Ptr_lnf_p_poc = 0.5 * log((abs(pqt_wc)).^2);
Ptr_lna_p_poc = hilbert(0.5*Ptr_lnf_p_poc);
Ptr_lnfae_p_poc = exp(Ptr_lna_p_poc);
tr_min_tmp_p_poc = real(ifft(Ptr_lnfae_p_poc));
for i = 1:3277
    tr_min_p_poc(i) = tr_min_tmp_p_poc(3278-i);
end
tr_min_p_poc(3277) = 0;

% Second to transform the TR without free surface multiples into min-phased TR
% FSMR can be done directly on the min-phased TR 
pqt_w = fft(pqt);
Ptr_lnf_p_po = 0.5 * log((abs(pqt_w)).^2);
Ptr_lna_p_po = hilbert(0.5*Ptr_lnf_p_po);
Ptr_lnfae_p_po = exp(Ptr_lna_p_po);
tr_min_tmp_p_po = real(ifft(Ptr_lnfae_p_po));
for i = 1:3277
    tr_min_p_po(i) = tr_min_tmp_p_po(3278-i);
end
tr_min_p_po(3277) = 0;

% Calculate the free surface multiples of transmission response
fsbsw = T_n_conv - pqt;

figure
subplot(2,1,1)
plot(x, Rs_conv, x, pq,'r'); axis tight
legend('RR before FSMR', 'RR after FSMR')
title('reflection response before(blue) and after(red) FSMR')
subplot(2,1,2)
plot(x, -fsbsw); axis tight
title('free surface multiples of transmission response (with opposite sign of amplitude)')
echo on;
%
% This figure illustrates an alternative way to extract the reflection
% response from teleseismic data discovered by Chengliang Fan.  He
% noticed that the data removed from the transmission response by
% the free surface multiple removal algorithm was the reflection 
% response with a sign reversal.  This provides an alternative way of
% constructing the reflection response directly from the output of the
% multiple removal algorithm applied to the transmission data.  
% 
pause;
%
% This illustrates the application of the Kolmogorov algorithm to produce
% a minimum phase wavelet.  In this example it does little more than time shifts 
% the impulse response so the main P arrival pulse would arrive at zero lag.
% The reason is that these impulse response functions are minimum phase.
% Baig, Bostock, and Mercier use this approach in their statistical method
% based on the blind deconvolution conjecture.  We include it in this 
% demo mainly to show how it can be implemented in matlab. When the 
% minimum phase assumption can be justified this is may prove to be a useful
% component in future methodogies for estimating impulse response functions.
echo off;

figure
subplot(2,1,1)
plot(x, T_n_conv, x, pqt,'r'); axis tight
legend('TR before FSMR', 'TR after FSMR')
title('transmission response before(blue) and after(red) FSMR')
subplot(2,1,2)
plot(x, tr_min_p_poc, x, tr_min_p_po,'r'); axis tight
legend('min-phase TR before FSMR', 'min-phase TR after FSMR');axis tight
title('min-phase transmission response before(blue) and after(red) FSMR')

timeinv = toc; % calculate the total computing time

end