Dispersion compensation in Fourier domain optical coherence tomography using the fractional Fourier transform

Norman Lippok; Stéphane Coen; Poul Nielsen; Frédérique Vanholsbeeck

doi:10.1364/OE.20.023398

1. Introduction

Developed in the early 1990’s, optical coherence tomography (OCT) is an established and successful imaging technology that enables high resolution, cross-sectional imaging of biological tissues and materials [1, 2]. Nowadays, it is commonly implemented in the frequency domain with so-called Fourier-domain OCT (FD OCT), which has significantly improved sensitivity and imaging speed [3–5]. Frequency domain systems can be implemented either by using a low coherent light source and a spectrometer (SD OCT) [6–8], or by using a coherent swept-source laser (SS OCT) that is sampled by a photo diode [9–12].

OCT imaging has proven to be a powerful diagnostic tool in many medical fields. In some systems, OCT tomograms have image quality comparable to histology with resolutions down to the 1 μm level [13]. Obtaining high isotropic resolution over a large imaging depth range is however hindered by two factors. First, the high NA objectives traditionally used to obtain high lateral resolution limit the depth of field and hence the imaging depth. This can be circumvented by using axicons [14–16], binary-phase filters [17], or multi-beam OCT [18, 19]. The second problem arises from the broadband nature of the light source. Increasing the axial resolution requires larger optical bandwidths, which make OCT systems more sensitive to chromatic dispersion, especially for a long sensing range. For example, using a center wavelength around 800nm, imaging the human macula with an axial resolution below 15 μm makes dispersion compensation necessary in order to account for the dispersion produced by the vitreous body of the human eye [20, 21]. In ophthalmology, it is therefore advantageous to operate close to the zero dispersion wavelength of water, λ₀ = 1 μm [22].

With this in mind, it is no surprise that dispersion compensation methods have become increasingly important for OCT. Traditional methods rely on placing the right amount of dispersion balancing material in one interferometer arm of the OCT setup [23–25], but this is usually only practical for 2nd-order dispersion. Grating-based phase delay scanners [26] and dual optical fiber stretchers [27] can also be used for 2nd order dispersion compensation and present some degree of tunability. However, these approaches require bulky components and the latter is difficult to adapt for depth-dependent compensation of sample dispersion. We note that a fiber-stretching-based dispersion compensator has been recently combined with a grating-based, scanning free, time domain OCT system to compensate for both 2nd and 3rd-order dispersion [28], but to the detriment of added complexity.

Because of these drawbacks, OCT systems increasingly rely on numerical dispersion compensation techniques, which offer continuous adjustment capabilities and, in theory, can be optimized for any amount of dispersion. This was first demonstrated in time-domain systems [21, 29–31]. It is however more readily implemented with FD OCT systems which have direct access to the phase information of the signal as was originally demonstrated almost simultaneously by two groups, Cense et al. [32] and Wojtkowski et al. [33]. Their method is based on the complex conjugate of the dispersive spectral phase. It compensates for 2nd and higher-order dispersion by adding an optimized phase term to the analytical expression of the measured spectral fringe signal. The phase term is found using a sharpness function that searches for maximum signal magnitude in a given depth range. Note that FD OCT systems still require nearly dispersion-matched interferometer arms: too much dispersion spreads the point-spread-function (PSF) over a wide range of depths, leading to sensitivity decay and information loss, even after numerical dispersion compensation.

In our work, we revisit the problem of numerical dispersion compensation based on the use of the fractional Fourier transform (FrFT). In doing so, our aim is not to provide a replacement to the dispersive spectral phase compensation technique discussed above and which is rapidly becoming the de-facto standard of numerical dispersion compensation. Rather, we wish to illustrate how the FrFT provides a new fundamental perspective on the nature and role of group-velocity dispersion in Fourier domain OCT and how it highlights in a visual manner the physics behind dispersion compensation. Our work also provides new general insights into advanced problems associated with dispersion compensation such as defining a sharpness function that is not dependent on the presence of an isolated single back-scatterer, the issue of depth-dependent dispersion, or the prospect of differentiating materials by dispersion mapping [34, 35]. In the following, we first briefly present the theory of the FrFT and then demonstrate numerical dispersion compensation using FrFT experimentally. Furthermore, we show theoretically how we can extend our method, and use the FrFT for depth-dependent sample dispersion compensation. Finally, we show how we can extract the 2nd-order group-velocity dispersion coefficient β₂ of a sample from the order parameter of the fractional Fourier transform. All of our techniques are general and can be applied to SD as well as SS OCT, but also more generally to interferometry or optical coherence domain reflectometry.

2. Theory

The fractional Fourier transform (FrFT) is a generalization of the traditional Fourier transform (FT), with the addition of an order parameter a that can be interpreted as the FT to the power of a. The traditional FT and the associated frequency domain are special cases of the FrFT and its continuum of fractional Fourier domains, which are elegantly related to phase-space distributions as will be discussed below [36]. The FrFT was initially defined by Namias in the context of quantum mechanics [37] but has then been used in optics [38] and signal processing [39]. For various problems in which the Fourier transform (FT) is used, there is a potential for generalization and improvement by using the FrFT. By replacing the traditional Fourier transform (FT) with the FrFT, we gain an additional degree of freedom to a problem within a phase-space distribution (time-frequency distribution).

The a-th order FrFT is a linear transform. For the range 0 < |a| < 2 it is defined as [40]

f_{a} (u_{a}) = F^{a} [f (u)] = \int_{- \infty}^{\infty} A_{ϕ} exp [i π [u_{a}^{2} cot (ϕ) - 2 u u_{a} csc (ϕ) + u^{2} cot (ϕ)]] f (u) d u,

where

A_{ϕ} = \frac{exp {- i π \frac{sign [sin (ϕ)]}{4} + \frac{i ϕ}{2}}}{{| sin (ϕ) |}^{1 / 2}} and ϕ = \frac{a π}{2} .

Here u is the independent variable of the transform input function while u_a is that of the transform output, with corresponding fractional transform order a. These independent variables are assumed to be dimensionless (see below).

When a = 1, A_ϕ becomes unity, the first and third terms in the argument of the exponential in Eq. (1) vanish, and the FrFT reduces to the traditional FT (i.e., the FrFT to the power of 1). Therefore, the u₁ and u axes represent the (normalized) time τ and optical frequency ν axes, respectively (the normalization is such that uu₁ = ντ). Note that, in our work, these domains are swapped compared to traditional terminology because we deal here with FD OCT, i.e., the measurements f (u) are performed in the frequency domain.

For a → 0 and a → ±2 the integral kernel approaches δ(u_a − u) and δ(u_a + u), sampling f(u) as the identity and parity operators, respectively [40]. Only for a = 1 does the transform of a real function yield mirrored counterparts in the two halves of the Fourier space. Otherwise the transform yields different energy distributions in the two (fractional) Fourier half spaces. Except for the special cases a = 0 and a = 1, the transform output lies in neither the traditional frequency nor time domain.

The FrFT is related to a more general time-frequency representation of the function f(u), namely its Wigner distribution W_f (u, u₁) [36]. Recall that time-frequency distributions aim to represent how the spectral content of a signal change with time, or, in other words, how the energy of a signal is distributed both in time and frequency. Compare this with the traditional FT, which provides no information as to when specific frequency components occur in the signal. The Wigner distribution is one of the most commonly used time-frequency distributions and it is found that the Wigner distribution of f_a is simply a rotated version of that of f [40]. More specifically, W_{f_a} is W_f as seen in a reference frame (u_a, u_a+1) rotated by an angle ϕ= aπ/2 with respect to the time-frequency plane (u, u₁),

W_{f_{a}} (u_{a}, u_{a + 1}) = W_{f} (u_{a} cos ϕ - u_{a + 1} sin ϕ, u_{a} sin ϕ + u_{a + 1} cos ϕ) .

From the integral projection properties of the Wigner distribution,

\int_{- \infty}^{\infty} W_{f} (u, u_{1}) d u_{1} = {| f (u) |}^{2}

, it follows that the integral projection of the Wigner distribution W_f onto an axis u_a making an angle ϕ = aπ/2 with the u axis, yields the squared amplitude of the fractional Fourier transform of order a,

\int_{- \infty}^{\infty} W_{f_{a}} (u_{a}, u_{a + 1}) d u_{a + 1} = {| f_{a} (u_{a}) |}^{2} .

A projection at 90° (a = 1) corresponds to the traditional FT, and is consistent with the other marginal of the Wigner distribution,

\int_{- \infty}^{\infty} W_{f} (u, u_{1}) d u_{1} = {| f_{1} (u_{1}) |}^{2}

.

To illustrate these concepts, we consider a simple linear chirp signal

S (ω) = Re {I (ω) exp [j (ω τ + β_{2} l {(ω - ω_{c})}^{2})]} .

In an FD OCT setup, such a spectral signal would be observed in the presence of residual chromatic dispersion from a single reflection occurring at a delay τ from the point of zero-path difference of the interferometer, with β₂ the group-velocity dispersion coefficient of the dispersive element of optical thickness l. Here ω = 2πν is the angular optical frequency, and ω_c the center frequency of the source spectrum I(ω). After normalization, u = ν/η, u₁ = ητ, Eq. (5) reads

S (u) = Re {I (u) exp [j (2 π u u_{1} + {(2 π η)}^{2} β_{2} l {(u - u_{c})}^{2})]} = Re {I (u) exp [j φ]} .

For digital computation, the normalization factor η is set as [40],

η = \sqrt{\frac{Δ ν}{Δ τ}} = d ν \sqrt{N},

where Δν and Δτ are the width of the spectral and temporal intervals over which our signal is represented, respectively, dν = Δν/N = 1/Δτ is the spectral resolution, and N is the number of sampling points. With this scaling, the normalized length of both intervals are equal to the dimensionless quantity

\sqrt{Δ ν Δ τ} = \sqrt{N}

, the Wigner distribution is confined to a circle, and the samples in both domains are spaced

1 / \sqrt{N}

apart.

The Wigner distribution of the chirped signal S(u), Eq. (6), is shown schematically in Fig. 1, which can be understood as follows. The second-order dispersion term (β₂) leads to a linear frequency modulation of the detected spectral fringes S(u) as a function of optical frequency, which can be related to the inclination of the Wigner distribution in the time-frequency plane. Note that, because we are dealing with spectral fringes, the instantaneous “frequency” of the fringes has unit of inverse optical frequency, i.e., time, hence appears as the u₁ (i.e., temporal) domain. This instantaneous frequency can be obtained by taking the derivative of the phase of the complex spectral fringes with respect to angular optical frequency, hence is given by

\frac{1}{2 π} \frac{\partial φ}{\partial u} = u_{1} + {(2 π η)}^{2} \frac{β_{2} l}{π} (u - u_{c}),

which highlights the linear frequency modulation. As explained above, Eq. (4), the squared amplitude of the a-th order FrFT can be obtained as the integral projection of the Wigner distribution onto the u_a axis making an angle ϕ = aπ/2 with the u axis. In Fig. 1, we have chosen the order a and the u_a axis orientation in a way that the integral projection of our chirped signal leads to a maximally narrowed response (see green shaded area along the u_a axis). This illustrates that the FrFT, with an optimized order parameter a, can effectively lead to dispersion compensation of an FD OCT signal. In contrast, projection onto the u₁ axis, i.e., the traditional FT of the signal, leads to a much broader response (red shaded area along the u₁ axis), because dispersion is uncompensated.

Fig. 1 Wigner distribution of a chirped signal in the time-frequency plane (u₁, u). It illustrates how 2nd-order dispersion affects a spectral interferogram and how the continuum of domains of the FrFT can be interpreted geometrically by the rotated frame (u_a, u_a₊₁). The order parameter a chosen for the plot is optimal for the signal considered here, i.e., it leads to a dispersion-compensated OCT depth signal after FrFT.

Download Full Size | PDF

In essence, the FrFT order a can be chosen to adjust the chirp rate of a dispersed OCT signal. This amounts to a rotation of the time-frequency plane. In the optimally rotated frame, the instantaneous frequency of the spectral fringes remains constant. The fractional Fourier transformed OCT signal only results from a different projection angle in the time-frequency distribution compared to the standard FT, while entirely preserving the energy of the original spectrum. With this in mind, one could state that group-velocity dispersion does not degrade the axial resolution in FD OCT but only causes one to observe an OCT depth signal from a perspective in the time-frequency plane that is not suitable for imaging. The FrFT allows one to correct for that perspective.

The optimized FrFT order parameter a that leads to a dispersion-compensated OCT depth signal can be interpreted as an intuitive measure of the amount of chromatic dispersion present in an FD OCT system. If the optimized order is a = 1, the OCT A-scan is dispersion free. Normal and anomalous dispersion can easily be distinguished by a > 1 and a < 1, respectively.

3. Experimental setup

To experimentally demonstrate numerical dispersion compensation using the FrFT, we have used an SD OCT configuration as shown in Fig. 2. It is based on a Michelson interferometer built around a 50/50 fiber coupler made up of SM800 single-mode optical fiber (Thorlabs). The sample arm incorporates a galvanometric mirror for transverse scanning so that we can obtain 2D images. On the detector side, we use a custom-built spectrometer with a 30 mm focal length achromatic collimating lens at the input. The light is then spectrally dispersed with a transmission volume phase holographic grating (1200 lines/mm, Wasatch Photonics Inc.) before being imaged onto a CMOS line scan camera (BASLER Sprint spl2048–70km) using another achromatic doublet lens with 75 mm focal length. The camera offers 2048 pixels, 12 bit resolution, and a maximum acquisition rate of 70 kHz. The spectrometer provides a spectral bandwidth of 220 nm with a spectral resolution of 0.1 nm. It has been carefully calibrated to avoid any coupling with dispersion compensation [41, 42].

Fig. 2 SD OCT configuration used to demonstrate FrFT-based dispersion compensation. A 26 mm water cell induces a dispersion imbalance between the interferometer arms. M: mirror, L: lens, PC: polarization controller, FC: fiber collimator, GM: galvo mirror, S: sample, SLD: superluminescent diode, Disp: dispersion matching material, G: grating, LSC: line scan camera.

Download Full Size | PDF

The light source of our OCT system is a superluminescent diode with a center wavelength of 843 nm and an optical bandwidth (full width at half maximum, FWHM) of Δλ_FWHM = 110 nm. The theoretical axial resolution, calculated by inverse Fourier transforming the source spectrum (Wiener-Khinchin theorem), is 3.6 μm FWHM. Although care was taken to equalize the fiber lengths between the two arms of the interferometer, dispersion was not completely balanced in the initial setup. Coarse dispersion compensation was done physically by inserting two BK7 microscope slides into the reference arm, which led to an actual axial resolution of 3.8 μm that was close to the theoretical minimum. Dispersion may be matched more precisely e.g. by dispersion compensating prisms of the correct material. The corresponding point-spread-function (PSF), obtained by using a mirror in place of the sample, is plotted as the dashed red curve in Fig. 3. Note that here the measured spectrum has been processed with the traditional FT to obtain the PSF. Figure 3 also shows the theoretical PSF (solid blue curve) for comparison. The difference can be explained by the fact that the two microscope slides do not exactly compensate the residual setup dispersion.

Fig. 3 PSF of our OCT system after coarse physical dispersion compensation (dotted red) compared with the theoretical one (solid blue). Inset: corresponding source spectrum.

Download Full Size | PDF

For all our measurements, the signal coming from the reference arm, measured by blocking the sample arm, and averaged over 100 spectra, is subtracted to obtain the spectrally modulated signal only. The modulated signal is then re-sampled to account for the hyperbolic dependence between wavelength and frequency. Finally, zero padding is performed in order to improve the digital sampling resolution in the transformed domain. The sensitivity of our OCT system was measured as 98 dB at 50 μs exposure time and 100 μm depth. A sensitivity fall off of 11.6 dB was measured at a depth of 1 mm.

4. Results

4.1. Point-spread function measurements

To test the efficiency of the FrFT for numerical dispersion compensation, the interferometer was first purposely dispersion mismatched by placing a 26 mm long water cell in the sample arm. Applying the traditional FT to the acquired spectrum resulted in a dispersed PSF 49 μm wide (FWHM). This is shown as the solid blue curve in Fig. 4(a) and compared with the PSF obtained without the dispersive water cell (dotted red curve), which is the same as that plotted in Fig. 3. Using the FrFT with an optimized order parameter a_opt = 1.0555, one obtains the OCT depth scan while simultaneously fully compensating for group velocity dispersion, as revealed by the dashed black curve in Fig. 4(a). Using the FrFT leads to an axial resolution of 3.65 μm, closer to the theoretical minimum than that observed with coarse physical dispersion compensation before introducing the water cell (3.8 μm). This clearly highlights that the FrFT can efficiently compensate dispersion with continuous adjustment capabilities.

Fig. 4 (a) The dotted red curve is the PSF of our OCT system (3.8 μm FWHM) using only BK7 coarse physical dispersion compensation and the traditional FT (same as in Fig.3). The two other curves are obtained with an additional 26 mm long dispersive water cell. Using the traditional FT (solid blue curve), the PSF broadens to 49 μm FWHM, while the FrFT with optimized order a_opt = 1.0555 numerically compensates dispersion, leading to a 3.65 μm FWHM PSF (dashed black curve). (b) Sharpness function as a function of fractional Fourier transform order, a.

Download Full Size | PDF

The optimized FrFT order a_opt was found by looking for the FrFT order a for which the peak intensity of the PSF is maximized. Figure 4(b) presents a graph of the PSF peak intensity versus a, which can be interpreted as a sharpness function. The peak intensity obtained with the optimized FrFT order is approximately 3.3 times higher than that obtained with the traditional Fourier transform (i.e., without dispersion compensation). This optimization only needs to be done occasionally and subsequent images can be analyzed using the same order value. Note that during in-vivo imaging it can be challenging to find a reference-interface that is suitable for optimization. For retinal imaging, it has been suggested to use the center of the fovea (foveal umbo) for this purpose [32]. In the next subsection, we will introduce another more systematic method.

To provide additional physical insights into FrFT-based dispersion compensation, Fig. 5(a) compares the measured unwrapped spectral phase and Fig. 5(b) the instantaneous spectral fringe frequency [i.e., 1st order derivative of the spectral phase, see Eq. (8)] without the water cell (dotted red), with the water cell (solid blue), and with the water cell plus FrFT dispersion compensation (dashed black). As can be seen, the presence of the water cell induces a strong spectral frequency modulation visible in the inclination of the solid blue curve in Fig. 5(b) compared with the dotted red curve, as was already discussed schematically in Fig. 1. After FrFT dispersion compensation, the curve becomes horizontal. It is even flatter than before the water cell is introduced, showing that the FrFT numerically compensates for both the dispersion of the water cell and the residual dispersion of the setup, i.e., what was left uncompensated by the physical introduction of the BK7 glass. Figure 5(c) and Fig. 5(d) show the Wigner distribution of the signal, respectively before and after FrFT dispersion compensation. It is easily seen that the signal energy is entirely preserved and only rotated about the domains after FrFT.

Fig. 5 (a) Spectral phase and (b) instantaneous spectral fringe frequency without (dotted red) and with the water cell (solid blue), and with the water cell plus FrFT dispersion compensation (dashed black). (c) and (d) shows the Wigner distribution of the spectrum acquired with the water cell before and after FrFT dispersion compensation.

Download Full Size | PDF

For completeness, let us point out that the spectral phase of the FrFT dispersion compensated signal [dashed black curve in Fig. 5(a) and Fig. 5(b)] was obtained from the inverse Fourier transform of the FrFT-compensated complex spectrum (i.e. analytic signal),

φ_{FrFT} = arg {F^{- 1} {F^{a_{opt}} {S (u) + i \hat{S} (u)}}} .

Here Ŝ(u) is the Hilbert transform of the acquired spectrum S(u), and F⁻¹ represents the traditional inverse Fourier transform. To understand this expression, we need to recall that the Wigner transform of a real signal S(u) is symmetric with respect to the transformed axis, W_S(u, u₁) = W_S(u, − u₁). Upon projection on the optimal rotated u_a axis [which amounts to perform the FrFT, see Eq. (4)], out of the two mirror parts of the Wigner transform only one leads to a dispersion compensated response while the other leads to a broadened response in the other half of the transformed domain (not shown in Fig. 1 for simplicity). The broadened part of the response is in essence affected by twice the amount of dispersion present. It is the equivalent of the complex conjugate term encountered while using the traditional FT on a real signal. We will refer to it as such, although strictly speaking it is not, in the general case (a ≠ 1), the mirror image of the other half of the fractional Fourier domain. This “complex conjuguate” term needs to be eliminated or it will distort the retrieved spectral phase after inverse Fourier transform. We solve this problem by starting from the analytic signal.

Finally, the width of the PSF obtained with the water cell after FrFT numerical dispersion compensation was measured for different axial delays. These measurements are plotted in Fig. 6 and compared with that obtained without water cell (and using the traditional FT). The same FrFT order was used for all depths. The resolution observed in the two cases are in very good agreement. Higher order dispersion terms did not disturb our measurements. However, we need to point out a slight increase (7 %) of asymmetric side lobes of the PSF after numerical dispersion compensation, which we attribute to the third order dispersion of the water cell. A similar effect may also arise if higher bandwidths were used. In our current state of knowledge, we cannot readily compensate higher order dispersion using the FrFT but we do not preclude that this may be possible using transformations of the time-frequency plane more complex than rotations (i.e. linear projections).

Fig. 6 Measured PSF for different delays for matched dispersion (dots) and FrFT compensated dispersion (cross).

Download Full Size | PDF

4.2. Imaging

To further validate our method, we have performed FrFT numerical dispersion compensation on two particular samples. The first consists of two stacked 100 μm thick microscope cover slides. The poor surface flatness of the slides produced an air gap of approximately 20 μm between them. Figure 7(a) was obtained using the traditional FT. The gap is blurred due to the poor axial resolution caused by the presence of the dispersive water cell. Using the optimized FrFT, the air gap between the two microscope slides can be resolved as seen in Fig. 7(b).

Fig. 7 OCT depth scan and intensity image of two cover slides (a) using the traditional FT and (b) using FrFT numerical dispersion compensation with the optimized order parameter. The 26 mm-long dispersive water cell is present in the setup in both cases.

Download Full Size | PDF

Similarly, Fig. 8(b) shows the cross sectional image of a grape using the traditional FT and Fig. 8(c) using the FrFT with optimized order, which offered an instantaneous dispersion compensated tomogram. Here, because of the complexity of the image, it was more challenging to use the same sharpness function as described in the previous Section to find the optimized order parameter a_opt. Instead, we have looked for the maximum of the Radon transform of the spectrogram of the full complex depth signal of one random A-scan. This Radon transform plot is shown in Fig. 8(a). Recall that the Radon transform of a two-dimensional distribution consists of a set of integral projections for various projection angles δ [43]. In essence, this procedure is therefore similar to taking the FrFT of the signal for a range of order parameter a [which would correspond to the Radon transform of the Wigner distribution, see Eq. (4)] where we can relate the projection angle δ and the order parameter a through δ = aπ/2. The important difference is that here we used the spectrogram of the signal rather than the Wigner distribution. The spectrogram is another time-frequency distribution that employs a windowed FT of the complex signal, i.e., the short-time FT of the complex signal. It has the advantage of not exhibiting cross-terms, in contrast to the Wigner distribution, because it is phase insensitive. Consequently, the point of the Radon transform where the energy converges determines the optimized order parameter. For our grape image, a quick examination of Fig. 8(a) suggests an average optimized FrFT order a_opt = 1.04 highlighted by the green line. It was used for all A-scans of the tomographic B-scan shown in Fig. 8(c). In contrast, the orange line in Fig. 8(a) represents the traditional FT which leads to Fig. 8(b). Note that the grape was imaged from the bottom through a microscope cover slide and that the front interface of that slide was positioned ahead of the point of zero path difference of the interferometer in order to minimize the depth dependent sensitivity fall-off. The intensity maximum at approximately 2000 Radon bins and FrFT order a = 0.96 therefore corresponds to the complex conjugate term of the front interface of the glass slide. Its intensity maximum is at FrFT order 2 − a_opt as it has opposite dispersion and therefore experiences twice the amount of dispersion after compensation compared to the traditional FT.

Fig. 8 (a) Radon transform of the spectrogram of one A-scan of a grape. (b) B-scan of the grape using the traditional FT (orange line). (c) B-scan of the grape using the optimized FrFT (a_opt = 1.04, green line). The same FrFT order was used for all A-scans. The white bars correspond to 500 μm. Note that the x-axis of (a) only corresponds strictly to depth for the optimized FrFT so is labeled in terms of the pixels of the projections (or Radon bins).

Download Full Size | PDF

To complete this Section, we must note that the algorithm used for numerical computation of the FrFT, and which was introduced by Ozaktas et al. [40], is currently one order of magnitude slower than the Fast Fourier Transform (FFT) algorithm used to compute traditional FTs. With current computing technology, this would however not preclude the use of the FrFT for real time processing.

5. Depth-dependent sample dispersion compensation

An axial resolution below 3 μm is generally sensitive to sample dispersion, so that even thin sample layers cause a broadening of the PSF during imaging. In such situations, the depth-dependent dispersion of the sample is not negligible and must be taken into account in order to get the sharpest image at all depths. For a Gaussian source spectrum, the factor of axial resolution broadening σ_PSF due to a dispersive element of thickness l and group velocity dispersion coefficient β₂ is given by

σ_{PSF} = \sqrt{1 + {(\frac{π^{2} c^{2}}{2 ln (2)} \frac{l β_{2} Δ λ_{FWHM}^{2}}{λ_{c}^{4}})}^{2}} .

Applying this relation to the parameters of the source used in our experiments (see Section 3) reveals that sample dispersion is negligible in our case: a 1 mm thick slide of BK7 glass (β₂ = 41 ps²/km) would only broaden the PSF by a factor of 1.2. Given the limitations of our equipment, we have therefore been unable to investigate experimentally depth-dependent dispersion compensation at this stage. However, in order to demonstrate that the FrFT can be used to handle this problem, and for the completeness of this article, we present below a proof-of-principle demonstration of such capability based on numerical simulations. In all these simulations, we consider a source bandwidth Δλ_FWHM = 210 nm at λ₀ = 710 nm center wavelength.

When the sample dispersion is not negligible, the optimal FrFT order required for dispersion compensation becomes a function of imaging depth. Figure 9 highlights this issue for a simulated sample made up of five identical 100 μm-thick microscope cover slides, stacked against each other, and for which we assume an average sample group velocity dispersion of β₂ = 54 ps²/km. In Fig. 9(a), we have plotted the sharpness function [refer to Fig. 4(b)] at each interface, i.e., the intensity of the OCT signal at those depths, versus the FrFT order. Because dispersion gradually increases across the sample, choosing the FrFT order to compensate dispersion and maximize the PSF at a certain depth means the signals simultaneously obtained at the other depths are not optimal. For completeness, Fig. 9(b) shows the optimum FrFT order required for dispersion compensation at each depth.

Fig. 9 (a) Sharpness function at each interface of a simulated sample made up five identical 100 μm-thick cover slides as a function of FrFT order. (b) Corresponding optimal FrFT order as a function of imaging depth.

Download Full Size | PDF

To perform depth-dependent dispersion compensation, it is possible to generalize the approach used to sharpen the image of the grape (Fig. 8) and that was based on looking for a maximum of the Radon transform of a time-frequency distribution of the OCT signal. To illustrate this, we show in Fig. 10(a) a two-dimensional plot of the FrFT of the OCT signal of our simulated sample for a range of orders, which is equivalent to the Radon transform of the Wigner distribution of the signal. We can immediately see that taking a cross-section of that plot along the diagonal dashed yellow line that links the maxima, and projecting it along the depth axis, leads to the depth-dependent dispersion compensated OCT depth scan shown in Fig. 10(d). In essence, this procedure can be interpreted as performing a “short time” FrFT of the signal and using a different optimal FrFT order for each depth. For comparison, we also show in Fig. 10(b) and Fig. 10(c) the depth scans obtained, respectively, when performing a traditional FT (corresponding to the solid green line in the Radon plot) and when using the same FrFT order at all depths so as to compensate for the average dispersion (dotted red line in the Radon plot). In both cases, dispersion is compensated only for one optimized depth region, whereas other regions appear blurred.

Fig. 10 (a) Color plot of OCT depth scans obtained for a range of FrFT orders for a simulated sample made up five identical 100 μm-thick cover slides. Three cross-sections are highlighted: (b) traditional FT (a = 1, solid green line), (c) average sample dispersion compensation (dotted red line), and (d) depth-dependent dispersion compensation (dashed yellow line).

Download Full Size | PDF

Performing depth-dependent dispersion compensation by finding an optimized cross-section of the Radon plot of the signal can be generalized to more complex samples. In Fig. 11, we show similar plots to Fig. 10 but for a simulated sample which presents normal dispersion in a first 150 μm-thick layer, with a group-velocity dispersion coefficient β₂ linearly increasing with depth from 54 ps²/km to 97 ps²/km, followed by a second layer with uniform anomalous dispersion, β₂ = −38 ps²/km. Here the optimized cross-section is made up of a quadratic and a linear section, which leads to the optimized A-scan plotted in Fig. 10(c). For comparison, Fig. 10(b) is the less satisfactory result obtained with the traditional FT.

Fig. 11 Same as in Fig. 10 but for a simulated sample which presents linearly increasing normal dispersion in a first 150 μm-thick layer, followed by a second layer with uniform anomalous dispersion. (b) is obtained with the traditional FT (solid green line), while (c) comes from the optimized cross-section of the Radon plot (dashed yellow line).

Download Full Size | PDF

OCT signals from real samples may be more challenging to process, but the problem of finding an optimal cross-section passing through maxima of the Radon plot as defined above seems amenable to an appropriate algorithm and may even be applicable to samples with no isolated scatterers. Indeed, a quick look at Fig. 8(a) obtained for our grape image shows that this procedure would have led to the horizontal dashed green line that corresponds to the optimal FrFT order, a = 1.04. Only one maximum, on the left, is excluded, but we have identified this as an artefact resulting from a spurious reflection in our OCT setup. It explains the thin horizontal line visible in both Fig. 8(b) and Fig. 8(c), which always appear at the same depth. We must also point out that third-order dispersion is too small to have an impact on depth dependent dispersion compensation and can therefore be neglected up to an axial resolution of 2 μm, for a center wavelength in the normal dispersion regime of water (λ₀ < 1 μm) [44]. If water is the dominant medium in the sample, only if a better resolution is needed would the approach described here break down, as compensation of higher order terms may become necessary.

6. Group velocity dispersion measurement using FrFT

While dispersion causes broadening of the PSF and blurs images, it has the potential to provide additional functional information. For example, in tissues which exhibit regions with different group-velocity dispersion coefficients, dispersion information can be used for tissue differentiation. Liu et al. proposed to use material dispersion in order to differentiate water and lipid as a diagnosis of vulnerable plaques in the coronary arteries [34]. Obviously, this requires a mean to extract absolute or relative information about the group-velocity dispersion coefficient β₂. We show below that our FrFT-based dispersion compensation routine makes such information readily available. Note that the experimental estimation of the dispersion coefficient β₂ of a material is not only of interest for imaging, but has also applications in many other fields such as in fiber optics or in non-linear optics.

From the analysis of a chirp signal, Eq. (5), we have shown in Fig. 1 how the time-frequency representation provides a simple geometric relationship between the group-velocity dispersion coefficient β₂ and the FrFT order required for optimal dispersion compensation, a_opt. In particular, the modulation rate of the spectral fringes of the OCT signal [see Eq. (8)] can be readily extracted from Fig. 1 without the need of additional phase analysis:

tan (π - ϕ) = \frac{π}{l β_{2} {(2 π η)}^{2}},

where ϕ = a_optπ/2 defines the orientation of the optimal fractional Fourier domain. Solving for the group-velocity dispersion coefficient, and with

η = d ν / \sqrt{N}

the normalization parameter used in our digital implementation [Eq. (7)], yields

β_{2} = \frac{- π}{l N d ω^{2} tan (a_{opt} π / 2)},

where dω = 2πdν is the angular frequency spacing between the samples of the OCT spectral signal. The above equation clearly shows that the absolute value of β₂, including its sign, can be extracted once the optimal FrFT order for dispersion compensation has been obtained.

To demonstrate this capability, we have measured with this technique the group velocity dispersion of a short length of single mode fiber and of distilled water, both at a center wavelength of 843 nm. The fiber sample was a 244 mm long SM800 (Thorlabs) fiber, which has a mode field diameter of 5.6 μm. The water sample was the same as used for the dispersion compensation experiments, i.e., 26 mm thick. In both cases, we first determined the residual dispersion of the setup (lβ₂)^(setup), using Eq. (12) without the sample-under-test, by optimizing the point spread function with only a mirror in the sample arm. This was then subtracted from the corresponding measurement with the sample present (lβ₂)^(total) and normalized to the sample length to yield the dispersion coefficient β₂ of the sample itself,

β_{2}^{(sample)} = \frac{1}{l_{sample}} [{(l β_{2})}^{(total)} - {(l β_{2})}^{(setup)}] .

The results are shown in Table 1, and are in excellent agreement with reference values. The measured fiber dispersion agrees with a measurement that was obtained by traditional white light interferometry [27] while that of water is in good agreement with Van Engen et al., who measured a group velocity dispersion of water of approximately 22 ps²/km at 840 nm. The table also shows the optimal FrFT orders that were obtained first without then with the sample, the former value being reminiscent of the residual uncompensated dispersion of the setup. Note how the residual setup dispersion is anomalous,

a_{opt}^{(setup)} < 1

, for the fiber sample measurement but normal,

a_{opt}^{(setup)} > 1

, for the water measurement. This is due to the empty water cell being included in the setup dispersion in the latter case.

Table 1. FrFT-based measurements of the group velocity dispersion coefficient β₂ of a single mode fiber and of distilled water. All the measurements have been done at 843 nm but the reference value for water is for a 840 nm wavelength. N = 1001, dω = 34.5×10¹⁰ rad/s.

View Table

In OCT one generally has no information about the absolute thickness of depth layers in a cross sectional image as the refractive index of the material is unknown. However, using Eq. (12), one can readily obtain relative dispersion information by using the optimized FrFT order and the relative thickness of a layer. The relative thickness corresponds to l = l̄/n, where n may be assumed as 1.33 and l̄ is the layer thickness obtained from the OCT depth scan.

7. Conclusion

The fractional Fourier transform (FrFT) is a generalization of the traditional Fourier transform. The FrFT provides a new fundamental perspective on the nature and role of group-velocity dispersion in Fourier domain OCT and can be seen as a visual tool to highlight the physics behind dispersion compensation. Using the FrFT one obtains depth information in FD OCT while simultaneously compensating for group velocity dispersion. Our theoretical axial resolution of 3.6 μm was recovered by optimizing the order parameter of the FrFT to compensate the group velocity dispersion induced by a 26 mm long water cell for a source spectral bandwidth of 110 nm. The technique was successfully demonstrated on a biological sample. Furthermore we provided new insights on the issue of depth-dependent dispersion. Simulations have shown that both normal and anomalous sample dispersion can be addressed by dynamically adapting the order parameter as a function of depth. This method can be seen as analogous to a “short time” FrFT but is more efficient and intuitive and may even be applicable to samples without isolated scatterers. From the optimized FrFT order parameter one also readily obtains some quantitative information about the amount of dispersion in an OCT configuration and sample. We have derived the relationship between the FrFT order parameter and group velocity dispersion and used it to successfully measure the group velocity dispersion coefficient of distilled water and some single mode fiber (at 840 nm).

Acknowledgments

This work was supported by a NERF grant from the Foundation for Research Science and Technology from the New Zealand government.

References and links

1. D. Huang, E. A. Swanson, C. P. Lin, J. S. Schuman, W. G. Stinson, W. Chang, M. R. Hee, T. Flotte, K. Gregory, C. A. Puliafito, and J. G. Fujimoto, “Optical coherence tomography,” Science 254, 1178–1181 (1991). [CrossRef] [PubMed]

2. A. F. Fercher, C. K. Hitzenberger, G. Kamp, and S. Y. El-Zaiat, “Measurement of intraocular distances by backscattering spectral interferometry,” Opt. Commun. 117, 43–48 (1995). [CrossRef]

3. R. Leitgeb, C. K. Hitzenberger, and A. F. Fercher, “Performance of Fourier domain vs. time domain optical coherence tomography,” Opt. Express 11, 889–894 (2003). [CrossRef] [PubMed]

4. J. F. de Boer, B. Cense, B. H. Park, M. C. Pierce, G. J. Tearney, and B. E. Bouma, “Improved signal-to-noise ratio in spectral-domain compared with time-domain optical coherence tomography,” Opt. Lett. 28, 2067–2069 (2003). [CrossRef] [PubMed]

5. M. A. Choma, M. V. Sarunic, C. Yang, and J. A. Izatt, “Sensitivity advantage of swept source and Fourier domain optical coherence tomography,” Opt. Express 11, 2183–2189 (2003). [CrossRef] [PubMed]

6. G. Häusler and M. W. Lindner, “Coherence radar and spectral radar — New tools for dermatological diagnosis,” J. Biomed. Opt. 3, 21–31 (1998). [CrossRef]

7. M. Wojtkowski, R. Leitgeb, A. Kowalczyk, T. Bajraszewski, and A. F. Fercher, “In vivo human retinal imaging by Fourier domain optical coherence tomography”, Journal of Biomedical Optics 7, 457–463 (2002). [CrossRef] [PubMed]

8. N. A. Nassif, B. Cense, B. H. Park, M. C. Pierce, S. H. Yun, B. E. Bouma, G. J. Tearney, T. C. Chen, and J. F. de Boer, “In vivo high-resolution video-rate spectral-domain optical coherence tomography of the human retina and optic nerve,” Opt. Express 12, 367–376 (2004). [CrossRef] [PubMed]

9. S. R. Chinn, E. A. Swanson, and J. G. Fujimoto, “Optical coherence tomography using a frequency-tunable optical source,” Opt. Lett. 22, 340–342 (1997). [CrossRef] [PubMed]

10. B. Golubovic, B. E. Bouma, G. J. Tearney, and J. G. Fujimoto, “Optical frequency-domain reflectometry using rapid wavelength tuning of a Cr⁴⁺: forsterite laser,” Opt. Lett. 22, 1704–1706 (1997). [CrossRef]

11. S. H. Yun, G. J. Tearney, J. F. de Boer, N. Iftimia, and B. E. Bouma, “High-speed optical frequency-domain imaging,” Opt. Express 11, 2953–2963 (2003). [CrossRef] [PubMed]

12. R. Huber, M. Wojtkowski, and J. G. Fujimoto, “Fourier Domain Mode Locking (FDML): A new laser operating regime and applications for optical coherence tomography,” Opt. Express 14, 3225–3237 (2006). [CrossRef] [PubMed]

13. L. Liu, J. A. Gardecki, S. K. Nadkarni, J. D. Toussaint, Y. Yagi, B. E. Bouma, and G. J. Tearney, “Imaging the subcellular structure of human coronary atherosclerosis using microoptical coherence tomography,” Nature Med. 17, 1010–1014 (2011). [CrossRef] [PubMed]

14. R. A. Leitgeb, M. Villiger, A. H. Bachmann, L. Steinmann, and T. Lasser, “Extended focus depth for Fourier domain optical coherence microscopy” Opt. Lett. 31, 2450–2452 (2006). [CrossRef] [PubMed]

15. K.-S. Lee and J. P. Rolland, “Bessel beam spectral-domain high-resolution optical coherence tomography with micro-optic axicon providing extended focusing range,” Opt. Lett. 33, 1696–1698 (2008). [CrossRef] [PubMed]

16. C. Blatter, B. Grajciar, C. M. Eigenwillig, W. Wieser, B. R. Biedermann, R. Huber, and R. A. Leitgeb, “Extended focus high-speed swept source OCT with self-reconstructive illumination,” Opt. Express 19, 12141–12155 (2011). [CrossRef] [PubMed]

17. L. Liu, F. Diaz, L. Wang, B. Loiseaux, J.-P. Huignard, C. J. R. Sheppard, and N. Chen, “Superresolution along extended depth of focus with binary-phase filters for the Gaussian beam,” J. Opt. Soc. Am. A 25, 2095–2101 (2008). [CrossRef]

18. J. Holmes, S. Hattersley, N. Stone, F. Bazant-Hegemark, and H. Barr, “Multi-channel Fourier domain OCT system with superior lateral resolution for biomedical applications,” Proc. of SPIE 684768470O (2008). [CrossRef]

19. B. A. Standish, K. K. C. Lee, A. Mariampillai, N. R. Munce, M. K. K. Leung, V. X. D. Yang, and I. A. Vitkin, “In vivo endoscopic multi-beam optical coherence tomography,” Phys. Med. Biol. 55, 615–622 (2010). [CrossRef] [PubMed]

20. C. K. Hitzenberger, A. Baumgartner, W. Drexler, and A. F. Fercher, “Dispersion effects in partial coherence interferometry: implications for intraocular ranging,” J. of Biomed. Opt. , 144–151 (1999). [CrossRef]

21. A. F. Fercher, C. K. Hitzenberger, M. Sticker, R. Zawadzki, B. Karamata, and T. Lasser, “Dispersion compensation for optical coherence tomography depth-scan signals by a numerical technique,” Opt. Comm. 204, 67–74 (2002). [CrossRef]

22. A. G. Van Engen, S. A. Diddams, and T. S. Clement, “Dispersion measurements of water with white-light interferometry,” Appl. Opt. 37, 5679–5686 (1998). [CrossRef]

23. B. Bouma, G. J. Tearney, S. A. Boppart, M. R. Hee, M. E. Brezinski, and J. G. Fujimoto, “High-resolution optical coherence tomographic imaging using a mode-locked Ti:Al₂O₃ laser source,” Opt. Lett. 20, 1486–1488 (1995). [CrossRef] [PubMed]

24. W. Drexler, U. Morgner, F. X. Kartner, C. Pitris, S. A. Boppart, X. D. Li, E. P. Ippen, and J. G. Fujimoto, “In vivo ultrahigh-resolution optical coherence tomography,” Opt. Lett. 24, 1221–1223 (1999). [CrossRef]

25. W. Drexler, U. Morgner, R. K. Ghanta, F. X. Krtner, J. S. Schuman, and J. G. Fujimoto, “Ultrahighresolution ophthalmic optical coherence tomography,” Nature Medicine 7, 502–507 (2001). [CrossRef] [PubMed]

26. G. J. Tearney, B. E. Bouma, and J. G. Fujimoto, “High-speed phase- and group-delay scanning with a grating-based phase control delay line,” Opt. Lett. 22, 1811–1813 (1997). [CrossRef]

27. S. Iyer, S. Coen, and F. Vanholsbeeck, “Dual-fiber stretcher as a tunable dispersion compensator for an all-fiber optical coherence tomography system,” Opt. Lett. 34, 2903–2905 (2009). [CrossRef] [PubMed]

28. L. Froehly, S. Iyer, and F. Vanholsbeeck, “Dual-fibre stretcher and coma as tools for independent 2nd and 3rd order tunable dispersion compensation in a fibre-based ‘scan-free’ time domain optical coherence tomography system,” Opt. Commun. 284, 4099–4106 (2011). [CrossRef]

29. A. F. Fercher, C. K. Hitzenberger, M. Sticker, R. Zawadzki, B. Karamata, and T. Lasser, “Numerical dispersion compensation for Partial Coherence Interferometry and Optical Coherence Tomography,” Opt. Express 9, 610–615 (2001). [CrossRef] [PubMed]

30. J. F. de Boer, C. E. Saxer, and J. S. Nelson, “Stable carrier generation and phase-resolved digital data processing in optical coherence tomography,” Appl. Opt. 40, 5787–5790 (2001). [CrossRef]

31. D. L. Marks, A. L. Oldenburg, J. J. Reynolds, and S. A. Boppart, “Digital algorithm for dispersion correction in optical coherence tomography for homogeneous and stratified media,” Appl. Opt. 42, 204–217 (2003). [CrossRef] [PubMed]

32. B. Cense, N. A. Nassif, T. C. Chen, M. C. Pierce, S.-H. Yun, B. H. Park, B. E. Bouma, G. J. Tearney, and J. F. de Boer, “Ultrahigh-resolution high-speed retinal imaging using spectral-domain optical coherence tomography,” Opt. Express 12, 2435–2447 (2004). [CrossRef] [PubMed]

33. M. Wojtkowski, V. J. Srinivasan, T. H. Ko, J. G. Fujimoto, A. Kowalczyk, and J. S. Duker, “Ultrahigh-resolution, high-speed, Fourier domain optical coherence tomography and methods for dispersion compensation,” Opt. Express 12, 2404–2422 (2004). [CrossRef] [PubMed]

34. 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 system,” Phys. Med. Biol. 49, 923–930 (2004). [CrossRef] [PubMed]

35. J. Liebermann, C. Brckner, B. Grajciar, J. Haueisen, and A. F. Fercher, “Dual-band refractive Low Coherence Interferometry in the spectral domain for dispersion measurements,” Proc. of SPIE 7889, 788922 (2011). [CrossRef]

36. L. Cohen, “Time-frequency distributions — A review,” Proc. of the IEEE 77, 941–981 (1989). [CrossRef]

37. V. Namias, “The fractional order Fourier transform and its application to quantum mechanics,” IMA J. Appl. Math. 25, 241–265 (1980). [CrossRef]

38. D. Mendlovic, H. M. Ozaktas, and A. W. Lohmann, “Graded-index fibers, Wigner-distribution functions, and the fractional Fourier transform,” Appl. Opt. 33, 6188–6193 (1994). [CrossRef] [PubMed]

39. L. Durak and S. Aldirmaz, “Adaptive fractional Fourier domain filtering,” Sig. Proc. 90, 1188–1196 (2010). [CrossRef]

40. H. M. Ozaktas, O. Arıkan, M. A. Kutay, and G. Bozdağı, “Digital computation of the fractional Fourier transform,” IEEE Trans. Sig. Proc. 44, 2141–2150 (1996). [CrossRef]

41. C. Dorrer, N. Belabas, J. P. Likforman, and M. Joffre, “Spectral resolution and sampling issues in Fourier-transform spectral interferometry,” J. Opt. Soc. Am. B 17, 1795–1802 (2000). [CrossRef]

42. M. Mujat, B. H. Park, B. Cense, T. C. Chen, and J. F. de Boer, “Autocalibration of spectral-domain optical coherence tomography spectrometers for in vivo quantitative retinal nerve fiber layer birefringence determination,” J. Biomed. Opt. 12(4), 041205 (2007). [CrossRef] [PubMed]

43. A. C. Kak and M. Slaney, Principles of Computerized Tomographic Imaging (IEEE Press, 1988).

44. T. R. Hillman and D. D. Sampson, “The effect of water dispersion and absorption on axial resolution in ultrahigh-resolution optical coherence tomography,” Opt. Express 13, 1860–1874 (2005). [CrossRef] [PubMed]

Sample	l	$a_{opt}^{(setup)}$	$a_{opt}^{(total)}$	$β_{2}^{(sample)}$	Reference value

SM800	244 mm	0.9986	1.2161	38.3 ± 0.9 ps²/km	38 ps²/km [27]
H₂O	26 mm	1.0059	1.0193	21.4 ± 0.4 ps²/km	22 ps²/km [22]

Dispersion compensation in Fourier domain optical coherence tomography using the fractional Fourier transform

Abstract

1. Introduction

2. Theory

3. Experimental setup

4. Results

4.1. Point-spread function measurements

4.2. Imaging

5. Depth-dependent sample dispersion compensation

6. Group velocity dispersion measurement using FrFT

7. Conclusion

Acknowledgments

References and links

Cited By

Figures (11)

Tables (1)

Equations (13)

Optics Express