% time2freq.m
%	Illustrates the conversion of time-resolved data
% into frequency-domain data using the fast Fourier transform
% function, fft(), in MATLAB.
%	Uses Diffusion Theory to generate time-resolved data F(r,t).
%	Uses fft() to convert F(r,t) --> F(r,f)
% where     F(r,f) = Fre + i*Fim
%           Magn   = Fre^2 + Fim^2 = F*conj(F)
%           Phase  = atan2(Fre, Fim)

clear

%%%%%%%%%%%%%%%%%%%%%%
% USER CHOICES 
%%%
%%% Introduce test frequency to test calibration of frequency domain
calib = 0;             % 1 = ON, 0 = OFF
testfreq = 1e10;       % [Hz]  CHOOSE A TEST FREQUENCY FOR CALIBRATION
testamplitude = 20e-4; % CHOOSE an amplitude for freq. calibration
musp = 10;             % reduced scattering coefficient [cm^-1]
%%%
%%%%%%%%%%%%%%%%%%%%%%


figure(1);clf
figure(2);clf
sym = 'rgbrgbmckmck';

% Run computation for 10 values of absorption (mua)
for j=1:10

mua = (j-1)	*0.1; % [cm^-1]
	
%%%%%
% Setup time-resolved diffusion theory problem
%%%%
Uo = 1;         % [J] impulse energy of source at time zero at origin
nt = 1.33;      % refractive index
c = 2.98e10/nt; % light speed [cm/s]
mfpp = 1/(mua + musp); % transport mean free path [cm]
D = mfpp/3;     % diffusion length [cm]
alpha = c*D;    % diffusivity for light [cm^2/s]
r = 1;          % distance of observation point from source

tmax = 8000e-12;          % [s] max time of interest
N = 2^12;                 % number of time points.
dt = tmax/(N-1);          % [s] time step
t  = (0:dt:tmax)' + dt/2; % [s] time

%%%% Time-resolved diffusion theory
Frt = Uo*c*exp(-r^2./(4*alpha*t))./(4*pi*alpha*t).^(3/2).*exp(-mua*c*t);
%%%%

%%% Introduce test frequency to test calibration of frequency domain
if calib  
	Frt = Frt + testamplitude*sin(2*pi*testfreq*t);
end
%%%

%%%%%% FIG 1a TIME-RESOLVED %%%%%%
figure(1)
subplot(2,1,1)
semilogx(t,Frt,[sym(j) '-'])
hold on
set(gca,'fontsize',18)
title(sprintf('@ r = %1.0f cm', r),'fontsize',18)
text(1e-11,1.9e9,'1-J impulse response','fontsize',18)
axis([1e-12, 1e-8, 0, 2e9])
xlabel('t [s]')
ylabel('F(r,t) [W/cm^2]')

drawnow

%%% take fft()  %%%%%%%
u  = fft(Frt,N);
FD = u(1:N/2);
P  = FD.*conj(FD);
P2 = real(FD).^2 + imag(FD).^2; % alternative method of calculating P
f  = 1/dt*(0:N/2-1)'/N; % only use 1st half of FD = fft(). Remainder is symmetric.
df = f(2) - f(1);
P = P/sum(P*df);
unity = sum(P*df)

%%%%%% FIG 1b FREQUENCY-DOMAIN %%%%%%
subplot(2,1,2)
loglog(f,P,[sym(j) '-'])
hold on
%loglog(f,P2,'ro')

if calib
	plot([testfreq testfreq],[1e8 1e10],'k--')
end

set(gca,'fontsize',18)
xlabel('frequency [Hz]')
ylabel('Power spectrum [W/Hz]')
axis([1e8 1e12 1e-20 1e-5])


%%%%%%%% FIGURE 2 Phase %%%%%%%%%%%%%%%
figure(2)
phase = atan2(imag(FD(1:N/2)), real(FD(1:N/2)));
semilogx(f,phase*180/pi,[sym(j) '-'])
hold on
set(gca,'fontsize',18)
ylabel('Phase [degrees]')
xlabel('frequency [Hz]')

end % j


%%% 
% extra labels
%%%%

figure(1)

subplot(2,1,1)
tt = r/c;
plot([tt tt],[0 2e9],'k--')
plot([tt/2 tt],[1e9 1e9],'k-')
text(tt/2,1.05e9,'r/c','fontsize',18)
% x-axis tic marks
for jj=-12:-9
	for ii=1:10
		x = 10^jj*ii;
		plot([x x],[0 0.002],'k-')
	end
end

subplot(2,1,2)
plot([1/tt 1/tt],[1e-20 1e-5],'k--')
plot([2/tt 1/tt],[1e-12 1e-12],'k-')
text(2/tt,1e-12,'c/r','fontsize',18)
% x-axis tic marks
for jj=8:11
	for ii=1:10
		x = 10^jj*ii;
		plot([x x],[1e-20 3e-20],'k-')
	end
end
% y-axis tic marks
for jj=-20:5
	y = 10^jj;
	plot([1e8 1.3e8],[y y],'k-')
	for ii=1:10
		y = 10^jj*ii;
		plot([1e8 1.1e8],[y y],'k-')
	end
end
	y = 10^0;
	plot([1e8 1.5e8],[y y],'k-','linewidth',2)
	y = 10^-5;
	plot([1e8 1.5e8],[y y],'k-','linewidth',2)
	y = 10^-10;
	plot([1e8 1.5e8],[y y],'k-','linewidth',2)

figure(2)
% x-axis tic marks
for jj=8:11
	for ii=1:10
		x = 10^jj*ii;
		plot([x x],[-200 -190],'k-')
	end
end

