SPIDER_simulation.m

%	SPIDER_simulation.m
%	Spectral phase interferometry for direct electric-field reconstruction
%
%	This script simulates SPIDER acquisition and reconstruction.
%	A pulse with chirp according to the dispersion in BK7 is generated and a SPIDER
%	interferogram as recorded by a spectrometer (i.e. with noise and quantisation) is 
%	simulated. The interferogram is then inverted using the SPIDER reconstruction 
%	algorithm. The reconstructed & simulated pulses are then compared.

%%	Some useful functions / constants
c = 299.792458;					%	Speed of light (nm/fs)
chng_rng = @(x) 2*pi*c./x;		%	Change nm <--> rad/fs

%	Gaussian function:
%		x	=	Co-ordinate
%		x0	=	Centre
%		dx	=	FWHM
%			= sqrt(2*log(2))*sigma_x, where sigma_x is 1/e^2 half width
%		n	=	Order (normal = 1)
gauss1D = @(x, x0, dx, n) exp(-log(2)*(2*(x-x0)/dx).^(2*abs(floor(n))));

%%	Data ranges
Nw = 1024;						%	Number of freq/wavelength data points

%	Spectrometer range
%	Note:
%		The data acquired in an experiment is taken with a spectrometer
%		which is approximately linear in wavelength, thus all calculations
%		will be used with a linear wavelength grid (lambda) and it's
%		corresponding nonlinear angular frequency grid (omega). Later on, I
%		shall create a linear angular frequency grid (w). From now on, I
%		refer to angular frequency as just frequency (or freq.)

%	Wavelength
lambda_min = 650;
lambda_max = 950;
dlambda = (lambda_max-lambda_min)/(Nw-1);
lambda = (0:Nw-1)'*dlambda + lambda_min;

%	Angular frequency
omega = chng_rng(lambda);
omega_max = max(omega);
omega_min = min(omega);
domega = (omega_max-omega_min)/(Nw-1);

%%	Dispersion of BK7
%	Sellmeier co-efficients (from Wiki)
B = [1.03961212,	2.31792344e-1,	1.01046945];
C = [6.00069867e-3,	2.00179144e-2,	1.03560653e2];

%	Refractive index
%	NOTE: lambda is assumed to be in micrometers!!
n = @(lambda) sqrt(1 + sum(((lambda.^2)*B) ./ ...
	((lambda.^2)*ones(1,length(C)) - ones(Nw,1)*C), 2));

%	Phase accumulated (rad)
%	NOTE: lambda & z are assumed to be in nm
phi_z = @(lambda, z) n(lambda*1e-3).*(2*pi./lambda)*z;

%%	Generate pulse
Dt = 30;						%	Fourier transform limited (FTL) FWHM intensity pulse duration (fs)
Dw = 4*log(2)/Dt;				%	FWHM intensity bandwidth (rad/fs) for FTL duration Dt
lambda_0 = 800;					%	Central wavelength (nm)
w0 = chng_rng(lambda_0);		%	Central frequency (rad/fs)


%	Notes:
%		The delay cannot be found, only relative to some other pulse,
%		therefore the absolute delay needs to be removed. Zero delay is at
%		the entrance of our imaginary glass (BK7), so the delay at w0 is
%		removed (equal to the gradient of the phase at w0). The absolute
%		phase is also removed. Technically, this does not need to be done
%		as SPIDER (or any other self-referencing technique) cannot measure
%		these quantities, and so this is only for a comparison between the
%		reconstructed and 'actual' phase.
%
%		The full width half max (FWHM) values are quoted for the intensity
%		(in boh the time and (angular) frequency domains, so this value
%		needs to be scaled fo the electric field (i.e. by sqrt(2) for
%		Gaussian pulses). For non-Gaussian pulses, the sqrt of the
%		intensity can be taken.
phi = phi_z(lambda, 2e7);						%	Phase accumulated from 2cm BK7
t0 = (spline(omega, phi, w0+w0/1e3) ...
	- spline(omega, phi, w0-w0/1e3))/(w0/.5e3);	%	Delay from BK7
phi = phi-omega*t0;								%	Remove delay
phi = phi - spline(omega, phi, w0);				%	Remove absolute phase
Ew = gauss1D(omega, w0, sqrt(2)*Dw, 1).*exp(i*phi);			%	Pulse in freq domain

%%	SPIDER Parameters
delay = 500;
shear = 2*pi/(30*Dt);		%	Maximum pulse window is 10*Dt

%	Note:
%		In a real experiment, the data will contain noise (additive,
%		multiplicative & quantisation). The data should always be taken
%		with a small background, which can be removed via software
%		processing. Hardware removal might remove real data which cannot be
%		recovered. I assume noise is noramally distributed.
alpha_a = .01;				%	Addative noise
alpha_m = .01;				%	Multiplicative noise
alpha_q = 16;				%	Bit depth (quantisation error)
err_delay = 0e-3;			%	Pulse-pulse fluctuations in the delay

%	Normally, the trace is acquired over many laser shots, with a
%	fluctuating delay between the pulses, thus an average over many pulses
%	is acquired.
N_acq = 1;					%	Number of pulses to integrate over

Iw = zeros(Nw, N_acq);
Iw_cal = Iw;

Ew1 = spline(omega, Ew, omega-shear).*(omega>=shear);
Iw_true = abs(Ew + Ew1.*exp(i*omega*delay)).^2;
Iwc_true = abs(Ew.*(1+exp(i*omega*delay))).^2;
for n=1:N_acq
	%	SPIDER Interferogram
	Iw(:,n) = abs(Ew + Ew1.*exp(i*omega*delay*(1+randn(1)*err_delay))).^2;
	%	Calibration Interferogram
	Iw_cal(:,n) = abs(Ew.*(1+exp(i*omega*delay*(1+randn(1)*err_delay)))).^2;
end

%	Add noise & quantise
Iw_n = sum(Iw, 2);									%	Integrate
Iw_n = (Iw_n + randn(Nw, 1)*alpha_a*max(Iw_n));		%	Additive noise
Iw_n = Iw_n.*(1 + randn(Nw, 1)*alpha_m);			%	Multiplicative noise
Iw_n = Iw_n/max(Iw_n);								%	Normalise
Iw_n = floor((2^(alpha_q-1)-1)*Iw_n);				%	Quantisation

Iwc_n = sum(Iw_cal, 2);								%	Integrate
Iwc_n = Iwc_n + randn(Nw, 1)*alpha_a*max(Iwc_n);	%	Additive noise
Iwc_n = Iwc_n.*(1 + randn(Nw, 1)*alpha_m);			%	Multiplicative noise
Iwc_n = Iwc_n/max(Iwc_n);							%	Normalise
Iwc_n = floor((2^(alpha_q-1)-1)*Iwc_n);				%	Quantisation

%%	SPIDER Reconstruction

%	Note:
%		For large bandwidths, the fringes in the calibration trace will not
%		appear to have the same period due to nonlinear frequency grid.
%		Thus, one should interpolate the data onto a linear frequency grid
%		1st (or perform a non-linear FFT routine), otherwise it may not be
%		possible to isolate the sidebands.
omega_lin = (0:Nw-1)'*domega + omega_min;			%	Generate linear freq grid
pxls = (1:Nw)';										%	pxl grid (for filter)

%	Fourier transform data.
%	Note:
%		It doesn't matter which Fourier transform you perform (forward or
%		backward). The only difference is the phase will have opposite
%		signs, so it may be necessary to flip the sign of the extracted
%		phase before integration/concatenation. Note that this is a
%		deterministic effect (do the maths with the correct Fourier
%		transform, shear and sidelobe.)
It = fftshift(fft(Iw_n));
Itc = fftshift(fft(Iwc_n));
It_lin = fftshift(fft(spline(omega, Iw_n, omega_lin)));
Itc_lin = fftshift(fft(spline(omega, Iwc_n, omega_lin)));

%	Filter
%	Note:
%		Both sidelobed contain the exact amount of data, thus it does not
%		matter which one you filter. However, as they are complex
%		conjugates, the sign of the phase is opposite for each. Thus, if
%		you find the reconstructed phase is inverted, simply flip the
%		sign of the phase before integration/concatenation. Note that this is a
%		deterministic effect (do the maths with the correct Fourier
%		transform, shear and sidelobe.)
ft = gauss1D(pxls, 583, 40, 4);
Iwf = ifft(ifftshift(It.*ft));
Iwcf = ifft(ifftshift(Itc.*ft));
Iwf_lin = ifft(ifftshift(It_lin.*ft));
Iwcf_lin = ifft(ifftshift(Itc_lin.*ft));

%	Extract phase & remove calibration
%	Note:
%		Because omega = 2*pi/lambda = (0:-1:-(Nw-1))'*dw +wmax, the phase
%		for the linearly interpolated data is inverted with respect to the
%		other.
%
%		If the shear is small enough to allow the pulse to be sampled
%		adequately, then the phase extracted should not go outside +/-pi,
%		therefore this is a good check of the routine. This also means
%		unwrapping is unecessary, although it could be used if one needed.
theta = unwrap(angle(Iwf.*conj(Iwcf)));
theta_lin = -unwrap(angle(Iwf_lin.*conj(Iwcf_lin)));

%	New frequency ranges - this is beacuse the reconstructed phase falls on
%	the frequency grid between the two sheared pulses.
omega_s = omega-shear/2;
omega_s_lin = omega_lin-shear/2;

%	Integrate phase via trapezium (mid-point) rule
phi_i = cumtrapz(omega, (theta-spline(omega_s, theta, w0))/shear);	%	Integrate (removing linear slope)
phi_i = phi_i - spline(omega_s, phi_i, w0);							%	Remove absolute phase

phi_i_lin = cumtrapz(omega_s_lin, (theta_lin-spline(omega_s_lin, theta_lin, w0))/shear);	%	Integrate (removing slope)
phi_i_lin = phi_i_lin - spline(omega_s_lin, phi_i_lin, w0);			%	Remove absolute phase

%	Concatenate phase
%	Note:
%		The gradient is taken half-way between the two sheared pulses (i.e.
%		at omega-shear/2 = omega_s;
omega_c = (omega_min:abs(shear):omega_max)';

phi_c = cumsum(spline(omega, theta, omega_c)...
	-spline(omega_s, theta, w0));				%	Concatenate, removing linear phase
phi_c = phi_c - spline(omega_c, phi_c, w0);		%	Remove absolute phase

phi_c_lin = cumsum(spline(omega_lin, theta_lin, omega_c)...
	-spline(omega_s_lin, theta_lin, w0));			%	Concatenate, removing linear phase
phi_c_lin = phi_c_lin - spline(omega_c, phi_c_lin, w0);		%	Remove absolute phase

% plot(omega, phi, '.', omega_c, [phi_c phi_c_lin], omega_s, phi_i, omega_s_lin, phi_i_lin);
%%	Obtain spectral intensity
%	Assuming an accurate, independent spectral measurement can be obtained
Ew_sp = abs(Ew)/max(abs(Ew));

%	Obtain |Ew| = sqrt(Iw) from the dc term in the interferogram, allowing
%	online retrieval. Note that this only works when the shear is small,
%	and the spectrum is relatively smooth.
ft0 = gauss1D(pxls, 513, 40, 4);
Iwf0 = ifft(ifftshift(It.*ft0));
Iwf0_lin = ifft(ifftshift(It_lin.*ft0));

Ew_int = sqrt(abs(Iwf0)/max(abs(Iwf0)));
Ew_int_lin = sqrt(abs(Iwf0_lin)/max(abs(Iwf0_lin)));

%%	Temporal info
%	Linear range in frequency
dw = shear;
w = (0:Nw-1)'*dw;
Nt = 4096;
dt = 2*pi/(Nt*dw);
t = ((0:Nt-1)'-floor(Nt/2))*dt;

%	Interpolate data onto linear freq grid
%	Naming convention:
%		sp	=	Used actual spectrum
%		int =	Used spectrum from interferogram
%		i	=	Used integrated phase
%		c	=	Used concatenated phase
%		lin =	Used interferogram which was interpolated on linear grid
fw = (w>=omega_min).*(w<=omega_max);
Ew_true = spline(omega, Ew, w).*fw;

Ew_sp_ = spline(omega, Ew_sp, w).*fw;
Ew_int_ = spline(omega_s, Ew_int, w).*fw;
Ew_int_lin_ = spline(omega_s_lin, Ew_int_lin, w).*fw;

phi_i_ = spline(omega_s, phi_i, w).*fw;
phi_c_ = spline(omega_c, phi_c, w).*fw;
phi_i_lin_ = spline(omega_s_lin, phi_i_lin, w).*fw;
phi_c_lin_ = spline(omega_c, phi_c_lin, w).*fw;

Ew_sp_i = Ew_sp_.*exp(i*phi_i_);
Ew_sp_c = Ew_sp_.*exp(i*phi_c_);
Ew_int_i = Ew_int_.*exp(i*phi_i_);
Ew_int_c = Ew_int_.*exp(i*phi_c_);

Ew_sp_i_lin = Ew_sp_.*exp(i*phi_i_lin_);
Ew_sp_c_lin = Ew_sp_.*exp(i*phi_c_lin_);
Ew_int_i_lin = Ew_int_.*exp(i*phi_i_lin_);
Ew_int_c_lin = Ew_int_.*exp(i*phi_c_lin_);

%	Calculate temporal profile
%	The phase term is to shift the carrier frequency to the true carrier
%	frequency. Padding is used to obtain a higher temporal resolution
Et_true = fftshift(fft(Ew_true, Nt));%.*exp(-i*min(w)*t);

Et_sp_i = fftshift(fft(Ew_sp_i, Nt));%.*exp(-i*min(w)*t);
Et_sp_c = fftshift(fft(Ew_sp_c, Nt));%.*exp(-i*min(w)*t);
Et_int_i = fftshift(fft(Ew_int_i, Nt));%.*exp(-i*min(w)*t);
Et_int_c = fftshift(fft(Ew_int_c, Nt));%.*exp(-i*min(w)*t);

Et_sp_i_lin = fftshift(fft(Ew_sp_i_lin, Nt));%.*exp(-i*min(w)*t);
Et_sp_c_lin = fftshift(fft(Ew_sp_c_lin, Nt));%.*exp(-i*min(w)*t);
Et_int_i_lin = fftshift(fft(Ew_int_i_lin, Nt));%.*exp(-i*min(w)*t);
Et_int_c_lin = fftshift(fft(Ew_int_c_lin, Nt));%.*exp(-i*min(w)*t);


%%	Plot results
figure(1)
plot(lambda, [Iw_true/4 Iw_n/max(Iw_n)], lambda, .25*abs([Ew Ew1]).^2, '.-')
title('SPIDER interferograms and original sheared spectra');
xlabel('Wavelength /nm');
ylabel('Intensity /arb.u.')
legend('Perfect interferogram', 'Measured interferogram', 'E(\omega)', 'E(\omega+\Omega)');
figure(2)
plot(t, real(Et_true), '.-', t, real([Et_sp_i Et_sp_c]));
title('True and reconstructed temporal field: Independent spectrum');
xlabel('Time /fs');
ylabel('Electric field /arb.u.');
legend('True', 'Integrated', 'Concatenated');

figure(3)
plot(t, real(Et_true), '.-', t, real([Et_sp_i_lin Et_sp_c_lin]));
title('True and reconstructed temporal field: Independent spectrum (linear interpolation)');
xlabel('Time /fs');
ylabel('Electric field /arb.u.');
legend('True', 'Integrated', 'Concatenated');

figure(4)
plot(t, real(Et_true), '.-', t, real([Et_int_i Et_int_c]));
title('True and reconstructed temporal field: Extracted spectrum');
xlabel('Time /fs');
ylabel('Electric field /arb.u.');
legend('True', 'Integrated', 'Concatenated');

figure(5)
plot(t, real(Et_true), '.-', t, real([Et_int_i_lin Et_int_c_lin]));
title('True and reconstructed temporal field: Extracted spectrum (linear interpolation)');
xlabel('Time /fs');
ylabel('Electric field /arb.u.');
legend('True', 'Integrated', 'Concatenated');