%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Put data files in the same folder of this code: The first one is the 
% time-delayed spectral interferograms (SI) and the other is the spectrum 
% of the local oscillator. 
% Note: The wavelength mesh of the two files must be identical. 
%
% The analyzed results include
% (1) the figure of the input SI, and the corresponding temporal profile
% (2) the retrieved modulus and phase by wavelet-transform,
% (3) the magnitude and real of the echo field in the 2D plot (absorption 
% frequency (THz: y-axis) vs emission frequency: x-axis). 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;
clc;
%%
% Names of the input files and some flag parameters
SI_data_input='cg_tauScan2.txt';
LO_spectrum_input='LOSpec.txt' ;
% Choose delay_flag from 1: Echo-scan; 2: T-scan; 3: LO-scan
delay_flag=1;
% Time-zero position (mm) of the selected delay-line
time_zero=0.008 ;
%%
% Select the CCD_flag for camera= 1: NIRCCD; 2: VISCCD
CCD_flag=1;
if CCD_flag==1
    widt=256;    % the number of pixels of the NIR CCD
    Nw = 512;    % size of the linearized frequency mesh  
elseif CCD_flag==2
    widt=1340;   % the number of pixels of the Visible light CCD
    Nw = 1024;   % size of the linearized frequency mesh 
end
%% data input
% Input spectrum of local oscillator (LO): function of lambda [nm]
fid=fopen(LO_spectrum_input);
lambda=textscan(fid, '%f ','delimiter','\n','HeaderLines', 1);
lambda=lambda{1,1};
lambda(end,:)=[];         % a column vector array of [widt,1]
fid=fopen(LO_spectrum_input);
    t_delay=[];
    t_delay=textscan(fid, '%f ','delimiter','\n','HeaderLines', 4);
    t_delay=t_delay{1,1};
    t_delay(end,:)=[];
    time_delay=[];
    time_delay=transpose(t_delay); % time_delay: a row vector
fid=fopen(LO_spectrum_input);
    t_delay=[];
    t_delay=textscan(fid, '%f ','delimiter','\n','HeaderLines', 7);
fid=fopen(LO_spectrum_input);
    t_delay=[];
    t_delay=textscan(fid, '%f ','delimiter','\n','HeaderLines', 10);
    t_delay=[];
fid=fopen(LO_spectrum_input);
    intensity_temp=textscan(fid, '%f ','delimiter','\n','HeaderLines', 13);
    intensity_temp=intensity_temp{1,1};
fclose(fid);
%
nt=length(time_delay);
nlam=length(lambda);
intensity=[];
for i=1:nt
    intensity(:,i)=intensity_temp(1+nlam*(i-1):nlam*i,1);
end
ref=zeros(length(lambda),1);
for it=1:nt
    ref(:,1)=ref(:,1)+intensity(:,it);
end
ref(:,1)=ref(:,1)./nt;
for k=1:length(ref)
    if(ref(k,1)<1)
        ref(k,1)=1.0;
    end
end
% Input SI at the specified delay times
lambda=[];
fid=fopen(SI_data_input);
lambda=textscan(fid,'%f ','delimiter','\n','HeaderLines', 1);
lambda=lambda{1,1};
lambda(end,:)=[];     % a column vector array of [widt,1]
nlam=length(lambda);
%
if delay_flag==1
    fid=fopen(SI_data_input);
    t_delay=[];
    t_delay=textscan(fid,'%f ','delimiter','\n','HeaderLines', 4);
    t_delay=t_delay{1,1};
    t_delay(end,:)=[];
    time_delay=[];
    time_delay=transpose(t_delay); % time_delay: a row vector [1,nt]
    fid=fopen(SI_data_input);
    t_delay=[];
    t_delay=textscan(fid,'%f ','delimiter','\n','HeaderLines', 7);
    fid=fopen(SI_data_input);
    t_delay=[];
    t_delay=textscan(fid,'%f ','delimiter','\n','HeaderLines', 10);
    t_delay=[];
elseif delay_flag==2
    fid=fopen(SI_data_input);
    t_delay=[];
    t_delay=textscan(fid,'%f ','delimiter','\n','HeaderLines', 4);
    t_delay=t_delay{1,1};
    t_delay(end,:)=[];
    time_delay=[];
    time_delay=transpose(t_delay);
    fid=fopen(SI_data_input);
    t_delay=[];
    t_delay=textscan(fid,'%f','delimiter','\n','HeaderLines', 7);
    t_delay=[];
    fid=fopen(SI_data_input);
    t_delay=[];
    t_delay=textscan(fid,'%f ','delimiter','\n','HeaderLines', 10);
    t_delay=[];
elseif delay_flag==3
    fid=fopen(SI_data_input);
    t_delay=[];
    t_delay=textscan(fid,'%f ','delimiter','\n','HeaderLines', 4);
    t_delay=t_delay{1,1};
    t_delay(end,:)=[];
    time_delay=[];
    time_delay=transpose(t_delay);
    fid=fopen(SI_data_input);
    t_delay=[];
    t_delay=textscan(fid,'%f ','delimiter','\n','HeaderLines', 7);
    t_delay=[];
    fid=fopen(SI_data_input);
    t_delay=textscan(fid,'%f ','delimiter','\n','HeaderLines', 10);
    t_delay=[];
end
fid=fopen(SI_data_input) ;
intensity_temp=[];
intensity_temp=textscan(fid,'%f ','delimiter','\n','HeaderLines', 13);
intensity_temp=intensity_temp{1,1};  % a row array of [nt*widt,1]
fclose(fid);
%% Preprocessing data file
intensity=[];
nt=length(time_delay);
for i=1:nt
    intensity(:,i)=intensity_temp(1+nlam*(i-1):nlam*i,1);
end
for it=1:nt
   for j=1:length(intensity(:,it))
      if(intensity(j,it)<1)
         intensity(j,it) =1.0;
      end
   end
end
%% Some scalings
% Lambda: wavelength array [in nm] and the corresponding freq [in THz]
% c=299.792458 [nm/fs], freq [THz]=1000*c[nm/fs]/lambda[nm].
% omega [rad/fs] = 2000*pi*freq [THz]
c = 299.792458;              % [nm/fs]
chng_rng = @(x) 2*pi*c./x;   % Change from [nm] to [rad/fs]
omega = chng_rng(lambda);    % angular frequency mesh [rad/fs]
freq = omega.*159.155 ;      % frequency mesh [THz] for SI
freq_max = freq(1);
freq_min = freq(widt);
dfreq = (freq_max-freq_min)/(Nw-1);
freq_lin = freq_max-(0:Nw-1)'*dfreq ;  % Linear freq grid (1-column array)
time_delay = (time_delay-time_zero)*(2*10^6/c); % in [fs](1-row array)
tdelay_max = max(time_delay);
tdelay_min = min(time_delay);
% Arrange the time_delay to be increasing from negative to positive
if (time_delay(1)>time_delay(end))
   int_rev=[];
   dt_rev=[];
   for j=1:nt
      int_rev(:,nt-j+1)=intensity(:,j) ;
      dt_rev(nt-j+1)=time_delay(j) ;
   end 
   for j=1:nt
      intensity(:,j)=int_rev(:,j) ;
      time_delay(j)=dt_rev(j) ;
   end 
   int_rev=[];
   dt_rev=[];
end
%% Plotting time-delayed SI data
figure(1)
% intensity(y=freq,x=td), intensity'(y=td,x=freq)
imagesc(freq,time_delay,transpose(intensity));
set(gca,'FontSize',16);
xlabel('\fontsize{18}Emission Frequency (THz)');
ylabel('\fontsize{18}Delay Time (fs)');
title('\fontsize{20}\bfSpectral Interferograms');
colormap(hot) ;
cb=colorbar;
set(cb,'FontSize',16);
%% Plotting temporal profile of SI
bkg=zeros(nlam,1);  %1-column array
tdprf=zeros(nt,1);  %1-column array
for it=1:10
    bkg(:,1)=bkg(:,1)+intensity(:,it);
end
bkg(:,1)=bkg(:,1)./10 ;
for it=1:nt
    sum=0.;
    for jw=1:length(intensity(:,it))
      sum=sum+intensity(jw,it)-bkg(jw,1) ;
    end
    tdprf(it,1)=sum;
end
tdprf_min=min(tdprf);
tdprf(:,1)=tdprf(:,1)-tdprf_min;
% Find the peak position
[tdprf_max,imax]=max(tdprf);
tdprf(:,1)=tdprf(:,1)/tdprf_max;
origin=time_delay(imax)*0.5*c*10^-6+time_zero
fid = fopen('tdprf.txt', 'w');
for it=1:nt
    fprintf(fid, '%12.4g %12.4g\r\n', time_delay(it),tdprf(it,1));
end
fclose(fid);
figure(2)
plot(time_delay',tdprf,'-r','LineWidth',2);
axis([tdelay_min tdelay_max 0. 1.]);
set(gca,'FontSize',16);
xlabel('\fontsize{18}Time Delay (fs)');
ylabel('\fontsize{18}Intensity (arb. units)');
title('\fontsize{18}\bfTime-Delayed Integrated SI');
%
%% Mapping LO spectrum to the linearized mesh freq_lin
ref_lin = spline(freq, ref, freq_lin) ;
for k=1:Nw
    if(ref_lin(k,1)<1)
        ref_lin(k,1) =1.0;
    end
end
ref=[];
for k=1:Nw
    ref_lin(k,1) =sqrt(ref_lin(k,1)) ;
end
%% Retrieving phase information from SI with Wavelet transform
%
% Note: tau_w is the time separation between echo and LO, which
% can vary from 1-ps (0.33 mm) to 10.0 ps (3.3 mm). The corresponding angular
% frequency spacing between two neighboring SI peaks varies from 6 [THz]
% for 1-ps delay to 0.6 [THz] for 10-ps delay. 
%
% Set up the scales of the Morlet wavelets.
% Low scale values compress the wavelet and correlate better with high 
% frequencies and the fine-scale features in the input signal are better 
% revealed.
amin=0.1 ;   
% High scale values stretch the wavelet and correlate better with the low 
% frequency content of the signal.
amax=30;        
step=0.2;                  
a=[amin:step:amax];
na=length(a) ;
wphy=zeros(Nw,nt);
phase=zeros(Nw,nt);
echo_amp=zeros(Nw,nt);
sig_lin=zeros(Nw,1);
for ii=1:nt
    wfab=[];
    mwf=[];
    pwf=[];
    sig_lin(:,1) = spline(freq, intensity(:,ii), freq_lin) ;
    wfab=cwt(sig_lin,a,'cmor1-1');
                    % coefs = cwt(x,scales,'wname'):
                    % Scales  determine the degree to which the wavelet
                    % is compressed or stretched.
                    %
                    % wname='cmorlFb-Fc' selects the complex Morlet wavelet 
                    % function of [1/sqrt(pi*Fb)]*exp[2*pi*i*Fc*x]*exp[-x^2/Fb],
                    % where Fb is a bandwidth parameter, and Fc the central
                    % frequency. The pseudo-frequency[Hz]=Fc/(a*delta)
                    % with a=scale factor and delta the sampling period. 
    mwf=abs(wfab);   
    pwf=angle(wfab);  % angle(Z) returns the phase angles, in radians, for each 
                      % element of complex array Z. 
                      % The phase angles lie between [-pi, +pi].
% For each spectral components {m}, search for the position of the maximum 
% modulus 
    for m=1:Nw        
        maxwtm=mwf(1,m);
        k=1;
        for n=1:na
            if mwf(n,m)>maxwtm
                maxwtm=mwf(n,m);
                k=n;
            end
        end
        wphy(m,ii)=pwf(k,m); % Retrieve the phase at k (the position of 
                             % maximum modulus.
        ml(m,ii)=k;
        echo_amp(m,ii)=mwf(k,m)*(2.0/sqrt((amin+k*step))); 
                                  % The amplitude of the oscillating
                                  % component can be deduced from
                                  % abs(C(a_ridge,b))*2/sqrt(a_ridge),
                                  % giving echo_amp(m,ii)=0.505*mwf(k,m);
    end
    phase(:,ii)=unwrap(wphy(:,ii)); % unwrap the phase profle.
end
%% Removing linear phase from the phase profiles
fmfit=zeros(11,1);
prf=zeros(11,1);
for k=1:11
    prf(k)=phase(Nw/2+k-6,imax);
    fmfit(k)=freq_lin(Nw/2+k-6);
end
p = polyfit(fmfit,prf,1);
for ii=1:nt
    for j=1:Nw
        phase(j,ii)=phase(j,ii)-polyval(p,freq_lin(j));
    end
end
%% Removing the local oscillator amplitude 
for ii=1:nt
    echo_amp(:,ii)=echo_amp(:,ii)./ref_lin(:,1); 
end
%% Plotting the amplitude and phase of echo field Es'(y=tau, x=freq_lin)
figure(3)
imagesc(freq_lin,time_delay,phase');
set(gca,'FontSize',16);
title('\fontsize{18}\bfSpectral Phase (rad)');
xlabel('\fontsize{18}Emission Frequency (THz)');
ylabel('\fontsize{18}Delay Time (fs)');
colormap(hot) ;
cb=colorbar;
set(cb,'FontSize',16);

figure(4)
imagesc(freq_lin,time_delay,echo_amp');
set(gca,'FontSize',16);
title('\fontsize{18}\bfMagnitude');
xlabel('\fontsize{18}Emission Frequency (THz)');
ylabel('\fontsize{18}Delay Time (fs)');
colormap(hot) ;
cb=colorbar;
set(cb,'FontSize',16);
%% tau-FFT of the echo field
delt=abs(time_delay(2)-time_delay(1));  % [fs]
N_ex=512 ;
tmax=delt*N_ex/2;
t_ex = -tmax+(0:N_ex-1)*delt ;
teps=0.2*delt;
td_1=time_delay(1);
k=1;
for n=1:N_ex
    if abs(t_ex(n)-td_1)<teps
       k=n;
    end
end
k_st=k;
echo_field=zeros(Nw,N_ex);
echo_a=zeros(N_ex,Nw);
% echo_p=zeros(N_ex,Nw);
%
% Assemble the time-dependent echo field spectrum E(tau,freq), in which 
% tau is the time_delay between the first two pulses.
% According to Nyquist theorem, at a sampling rate 1/delt, the highest 
% frequency that can be correctly sampled is f_max=500(THz)/delt(fs).
% For delt= 3.34-fs, f_max=150 THz. Any signal with frequency f>fmax will 
% be mapped to f-2*fmax.
fmax=500./delt
delf=fmax/(N_ex/2)
w_ex = -fmax+(0:N_ex-1)'*delf ;
% Note: tau shall be positive for rephasing process, so here we insert 
% negative sign in tk.
for k=1:nt
    echo_field(:,k-1+k_st)=echo_amp(:,k).*exp(sqrt(-1)*phase(:,k)) ;
end
for n=1:Nw
    echo_a(:,n)=fftshift(fft(transpose(echo_field(n,:)),N_ex));
%    echo_p(:,n)= unwrap(angle(echo_a(:,n)));
end
%% Data Preparation for Plotting
f_max=max(freq_lin);
f_min=min(freq_lin);
df=(f_max-f_min)/99.;
dLine_x=f_min+(0:99)*df;
dLine_y=-dLine_x;
yaxis=zeros(N_ex,1);
yaxis(:,1)=w_ex(:)-2*fmax;
% Positive frequency w1 in Echo appears at -w1 in the w_ex mesh because of
% delta(w_ex+w1). However due to the aliasing effect from w1>fmax, the 
% spectral peak instead appears at 2*fmax-w_ex. So the real frequency -w1
% shall be presented with w_ex-2*fmax. 
feps=0.4*delf;
kmin=1;
for n=N_ex/2:N_ex
    if abs(yaxis(n)+f_max)<feps
       kmin=n ;
    end
end
kmax=1;
for n=N_ex/2:N_ex
    if abs(yaxis(n)+f_min)<feps
       kmax=n;
    end
end
knum=abs(kmax-kmin)+1
twoD_sigr=zeros(knum,Nw);
twoD_sigm=zeros(knum,Nw);
twoD_yaxis=zeros(knum,1);
for ii=kmin:kmax
    twoD_yaxis(ii-kmin+1,1)=yaxis(ii) ;
    twoD_sigr(ii-kmin+1,:)=real(echo_a(ii,:));
    twoD_sigm(ii-kmin+1,:)=abs(echo_a(ii,:));
end
%% Plotting the 2D data
figure(5)
contourf(freq_lin,twoD_yaxis,twoD_sigm,10);
set(gca,'YDir','normal');
axis([f_min f_max -f_max -f_min],'square');
hold on;
plot(dLine_x,dLine_y, '-y','LineWidth',1);
set(gca,'FontSize',16);
title('\fontsize{18}\bfFFT: Abs(Echo)');
xlabel('\fontsize{18}Emission Frequency (THz)');
ylabel('\fontsize{18}Excitation Frequency (THz)');
colormap(jet) ;
cb=colorbar;
set(cb,'FontSize',16);
%
figure(6)
contourf(freq_lin,twoD_yaxis,twoD_sigr,10);
hold on;
plot(dLine_x,dLine_y, '-w','LineWidth',1);
set(gca,'YDir','normal');
axis([f_min f_max -f_max -f_min],'square');
set(gca,'FontSize',16);
title('\fontsize{18}\bfFFT: Real(Echo)');
xlabel('\fontsize{18}Emission Frequency (THz)');
ylabel('\fontsize{18}Excitation Frequency (THz)');
colormap(jet) ;
cb=colorbar;
set(cb,'FontSize',16);



