1. Introduction

Light scattering spectroscopy (LSS) is a measurement of the intensity of light scattered from tissue samples as a function of wavelength and/or angle [1]. Cell nuclei, cellular membranes, and mitochondria are all strong scatterers of light in the visible to near IR wavelength range. LSS is extremely sensitive to small changes in tissue architecture and has application for general tissue morphology, but of particular interest for cancer screening is the light scattered from the cell nuclei. LSS can detect small changes in the size and refractive index of cell nuclei. It has been shown previously that nuclei can be modeled as spherical Mie scatterers with characteristic diameters of 4 to 7 μm for healthy samples and diameters as large as 20 μm for dysplastic or cancerous samples [1].

Optical coherence tomography (OCT) is a tissue imaging modality that is capable of generating high resolution in vivo images of tissue based on low-coherence interferometric measurements [2]. OCT can be implemented in a conventional time-domain configuration using a low-coherence light source or in a frequency-domain configuration using a swept-wavelength laser or spectrometer. It has been demonstrated that the image contrast of OCT measurements can be enhanced through spectroscopic OCT in which the spectral content of the light is also determined [3, 4]. More recently, it has been demonstrated that LSS could be combined with spectroscopic OCT [5]. This combination enables coherence gating of the backscattered light, eliminating the diffuse scattering background that complicates LSS measurements. It has been demonstrated that the size of the scatterer can be determined from low-coherence interferometric measurements of the intensity of the scattered light as a function of wavelength [5,6], as well as from measurements of the scattering as a function of angle [7–9].

Conventionally, spectroscopic OCT has focused on measurement and analysis of the intensity of scattered light as a function of wavelength. In this paper we demonstrate that OCT measurements of the dispersion of the scattered light as a function of wavelength can also provide valuable information that is independent from the intensity. Dispersion is often viewed in terms of its deleterious effect on the resolution of OCT images. However, it has recently been proposed that measurements of dispersion could be used for clinical diagnostics, for example, to characterize plaque morphology for heart disease screening [10]. Additionally, it has been demonstrated that the phase velocity in a random medium is strongly influenced by scatterer size [11]. We demonstrate a spectral-domain OCT measurement technique that can determine the dispersion of a sample with high sensitivity. This technique could be applicable to the in vivo measurements of the sizes and refractive indices of scattering particles in biological tissue.

Tissue dispersion has two components: material dispersion, which arises from the wavelength dependence of the refractive index, and scattering dispersion, which results from the wavelength dependence of the scattering that typically occurs in tissue. Dispersion can be described by the group delay, which is defined as dϕ/dω, where ϕ is the phase of the electric field and ω is the angular frequency. We have developed Mie scattering analysis to predict the scattering group delay as a function of wavelength, simple phantoms to emulate the nuclei of cells, and a spectral domain OCT system capable of characterizing group delay as a function of wavelength with high resolution. In this paper, we show that the group delay of the scattered field is also strongly affected by the diameter and refractive index of the scatterer.

When compared with OCT measurements of the scattered intensity as a function of wavelength, spectral phase-dispersion measurements are important for the following reasons:

2. Mie scattering calculations

The complete solution to the scattering of an electromagnetic plane wave by a homogeneous sphere (absorbing or nonabsorbing) embedded in a homogeneous nonabsorbing infinite medium was found by Mie in 1908 [14–16]. We center our coordinate system at the center of the spherical scatterer. Assuming an exp(-i2πft) time dependence of the fields, and an incident plane wave traveling in the positive z-direction, the scattered electric fields, E _∥s and E _⊥s, in the far field are related to the incident fields by

where E _∥i and E _⊥i are components of the incident electric field, r is the distance of the observation point from the center of the sphere, f is the frequency, and k = 2πn / λ ₀ is the propagation constant. The subscripts ∥ and ⊥ signify the components parallel and perpendicular to the plane containing the incident and scattered propagation directions, respectively. The complex matrix elements S ₁(θ,f) and S ₂(θ,f) contain both the amplitude and the phase of the scattered light as a function of the scattering angle θ. Calculation of S ₁(θ, f) and S ₂(θ, f) is straightforward; public-domain software is available 17]. In the case of scattering from a dielectric sphere, their functional forms depend in a nontrivial way on the sphere diameter and the refractive indices of the sphere and surrounding medium (D, n_s and n_m, respectively). It is precisely this dependence that enables the “inversion” of spectral measurements for these physical parameters. For both the experimental results and computational analysis that follows, we restrict our results to the “backscattering” or retro-reflection direction, in which case S ₁(π, f) = S ₂(π, f) = S(f).

 

Fig. 1. Relative group delay predicted from Mie theory for different sphere diameters. The group delay of 3 mm of water is also shown for comparison.

Download Full Size | PDF

The group delay is calculated from the complex scattering function S(f) as follows:

where arg(S(f)) is the phase of the scattering function. Care must be taken when implementing Eq. (2) numerically, because the complex argument exhibits discontinuities (often referred to as phase wraps).

Figure 1 shows the predicted group delay as we vary the sphere’s diameter, assuming a single sphere with complex refractive index of n_s = 1.576 + i0.003 in a homogeneous, lossless medium with refractive index n_m = 1.346 + i0.000. We chose these refractive index values to be consistent with our measurements described below. In Fig. 1 we added an arbitrary constant to each of the group delay curves so that they are shifted with respect to each other to enhance the clarity of the plot. This is consistent with our measurements results, which include an arbitrary constant and are described as relative group delay (RGD). From Fig. 1 it is clear that the sphere size and the periodicity of the group delay are directly related and that the periodicity could potentially be applied to determine sphere size with high sensitivity. Also shown in Fig. 1 is the material dispersion of a 3 mm thick sample of water calculated from Schiebener’s formula [18]. The material dispersion of water is featureless compared with the scattering signatures, indicating that the material dispersion is not expected to have a significant effect on the results in measurements of water-based scattering samples, such as tissue.

3. Minimum and non-minimum phase response

In this section we closely examine the analytic character of the phase response of the scattered fields. We will demonstrate that the phase and the magnitude of the light scattered from a single sphere are independent. Using the language of linear response theory, we show that phase of the frequency response function of a spherical scatterer can exhibit both minimum and non-minimum phase characteristics. A minimum-phase function has the property that, in principle, its phase can be computed from measurements of its magnitude over a broad bandwidth. Therefore, when a linear response function is non-minimum phase, its phase and magnitude are independent quantities and both must be measured for a more complete characterization of the sample.

Through the detailed computations below, we show that even the simplest case of potential biological interest, back-scattering off an isolated sphere, exhibits multiple resonances in its group delay as a function of wavelength. Some of these resonances can be explained entirely by minimum phase computations. Others cannot. The underlying physics generating the two classes of response are not entirely understood. We intend to pursue this in the future as it is indicative of a fundamental distinction in the classification of electromagnetic resonance phenomena that would transcend the present application of OCT measurements of biological tissue. This question aside, the results below establish the independence of phase and intensity.

3.1 Linear system theory

The theory of linear systems is large and well-developed. We briefly present definitions and results below and refer to [19] for more detailed discussions. In the time domain, the response of a linear, time-invariant system to an arbitrary input is given by a convolution

where x(t) represents the system input, y(t) is the system output, and h(t) is the time-domain impulse response function. We limit our consideration to causal systems where, by definition, h(t) = 0 for t < 0, and the lower limit of the convolution integral above can be changed to 0. The system’s frequency response function H(f) is given by

In the frequency domain, the implications of causality are subtle but well known. The resulting statement is that a complex-valued function H(f) = H_r (f) + iH_i (f) may be continued as an analytic function for f in the upper half of the complex plane, C ⁺ = {z = f+iγ, γ > 0}. As a result of this analyticity, one derives the usual “Kramers-Kronig” (KK) relations between H_r (f) and H_i (f).

The derivation of the KK relations requires that the impulse response be linear, time-invariant, causal, and real. Although this may appear restrictive, in fact these properties are sufficiently generic that in essence they define the classical analysis of linear response functions. By contrast, minimum phase is both a more restrictive and subtle type of linear response function. Minimum phase is easiest to define in the frequency domain. We write H(f) = ρ(f)exp(iϕ(f)), where ρ(f) represents the magnitude and ϕ(f) is the phase of the electric field. Taking a logarithm, we have

If ln(H(f)) has the same analyticity properties as H(f) itself, then phase may be computed from magnitude from the appropriate KK relation [12,19],

The integral operator in Eq. (6) is singular at the “target” frequency, s = f, and is interpreted in a principal value sense. Functions whose phase and magnitude satisfy this relationship are called minimum phase.

3.2 Dielectric sphere

As in Sec. 2, we study the phase of the frequency response of a dielectric sphere in the back-scattering direction, S ₁(f, π) = S ₂(f, π) = S(f). The sphere is situated at the origin in a non-absorbing medium with refractive index n _m = 1.33. The diameter is D = 10 μm. In one case, the refractive index of the sphere is n _s = 1.39 + i0.003; in another case, it is n _s = 1.41 + i0.003. The imaginary part was chosen to match that of polystyrene microspheres.

Assuming an incident plane wave traveling in the positive z-direction, linearly polarized along the x-axis, the back-scattered electric field along the same axis was computed from Eq. (1) as discussed in Sec. 2. In Fig. 2 we show the logarithm of the magnitude of the response as a function of wavelength in nanometers, and we observe the characteristic resonances. The uniform shift of the resonances toward smaller wavelengths for the smaller refractive index is expected from the notion that a geometrically induced resonance should occur at invariant points of kD = 2πn_sD/λ.

We compute the minimum phase from the logarithm of the intensity using the singular integral operator in Eq. (6). The group delay of the true phase (calculated from the Mie analysis) and that of the minimum phase are plotted for both refractive indices in Fig. 3. Note that while the true and minimum phase group delays agree for the Re(n_s ) = 1.39 case, indicating that this case is indeed minimum phase, there are sign differences at λ = 1313 nm and 1373 nm for the Re(n_s ) = 1.41 case, indicating that its scattering function is non-minimum phase.

We found the appearance of non-minimum phase behavior shown in Fig. 3 to be somewhat surprising, so we went to great lengths to verify that this is indeed a non-minimum phase example and not a computational error. First, we ensured convergence of our computations to 10 digits, so that we were confident that we were not making an error in our numerical analysis. We were also careful to implement phase unwrapping correctly, but clearly phase wraps would appear as abrupt features, rather than the smooth curves of Fig. 3. Additionally, our most important verification came from a comparison of the deviations from minimum phase to a fundamental theorem that describes the structure of non-minimum phase functions [19]. This comparison is described in detail in Appendix A, but to summarize our findings, we found that the deviations from minimum phase shown in Fig. 3 are entirely consistent with the theoretical structure of non-minimum phase functions; therefore we conclude that this is truly an example of a non-minimum phase function rather than a computational error.

 

Fig. 2. Predicted magnitudes of the backscattered light from a single 10 μm sphere, calculated for two different refractive indices of the sphere.

Download Full Size | PDF

 

Fig. 3. Group delay of computed “true” phase (solid) and minimum phase (dashed). Note that for Re(n_s ) = 1.41, the true group delay resonances at λ = 1313 nm and 1373 nm have opposite sign compared to the resonances predicted assuming minimum phase.

Download Full Size | PDF

One important point comes from a comparison of the shape of the curves in Figs. 1 and 3. In Fig. 1, the RGD curves are approximately sinusoidal, while in Fig. 3 the curves are approximate delta functions. We found that the shape of the RGD curve is determined by the refractive index contrast between the sphere and the surrounding medium. The refractive indices (high index contrast) used to calculate Fig. 1 was chosen to match the materials used in our measurements of phantoms described below. The refractive indices (low index contrast) used to calculate Fig. 3 were chosen to be closer to what one might expect to see in actual tissue. The key point is that the light scattered from a single sphere is, in general, not minimum phase, and it is logical to expect that the more complex geometries present in tissue would be non-minimum phase as well.

Returning to the absorptive part of the refractive index of the sphere, Im(n_s ) = 0.003, we found that specific details of the phase response depend quite sensitively on this choice. In fact, it appears that a stable and sensitive procedure could be developed to invert direct measurements of scattered phase to yield this absorptive coefficient. It is not clear whether this would be of biological value, as most research to date involving the use of LSS-OCT as a medical diagnostic tool has focused on inversion of intensity scattering data for the cell nucleus size and real part of the refractive index. Performing several computations over a range of values from 0≤Im(n_s )≤0.006, we observed non-minimum phase resonances at some of the frequencies corresponding to resonances in the magnitude response. Where the non-minimum phase character occurs, and why, are topics for future study. Nevertheless, the persistence was clear. In general, given any arbitrary combination of one or more scatterers, our conclusion is that one must assume the sample is non-minimum phase unless specifically proven otherwise. Furthermore, if a sample is non-minimum phase, then both magnitude and phase carry independent information about the sample properties and both must be measured for a more full characterization of the sample.

4. Measurement

Our spectral domain OCT system is shown in Fig. 4. The system consists of two fiber-optic Michelson interferometers: a reference interferometer to track the wavelength of the tunable laser as it is swept, and a measurement interferometer. Our tunable laser has a nominal wavelength of 1300 nm, a 117 nm tuning range, and a maximum sweep speed of 80 nm/sec. The phantom (sample) is placed in one arm of the measurement interferometer and is coarsely aligned by monitoring the light scattered off the sample with the microscope’s IR camera. After the coarse alignment step, the sample alignment is optimized by monitoring the power of the backscattered light that is coupled back into the interferometer. After sample alignment, the laser is swept while the signal from the reference interferometer is sent to a zero-crossing circuit, which is used to trigger sampling of the measurement interferometer signal.

 

Fig. 4. Diagram of the OCT system used for dispersion metrology.

Download Full Size | PDF

Our measurement technique is similar to that used in previously demonstrated spectral domain measurements of the RGD of optical fiber components, particularly fiber Bragg gratings [20]. As we did in Sec. 3, we treat our sample as a linear filter, and we describe its complex spectral frequency response as H(f) = ρ(f)exp(iϕ(f)). We model the detected power at the output of the measurement interferometer as

where E₀ is the field incident on the interferometer, X₁ and X₂ are the optical path lengths of the two arms of the measurement interferometer, S(f) is the power spectrum of the source, and C₁ is a constant. We perform an inverse Fourier transform of this signal to obtain

where FT^-1 represents the inverse Fourier transform and h(t) is the analytic version of the impulse response of the sample [2]. We then use a window to remove the DC and autocorrelation terms [the first and second terms of Eq. (8)] and one of the complex conjugate impulse response terms. The one remaining term is the complex time-domain interferogram I(t). We can calculate the RGD from the following formula [21]:

where t_g is the relative group delay and Re represents the real part of the function. The wavelength resolution of the RGD result is determined by the width of the window. The optimal window width is a tradeoff between wavelength resolution, depth resolution, and RGD resolution; larger windows give better wavelength resolution but worse depth resolution. Also, larger windows include more noise, thereby degrading the RGD resolution.

From the shift theorem of Fourier transforms, we know that the function ϕ(f) includes a linear term, the slope of which is determined by the choice of the t = 0 point of the function I(t). For this reason, our group delay results include an additive arbitrary constant and are described as relative group delay.

Another important point to note is that we have configured our reference interferometer with a 1.5 m fiber length imbalance, which creates single-pass optical path difference (OPD) of approximately 2 m. This gives a high frequency reference signal, which requires a rapid data acquisition card (100 ksamples/s assuming the laser is swept at 40 nm/s). This results in significant oversampling of our frequency domain signal and a very large range of our time domain signal (corresponding to approximately 300,000 data points or 2 m of total depth range). Scattering and absorption limit the total penetration depth of light in tissue to less than 3 mm, so the additional range creates large regions of the time domain signal that are unused, creating the notion of “wasted data.” Practically, the disadvantage of oversampling is the computational overhead created by the large arrays. However, the advantages of oversampling include slight improvements in noise performance, reduced sensitivity to sampling errors, and reduced need for exact path matching of the measurement interferometer. In the future we will likely reduce our sample rate and array size in the interest of reducing computation time, as well as increase the laser sweep speed to be compatible with in vivo imaging; but for now our efforts are focused on proof of principle rather than measurement speed.

As a result of our large depth range, we can configure our measurement interferometer with a large OPD, eliminating any possibility of aliasing of the terms in Eq. (8) and simplifying our windowing process. However, we did approximately match the lengths of fiber in the sample and reference arms of our measurement interferometer so that the only dispersion in the interferometer was introduced by our sample.

A recent investigation has made the claim that, in a frequency-domain OCT measurement, the effective scattering function of tissue is either minimum phase or close to minimum phase, based on the intuitive notion that the time domain impulse response of any minimum phase function has a first or primary lobe that is dominant to all later responses and reflections [22]. In Eq. (8) the delta function could be thought of as a dominant first lobe. However, that delta function is windowed out in our analysis. Clearly the intent of the experiment is not to measure the reflectance of the mirror but rather the reflections generated by the sample. For scattering samples, such as tissue, these are inherently smaller than those generated by the mirror. Although it may be true that the combined response function (mirror and tissue sample reflections) is minimum phase, analysis of the combined response would relegate the signal of interest to small perturbations of the large and uninteresting reflection caused by the mirror. Although the analysis described in Ref. [22] may be useful in cases where the autocorrelation term of Eq. (8) aliases the other terms, it does not affect our derivation in Sec. 3 above, where we find that even a single spherical Mie scatterer is, in general, non-minimum phase.

We tested our measurement system by measuring the RGD of a 5 cm thick sample of well characterized commercial high-index glass. We compared our result with the change in RGD with wavelength predicted from the Sellmeier equation for that glass and found agreement better than 20 fs over a 45 nm wavelength range (representing the central 65 % of the laser tuning range, as measured using our older laser, which was capable of tuning only over a 70 nm range). We limited our wavelength range for this comparison because our RGD results are less accurate at the extreme ends of the tuning range, which may result from edge effects in the Fourier transform calculation or from uncertainties created by the laser tuning as it starts and finishes a wavelength sweep, combined with the limited frequency response of the zero-crossing trigger circuit.

We developed phantoms by embedding a low density (approximately 8 spheres per mm³) of small polystyrene spheres in a 20 % by weight glycerin/water solution. We placed a small drop of the glycerin/water/spheres mixture on a microscope slide and topped it with a cover slip. The glycerin/water solution matches the density of the polystyrene spheres, providing a neutral buoyancy for the spheres in solution [23]. The real part of the refractive index of polystyrene at a wavelength of 1300 nm is given by Re(n_s ) = 1.576 ± 0.004, while the absorptive part is given by Im(n_s ) = 0.003 ± 0.001 [24].

We both calculated and measured the refractive index of the glycerin/water solution. We calculated the refractive index from the following formula:

where n_m is the refractive index, M is the molecular weight, and ρ is the density of the composite solution, while A and B are the two composite materials with mole fractions f_A and f_B , molecular weights M_A and M_B , and densities ρ_A and ρ_B [25]. We used the refractive indices of water [26] and glycerin [27] at a wavelength of 1300 nm in the calculation. From this, we obtained a refractive index of 1.346, compared with the 1.357 published value for the glycerin/water solution [23]; the difference between the two values might be attributed to the change in refractive index with wavelength. We also determined the refractive index from an OCT measurement of the optical path difference of a 1 mm thick cuvette filled with the glycerin/water solution. We obtained a measured refractive index of 1.38 ± 0.03.

We measured our phantoms by first collimating the light from the sample arm of our OCT system to a diameter of 1.3 mm. We then focused on the sample using a 10X objective, giving a focal spot size of approximately 25 μm with an effective numerical aperture of 0.04. The microscope slide was tilted with respect to the incident beam to avoid specular reflections from the slide. We aligned the sample as described above. We then measured the frequency-domain interference and calculated the RGD of each sphere using the calculations described above. We applied a Hamming window to the time domain interferogram corresponding to a truncation length of 240 μm. This gives a wavelength resolution in the RGD of approximately 7 nm. Figure 5 shows the RGD determined from a single, unaveraged measurement of four different sizes of spheres in the glycerin/water solution.

 

Fig. 5. Comparison of measured and predicted group delay of phantoms constructed from polystyrene spheres embedded in a glycerin/water solution. We added arbitrary constants to each of the curves to separate them.

Download Full Size | PDF

From Fig. 5, we see good qualitative agreement between theory and experiment (i.e., the periodicity of the RGD oscillations matches that predicted by theory). We think that the main factor inhibiting exact agreement of our measured and predicted RGD is the uncertainty associated with independent measurements of sphere size and refractive indices, combined with the fact that the Mie prediction is very sensitive to small changes in these numbers. For example, changing the surrounding refractive index from our predicted value of 1.346 to the measured value of 1.38 dramatically affects the locations of the RGD peaks. The Mie scattering predictions are also strongly affected by small changes in sphere diameter. We attempted to obtain accurate and independent measurements of the diameter of each sphere using a microscope with a reticule eyepiece, but we discovered that the uncertainty of that measurement is no better than the uncertainty of sphere size specifications provided by the sphere manufacturer (5.0 ± 0.5 μm, 8.0 ± 0.8 μm, 10.0 ± 1.0 μm, and 14.5 ± 1.0 μm). Intensity-scattering measurements exhibit a similar sensitivity to small changes in diameter and refractive index, and researchers have addressed this by determining sphere diameter and/or refractive index directly from the Mie scattering measurements. This can either be done by performing a Fourier transform on the measured data to determine the principal frequency [5], or by calculating the χ ² difference between measurement and a database of Mie predictions as a function of diameter and refractive index [9]. Other factors affecting our measurement uncertainty include the wavelength scale of our RGD results, which may include an uncertainty as large as 1 nm. Also, the Mie predictions were calculated assuming plane wave illumination and direct backscattering, whereas our measurements were made with a focused beam to increase the total backscattered power collected. However, for the case of intensity-scattering OCT measurements, it has been demonstrated that the Mie predictions are not strongly affected by focusing, provided that the beam waist is large compared with the wavelength [6], as is the case in our measurements.

A detailed analysis of the uncertainty of the scatterer sizes determined from phase dispersion measurements is beyond the scope of this paper. However, we can say that, when calculating the diameter of the scatterer from a Fourier transform of the group delay, a key factor in the uncertainty is the resolution. Given that our laser can only tune over a 117 nm wavelength range, the resolution of a Fourier-transform sizing calculation is 5 μm (assuming that the scatterer has a refractive index of 1.4). We are currently developing a system with a broader wavelength range to overcome this resolution limit.

As we mentioned above, the group delay can appear sinusoidal or as a combination of delta functions, depending on the ratio of the refractive index of the sphere to that of the surrounding medium. In Fig. 6, we show an experimental result where the group delay takes the form of a combination of approximate delta functions. This measurement was obtained from an 18 μm polystyrene sphere embedded in a porcine gelatin (10 % by weight gelatin in water), under the same focusing, alignment, windowing, and calculation conditions as described above. In Fig. 6 there is an obvious sign reversal between the peaks, possibly indicating that this is an example of a non-minimum phase function, as described in Sec. 3. We were unable to find published data for the refractive index of 10 % by weight gelatin, but given that a 1 % by weight gelatin sample has a refractive index of 1.5118 [28], we believe that our polystyrene sphere in a 10 % by weight gelatin sample would likely be a low-index contrast example.

 

Fig. 6. Measured RGD of an 18 μm sphere embedded in porcine gelatin. The curve has the shape of approximate delta functions, and there is a clear sign reversal between peaks, qualitatively similar to that predicted by the non-minimum phase analysis.

Download Full Size | PDF

We limited our considerations to single scatterers in the above examples to better facilitate comparisons of predicted and measured data. However, given the multitude of scatterers that exist within a biological tissue sample, it is critical to consider situations under which more than one scatterer is simultaneously illuminated. In Fig. 7 we show measured group delay from a single 15 μm polystyrene sphere on a microscope slide surrounded by air. This sample was aligned and measured as described above. Also shown in Fig. 7 is a measurement of the group delay from a monolayer of 15 μm polystyrene spheres on a slide. The monolayer was illuminated with a collimated beam (1.3 mm diameter) normally incident on the rear side of the microscope slide. The RGD was calculated from the interference of the light scattered from the monolayer with of the rear surface reflection from the microscope slide. These results are preliminary, but the important conclusion of Fig. 7 is that, although the group delay amplitude is reduced in the case of multiple scatterers, the periodicity of the group delay is unchanged.

 

Fig. 7. Comparison of the RGD measured from a single 15 μm sphere measured with a focused beam and the RGD from a monolayer of 15 μm spheres measured with a 1.3 mm diameter collimated beam. The sphere monolayer data was multiplied by a factor of 5 so that both curves could be shown on the same graph. Although amplitudes are different and the curves are phase shifted with respect to each other, both curves appear to have the same period.

Download Full Size | PDF

5. Conclusions

To measure the size of optical Mie scatterers, we have demonstrated a new technique that is based on OCT measurements of the group delay of the scattered light. This technique is potentially very sensitive to small differences in the size and refractive index of scatterers. Additionally, we have demonstrated that a single spherical scatterer is a non-minimum phase function, which means that the phase spectrum is independent of the intensity spectrum of the scattered light. Therefore, the group delay provides information about the scatterer that is not available from a measurement of intensity alone. Given that single sphere scattering is non-minimum phase, we believe it is likely that the more complex structures present in tissue would be non-minimum phase as well. Measurements of group delay, either alone or in combination with intensity measurements, may have application to general tissue morphology measurements and the in vivo detection of precancerous dysplasia. Our future efforts will focus on extending the scattering predictions, minimum phase analysis, and our measurements to include the interactions between multiple scatterers. Additionally, we plan to perform measurements on more complex tissue phantoms and real tissue samples.

Appendix A

In this section, we closely examine the properties of minimum phase functions with the goal of verifying our conclusion that our example of single sphere scattering shown in Fig. 3 is indeed non-minimum phase. We start with Eq. (6), which follows from analyticity, thus minimum phase is equivalent to the analyticity of ln[H(z)] in C ⁺. As H(z) is analytic in this half-plane, one obstruction to H(z) being minimum phase would be a point z_n ∈ C ⁺ such that H(z_n ) = 0, as ln[H(z)] would then have a branch point at z = z_n = f_n + iγ_n .

One example of a non-minimum phase function is the all-pass, or Blaschke factor, defined as

(A1)$$B\left(f\right)=\frac{f-z_n}{f-z_n^{*}},$$

where z_n ∈ C ⁺. We readily verify that B is analytic in the upper half plane, it has a zero there, and that |B(f)| = 1 for all real f. Forming the logarithm we find (see below for choice of branch for the inverse trig function)

(A2)$$\ln \left(B\left(f\right)\right)=\ln \left(\mid B\left(f\right)\mid \right)+iArg\left(B\left(f\right)\right)=i2{\cos }^{-1}\left(\frac{f-f_n}{\sqrt{{\left(f-f_n\right)}^2+{\gamma }_n^2}}\right).$$

We find that the real part of ln(B(f)) is zero, yet the function has a nontrivial phase; therefore, as expected, B(f) is not a minimum phase function.

A fundamental theorem states that, assuming causality, any arbitrary frequency response function can be expressed as a product of all-pass factors and a minimum phase function as follows [19]:

(A3)$$H\left(f\right)=H_{\mathit{mp}}\left(f\right)\exp \left(i2\pi f\tau \right)\prod _n\frac{f-z_n}{f-z_n^{*}}.$$

In Eq. (A3), H_mp (f) is minimum phase, τ is real, and the product of all-pass terms is over all of the upper half plane zeros of H(z). We note that |H_mp (f)| = |H(f)| regardless of whether H(f) is or is not minimum phase. If H(f) itself is minimum phase, then the exponential and all-pass terms in Eq. (A3) are identically equal to 1, i.e., τ = 0 and {z_n } = ∅. The point of this discussion is to show that, using Eq. (6), one may compute any system’s minimum phase response ϕ_mp (f) from a measurement of the magnitude of the system frequency response. If the true phase of the frequency response,ϕ_E (f), is known independently, then the difference (if any) necessarily can be expressed as

(A4)$$\Delta \left(f\right)={\varphi }_E(f)-{\varphi }_{mp}\left(f\right)=2\pi f\tau +\sum _{n}2{\cos }^{-1}\left(\frac{f-f_n}{\sqrt{{\left(f-f_n\right)}^2+{\gamma }_n^2}}\right).$$

Thus the factorization given in Eq. (A3) restricts the theoretical structure of deviations from minimum phase.

One key technical point: the time-domain impulse response functions of interest are always real. Parity considerations of standard Fourier analysis then imply that nontrivial all-pass terms occur in pairs; i.e., if the product in Eq. (A3) contains a term with a zero at z_n = f_n + iγ_n , then it also has a term with a zero at z_n′ = -$z_n^{*}$ = -f_n + iγ_n. In this analysis, the symmetry is accounted for implicitly, and we discuss only the terms with resonances at positive real parts. Additionally, the choice of cos^-1 as the inverse trig function in Eq. (A2) has the virtue that it unwraps the phase as the frequency f sweeps over the resonances f_n , thus it may be differentiated with respect to f[see Eq. (A5)]. However, when summing over multiple terms corresponding to the individual factors of the product in Eq. (A3), we should subtract 2πm from the result where m is the number of paired factors to maintain the Fourier parity. This constant term shifts the overall branch but does not affect the group delay results.

We can now test the Re(n_s ) = 1.41 case of Fig. 3 to verify our conclusion that it is a non-minimum phase function as follows. We found that it is easier to express the deviations from minimum phase in terms of group delay rather than phase. Differentiating Eq. (A2), we find that the group delay of an all-pass term is given by

(A5)$$t_g=\frac{1}{2\pi }\frac{d{\varphi }_B}{df}=\frac{{\gamma }_n}{\pi \left({\left(f-f_n\right)}^2+{\gamma }_n^2\right)}.$$

Equation (A5) is recognizable as an “approximate delta function” centered at f_n with a height 1/(πγ_n ). We differentiate Eq. (A4), so that we can express the deviations from minimum phase in terms of group delay rather than phase, giving

(A6)$$\Delta t_g\left(f\right)=\frac{1}{2\pi }\left(\frac{d{\varphi }_E}{df}-\frac{d{\varphi }_{\mathit{mp}}}{df}\right)=\tau +\frac{1}{\pi }\sum _n\frac{{\gamma }_n}{\left({\left(f-f_n\right)}^2+{\gamma }_n^2\right)}.$$

For the Re(n_s ) = 1.41 case of Fig. 3, we performed a nonlinear curve fit to the difference between the RGD computed from the Mie calculations and the minimum phase RGD for two “all-pass” zeros as well as a constant offset, i.e., two frequencies f_n and heights γ_n , and a non-trivial τ. For the Re(n_s ) = 1.39 case, we fit the same difference for τ alone. In Fig. 8 we plot the difference between the true RGD and the minimum phase RGD after subtracting out these curve fits. Note that the resonances of Fig. 3 are removed entirely, and the residual is smooth. The deviations at the minimum and maximum wavelengths are created by the truncation of the integral in Eq. (6) [29]. We conclude that the discrepancy between the Mie calculation and minimum phase results is entirely explained by a fundamental theorem that describes the shape of non-minimum phase functions [Eq. (A3)], therefore the Re(n_s ) = 1.41 case truly is non-minimum phase.

 

Fig. 8. Plot of the difference between true and minimum phase group delays after subtracting the all-pass and linear terms determined from a nonlinear fit. This is consistent with Eq. (A3), and indicates that the Re(n_s ) = 1.41 case is truly an example of a non-minimum phase scattering function. The dashed vertical lines show the positions of the resonances for reference. The small features at λ = 1237 nm and 1294 nm are numerical artifacts.

Download Full Size | PDF

Acknowledgment

The authors would like to thank Tom Milner of the University of Texas at Austin for many helpful discussions.

References and links

1 . L.T. Perelman , V. Backman , M. Wallace , G. Zonios , R. Manoharan , A. Nusrat , S. Shields , M. Seiler , C. Lima , T. Hamano , I. Itzkan , J. Van Dam , J. M. Crawford , and M. S. Feld , “ Observation of periodic fine structure in reflectance from biological tissue: a new technique for measuring nuclear size distribution ,” Phys. Rev. Lett. 80 , 627 – 630 ( 1998 ). [CrossRef]  

2 . A. F. Fercher , W. Drexler , C. K. Hitzenberger , and T. Lasser , “ Optical coherence tomography - principles and applications ,” Rep. Prog. Phys. 66 , 239 – 303 ( 2003 ). [CrossRef]  

3 . U. Morgner , W. Drexler , F. X. Kärtner , X. D. Li , C. Pitris , E. P. Ippen , and J.G. Fujimoto , “ Spectroscopic optical coherence tomography ,” Opt. Lett. 25 , 111 – 113 ( 2000 ). [CrossRef]  

4 . D. C. Adler , T. H. Ko , P. R. Herz , and J. G. Fujimoto , “ Optical coherence tomography contrast enhancement using spectroscopic analysis with spectral autocorrelation ,” Opt. Express 12 , 5487 – 5501 ( 2004 ). [CrossRef]   [PubMed]  

5 . A. Wax , C. Yang , and J. A. Izatt , “ Fourier-domain low-coherence interferometry for light-scattering spectroscopy ,” Opt. Lett. 28 , 1230 – 1232 ( 2003 ). [CrossRef]   [PubMed]  

6 . C. Xu , P. S. Carney , and S. A. Boppart , “ Wavelength-dependent scattering in spectroscopic optical coherence tomography ,” Opt. Express 13 , 5450 – 5462 ( 2005 ). [CrossRef]   [PubMed]  

7 . A. Wax , C. H. Yang , R. R. Dasari , and M. S. Feld , “ Measurement of angular distributions by use of low-coherence interferometry for light-scattering spectroscopy ,” Opt. Lett. 26 , 322 – 324 ( 2001 ). [CrossRef]  

8 . A. Wax , C. Yang , V. Backman , K. Badizadegan , C. W. Boone , R. R. Dasari , and M. S. Feld , “ Cellular organization and substructure measured using angle-resolved low-coherence interferometry ,” Biophys. J. 82 , 2256 – 2264 ( 2002 ). [CrossRef]   [PubMed]  

9 . J. W. Pyhtila , R. N. Graf , and A. Wax , “ Determining nuclear morphology using an improved angle-resolved low coherence interferometry system ,” Opt. Express 11 , 3473 – 3484 ( 2003 ). [CrossRef]   [PubMed]  

10 . B. Liu , E. A. Macdonald , D. L. Stamper , and M. E. Brezinski , “ Group velocity dispersion effects with water and lipid in 1.3 μm optical coherence tomography ,” Phys. Med. Biol. 49 , 923 – 930 ( 2004 ). [CrossRef]   [PubMed]  

11 . C. Yang , A. Wax , and M. S. Feld , “ Measurement of the anomalous phase velocity of ballistic light in a random medium by use of a novel interferometer ,” Opt. Lett. 26 , 235 – 237 ( 2001 ). [CrossRef]  

12 . M. Beck , I. A. Walmsley , and J. D. Kafka , “ Group delay measurements of optical components near 800 nm ,” IEEE J. Quantum Electron. 27 , 2074 – 2081 ( 1991 ). [CrossRef]  

13 . G. Lenz , B. J. Eggleton , C. R. Giles , C. K. Madsen , and R. E. Slusher , “ Dispersive properties of optical filters for WDM systems ,” IEEE J. Quantum Electron. 34 , 1390 – 1402 ( 1998 ). [CrossRef]  

14 . G. Mie , “ Beitrage zur Optik trüber Medien speziell kolloidaler Metallösungen ,” Ann. Phys. 25 , 377 – 445 ( 1908 ). [CrossRef]  

15 . C. F. Bohren and D. R. Huffman , Absorption and Scattering of Light by Small Particles ( John Wiley & Sons, New York , 1983 ), Chap. 4.

16 . H. C. van de Hulst , Light Scattering by Small Particles ( Dover, New York , 1981 ), Chap. 9.

17 . T. A. Germer , SCATMECH: Polarized Light Scattering C+ + Class Library , available at http://physics.nist.gov/scatmech .

18 . P. Schiebener , J. Straub , J. M. H. L. Sengers , and J. S. Gallagher , “ Refractive-index of water and steam as function of wavelength, temperature, and density ,” J. Phys. Chem. Reference Data 19 , 677 – 717 ( 1990 ). [CrossRef]  

19 . H. Dym and H. P. McKean , Fourier Series and Integrals ( Academic Press, New York , 1972 ), Chap. 3.5.

20 . M. E. Froggatt , T. Erdogan , J. Moore , and S. Shenk , “ Optical frequency domain characterization (OFDC) of dispersion in optical fiber Bragg gratings ,” in Bragg Gratings, Photosensitivity, and Poling in Glass Waveguides , OSA Technical Digest ( Optical Society of America, Washington, DC , 1999 ), pp. 227 – 229 .

21 . V. Laude , “ Noise analysis of the measurement of group delay in Fourier white-light interferometric cross correlation ,” J. Opt. Soc. Am. B 19 , 1001 – 1008 ( 2002 ). [CrossRef]  

22 . A. Ozcan , M. J. F. Digonnet , and G. S. Kino , “ Frequency-domain optical coherence tomography based on minimum-phase functions ,” in Coherence Domain Optical Methods and Optical Coherence Tomography in Biomedicine X , V. V. Tuchin , J. A. Izatt , and J. G. Fujimoto , eds., Proc. SPIE 6079 , 607912 ( 2006 ). [CrossRef]  

23 . http://www.dow.com/glycerine/resources/dwnlit.htm http://www.2spi.com/catalog/standards/microspheres.shtml .

24 . X. Ma , J. A. Lu , R. S. Brock , K. M. Jacobs , P. Yang , and X.-H. Hu , “ Determination of complex refractive index of polystyrene microspheres from 370 to 1610 nm ,” Phys. Med. Biol. 48 , 4165 – 4172 ( 2003 ). [CrossRef]  

25 . M. Born and E. Wolf , Principles of Optics ( Pergamon Press , 1975 ), Chap. 2.3.3.

26 . J. E. Bertie and Z. Lan , “ Infrared intensities of liquids XX: the intensity of the OH stretching band of liquid water revisited, and the best current values of the optical constants of H ₂ 0(l) at 25 °C between 15,000 and 1 cm ^-1 ,” Appl. Spectrosc. 50 , 1047 – 1057 ( 1996 ). [CrossRef]  

27 . P. R. Cooper , “ Refractive-index measurements of liquids used in conjunction with optical fibers ,” Appl. Opt. 22 , 3070 – 3072 ( 1983 ). [CrossRef]   [PubMed]  

28 . K. Nikolova , I. Panchev , and S. Sainov , “ Optical characteristics of biopolymer films from pectin and gelatin ,” J. Optoelectronics Adv. Mater. 7 , 1439 – 1444 ( 2005 ).

29 . A. Dienstfrey , P. D. Hale , D. A. Keenan , T. S. Clement , and D. F. Williams , “ Minimum-phase calibration of sampling oscilloscopes ,” IEEE Trans. Microwave Theory Tech. (to be published).