Numerical compensation of system polarization mode dispersion in polarization-sensitive optical coherence tomography

Ellen Ziyi Zhang; Wang-Yuhl Oh; Martin L. Villiger; Liang Chen; Brett E. Bouma; Benjamin J. Vakoc

doi:10.1364/OE.21.001163

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.

Fig. 1 (a) A simplified fiber-based PS-OCT system based on a wavelength-swept source is presented. The system is constructed with a polarization-diverse receiver (unbalanced for simplicity) which in combination with heterodyne detection enables wavelength-by-wavelength polarimetry measurements of the sample arm light at D. (b) An optically equivalent description of the system in (a) where the distributed sample arm properties have been combined into three blocks INST1, SMPL, and INST2.

Download Full Size | PDF

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) = T_{out} (k) \cdot R (k) \cdot T_{in} (k) \cdot [\begin{matrix} 1 \\ 0 \end{matrix}],

where 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

M (k) = [C_{out} (k) \cdot T_{out} (k)] \cdot R (k) \cdot [T_{in} (k) \cdot C_{in} (k)] \cdot [\begin{matrix} 1 \\ 0 \end{matrix}],

where the first compensator (COMP1) is described by 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

M (k) = R (k) \cdot [\begin{matrix} 1 \\ 0 \end{matrix}] .

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.

Fig. 2 The influence of PMD induced by the instrument (elements INST1 and INST2) can be eliminated by inserting two PMD compensators (COMP1 and COMP2) configured to cancel the affects of INST1 and INST2 respectively.

Download Full Size | PDF

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.

Fig. 3 A PS-OCT system with simultaneous launching and independent detection of two polarization states. The source is split into two orthogonal polarizations by a first PBS, and each is encoded with a unique acousto-optic frequency shifter (AOFS). The signals are then recombined by a second PBS. System design simplified from ref [11].

Download Full Size | PDF

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:

M (k) = T_{out} (k) \cdot R (k) \cdot T_{in} (k) \cdot S (k),

where 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_{out} (k) \cdot M (k) \cdot C_{in} (k),

where C_in(k) and C_out(k) are our calibration matrices:

\begin{array}{l} C_{in} (k) = T_{in}^{- 1} (k) \\ C_{out} (k) = T_{out}^{- 1} (k) . \end{array}

Then we are able to remove the influence of T_in(k) and T_out(k) from the measurement system

C_{out} (k) \cdot M (k) \cdot C_{in} (k) = R (k) .

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₁(k), M₂(k), and M₃(k) can be expressed as

M_{1} (k) = T_{out} (k) \cdot R_{1} (k) \cdot T_{in} (k)

M_{2} (k) = T_{out} (k) \cdot R_{2} (k) \cdot T_{in} (k)

M_{3} (k) = T_{out} (k) \cdot R_{3} (k) \cdot T_{i n} (k),

where we have here explicitly allowed the Jones matrices R₁(k), R₂(k), and R₃(k) to vary with wavelength. Because these are calibrating signals from a known optical path, we assume that the Jones matrices R₁(k), R₂(k), and R₃(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):

[M_{2} (k) \cdot M_{1}^{- 1} (k)] = T_{out} (k) \cdot [R_{2} (k) \cdot R_{1}^{- 1} (k)] \cdot T_{out}^{- 1} (k)

[M_{3} (k) \cdot M_{1}^{- 1} (k)] = T_{out} (k) \cdot [R_{3} (k) \cdot R_{1}^{- 1} (k)] \cdot T_{out}^{- 1} (k),

where [R₂(k) · R₁⁻¹(k)] and [R₃(k) · R₁⁻¹(k)] are known a priori, and [M₂(k) · M₁⁻¹(k)] and [M₃(k) · M₁⁻¹(k)] are measured. If the eigen decomposition of [R₂(k) · R₁⁻¹(k)] and [R₃(k) · R₁⁻¹(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.

Fig. 4 The optical block design of an OCT system providing three independent signal paths/Jones matrices. From these signals, the required compensating matrices C_in(k) and C_out(k) can be derived analytically.

Download Full Size | PDF

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 60kHz over 140nm sweeping range centered at 1270nm and yielded 6 μm axial resolution in tissue. The carrier frequencies of two simultaneous polarizations were centered at 25MHz and 75MHz by placing a 50MHz frequency shifter in the reference arm and two identical 25MHz 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 10mW, and the back-reflected light was directed to a polarization-diverse receiver.

Fig. 5 Experimental setup of frequency multiplexing PS-OFDI system (FS: frequency shifter; FS1: 50MHz, FS2: 25MHz upshift, FS3: 25MHz downshift; LP: linear polarizer; PC: polarization controller; BBS: broad beamsplitter; PBS: polarization beamsplitter; BR: balanced receiver; A/D: analog-to-digital converter).

Download Full Size | PDF

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 100MHz (the highest RF frequency in our detection bandwidth), a 1cm mismatch induces 0.02radians 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₁(k), R₂(k), and R₃(k) are given by

\begin{array}{l} R_{1} (k) = [\begin{matrix} 1 & 0 \\ 0 & 1 \end{matrix}] \\ R_{2} (k) = [\begin{matrix} 1 & 0 \\ 0 & e^{i \cdot 2 Γ (k)} \end{matrix}] \\ R_{3} (k) = \frac{1}{2} \cdot [\begin{matrix} 1 + e^{i \cdot 2 Γ (k)} & 1 - e^{i \cdot 2 Γ (k)} \\ 1 - e^{i \cdot 2 Γ (k)} & 1 + e^{i \cdot 2 Γ (k)} \end{matrix}], \end{array}

where

Γ (k) = \frac{π}{2} \cdot \frac{k}{k_{°}}

is the wavelength dependent retardance. Here the central wavelength of the waveplate is 1300nm, and with the laser swept range we can calculate Γ(k).

Fig. 6 The microscope used to provide three calibrating signals (BBS: broad beamsplitter; W1: quarter waveplate (QWP) oriented at 0°; W2: QWP oriented at 45°)

Download Full Size | PDF

5.3. Extracting depth-encoded signals M₁(k), M₂(k), and M₃(k) from a single A-line

The digitized fringe data were first demodulated about their 25MHz and 75MHz 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₁(k), M₂(k), and M₃(k) were constructed.

5.4. Normalizing the measured signals M₁(k), M₂(k), and M₃(k)

By constructing the measured signals M₁(k), M₂(k), and M₃(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:

\begin{array}{l} M_{1} (k) = c_{1} (k) \cdot T_{out} (k) \cdot R_{1} (k) \cdot T_{in} (k) \cdot K (k) \\ M_{2} (k) = c_{2} (k) \cdot T_{out} (k) \cdot R_{2} (k) \cdot T_{in} (k) \cdot K (k) \\ M_{3} (k) = c_{3} (k) \cdot T_{out} (k) \cdot R_{3} (k) \cdot T_{in} (k) \cdot K (k) . \end{array}

Here, the scalar factors c₁(k), c₂(k), and c₃(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₁(k), c₂(k), and c₃(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₁(k), c₂(k), c₃(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

(\frac{c_{1} (k)}{c_{2} (k)}) \cdot [M_{2} (k) \cdot M_{1}^{- 1} (k)] = T_{out} (k) \cdot [R_{2} (k) \cdot R_{1}^{- 1} (k)] \cdot T_{out}^{- 1} (k)

(\frac{c_{1} (k)}{c_{3} (k)}) \cdot [M_{3} (k) \cdot M_{1}^{- 1} (k)] = T_{out} (k) \cdot [R_{3} (k) \cdot R_{1}^{- 1} (k)] \cdot T_{out}^{- 1} (k),

where the matrix K(k) has been cancelled. Here, [R₂(k) · R₁⁻¹(k)] and [R₃(k) · R₁⁻¹(k)] are known a priori, and [M₂(k) · M₁⁻¹(k)] and [M₃(k) · M₁⁻¹(k)] are measured, but the scalars c₁(k), c₂(k), and c₃(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

(\frac{c_{1} (k)}{c_{2} (k)}) \cdot [M_{2} (k) \cdot {M_{1}}^{- 1} (k)]

is a similarity transformation of [R₂(k) · R₁⁻¹(k)]. Therefore, each has the same eigenvalues upon eigen decomposition [22]. Stated alternatively, if the eigen decomposition of [R₂(k) · R₁⁻¹(k)] is written by

R_{2} (k) \cdot R_{1}^{- 1} (k) = V_{21} (k) \cdot D_{21} (k) \cdot V_{21}^{- 1} (k),

where V₂₁(k) contains the eigenvectors and D₂₁(k) is a diagonal matrix containing the eigenvalues, then the eigendecomposition of [M₂(k) · M₁⁻¹(k)] (absent the scalar factors) can be written by

M_{2} (k) \cdot M_{1}^{- 1} (k) = U_{21} (k) \cdot d_{21} (k) \cdot U_{21}^{- 1} (k),

where again U₂₁(k) contains the eigenvectors, d₂₁(k) is a diagonal matrix containing the eigenvalues, and D₂₁(k) and d₂₁(k) differ only because of the presence of the scalar factors c₁(k), c₂(k), and c₃(k). We can therefore calculate the lhs (left-hand side) of Eq. (12a) by first calculating the eigendecomposition of [M₂(k) · M₁⁻¹(k)] to get U₂₁(k), and then using the eigenvalues matrix D₂₁(k):

(\frac{c_{1} (k)}{c_{2} (k)}) \cdot M_{2} (k) \cdot M_{1}^{- 1} (k) = U_{21} (k) \cdot D_{21} (k) \cdot U_{21}^{- 1} (k) .

A similar process provides the lhs of Eq. (12b). With values for the lhs of Eq. (12a) and (12b), we can now solve this set of equations for 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

\begin{array}{l} C_{in} (k) = {[T_{in} (k) \cdot K (k)]}^{- 1} \\ C_{out} (k) = T_{out}^{- 1} (k) . \end{array}

It is worth noting that while this method assumes a known functional form for R₁(k), R₂(k), and R₃(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₂₁(k) and D₃₁(k)) or the eigenvectors (V₂₁(k) and V₃₁(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₂(k) · R₁⁻¹ (k)] and [R₃(k) · R₁⁻¹(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₁(k), R₂(k), and R₃(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

The normalizing procedure of Section 5.4 describes a method for calculating T_out(k) which includes an eigendecomposition of the measured data [M₂(k) · M₁⁻¹(k)] and [M₃(k) · M₁⁻¹(k)]. This eigendecomposition is applied for each wavelength, resulting in a set of eigenvectors U₂₁(k) and U₃₁(k). There is an order ambiguity, however, in how these eigenvector matrices are constructed; their columns (and associated eigenvalues in the matrices d₂₁(k) and d₃₁(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₂₁(k) and U₃₁(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₂₁(k) (or U₃₁(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.04ps per path) and high-PMD (0.08ps 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.08ps 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.

Fig. 7 Validation of compensation method by imaging a mirror sample, normalized Stokes vectors are plotted along swept wavelengths on Poincaré Sphere: (a) Mirror signal without and with compensation of both polarization states; (b) With another configuration of input and output polarization controllers, mirror signal without and with compensation.

Download Full Size | PDF

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.08ps (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.08ps patchcord and without PMD compensation (Fig. 8(c)), 0.0076 deg/μm with the 0.08ps patchcord and with PMD compensation (Fig. 8(d)), and 0.0090 deg/μm without the patchcord and without PMD compensation (Fig. 8(b)).

Fig. 8 Numerical PMD compensation reduces noise in PS-OCT images of local birefringence in intralipid phantom. (a) Intensity image from one input polarization state; (b) local birefringence image without PMF patchcord and without numerical PMD compensation; (c) local birefringence image with 0.08ps PMF patchcord and without numerical PMD compensation; (d) local birefringence image with 0.08ps PMF patchcord and with numerical PMD compensation. Local birefringence was calculated across an axial extent of 6 pixels. Colorbar: 0 − 1.3 deg/μm

Download Full Size | PDF

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.04ps and 0.08ps patchcord respectively and without PMD compensation are presented. In Fig. 9(d) and 9(f), the images with the 0.04ps and 0.08ps 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.

Fig. 9 Numerical PMD compensation reduces noise in PS-OCT images of local birefringence of muscle/tendon. (a) Intensity image from one input polarization state; (b) baseline local birefringence without PMF patchcords and without PMD compensation; (c) local birefringence image with 0.04ps PMF patchcord and without numerical PMD compensation; (d) local birefringence image with 0.04ps PMF patchcord and with numerical PMD compensation; (e) local birefringence image with 0.08ps PMF patchcord and without numerical PMD compensation; (f) local birefringence image with 0.08ps PMF patchcord and with numerical PMD compensation. Birefringence values of tendon were calculated to be 1.07 deg/μm, 1.10 deg/μm and 1.10 deg/μm from ROI(t) in Fig. 9(b), Fig. 9(d) and Fig. 9(f), respectively. Birefringence values of muscle were calculated to be 0.38 deg/μm, 0.39 deg/μm and 0.39 deg/μm from ROI(m) in Fig. 9(b), Fig. 9(d) and Fig. 9(f), respectively. Local birefringence was calculated across an axial extent of 6 pixels. ROI: 40*40 pixels. Colorbar: 0 − 2 deg/μm

Download Full Size | PDF

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

[M_{2} (k) \cdot M_{1}^{- 1} (k)] = T_{out} (k) \cdot [R_{2} (k) \cdot R_{1}^{- 1} (k)] \cdot T_{out}^{- 1} (k),

where [R₂(k) · R₁⁻¹(k)] is known a priori from calibration signals, and [M₂(k) · M₁⁻¹(k)] is measured. Note that [M₂(k) · M₁⁻¹(k)] is a similarity transformation of [R₂(k) · R₁⁻¹(k)], their eigenvalues are equal upon eigendecomposition [22]:

[U_{21} \cdot D \cdot U_{21}^{- 1}] = T_{out} \cdot [V_{21} \cdot D \cdot V_{21}^{- 1}] \cdot T_{out}^{- 1},

where U₂₁ and V₂₁ are eigenvectors of [M₂(k) · M₁⁻¹(k)] and [R₂(k) · R₁⁻¹(k)], respectively, while D is their eigenvalue matrix. Equation (17) can be rearranged to write

\begin{array}{l} D & = U_{21}^{- 1} \cdot T_{out} \cdot [V_{21} \cdot D \cdot V_{21}^{- 1}] \cdot T_{out}^{- 1} \cdot U_{21} \\ = [U_{21}^{- 1} \cdot T_{out} \cdot V_{21}] \cdot D \cdot [V_{21}^{- 1} \cdot T_{out}^{- 1} \cdot U_{21}] \\ = [U_{21}^{- 1} \cdot T_{out} \cdot V_{21}] \cdot D \cdot {[U_{21}^{- 1} \cdot T_{out} \cdot V_{21}]}^{- 1} . \end{array}

Since D is a diagonal matrix, it can be shown that a necessary and sufficient condition for this equation to be solved is that [U₂₁⁻¹ · T_out · V₂₁] is diagonal. With this requirement, we can write a general form for [U₂₁⁻¹ · T_out · V₂₁] as

[U_{21}^{- 1} \cdot T_{out} \cdot V_{21}] = [\begin{matrix} α_{11} & 0 \\ 0 & α_{22} \end{matrix}],

with α₁₁ and α₂₂ can be any complex value. The analytic solution for T_out becomes any linear combination of two orthogonal bases T_out = α₁₁A₁ + α₂₂A₂, where A₁ and A₂ can be calculated by

{\begin{matrix} [U_{21}^{- 1} \cdot A_{1} \cdot V_{21}] = [\begin{matrix} 1 & 0 \\ 0 & 0 \end{matrix}] \\ [U_{21}^{- 1} \cdot A_{2} \cdot V_{21}] = [\begin{matrix} 0 & 0 \\ 0 & 1 \end{matrix}] . \end{matrix}

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₃(k) must be added in. We can obtain another equation about T_out(k) as Eq. (9b)

[M_{3} (k) \cdot M_{1}^{- 1} (k)] = T_{out} (k) \cdot [R_{3} (k) \cdot R_{1}^{- 1} (k)] \cdot T_{out}^{- 1} (k) .

Following the same steps, we arrive a similar form for T_out as Eq. (19)

[U_{31}^{- 1} \cdot T_{out} \cdot V_{31}] = [\begin{matrix} β_{11} & 0 \\ 0 & β_{22} \end{matrix}],

where U₃₁ and V₃₁ are eigenvectors of [M₃(k) · M₁⁻¹(k)] and [R₃(k) · R₁⁻¹ (k)], respectively, and β₁₁ and β₂₂ an be any complex value. Again, T_out becomes any linear combination of two orthogonal bases T_out = β₁₁B₁ +β₂₂B₂, where B₁ and B₂ can be calculated by

{\begin{matrix} [U_{31}^{- 1} \cdot B_{1} \cdot V_{31}] = [\begin{matrix} 1 & 0 \\ 0 & 0 \end{matrix}] \\ [U_{31}^{- 1} \cdot B_{2} \cdot V_{31}] = [\begin{matrix} 0 & 0 \\ 0 & 1 \end{matrix}] . \end{matrix}

Combing Eq. (19) and Eq. (21), it turns out to solve T_out = α₁₁A₁ + α₂₂A₂ = β₁₁B₁ + β₂₂B₂ by linear algebra

[\begin{matrix} A_{1}^{(11)} & A_{2}^{(11)} & - B_{1}^{(11)} & - B_{2}^{(11)} \\ A_{1}^{(12)} & A_{2}^{(12)} & - B_{1}^{(12)} & - B_{2}^{(12)} \\ A_{1}^{(21)} & A_{2}^{(21)} & - B_{1}^{(21)} & - B_{2}^{(21)} \\ A_{1}^{(22)} & A_{2}^{(22)} & - B_{1}^{(22)} & - B_{2}^{(22)} \end{matrix}] \cdot [\begin{matrix} α_{11} \\ α_{22} \\ β_{11} \\ β_{22} \end{matrix}] = 0,

where the superscript (ij) denotes the ith row, jth column element of the 2 × 2 matrix. By solving this equation for the constants α₁₁, α₂₂, β₁₁ and β₂₂, T_out can then be constructed through either T_out = α₁₁A₁ + α₂₂A₂ or T_out = β₁₁B₁ + β₂₂B₂. With 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]

Numerical compensation of system polarization mode dispersion in polarization-sensitive optical coherence tomography

Abstract

1. Introduction

2. PS-OCT system model

3. Physical PMD compensation