% Importing intensity data for white light and HeNe laser
W_intensity = importdata('C:\Users\Asus G55\Desktop\FTS\Dec 17th\white_w.txt');
L_intensity = importdata('C:\Users\Asus G55\Desktop\FTS\Dec 17th\hene_w.txt');

% Subtract the average of the intensities to make the oscillations
% around 0
W_avg = mean(W_intensity);
W_intensity = W_intensity - mean(W_intensity);

L_avg = mean(L_intensity);
L_intensity = L_intensity - mean(L_intensity);

%Search for zero point crossings

% first look for exact zeros (Laser Light)
ind0_L = find( L_intensity == 0 );
% Now ind_0_L contains the indices of value 0

% then look for zero crossings between data points (Laser Light)
L1 = L_intensity(1:end-1) .* L_intensity(2:end);
ind1_L = find( L1 < 0 );
% Now ind_1_L contains the indices of negative values that are followed by
% positive ones or vice versa

% bring exact zeros and "in-between" zeros together (Laser Light)
ind_L = sort([ind0_L ind1_L]);

% Get "fractional" numbers of sample points at which the zeroes are
% estimated to occur through linear interpolation (Laser Light)
for (i=1:length(ind_L))
if (L_intensity(ind_L(i))==0)
ind_L(i)=ind_L(i);
elseif (L_intensity(ind_L(i))<0 && L_intensity(ind_L(i)+1)>0)
% linear interpolation to get fractional number of sample points
% corresponding to zero crossing
ind_L(i)=ind_L(i)+(abs(L_intensity(ind_L(i))))/(abs(L_intensity(ind_L(i)))+abs(L_intensity(ind_L(i)+1)));
elseif (L_intensity(ind_L(i))>0 && L_intensity(ind_L(i)+1)<0)
ind_L(i)=ind_L(i)+(abs(L_intensity(ind_L(i))))/(abs(L_intensity(ind_L(i)))+abs(L_intensity(ind_L(i)+1)));
elseif (L_intensity(ind_L(i))<0 && L_intensity(ind_L(i)-1)>0)
ind_L(i)=ind_L(i)-(abs(L_intensity(ind_L(i))))/(abs(L_intensity(ind_L(i)))+abs(L_intensity(ind_L(i)-1)));
elseif (L_intensity(ind_L(i))>0 && L_intensity(ind_L(i)-1)<0)
ind_L(i)=ind_L(i)-(abs(L_intensity(ind_L(i))))/(abs(L_intensity(ind_L(i)))+abs(L_intensity(ind_L(i)-1)));

end
end

% Now ind_L is expected to contain "mainly" non-integer values of sample
% points corresponding to zero point crossings

% Assign phases to all points

%Phases will be assigned to points lying between the first zero crossing
%and the last zero crossing
first=ceil(ind_L(1));
last=ceil(ind_L(end))-1;

% Initiate a vector of associated phases (Laser Light)
L_phase = zeros([1,last-first+1]);

for (j=2:length(ind_L))
for (i=first:last)
if (i>ind_L(j))
j=j+1;
break
end
if (i>=ind_L(j-1) && ind_L(j)>=i)
if (ind_L(j-1)>i-1 && i ~= first)
L_phase(i-first+1)= L_phase(i-first)+((i-ind_L(j-1))/(ind_L(j)-ind_L(j-1)))*pi+((ind_L(j-1)-(i-1))/(ind_L(j-1)-ind_L(j-2)))*pi;
elseif (ind_L(j-1)>i-1 && i == first)
L_phase(i-first+1)= L_phase(i-first+1)+((i-ind_L(j-1))/(ind_L(j)-ind_L(j-1)))*pi;
elseif (ind_L(j-1)<=i-1)
L_phase(i-first+1)= L_phase(i-first)+(1/(ind_L(j)-ind_L(j-1)))*pi;
end
end
end
end
% Now L_phase should contain the phases of all data points

% Define segments of equal phase
% Phase vector is the same length as the actual intensity vector to prevent
% FFT problems.
phases=linspace(L_phase(1), L_phase(end),length( W_intensity (first:last)));

% Fit points using cubic spline interpolation (White Light)
New_W_intensity = W_intensity (first:last);
W_proper=spline(L_phase,New_W_intensity,phases);

%%FTS part

I = W_proper;

%subtract avg of the intensity to make its avg = 0V
I_avg = mean(I);
I = I - mean(I);

%need to window out 0Thz to 100Thz frequency range

fs = 5000;
T = 1/fs;
V = 0.0499942e-3;
D = T*V;
Tao = D/3e8;
fs_tao = 1/Tao;
nyquist = fs_tao/(2*1e12) %nyquist frequency (THz)

m = length(I); % Window length
n = pow2(1.5*nextpow2(m)) % Transform length
y = fft(I,n); % DFT
f = (-n/2:n/2-1)*(fs_tao/n)/1e12; % Frequency range (THz)
nm = (1./(f*1e12))*3e8*1e9; %wavelength range (nm)
wl = ((1./f)*3e8)*1e-9; %wavelength (nm)
power = abs(y)/n; % Power of the DFT
power = power/max(power);

subplot(2,1,1);
plot(I);
xlabel('Bin');
ylabel('Intensity(V)');
title('Intensity(V) vs Bin: White light filter, Fs = 10kHz');
subplot(2,1,2);
plot(nm,fftshift(power));
axis([200 900 0 1]);
xlabel('Wavelength (nm)');
ylabel('Transmission');
title('Transmission vs Wavelength (nm)');