Figure 1: (Top) Time-resolved Frt(t). (Bottom) Power spectrum P(f).
The MATLAB program "time2freq.m" illustrates the use of MATLAB's fast Fourier transform function, fft(), to convert time-resolved data into frequency-domain data.
An array Frt(N,1) holds N values of time-resolved fluence rate [W/cm2] generated by diffusion theory. An array t(N,1) holds N evenly spaced time points [s] associated with F(N,1). The observation point is located a distance r from the source. The source is an impulse of Uo [J] deposited at r = 0 at t = 0. The absorption coefficient (mua) and the diffusivity for light (alpha = c*D = c/(3*(mua + musp)) are specified by the user.
t = (0:dt:tmax)' + dt/2;
Frt = Uo*exp(-r^2./(4*alpha*t))./(4*pi*alpha*t).^(3/2).*exp(-mua*c*t);
N = length(Frt);
The array FD holds the fast fourier transform of Frt. Initially, a dummy array u holds the values fft(Frt), but only the lower half of this array is used. The upper half is symmetric with the lower half. Therefore, FD is calculated:
u = fft(Frt)
FD = u(1:N/2)
FD holds N/2 values that describe the fluence rate in the frequency domain. FD is an array with N/2 rows and 2 columns. The first column is the real component of FD, and the second column is the imaginary component of FD:
FD(j) = [Real(FD(j)), Imag(FD(j))]
The frequency f [Hz] associated with the jth row in FD(j,:) is given by the time step dt [s] and the original number of bins N:
f = 1/dt*(0:N/2-1)'/N;
The power spectrum [W/Hz] is obtained from FD by
P = (real(FD).^2 + imag(FD).^2 )/sum(P*df)
or alternatively by
P = (FD.*conj(FD))/sum(P*df)
where df is 1/dt/N or simply df = f(3) - f(2). The term sum(P*df) is a normalization factor such that after normalization has produced the final P, the sum(P*df) equals unity which is appropriate for a proper probability density function for power versus frequency.
Plots of Frt(t) and P(f) are shown in the figure below. The series of curves in different colors (r,g,b,r,g,b,m,c,k,m) denote increasing optical absorption coefficients (0.1:0.1:1) [cm-1]. The reduced scattering coefficient is 10 [cm-1].
Figure 1: (Top) Time-resolved Frt(t). (Bottom) Power spectrum P(f).
The above power spectrum is also called the "Magnitude" of FD. Although not obvious on a log-log plot, the area under each power spectrum curve is unity as is proper for a probability density function.
The "phase" of FD is calculated by
phase = atan2(imag(FD), real(FD))
A plot of phase(f) is shown in the figure below. The color scheme is for increasing mua as in Fig. 1.
Figure 2: Phase(f)
In both Figs. 1 and 2, the lowest absorption value (first red line) allows the Frt(t) to persist too long and when the array ends suddenly the Frt(t) has not yet dropped toward zero. Consequently, the discontinuity at this end of the array acts like an impulse that generates broad spectrum noise throughout the frequency spectrum and this noise is apparent as an obvious deviation at high frequencies. The other higher absorption values were sufficient to bring Frt(t) close enough to zero at the end of the array to avoid this artifact.