function [rfz,rf]=apply_spectral_whitening_deconvolution(z,r,W,a_large_gauss)
%--------------------------------------------------------------------------
%   Matlab function to compute spectrally whitened auto-correlation
%           and cross-correlation of two signals
%   (eq. 18-21 from Tauzin et al. GJI 2019)
%
%   input
%           z, structure containing trace for auto-correlation
%               z.trace: trace
%               z.npts: number of points in trace
%               z.delta: time step
%               z.o: delay from 1st sample
%           r, structure containing trace cross-correlated with z-trace
%               see z above for fields in structure
%           W, width of window for spectral whitening (in Hz)
%                   0.1-0.4 for Moho
%                   1.0-4.0 for shallow crust
%           a_large_gauss, define the cut-off of the gaussian low-pass
%               filter; for instance see here,
%               http://eqseis.geosc.psu.edu/~cammon/HTML/RftnDocs/seq01.html
%
%   output
%           rfz, structure containing auto-correlation of 1st signal (z)
%           rf, structure containing cross-correlation (r z*)
%               format similar to z, r above
%
%   Notes:
%
%     * If z is vertical component, r radial, then
%       cross-correlation gives radial receiver function 
%       (eq. 21; Tauzin et al., 2019)
%
%     * Additional parameters hard-coded
%           t_delay, delay before zero-lag in output correlations
%           centre_gauss, center of gaussian filter
%
%     * Depending on z.npts, smoothing at low frequency and spectral
%       division may introduce a long-period noise. In that case, I
%       recommand to high-pass filter, e.g.
%               fn=1/(2*z.delta); [bf,af]=butter(2,.1/fn,'high');
%               rfz.trace=filtfilt(bf,af,rfz.trace);
%               rf.trace=filtfilt(bf,af,rf.trace);
%
%     * A minus sign is applied in the output cross-correlation. 
%       It depends on the convention for the orientation of radial
%       component
%
%     * Users are requested to cite (as appropriate):
%
%           Tauzin, B., Pham, T. S., & Tkalčić, H. (2019).
%           Receiver functions from seismic interferometry:
%           a practical guide.
%           Geophysical Journal International, 217(1), 1-24.
%
%     * Modified from Matlab functions from
%           J. Vergne,
%           Ecole Et Observatoire des Sciences de la Terre, Strasbourg, FR
%
% Benoit TAUZIN (PhD)
% benoit.tauzin@univ-lyon1.fr / benoit.tauzin@anu.edu.au
%
% Research School of Earth Science
% Australian National University, Canberra, Australia
%
% Laboratoire de Geologie de Lyon
% Universite de Lyon, Villeurbanne, France
%
% last updated 08.02.2019
% ------------------------------------------------------------------------

switch nargin
    case 2
        W=0.1;                      % Default value (target Moho)
        a_large_gauss=5.0;          % Default value (target Moho)
end

% ADDITIONAL PARAMETERS ---------------------------------------------------
centre_gauss=0;                     % CENTER OF THE GAUSSIAN FILTER
t_delay=20;                         % DELAY BEFORE ZERO-LAG
%--------------------------------------------------------------------------

Data=[z.trace r.trace]';
m=z.npts; 
n=2^nextpow2(m); n=2*n;             % MODIFIED LENGTH FOR FFT
Data=fft(Data,n,2);                 % FFT
conz1=conj(Data(1,:));              % CONJUGATE OF VERTICAL COMPONENT
Data=Data.*[conz1;conz1];           % NUMERATORS EQ. 18, EQ. 21

df=1/(n*z.delta);                   % FREQUENCY STEP
dw=2*pi*df;                         % ANGULAR FREQUENCY STEP
w=[0:n-1]*dw;                       % ANGULAR FREQUENCY VECTOR
a2=4*(a_large_gauss).^2;            % GAUSSIAN FILTER PARAMETER

% APPLY PHASE DELAY AND GAUSSIAN FILTER (EQ. 18, 21)
gw=exp(-((w-centre_gauss*pi*2).*(w-centre_gauss*pi*2))/a2-w*t_delay*1i);	
Data2=Data.*[gw;gw];

% NUMBER OF SAMPLES OVER WHICH TO SMOOTH THE SPECTRUM
N=floor(W/(2*df));

% SMOOTH SPECTRUM OF VERTICAL COMPONENT AND DECONVOLVE
% I.E. EQ. 20; DENOMINATOR EQ. 18, 21
d=[];
for i=1:n % this loop could be avoided
    N0=length(max([1 i-N]):min([n i+N]));
    d(i)=1/N0*sum(abs(Data(1,max([1 i-N]):min([n i+N]))));
end

% SPECTRAL DIVISION (EQ. 18, 21)
Data_dec=Data2./[d;d];

% BACK TO TIME DOMAIN
Data_time=real(ifft(Data_dec,[],2));
Data_time=Data_time(:,1:m)/max(Data_time(1,:));

% OUTPUT AUTO-CORRELATION
rfz=z;
rfz.trace=Data_time(1,:)';
rfz.npts=length(rfz.trace);
rfz.o=t_delay;

% OUTPUT CROSS-CORRELATION
rf=z;
rf.trace=-Data_time(2,:)';          % MINUS SIGN DEPENDS ON ROTATION CONV.
rf.npts=length(rf.trace);
rf.o=t_delay;