## Abstract

Polarization mode dispersion (PMD), which can be induced by circulators or even moderate lengths of optical fiber, is known to be a dominant source of instrumentation noise in fiber-based PS-OCT systems. In this paper we propose a novel PMD compensation method that measures system PMD using three fixed calibration signals, numerically corrects for these instrument effects and reconstructs an improved sample image. Using a frequency multiplexed PS-OFDI setup, we validate the proposed method by comparing birefringence noise in images of intralipid, muscle, and tendon with and without PMD compensation.

© 2013 Optical Society of America

## 1. Introduction

Polarization-sensitive optical coherence tomography (PS-OCT) [1, 2] combines the structural imaging of conventional OCT [3, 4] with the ability to map birefringence, diattenuation, and polarization-dependent scattering. For biological tissues, this capability can potentially add compositional information to images of anatomy without requiring exogenous contrast agents [5–7]. Early PS-OCT systems were constructed using free-space interferometers that allowed for control of the light polarization state throughout the system. Fiber-optic implementations, however, are now needed for many clinical applications. The inclusion of optical fiber and the use of optical devices such as circulators invariably introduces polarization-mode dispersion (PMD) and polarization-dependent loss (PDL) into the system. It was recognized in early PS-OCT development that PMD can degrade system performance [8]. Recently, a framework for quantifying the polarimetry noise induced by PMD revealed the significance of this noise; it was shown that a few meters of single-mode fiber can induce sufficient PMD to influence the polarimetry noise in high-resolution PS-OCT systems [9].

Given that PMD has been shown to be a major source of noise in PS-OCT, one strategy for minimizing effects of PMD is to replace circulators inside the interferometer with fused fiber couplers and to minimize fiber lengths. The removal of circulators, however, results in a signal-to-noise ratio (SNR) penalty at least 4.5 *dB*, significantly degrading image quality. And in many systems, reducing the fiber lengths is inconsistent with practical clinical requirements. In this work, we present a method for compensating system PMD numerically using measurements from calibration signals. Therefore highly sensitive imaging in the presence of significant system PMD can be achieved without compromising SNR or limiting fiber lengths. The approach we present is fully general in the order of PMD that can be compensated and in the location of the PMD within the sample arm. This manuscript is organized as follows. We first show system PMD can be removed by multiplying the detected fringe data by two compensation matrices. Next, we demonstrate how these compensation matrices can be derived directly from measured fringe data. Finally, we demonstrate the application of this method to reduce PMD-induced noise in PS-OCT images of intralipid, tendon, and muscle.

## 2. PS-OCT system model

We start by defining a basic fiber-based PS-OCT system (Fig. 1(a)) alongside its simplified optical equivalent (Fig. 1(b)). The optical equivalent system provides a more straightforward model to use during the mathematical analyses of this work. In the optical equivalent system, the distributed transfer properties of the sample arm have been accumulated into discrete blocks (INST1, SMPL, and INST2), and the bi-directional propagation within a portion of the sample arm has been hidden by including its effect within the INST1 and INST2 blocks. The SMPL block describes only the tissue to be imaged. We will associate each block INST1, SMPL, and INST2 with a wavelength-dependent Jones matrix **T _{in}**(

*k*),

**R**(

*k*), and

**T**(

_{out}*k*) respectively. The reference arm in the fiber system includes a polarizer immediately prior to the combining beam-splitter. Because this polarizer prevents PMD in the reference arm from inducing polarimetry noise [9], we have not included a birefringence block element within the reference arm of the optical equivalent system.

PS-OCT systems are used to obtain information about the optical transfer properties of the sample (SMPL). Ideally, the system measures the properties of the sample (described by **R**(*k*)) without noise and confounding artifacts from the instrument transfer functions (described by **T _{in}**(

*k*) and

**T**(

_{out}*k*)). In practice, however, the detected signal is associated with both the sample properties and the instrument transmission properties according to

**M**(

*k*) denotes the Jones vector describing the light immediately prior to the receiver (location D in Fig. 1(b)), and where a single launched polarization state (linear in X) is assumed. Here all transfer functions (Jones matrices) are expressed in normalized form, ignoring an overall amplitude and phase scaling factor.

A brief note about the notations used in this work is appropriate here. For simplicity, we will associate the Jones vector **M**(*k*) with direct measurements from the polarization-diverse receiver. Because the receiver measures the real component of the fields, and not the complex signal as needed to construct the Jones vector, there is an additional processing step required to translate the direct measurements to the Jones vector **M**(*k*). Multiple approaches have demonstrated how to achieve this, and it is commonplace in OCT processing. In a system without complex-conjugate separation, if it is assumed that there are no signals in either the positive-delay or negative-delay space, then the complex fringe signals can be recovered from the measured real fringes through a mathematical Hilbert transformation. In systems that include complex-conjugate separation by including an acousto-optic frequency shifter, the complex fringes are obtained by simply demodulating the signal about its carrier frequency. This later approach [10,11] is the one used in the empirical studies of this work. For simplicity of the mathematical development, we do not explicitly include this step, and simply associate **M**(*k*) with the measured signal.

## 3. Physical PMD compensation

Because PMD acts upon optical signals in a linear manner, its effect can be undone. In principle, PMD compensation can be achieved physically (by inserting optical elements into the system) or numerically (by performing calculations on the measured signals). Physical PMD compensators are difficult to construct, and those that exist compensate only a limited magnitude of PMD [12]. By comparison, a numerical compensation strategy could offer arbitrary PMD compensation and would be implemented easily as an additional stage within the existing OCT processing framework. Thus, we focus on the numerical approach in this work. However, in the interest of communicating the basis for numerical compensation, it is helpful to first present the conceptually more intuitive physical approach, and to build upon this in describing the numerical approach.

Physical PMD compensation involves insertion of an optical element with PMD that cancels that of another element’s PMD. In an OCT system, full PMD compensation would require inserting a first compensator COMP1 before the sample and a second compensator COMP2 after the sample (Fig. 2). Mathematically, the Jones vector **M**(*k*) describing the optical field immediately prior to the detector is modified as

**C**(

_{in}*k*) and the second compensator (COMP2) is described by

**C**(

_{out}*k*). If the compensators are configured such that

**C**(

_{in}*k*) =

*inv*[

**T**(

_{in}*k*)] and

**C**(

_{out}*k*) =

*inv*[

**T**(

_{out}*k*)], the measured signal becomes associated only with the sample birefringence

**R**(

*k*) as It is important to note that for arbitrary

**R**(

*k*), as is the case for OCT imaging applications, the elements INST1 and INST2 must be compensated independently as described by Eq. (2), and each compensation must take place at the location of its associated optical element. The implications of this requirement on the design of a numerically compensated PS-OCT system are described in the next section.

## 4. Numerical PMD compensation

#### 4.1. Applying numerical compensation

To understand how numerical compensation of INST2 can work, we note that the compensator COMP2 performs superposition of the polarization fields that would have been detected (at detector X and detector Y) in its absence. If the fields had been directly detected in a coherent manner, i.e., with both amplitude and phase information and with a meaningful phase relationship between them, then the same mixing could be performed numerically. This is the basis for numerical compensation of PMD located after the sample.

The PMD induced by the element INST1 can also be removed numerically. Because it is located before the sample, however, it requires a different approach. Rather than mixing measurements of the X and Y detector, the numerical compensation of INST1 requires mixing of launched fields. Thus, for numerical compensation, we must design the PS-OCT system to launch two polarization states simultaneously (Fig. 3). These polarization states should nominally be orthogonal, but are in general only required to be non-degenerate. If the signals associated with each of these launched polarization states are separable at the receiver, and if the phase relationship between these signals are meaningful, then it becomes possible to simulate the effect of an arbitrary PMD element located prior to the sample by mixing the launched fields after detection. This mixing of launched fields is the basis for numerical compensation of PMD located before the sample.

Multiple approaches have been demonstrated for performing simultaneous illumination in two polarization states in PS-OCT. Examples include acousto-optic frequency encoding of each polarization state [11, 13], rapid intra A-line modulation of the polarization state [14, 15], or path-length encoding using fixed differential group delay lines [16]. In this work, we will adopt the frequency encoding model which is illustrated in Fig. 3 and described in detail in Ref. [11].

To develop a mathematical description of numerical compensation, we first modify Eq. (1) to include simultaneous launching of two polarization states:

**S**(

*k*) is a 2×2 matrix and each column describes the launched state of polarization (SOP) of the two distinguishable launched signals, and

**M**(

*k*) is now a 2×2 matrix, with each column describing the Jones vector associated with the respective launched states of

**S**(

*k*). For simplicity, we start by assuming we launch orthogonal polarization states, aligned with the laboratory coordinates X and Y, yielding

**S**(

*k*) = [1 0; 0 1]. It is shown at the end of this section that the analysis holds for arbitrary (non-degenerate)

**S**(

*k*). Because we have assumed

**S**(

*k*) to be given by the unity matrix, it can be discarded from Eq. (4).

We now assume a compensated complex fringe signal is calculated from the measured fringe signal **M**(*k*) as

**C**(

_{in}*k*) and

**C**(

_{out}*k*) are our calibration matrices:

**T**(

_{in}*k*) and

**T**(

_{out}*k*) from the measurement system

The compensated result **C _{out}**(

*k*) ·

**M**(

*k*) ·

**C**(

_{in}*k*) may differ from actual

**R**(

*k*) by an arbitrary but wavelength-independent rotation [17, 18]. For many PS-OCT measurements, local retardance calculation is not affected because parameters are calculated as differentials across depth. A recent work [18] also demonstrated measurement of absolute fast axis of the sample by extra calibration signals.

In the derivation of this result, we assumed launching of polarization states according to an identity matrix, i.e., **S**(*k*) = [1 0; 0 1]. In practice, errors in the implementation of a simultaneous polarization probing system will cause the exact form of **S**(*k*) to be unknown and potentially time-varying. However, assumption of identity form can still be used; in this case, deviations from this assumption can be included into the **T _{in}**(

*k*) matrix, solved for, and calibrated out. Thus, it is not essential that strictly orthogonal polarization states be launched, or that their amplitudes, phases, or group delays be matched. Similarly, many imperfections in the polarization-diverse receiver can be included into

**T**(

_{out}*k*), and removed from the system by this calibration process.

In summary, we have shown that if two signals of differing polarization states are launched and their associated signals can be separately detected (e.g., by frequency encoding each at a unique carrier), and if a polarization-diverse receiver is used, then it is possible to numerically remove the transfer function of the system including PMD from OCT systems. This is achieved by multiplying the measured complex fringe signal **M**(*k*) by wavelength-dependent compensation matrices **C _{in}**(

*k*) and

**C**(

_{out}*k*) according to Eq. (5).

#### 4.2. Solving for the compensation matrices

In the previous subsection, we demonstrated that for an appropriately designed PS-OCT system, we can numerically compensate for PMD induced before and after the sample by multiplying the detected fringe signals by compensation matrices **C _{in}**(

*k*) and

**C**(

_{out}*k*) (see Eq. (5)). In this section, we describe a method for solving for these compensation matrices based on measurements of calibration samples. The basic principle behind this approach is relatively straightforward, but its implementation in practice becomes complicated by the presence of amplitude and phase perturbations dependent on signal depth and by the lack of synchronization between the laser source and data acquisition. In this subsection, we focus on the basic principle, and leave the detailed implementation of numerical PMD compensation including some of the steps necessitated by the PS-OCT empirical system for Section 5.

We assume that the system is configured to provide three measurements (Appendix), each with a unique sample retardance but seeing the same instrument optical elements INST1 and INST2 transmission. Figure 4 illustrates the optical system. The three measured signals **M _{1}**(

*k*),

**M**(

_{2}*k*), and

**M**(

_{3}*k*) can be expressed as

**R**(

_{1}*k*),

**R**(

_{2}*k*), and

**R**(

_{3}*k*) to vary with wavelength. Because these are calibrating signals from a known optical path, we assume that the Jones matrices

**R**(

_{1}*k*),

**R**(

_{2}*k*), and

**R**(

_{3}*k*) are known a priori. By combining Eqs. 8a and 8b, and Eqs. 8a and 8c, we arrive the following set of two equations for

**T**(

_{out}*k*):

**R**(

_{2}*k*) ·

**R**

_{1}^{−1}(

*k*)] and [

**R**(

_{3}*k*) ·

**R**

_{1}^{−1}(

*k*)] are known a priori, and [

**M**(

_{2}*k*) ·

**M**

_{1}^{−1}(

*k*)] and [

**M**(

_{3}*k*) ·

**M**

_{1}^{−1}(

*k*)] are measured. If the eigen decomposition of [

**R**(

_{2}*k*) ·

**R**

_{1}^{−1}(

*k*)] and [

**R**(

_{3}*k*) ·

**R**

_{1}^{−1}(

*k*)] generate different eigenvectors, it can be shown that Eq. (9) uniquely defines

**T**(

_{out}*k*). They can be solved using numerical solving routines such as Levenberg-Marquardt algorithm [19]. With

**T**(

_{out}*k*),

**T**(

_{in}*k*) can be solved for using any of Eqs. 8a–c.

## 5. Implementation in a frequency-multiplexed PS-OCT system

#### 5.1. PS-OCT system design

The proposed approach for PMD compensation was implemented into a PS-OCT system using a wavelength-swept laser and frequency-multiplexing of launched polarization states. The system is presented in Fig. 5. The laser source [20] was operated at a repetition rate of 60*kHz* over 140*nm* sweeping range centered at 1270*nm* and yielded 6 *μm* axial resolution in tissue. The carrier frequencies of two simultaneous polarizations were centered at 25*MHz* and 75*MHz* by placing a 50*MHz* frequency shifter in the reference arm and two identical 25*MHz* shifters in the sample arm, one oriented to upshift and the other to downshift. The drive signals to each frequency shifter were phase-locked to each other, but not to the digitizer clock. In the sample arm, a fiber coupler was used in place of a circulator to allow calibrated PMD to be placed into the system using polarization-maintaining fiber (PMF) patchcords [9]. The light power launched to the sample was measured to be 10*mW*, and the back-reflected light was directed to a polarization-diverse receiver.

In contrast to the frequency multiplexing system demonstrated previously which launched polarizations oriented at linear 0° and linear 45° polarizations [11], this system was configured to launch orthogonal signals. Because a Jones matrix formalism was used that includes phase relation between two launched states, orthogonal polarizations can be used [13].

For this approach to work correctly, it was essential that the optical delays seen by each of the four paths of the balanced receiver be equalized. Variations in optical path length between the interfering beamsplitter (BBS in Fig. 5) and the receiver, or in the length of RF cabling between the detector and the digitizer can induce an RF frequency-dependent phase shift between channels, which manifests in PS-OCT as a depth-dependent change in polarization. If present, this will frustrate the compensation algorithm. In principle, if any differential delay is present, it could be removed as a separate step in post-processing. However, we chose to physically match these lengths to a tolerance of a few centimeters. At 100*MHz* (the highest RF frequency in our detection bandwidth), a 1*cm* mismatch induces 0.02*radians* of phase shift, and a correspondingly low rotation of polarization.

#### 5.2. Microscope design

To generate the three calibrating signals needed to solve for **C _{in}**(

*k*) and

**C**(

_{out}*k*), we built a specialized microscope that provides three fixed-depth signals as illustrated in Fig. 6. A small fraction of the imaging beam was reflected almost directly backward using a glass beam sampler (approximately 4% reflection). This light was then directed toward beamsplitters, separating the light into three beams directed at three corresponding mirrors. The first beam (signal 1) was not modified by a waveplate, while beams 2 and 3 included quarter-wave plates oriented at 0° and 45°. The beam sampler was implemented at nearly normal incidence to minimize the polarization-dependent loss. Each mirror included z-translation to adjust its fringe frequency. Based on this setup, the assumed functional form for

**R**(

_{1}*k*),

**R**(

_{2}*k*), and

**R**(

_{3}*k*) are given by

*nm*, and with the laser swept range we can calculate Γ(

*k*).

#### 5.3. Extracting depth-encoded signals **M**_{1}(k), **M**_{2}(k), and **M**_{3}(k) from a single A-line

_{1}

_{2}

_{3}

The digitized fringe data were first demodulated about their 25*MHz* and 75*MHz* carriers to generate complex fringe signals. Next, the fringes were interpolated to remove laser chirp and multiplied by a dispersion correcting vector. This results in four complex fringes defining the combinations of two input polarizations and two detection polarizations. To separate each of the three calibrating signals, each fringe was Fourier transformed to give A-lines. Each mirror signal appeared as a distinct spectral peak within the A-line. Each spectral peak was isolated by applying a Hanning window of 100*μm* FWHM width at its center, and the windowed A-lines were inverse Fourier transformed to yield the complex fringe associated with that mirror signal. This was repeated for each of four complex A-lines (two launched polarization states and two detected polarization states), and for each of the three mirror signals. From these 12 complex fringe signals, the wavelength-dependent matrices **M _{1}**(

*k*),

**M**(

_{2}*k*), and

**M**(

_{3}*k*) were constructed.

#### 5.4. Normalizing the measured signals **M**_{1}(k), **M**_{2}(k), and **M**_{3}(k)

_{1}

_{2}

_{3}

By constructing the measured signals **M _{1}**(

*k*),

**M**(

_{2}*k*), and

**M**(

_{3}*k*) as described in Section 5.3 using the system presented in Section 5.1, we must modify our previous model (Eq. (8)) to include effects associated with the phase and amplitude stability of the fringes, and of the depth-dependent phase variation in each fringe. In this section, we present a methodology to normalize out these factors. We start by writing a more generalized model for the measured signals including a set of unknown factors:

Here, the scalar factors *c*_{1}(*k*), *c*_{2}(*k*), and *c*_{3}(*k*) include the amplitude and phase variations that are associated with each calibration signal and that originate from the physical properties of the calibration paths (rather than from laser synchronization, for example). Such variations could arise from the wavelength-varying splitting ratio of beamsplitters in each path. Variations induced by the frequency-multiplexing arrangement, or which differ between one frequency channel relative to the other, are included in the matrix **K**(*k*). We note that unlike the scalars *c*_{1}(*k*), *c*_{2}(*k*), and *c*_{3}(*k*) which are unique for each mirror signal, the matrix **K**(*k*) is assumed to be constant for all mirror signals. Physically, **K**(*k*) can result in part from the phase-instability that is present uniquely in each frequency-multiplexed polarization channel, and therefore cannot be included in the scalars factors. This phase instability results in large part from the lack of synchronization between the laser and acquisition clock as described in [21]. It is further complicated by the presence of frequency shifters which introduces a second phase-shift into each A-line, and this phase shift differs between the two frequency-multiplexed channels (since the frequency shifts differ). The aline-to-aline instability of the matrix K(k) will depend on the implementation of the imaging system, and would be modified if, for example, k-clocking were implemented into the system. Finally, a slowly drifting path-length in each arm of the frequency multiplexing set-up will vary the relative phase of each frequency-multiplexed channel. The amplitude variation between frequency-multiplexed channels and instability can originate from the differing transmission properties of each arm of the frequency-multiplexing set-up including. All the factors *c*_{1}(*k*), *c*_{2}(*k*), *c*_{3}(*k*), and **K**(*k*) are assumed to be A-line dependent, i.e., to vary from one A-line to the next.

Using Eq. (11) in place of Eq. (8), we form the set of two equations analogous to Eq. (9) as

**K**(

*k*) has been cancelled. Here, [

**R**(

_{2}*k*) ·

**R**

_{1}^{−1}(

*k*)] and [

**R**(

_{3}*k*) ·

**R**

_{1}^{−1}(

*k*)] are known a priori, and [

**M**(

_{2}*k*) ·

**M**

_{1}^{−1}(

*k*)] and [

**M**(

_{3}*k*) ·

**M**

_{1}^{−1}(

*k*)] are measured, but the scalars

*c*

_{1}(

*k*),

*c*

_{2}(

*k*), and

*c*

_{3}(

*k*) are not known. Our goal is to either calculate these scalars, or equivalently to derive a new left-hand side of Eq. (12) that is independent of the scalars and allows

**T**(

_{out}*k*) to be solved for. We note that $\left(\frac{{c}_{1}\left(k\right)}{{c}_{2}\left(k\right)}\right)\cdot \left[{\mathbf{M}}_{\mathbf{2}}\left(k\right)\cdot {{\mathbf{M}}_{\mathbf{1}}}^{-1}\left(k\right)\right]$ is a similarity transformation of [

**R**(

_{2}*k*) ·

**R**

_{1}^{−1}(

*k*)]. Therefore, each has the same eigenvalues upon eigen decomposition [22]. Stated alternatively, if the eigen decomposition of [

**R**(

_{2}*k*) ·

**R**

_{1}^{−1}(

*k*)] is written by

**V**(

_{21}*k*) contains the eigenvectors and

**D**(

_{21}*k*) is a diagonal matrix containing the eigenvalues, then the eigendecomposition of [

**M**(

_{2}*k*) ·

**M**

_{1}^{−1}(

*k*)] (absent the scalar factors) can be written by

**U**(

_{21}*k*) contains the eigenvectors,

**d**(

_{21}*k*) is a diagonal matrix containing the eigenvalues, and

**D**(

_{21}*k*) and

**d**(

_{21}*k*) differ only because of the presence of the scalar factors

*c*

_{1}(

*k*),

*c*

_{2}(

*k*), and

*c*

_{3}(

*k*). We can therefore calculate the lhs (left-hand side) of Eq. (12a) by first calculating the eigendecomposition of [

**M**(

_{2}**k**) ·

**M**

_{1}^{−}

**(**

^{1}**k**)] to get

**U**(

_{21}*k*), and then using the eigenvalues matrix

**D**(

_{21}*k*):

**T**(

_{out}*k*) in the manner described in Section 4.2.

Substituting the solution for **T _{out}**(

*k*) into any of Eqs. 11 allows us to solve for the product

**T**(

_{in}*k*) ·

**K**(

*k*). While it might appear that we need to uniquely resolve

**T**(

_{in}*k*) from

**K**(

*k*), it can be appreciated from the form of Eq. (11) that we can simply include

**K**(

*k*) into the definition of

**T**(

_{in}*k*). We therefore require the product

**T**(

_{in}*k*) ·

**K**(

*k*) since the signal fringes (e.g., from the sample) will be modified by this product in the same manner as the calibrating signals (from the mirrors). However, by defining

**T**(

_{in}*k*) to include

**K**(

*k*), we impose the A-line to A-line variability of

**K**(

*k*) onto the product

**T**(

_{in}*k*) ·

**K**(

*k*). Thus, even for a stable system where the optical transfer function of the block elements INST1 and INST2 remain fixed, the effective optical transfer function for INST1 when including the phase instability described by

**K**(

*k*) varies from one A-line to the next. By contrast, the transfer function of INST2 (denoted by

**T**(

_{out}*k*)) is A-line stable. Thus, the solution for

**T**(

_{out}*k*) provided by solving Eq. (12) can be used repeatedly across A-lines for as long as the optical transfer function from the sample to receiver part is stable. In our experiment,

**T**(

_{out}*k*) was calculated from one A-line measurement within a frame, and used for all A-lines within the frame. It was thus updated every 1024 a-lines. However, the product

**T**(

_{in}*k*) ·

**K**(

*k*) must be recalculated for each A-line using the solution of

**T**(

_{out}*k*) and one of Eqs. 11. As a result, three signals are required to solve for the combination of

**T**(

_{out}*k*) and

**T**(

_{in}*k*) ·

**K**(

*k*), but once

**T**(

_{out}*k*) is known, one of the three signals must remain to be used to recalculate

**T**(

_{in}*k*) ·

**K**(

*k*). Finally, the compensation matrices

**C**(

_{in}*k*) and

**C**(

_{out}*k*) are calculated straightforwardly as

It is worth noting that while this method assumes a known functional form for **R _{1}**(

*k*),

**R**(

_{2}*k*), and

**R**(

_{3}*k*), the sensitivity of the approach to errors in the assumed and actual forms depend critically on whether these errors affect the resulting eigenvalues (

**D**(

_{21}*k*) and

**D**(

_{31}*k*)) or the eigenvectors (

**V**(

_{21}*k*) and

**V**(

_{31}*k*)). Errors that affect only the eigenvalues are effectively removed, since the measured data is normalized to the eigenvalues of the assumed form; it is critical therefore only that the eigenvalues between the lhs and rhs match, and not that they accurately describe the true transfer denoted by [

**R**(

_{2}*k*) ·

**R**

_{1}^{−1}(

*k*)] and [

**R**(

_{3}*k*) ·

**R**

_{1}^{−1}(

*k*)]. This flexibility however does not extend to the errors that affect the eigenvectors, as those are not normalized (set equal) on the lhs and rhs of Eq. (12). The implication of this is that in our system, for example, it is important that our expressions for

**R**(

_{1}*k*),

**R**(

_{2}*k*), and

**R**(

_{3}*k*) accurately describe the orientation of the waveplates, but need not accurately describe the magnitude of their retardance. This allows us, for example, to neglect any wavelength-dependence of the retardance of these waveplates. However, we note that the influence of the beamsampler and BBS are not included in the model of calibration signals (Eq. (8)), and it is possible that the performance could be improved if the model is extended to explicitly describe their transmission properties.

#### 5.5. Forcing **T**_{out}(k) to be continuous across wavelength

_{out}

The normalizing procedure of Section 5.4 describes a method for calculating **T _{out}**(

*k*) which includes an eigendecomposition of the measured data [

**M**(

_{2}*k*) ·

**M**

_{1}^{−1}(

*k*)] and [

**M**(

_{3}*k*) ·

**M**

_{1}^{−1}(

*k*)]. This eigendecomposition is applied for each wavelength, resulting in a set of eigenvectors

**U**(

_{21}*k*) and

**U**(

_{31}*k*). There is an order ambiguity, however, in how these eigenvector matrices are constructed; their columns (and associated eigenvalues in the matrices

**d**(

_{21}*k*) and

**d**(

_{31}*k*)) can be swapped. If the eigenvectors are not forced to be oriented within these matrices in a manner that is continuous across wavelength, then this discontinuity will be transferred to the solution of

**T**(

_{out}*k*) and finally to the compensated fringes.

It is necessary therefore to modify (i.e., swap columns within) the eigenvector matrices **U _{21}**(

*k*) and

**U**(

_{31}*k*) to ensure continuity before Eq. (15) is applied. To do this, we stepped across wavelengths and analyzed the change in the first column eigenvector of

**U**(

_{21}*k*) (or

**U**(

_{31}*k*)). Specifically, we calculated the Stokes vector associated with each column, and then looked for the angular displacement between it and that of an adjacent wavelength. A large angular shift indicates a discontinuity in eigenvectors, and triggers a swapping of columns to eliminate this discontinuity.

## 6. Results

All results in this section used the microscope described in Fig. 6 and the calibrating procedures described in Sections 4 and 5. In the experiment, two PM patchcords have been constructed to induce PMD consistent with low-PMD (0.04*ps* per path) and high-PMD (0.08*ps* per path) circulators.

#### 6.1. Mirror sample

We first imaged a mirror sample to illustrate that the influence of PMD was being canceled from the compensated interference fringe. In these measurements, a PMF patchcord was inserted into port 2 of the sample arm coupler. Thus, PMD was induced both before and after the sample. The patchcord induced 0.08*ps* of differential group delay (single-pass). A neutral density (ND) filter was placed before the sample mirror to prevent saturation from the mirror signal.

Compensated fringes were calculated according to Eq. eq5. In Fig. 7, the SOP of each column in the compensated fringe matrix is plotted in the Poincar*é* sphere representation [23], along with the SOP of the uncompensated fringes. An evolution of the SOP across wavelength is clearly evident without compensation, and is dramatically reduced with compensation. Similar results obtained for a second configuration of the polarization controller PC2 and PC3 in Fig. 5 confirm the robustness of the approach.

#### 6.2. Intralipid sample

To demonstrate that numerical PMD compensation leads to reduced noise in PS-OCT images of local birefringence, we imaged a 1% intralipid sample. Because intralipid has no intrinsic birefringence, it is well suited for characterizing noise in PS-OCT. Further, intralipid induces relatively little multiple scattering which can increase noise through an unrelated mechanisms [24]. In this section, local birefringence was calculated using a Jones matrix analysis with diagonalization [25, 26].

In these results, we compare local birefringence images obtained using numerical PMD compensation with local birefringence images obtained without PMD compensation. For images without numerical PMD compensation, we followed previous works [14, 27] and matched the path delay between the frequency-multiplexed channels by comparing the rate of linear phase accumulation in mirror signals. This differential path delay is induced by unequal free-space paths associated with frequency channel 1 and frequency channel 2. When performing full PMD compensation, this length mismatch is included in the instrument transfer functions and compensated along with other sources of PMD.

In Fig. 8(a), the intensity image is presented showing the three fixed calibrating signals located above the intralipid sample. Figure 8(c) and 8(d) present the local birefringence images without and with numerical PMD compensation respectively. In these images, a 0.08*ps* (single-pass) patchcord was inserted into the bi-directional portion of the sample arm. For comparison, Figure 8(b) presents the uncompensated local birefringence image without patchcord. These images demonstrate that numerical PMD compensation (Fig. 8(d)) provides an image qualitatively equivalent to that obtained by physically eliminating the dominant source of PMD (Fig. 8(b)). Quantitatively, the mean local birefringence measured across the topmost 500 *μ*m of the intralipid over 1024 A-lines was 0.35 *deg*/*μm* with the 0.08*ps* patchcord and without PMD compensation (Fig. 8(c)), 0.0076 *deg*/*μm* with the 0.08*ps* patchcord and with PMD compensation (Fig. 8(d)), and 0.0090 *deg*/*μm* without the patchcord and without PMD compensation (Fig. 8(b)).

#### 6.3. Tendon and muscle sample

To demonstrate reduced noise imaging of tissues containing intrinsic birefringence, we imaged chicken muscle and tendon. Figure 9(a) shows the intensity image again comprising the sample and three calibrating signals. In Fig. 9(b), the baseline local birefringence image is presented, obtained without insertion of any PMF patchcord and without PMD compensation. In Fig. 9(c) and 9(e), the local birefringence images obtained with a 0.04*ps* and 0.08*ps* patchcord respectively and without PMD compensation are presented. In Fig. 9(d) and 9(f), the images with the 0.04*ps* and 0.08*ps* patchcords and with numerical PMD compensation are presented. A similar noise reduction is achieved by numerical PMD compensation with these tissue samples as was demonstrated in intralipid.

## 7. Conclusions

We have demonstrated that instrument PMD can be removed numerically from PS-OCT systems. We also demonstrated a method for deriving the compensation matrices from a set of three signals measured simultaneously within a single A-line acquisition. Future work will focus on applying these techniques to promising applications of PS-OCT including intravascular imaging of coronary artery pathology and retinal imaging for Glaucoma diagnosis. Finally, while the discussion within this work focused on removal of PMD, we note that this approach measures and removes the generalized optical transfer function of the instrument including effects such as polarization-dependent loss. The elimination of PDL might have value for extracting, for example, undistorted diattenuation from tissues.

## Appendix: Why three calibration signals are required?

Here we demonstrate that three calibration signals are required in order to explicitly solve for instrument transfer functions **T _{in}**(

*k*) and

**T**(

_{out}*k*). In this appendix, analytic method is derived by ignoring amplitude and phase variation on measured calibration signals. We note, however, that given measured data contaminated by noise and by amplitude/phase variation, a more generalized model and numerical methods for solving for

**T**(

_{in}*k*) and

**T**(

_{out}*k*) are implemented in Section 5.4.

According to the compensation principle shown in Eq. (2), **T _{in}**(

*k*) and

**T**(

_{out}*k*) must be measured and compensated independently before recovering the sample property

**R**(

*k*). Therefore it is obvious to prove that only one calibration signal (for example, a mirror signal) is not enough to resolve both

**T**(

_{in}*k*) and

**T**(

_{out}*k*), as only the combined effect [

**T**(

_{out}*k*) ·

**T**(

_{in}*k*)] can be obtained.

Furthermore, we will demonstrate that two calibration signals are not able to uniquely derive **T _{in}**(

*k*) and

**T**(

_{out}*k*). We start with Eq. (9a) to solve for

**T**(

_{out}*k*) by combing Eqs. 8a and 8b

**R**(

_{2}*k*) ·

**R**

_{1}^{−1}(

*k*)] is known a priori from calibration signals, and [

**M**(

_{2}*k*) ·

**M**

_{1}^{−1}(

*k*)] is measured. Note that [

**M**(

_{2}*k*) ·

**M**

_{1}^{−1}(

*k*)] is a similarity transformation of [

**R**(

_{2}*k*) ·

**R**

_{1}^{−1}(

*k*)], their eigenvalues are equal upon eigendecomposition [22]:

**U**and

_{21}**V**are eigenvectors of [

_{21}**M**(

_{2}*k*) ·

**M**

_{1}^{−1}(

*k*)] and [

**R**(

_{2}*k*) ·

**R**

_{1}^{−1}(

*k*)], respectively, while

**D**is their eigenvalue matrix. Equation (17) can be rearranged to write

**D**is a diagonal matrix, it can be shown that a necessary and sufficient condition for this equation to be solved is that [

**U**

_{21}^{−1}·

**T**·

_{out}**V**] is diagonal. With this requirement, we can write a general form for [

_{21}**U**

_{21}^{−1}·

**T**·

_{out}**V**as

_{21}]*α*

_{11}and

*α*

_{22}can be any complex value. The analytic solution for

**T**becomes any linear combination of two orthogonal bases

_{out}**T**=

_{out}*α*

_{11}

**A**+

_{1}*α*

_{22}

**A**, where

_{2}**A**and

_{1}**A**can be calculated by

_{2}Hence, with two calibration signals, we find an infinity of linear combinations and not a unique solution for **T _{out}**(

*k*). In order to completely specify

**T**(

_{out}*k*), a third calibration signal

**R**(

_{3}*k*) must be added in. We can obtain another equation about

**T**(

_{out}*k*) as Eq. (9b)

**T**as Eq. (19)

_{out}**U**and

_{31}**V**are eigenvectors of [

_{31}**M**(

_{3}*k*) ·

**M**

_{1}^{−1}(

*k*)] and [

**R**(

_{3}*k*) ·

**R**

_{1}^{−1}(

*k*)], respectively, and

*β*

_{11}and

*β*

_{22}an be any complex value. Again,

**T**becomes any linear combination of two orthogonal bases

_{out}**T**=

_{out}*β*

_{11}

**B**+

_{1}*β*

_{22}

**B**, where

_{2}**B**and

_{1}**B**can be calculated by

_{2}**T**=

_{out}*α*

_{11}

**A**+

_{1}*α*

_{22}

**A**=

_{2}*β*

_{11}

**B**+

_{1}*β*

_{22}

**B**by linear algebra

_{2}*ij*) denotes the

*ith*row,

*jth*column element of the 2 × 2 matrix. By solving this equation for the constants

*α*

_{11},

*α*

_{22},

*β*

_{11}and

*β*

_{22},

**T**can then be constructed through either

_{out}**T**=

_{out}*α*

_{11}

**A**+

_{1}*α*

_{22}

**A**or

_{2}**T**=

_{out}*β*

_{11}

**B**+

_{1}*β*

_{22}

**B**. With

_{2}**T**(

_{out}*k*),

**T**(

_{in}*k*) can be easily solved using any calibration signal in Eq. (8).

## Acknowledgment

We would like to thank Dr. Wei Xu at Harvard Medical School for valuable discussion on analytic derivation part of this work. This project was supported by the Center for Biomedical OCT Research and Translation through Grant Number P41EB015903, awarded by the National Center for Research Resources and the National Institute of Biomedical Imaging and Bioengineering of the National Institutes of Health. This work was supported by NIH grant R01CA163528. This work was also supported by the National Research Foundation of Korea (NRF) grant 2012-0005633.

## References and links

**1. **M. R. Hee, D. Huang, E. A. Swanson, and J. G. Fujimoto, “Polarization-sensitive low-coherence reflectometer for birefringence characterization and ranging,” J. Opt. Soc. Am. B **9**, 903–908 (1992). [CrossRef]

**2. **J. F. de Boer, T. E. Milner, M. J. C. van Gemert, and J. S. Nelson, “Eye-length measurement by interferometry with partially coherent light,” Opt. Lett. **22**, 934–936 (1997). [CrossRef] [PubMed]

**3. **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]

**4. **A. F. Fercher, K. Mengedoht, and W. Werner, “Two-dimensional birefringence imaging in biological tissue using polarization sensitive optical coherence tomography,” Opt. Lett. **13**, 186–188 (1988). [CrossRef] [PubMed]

**5. **M. Pircher, C. K. Hitzenberger, and U. Schmidt-Erfurth, “Polarization sensitive optical coherence tomography in the human eye,” Prog. Retin. Eye Res. **30**, 431451 (2011). [CrossRef]

**6. **B.H. Park, C. Saxer, S.M. Srinivas, J.S. Nelson, and J.F. de Boer, “In vivo burn depth determination by high-speed fiber-based polarization sensitive optical coherence tomography,” J. Biomed. Opt. **6**, 474479 (2001). [CrossRef]

**7. **S. Nadkarni, M. C. Pierce, B. H. Park, J. F. de Boer, E. F. Halpern, S. L. Houser, B. E. Bouma, and G. J. Tearney, “Measurement of collagen and smooth muscle cell content in atherosclerotic plaques using polarization-sensitive optical coherence tomography,” J. Am. Coll. Cardiol. **49**, 1474–1481 (2007). [CrossRef] [PubMed]

**8. **C. E. Saxer, J.F. de Boer, B. H. Park, Y. Zhao, Z. Chen, and J.S. Nelson, “High-speed fiber-based polarization-sensitive optical coherence tomography of in vivo human skin,” Opt. Lett. **25**, 1355–1357 (2000). [CrossRef]

**9. **E. Z. Zhang and B. J. Vakoc, “Polarimetry noise in fiber-based optical coherence tomography instrumentation,” Opt. Express **19**, 16830–16842 (2011). http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-19-18-16830 [CrossRef] [PubMed]

**10. **S. Yun, G. Tearney, J. de Boer, and B. Bouma, “Removing the depth-degeneracy in optical frequency domain imaging with frequency shifting,” Opt. Express **12**, 4822–4828 (2004). http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-12-20-4822 [CrossRef] [PubMed]

**11. **W.Y. Oh, S. H. Yun, B. J. Vakoc, M. Shishkov, A. E. Desjardins, B. H. Park, J. F. de Boer, G. J. Tearney, and B. E. Bouma, “High-speed polarization sensitive optical frequency domain imaging with frequency multiplexing,” Opt. Express **16**, 1096–1103 (2008). http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-16-2-1096 [CrossRef] [PubMed]

**12. **H. Sunnerud, C. Xie, M. Karlsson, R. Samuelsson, and P. A. Andrekson, “A Comparison Between Different PMD Compensation Techniques,” J. Lightwave Technol. **20**, 368–378 (2002). [CrossRef]

**13. **K. H. Kim, B. H. Park, Y. Tu, T. Hasan, B. Lee, J. Li, and J. F. de Boer, “Polarization-sensitive optical frequency domain imaging based on unpolarized light,” Opt. Express **19**, 552–561 (2011). http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-19-2-552 [CrossRef] [PubMed]

**14. **M. Yamanari, S. Makita, and Y. Yasuno, “Polarization-sensitive swept-source optical coherence tomography with continuous source polarization modulation,” Opt. Express **16**, 5892–5906 (2008). http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-16-8-5892 [CrossRef] [PubMed]

**15. **W. Y. Oh, B. J. Vakoc, S. H. Yun, G. J. Tearney, and B. E. Bouma, “Single-detector polarization-sensitive optical frequency domain imaging using high-speed intra A-line polarization modulation,” Opt. Lett. **33**, 1330–1332 (2008). [CrossRef] [PubMed]

**16. **B. Baumann, W. Choi, B. Potsaid, D. Huang, J. S. Duker, and J. G. Fujimoto, “Swept source / Fourier domain polarization sensitive optical coherence tomography with a passive polarization delay unit,” Opt. Express **20**, 10229–10241 (2012). http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-20-9-10229 [CrossRef]

**17. **B. H. Park, M. C. Pierce, B. Cense, and J. F. de Boer, “Optic axis determination accuracy for fiber-based polarization-sensitive optical coherence tomography,” Opt. Lett. **30**, 2587–2589 (2005). [CrossRef] [PubMed]

**18. **Z. Lu and S. J. Matcher, “Absolute fast axis determination using non-polarization- maintaining fiber-based polarization-sensitive optical coherence tomography,” Opt. Lett. **37**, 1931–1933 (2012). [CrossRef] [PubMed]

**19. **http://www.mathworks.com/help/toolbox/optim/ug/fsolve.html

**20. **S. H. Yun, C. Boudoux, G. J. Tearney, and B. E. Bouma, “High-speed wavelength-swept semiconductor laser with a polygon-scanner-based wavelength filter,” Opt. Lett. **28**, 1981–1983 (2003). [CrossRef] [PubMed]

**21. **B. J. Vakoc, S. H. Yun, J. F. de Boer, G. J. Tearney, and B. E. Bouma, “Phase-resolved optical frequency domain imaging,” Opt. Express **13**, 5483–6593 (2005). http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-13-14-5483 [CrossRef] [PubMed]

**22. **J. Gentle, *Numerical Linear Algebra for Applications in Statistics* (Springer, 1998), chap. 2. [CrossRef]

**23. **J. F. de Boer, T. E. Milner, and J. S. Nelson, “Determination of the depth-resolved Stokes parameters of light backscattered from turbid media by use of polarization-sensitive optical coherence tomography,” Opt. Lett. **24**, 300–302 (1999). [CrossRef]

**24. **S. Jiao, G. Yao, and L. V. Wang, “Depth-resolved two-dimensional Stokes vectors of backscattered light and Mueller matrices of biological tissue measured with optical coherence tomography,” Appl. Opt. **39**, 6318–6324 (2000). [CrossRef]

**25. **B. H. Park, M. C. Pierce, B. Cense, and J. F. de Boer, “Jones matrix analysis for a polarization-sensitive optical coherence tomography system using fiber-optic components,” Opt. Lett. **29**, 2512–2514 (2004). [CrossRef] [PubMed]

**26. **S. Makita, M. Yamanari, and Y. Yasuno, “Generalized Jones matrix optical coherence tomography: performance and local birefringence imaging,” Opt. Express **18**, 854–876 (2010). http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-18-2-854 [CrossRef] [PubMed]

**27. **E. Götzinger, B. Baumann, M. Pircher, and C. K. Hitzenberger, “Polarization maintaining fiber based ultra-high resolution spectral domain polarization sensitive optical coherence tomography,” Opt. Express **17**, 22704–22717 (2009). http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-17-25-22704 [CrossRef]