%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % 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)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)