Fast recursive algorithm for broadband APFs using complex cepstrums

P. B. Phua; E. P. Ippen

doi:10.1364/OPEX.12.004896

1. Introduction

All-Pass Filters (APFs) integrated on optical integrated circuits are useful for optical equalization applications such as compensation of Chromatic Dispersion [1,2], Polarization Mode Dispersion (PMD) [3,4], and Polarization Dependent Loss (PDL) [5]. This approach is promising due to its compactness and potential economy of production, especially when the number of channels to be compensated is large. Inline broadband compensation over terahertz bandwidth is attractive since it can achieve simultaneous compensation for several channels without the need of wavelength demultiplexing cum multiplexing. For such broadband applications, a large number of APFs are involved to approximate the desired phase profiles. Thus the algorithm to find the optimum parameters of these APFs needs to be efficient and fast. This is especially true in the case of compensating for PMD and PDL, which are dynamically changing with time. Optimization algorithms, such as the nonlinear fit used in Ref. [3, 4], tend to get trapped in local optima when the number of search parameters is large. For this reason, they are efficient only when a small number of APFs is involved (<5).

Fast fitting algorithms that based on recursive equations are desired. Ref. [6–8] presented a simple and fast recursive method for designing a digital all-pass filter satisfying a given group delay specification using complex cepstrum [9]. The term cepstrum was introduced by Bogert et al. in 1963 [10] and has come to be accepted terminology for the inverse Fourier transform of the logarithm of the power spectrum of a signal. Cepstrum techniques have been applied extensively in speech analysis. The fundamental principle of these algorithms [6–8] is based on the relationship between the group delay of the minimum-phase denominator of a stable APF and the complex cepstral coefficients arising from its discrete Hilbert transform. The main difference among the three algorithms lies in how the cepstral coefficients are computed. Ref. [6] uses four Fast-Fourier Transforms (FFT) while [7] uses one FFT to approximate the cepstral coefficients in a least-square sense. Ref. [8] requires no FFT and uses a weighted-least-square approximation.

Fig. 1. A common building block in the broadband PMD compensator

Download Full Size | PDF

A broadband compensation scheme for PMD [3,4] and PDL [5] have been proposed. A common building block of these architectures is a frequency-dependent polarization rotator as shown in Fig. 1. The light is polarization beam-split into horizontal and vertical polarizations. Different frequency-dependent phase profiles for the two polarization arms, Φ_H(ω) for the horizontal and Φ_V(ω) for the vertical, are created using cascaded APFs. When the two polarizations are combined using a polarization beam combiner, this creates frequency-dependent rotation θ(ω) about {1,0,0} in Stokes space. This rotation profile θ(ω) corresponds to the differential phase profile of the two arms (i.e. θ(ω) = Φ_V(ω)-Φ_H(ω)). The first derivative of θ(ω) gives the Differential Group Delay (DGD) of the frequency-dependent polarization rotator. Fitting algorithms proposed in [6–8] are not optimized for fitting such a differential phase profile θ(ω) for two reasons. Firstly, these algorithms assume pairs of poles that are complex conjugate to one another. Thus, the APFs have to mirror-image one another about the center frequency resulting in the group delay response of even-symmetry about the center frequency. This reduces the working frequency range to half of the Free Spectral Range (FSR) of the filters and doubles the number of APFs required. One common implementation of APFs in an optical integrated circuit is a ring resonator coupled to a waveguide [1,2]. The FSR of the APF is inversely proportional to the ring radius. Since rings with small radii tend to have more bending loss and are more difficult to fabricate, we would like to utilize as much of the FSR as possible. Secondly, these algorithms concentrate on the higher order dispersion without accounting for the linear frequency dependence of the response. In applications such as chromatic dispersion compensation, this may not be an issue since this linear dependence corresponds to a group delay which does not distort the signal. However, in PMD and PDL compensation, the difference of the linear dependence in the two polarization arms, shown in Fig. 1, creates a first order PMD. This will introduce polarization-dependent signal distortion if it is unaccounted for. In this paper, we show how to rectify these two issues, and enhance the algorithms [7, 8] to efficiently fit the rotation angle profile of the frequency-dependent polarization rotator shown in Fig. 1. This algorithm is equally valid for fitting the differential phase profile required for the frequency-dependent variable attenuator proposed in [5].

2. Group delay of APFs, minimum-phase denominator polynomial and DGD

The transfer function, in the z-transform domain, of an all-pass filter of order N is given by

H (z) = z^{- N} \frac{N (z)}{D (z)} = z^{- N} \frac{\sum_{n = 0}^{N} a_{n}^{*} z^{n}}{\sum_{n = 0}^{N} a_{n} z^{- n}}

with a_o = 1. There are N poles and zeros. The zeros of H(z) are reciprocal complex conjugates of its poles. In [6–8], a_n is taken as a real number resulting in a pair of poles that are complex conjugates of one another. This limits the fitting algorithm to accept only group delay that has even symmetry about the center frequency. In general, a_n is complex. From Eq. (1), when z = e^jω ,

H (e^{j ω}) = e^{- j ω N} \frac{\sum_{n = 0}^{N} a_{n}^{*} e^{j ω n}}{\sum_{n = 0}^{N} a_{n} e^{- j ω n}}

where ω is the normalized frequency given by $\frac{2 π (f - f_{o})}{FSR}$ and f_o is the center optical frequency of the band of interest. Since the numerator N(e^jω ) is the complex conjugate of the denominator D(e^jω ), H(e^jω ) has unity magnitude response for all frequencies. Moreover, the phase response of the overall filter H(e^jω ) is given by

ϕ_{H} (ω) = - N ω + ϕ_{N} (ω) - ϕ_{D} (ω) = - N ω - 2 ϕ_{D} (ω)

where ϕ_N (ω) and ϕ_D (ω) denote the phase response of the numerator N(e^jω ) and the phase response of the denominator D(e^jω ), respectively. The normalized group delay of the overall filter is given by

τ_{H} (ω) = \frac{- d ϕ_{H} (ω)}{d ω} = N + 2 ϕ_{D}^{'} (ω)

where $ϕ_{D}^{'} (ω) = \frac{d ϕ_{D} (ω)}{d ω} = - τ_{D} (ω)$ is the negative group delay of the denominator D(e^jω ). Hence, the group delay response of the denominator can be obtained from the group delay response of the all-pass filter using the relation

τ_{D} (ω) = \frac{- τ_{H} (ω)}{2} + \frac{N}{2}

The coefficients a_n of the denominator D(z) and hence the filter H(z) are completely determined from τ_D (ω). From Eq. (1), for a stable APF, H(z) has all its poles within the unit circle in the complex z-plane. This implies that the denominator D(z) is a minimum phase function. A function is defined to be minimum phase when all its zeros and poles lies within the unit circle [9]. A minimum phase function is known to have zero average group delay over the FSR. From Eq. (5), this implies that the average group delay of a N ^th -order APFs over a FSR is fixed at N.

To approximate a desired rotation angle profile θ(ω) of the frequency-dependent polarization rotator shown in Fig. 1, we cascade N_Hor APFs in the horizontal polarization arm to generate the phase response Φ_Hor(ω) and N_Vert APFs in the vertical polarization arm to generate the phase response Φ_Vert(ω) so that

Φ_{Vert} (ω) - Φ_{Hor} (ω) = θ (ω)

Differentiating Eq. (6) with respect to ω, we obtain the corresponding normalized group delay as

τ_{Vert} (ω) - τ_{Hor} (ω) = τ_{DGD} (ω)

which is the Differential Group Delay (DGD) of the frequency-dependent polarization rotator. We assume all APFs have the same FSR. Since the average group delay of an N ^th order APF over the FSR is fixed at N, the average DGD over the FSR is fixed at N_Vert - N_Hor which is an integer. On the other hand, since PMD and PDL are stochastic processes, the required rotation angle profile θ(ω) for compensation is dynamically changing with an arbitrary non-integer value of average DGD τ̅_DGD over the frequency range of interest ∆f _int. To overcome this “incompatibility”, the FSR of the APFs must be larger than ∆f _int. We then express the required average DGD over ∆f _int as τ̅_DGD = N_diff + ∆τ̅_DGD where N_diff is an integer. Since the frequency range outside ∆f _int has no impact on the compensation, we appropriately extrapolate τ_DGD (ω) outside the frequency range ∆f _int to cancel ∆τ̅_DGD , thus rounding the average value of τ_DGD (ω) over the FSR to an integer value of N_diff . To account for this N_diff average DGD, we change the order of the APFs accordingly so that N_Vert - N_Hor = N_diff . One way to implement this is to have the same number of APFs in both polarization arms, and to “de-activate” N_diff APFs in the appropriate arms, depending on the sign of N_diff . From Eq. (7), if N_diff is a positive (negative) integer, we “de-activate” N_diff APFs from the horizontal (vertical) polarization arm, and vice versa. To “deactivate” an APF in an optical integrated circuit, we shift its resonant frequency out of ∆f _int and decrease the power coupling of the ring resonator to the waveguide so that it does not contribute any group delay within ∆f _int. The FSR of the APF is inversely proportional to the ring radius, and we would like to utilize the most of FSR as ∆f _int, in order to reduce the ring bending loss and ease its fabrication. In our simulations, we typically choose the FSR to be ~10% larger than ∆f _int. After accounting for the linear frequency dependence of the rotation angle, we now focus on fitting its higher order differential dispersion. We have the freedom to choose Φ_Hor(ω) and Φ_Vert(ω) so that the isotropic dispersion of the frequency dependent polarization rotator is a group delay.

τ_{Vert} (ω) + τ_{Hor} (ω) = N_{Vert} + N_{Hor}

Using Eq. (7) and Eq. (8), we get

τ_{Vert} (ω) = \frac{N_{Vert} + N_{Hor}}{2} + \frac{τ_{DGD} (ω)}{2}

and

τ_{Hor} (ω) = \frac{N_{Vert} + N_{Hor}}{2} - \frac{τ_{DGD} (ω)}{2}

which are the corresponding desired group delays for the respective polarization arms. We show in the following how to fit the desired group delay for the vertical polarization and the same procedures apply to the horizontal. For a desired group delay function τ_Vert (ω), we use Eq. (5) to determine the corresponding group delay response τ_DV (ω) of its minimum phase denominator as

τ_{DV} (ω) = \frac{- τ_{Vert} (ω)}{2} + \frac{N_{Vert}}{2}

We then decompose it into a sum of even function $τ_{DV}^{even}$ (ω) and odd function ro $τ_{DV}^{odd}$ (ω) of the form,

τ_{DV}^{even} (ω) = \frac{τ_{DV} (ω) + τ_{DV} (- ω)}{2}

τ_{DV}^{odd} (ω) = \frac{τ_{DV} (ω) - τ_{DV} (- ω)}{2}

3. Relationship between minimum-phase response and cepstral coefficients

Let D(ω) be a minimum phase frequency response. Since D(ω) is minimum phase with all its zeros and poles inside the unit circle, its natural logarithm D̂(ω) = ln D(ω) also has all its poles inside the unit circle. Thus D̂(ω) is a stable and causal frequency response. In addition, the inverse Fourier transform of D̂(ω) (also known as the cepstrum) yields a causal sequence of cepstral coefficients c(k). In terms of c(k), we can express

\hat{D} (ω) = \ln D (ω) = c (0) + \sum_{k = 1}^{\infty} c (k) e^{- j k ω}

From D(ω) = |D(ω)|e ^jϕ(ω), we see that the imaginary part of D̂(ω) gives the phase term as

ϕ (ω) + 2 υ π = - \sum_{k = 1}^{\infty} Re [c (k)] \sin k ω + \sum_{k = 1}^{\infty} Im [c (k)] \cos k ω

where v is an integer. Its corresponding group delay is

τ (ω) = \sum_{k = 1}^{\infty} k Re [c (k)] \cos k ω + \sum_{k = 1}^{\infty} k Im [c (k)] \sin k ω

, and the real part of D̂(ω) is given by

\ln ∣ D (ω) ∣ = c (0) + \sum_{k = 1}^{\infty} Re [c (k)] \cos k ω + \sum_{k = 1}^{\infty} Im [c (k)] \sin k ω

From Eq. (16), we can see that the real part of c(k) gives group delay of even symmetry and the imaginary part gives group delay of odd symmetry. In [6–8], c(k) is real since its group delay is only of even symmetry. For a given desired group delay τ_Vert (ω), we determine the corresponding group delay τ_DV (ω) of its minimum phase denominator response using Eq. (11). Then we decompose it into a sum of even group delay $τ_{DV}^{even}$ (ω) and odd group delay $τ_{DV}^{odd}$ (ω) using Eq. (12) and Eq. (13). Using Eq. (16), the even group delay function determines the real part

τ_{DV}^{even} (ω) = \sum_{k = 1}^{\infty} k Re [c (k)] \cos k ω

while the odd delay determines the imaginary part of the complex cepstral coefficients

τ_{DV}^{odd} (ω) = \sum_{k = 1}^{\infty} k Im [c (k)] \sin k ω

To compute the cepstral coefficients c(k) from Eq. (18) and (19), we can either use an inverse Fast Fourier transform (IFFT) [7] or compute them using a weighted-least-square approximation [8]. From the computed complex cepstral coefficients c(k) , we can now derive the coefficients a_n of the denominator D(z) using the recursive nonlinear difference equation [9] given by

a_{n} = \sum_{k = 0}^{n} (\frac{k}{n}) c (k) a_{n - k} n > 0

where c(0) = 0 and a ₀ = 1. To prove equation (20), we first differentiate D̂(z) = ln(D(z)) with respect to z so that

\frac{d \hat{D} (z)}{d z} = \frac{1}{D (z)} \frac{d D (z)}{d z}

Using the z-Transform pairs

n x [n] \overset{z - Transform}{\leftrightarrow} - z \frac{d X (z)}{d z}

and

x_{1} [n] * x_{2} [n] \overset{z - Transform}{\leftrightarrow} X_{1} (Z) . X_{2} (Z)

the inverse z-transform of Eq. (21) gives the recursive nonlinear difference equation in Eq. (20). Once the sequence a_k is known, the parameters of the APFs are determined by finding the roots of the denominator D(z) . The magnitude of each root gives the reflection coefficient of each APF and its argument gives the resonant frequency of the APF.

4. Summary of algorithm

The implementation procedure can be summarized as follows: let τ_DGD (ω) be the desired differential group delay function over a frequency range ∆f _int .

Compute average τ_DGD (ω) over ∆f _int and express it in terms of an integer N_diff and ∆τ̅_DGD . N_diff determines the required difference in orders of APFs used in the vertical and horizontal polarization arm. Since the frequency range outside does not affect the compensation, we extrapolate τ_DGD (ω) outside ∆f _int to cancel ∆τ̅_DGD.
From Eq. (9) and (10), compute the respective desired group delay for each polarization arm: τ_Vert (ω) and τ_Hor (ω).
For τ_Vert (ω) in the vertical polarization arm, use Eq. (11) to compute the corresponding group delay of the minimum phase denominator τ_DV (ω) .
Perform an even-odd decomposition of τ_DV (ω) using Eq. (12) and (13): $τ_{DV}^{even}$ (ω) and $τ_{DV}^{odd}$ (ω).
Compute the complex cepstral coefficients c(k) using Eq. (18) and (19). The real part of c(k) is computed using the even function $τ_{DV}^{even}$ (ω) while the imaginary part using the odd function $τ_{DV}^{odd}$ (ω) . This can be implemented using an IFFT [7] or computed in a weighted-least-square sense using a weighting function as in [8]
Compute the denominator coefficients a_n using the recursive relation in Eq. (20)
Compute the roots of the denominator D(z). The magnitude of each root gives the reflection coefficient of each APF and its argument gives the resonant frequency of the APFs.
Repeat Step 3 to 7 for the horizontal polarization arm.

Fig. 2. A random sample of the rotation angle profile θ(ω) required for one of the frequency-dependent polarization rotators used in the broadband PMD compensator. The desired rotation angle profile of θ(ω) is shown by the solid curve while the profile approximated by the APFs is shown by the dashed-curve. The number of APFs used for each polarization arms is 10 for (a), 20 for (b) and 30 for (c).

Download Full Size | PDF

5. Simulation

Figure 2 shows a sample of the rotation angle profile θ(ω) required for one of the frequency-dependent polarization rotators used in the broadband PMD compensator presented in Ref. [4]. The fiber to be compensated has a mean DGD of 6.75 ps. The desired rotation angle profile of θ(ω) is shown by the solid curve while the profile approximated by the APFs is shown by the dashed-curve. The cepstral coefficients were computed using IFFT as in [7]. The number of APFs used for each polarization arms is 10 for Fig. 2(a), 20 for Fig. 2(b) and 30 for Fig. 2(c). As expected, the approximation improves as the number of APFs increases. The FSR of the APFs is 1.1 THz while the frequency range of interest ∆f _int is 1 THz. This broadband compensator is capable of covering ten 100 GHz-spaced WDM channels simultaneously. Rings with such Tera-Hz-scale FSRs have been demonstrated in [11].

Since PMD is a stochastic process, the required rotation angle profile θ(ω) for compensation is dynamically changing. To gain statistics on the performance of the APFs in approximating the rotation angle profile, we generated an ensemble of 1000 random fibers by cascading 20 randomly oriented birefringence sections. The mean DGD of the ensemble is of 6.75 ps. For each fiber, the input signal is a 40 Gbit/s RZ pseudo-random bit sequence (2⁶-1) of Gaussian pulses of 10 ps (FHWM) pulse-width. PMD compensation is carried out at the fiber output. The broadband PMD compensator is comprised of three stages of frequency-dependent polarization rotators and the method to synthesize the required rotation angle profile θ(ω) for each stage is described in Ref. [4]. We then find the optimum parameters of the APFs to approximate these required rotation angle profiles using the recursive algorithm proposed in this paper. The number of APFs used is 40 for each arm of the polarization rotators. All APFs have the same FSR of 1 THz and the frequency range of interest ∆f _int is 0.9 THz. In order to show that PMD was compensated across the whole broadband frequency range, the center wavelength of the Gaussian pulses was randomly chosen within this 0.9 THz bandwidth. Figure 3 shows the cumulative probability distribution of the Bit Error Rate (BER) curve with and without compensation. The power was set at 3 dB above the power that provides a BER of 10^-9. The dashed curve is for the uncompensated case while the solid curve is for the compensated case. There is significant improvement in the signal quality after compensation, thereby illustrating the capability of the recursive algorithm in finding for the optimum fitting parameters of such large number of APFs.

Fig. 3. The cumulative probability distribution of the Bit Error Rate (BER) curve with and without PMD compensation.

Download Full Size | PDF

6. Conclusion

We present in this paper a fast recursive algorithm that uses the complex cepstrum to find the optimum fitting parameters of APFs in approximating arbitrary profile of differential group delay. It can handle a large number of APFs (>50) and thus it is suitable for broadband PMD and PDL compensation.

References

1. C. Madsen and J. Zhao, Optical Filter Design and Analysis: A Signal Processing Approach. New York, Wiley (1999).

2. C.K. Madsen, G. Lenz, A.J. Bruce, M.A. Capuzzo, L.T. Gomez, and R.E. Scotti, “Integrated All-Pass Filters for Tunable Dispersion and Dispersion Slope Compensation,” IEEE Photon. Technol. Lett. 11, 1623–1625 (1999). [CrossRef]

3. C.K. Madsen and P. Oswald, “Optical filter architecture for approximating any 2×2 unitary matrix,” Opt. Lett. 28, 534–536 (2003). [CrossRef] [PubMed]

4. P.B. Phua, H. A. Haus, and E. P. Ippen, “All-Frequency PMD Compensator In Feed-Forward Scheme,” J.of Lightwave Technol. 22, 1280–1289 (2004). [CrossRef]

5. P.B. Phua and E. P. Ippen, “A Deterministic Broadband Polarization-Dependent-Loss Compensator,” To appear in J. Lightwave Technol..

6. G.R. Reddy and M.N.S. Swamy, “Digital All-Pass Filter Design Through Discrete Hilbert Transform,” Proc. IEEE Int. Conf. Acoustics, Speech Signal Processing , 646–649 (1990)

7. K. Rajamani and Y.S. Lai, “A Novel Method for Designing Allpass Digital Filters,” IEEE Signal Processing Lett. 6, 207–209 (1999) [CrossRef]

8. G.J. Dolecek and J.D. Carmona, “Digital All-Pass filter design method based on Complex Cepstrums,” Electron. Lett. 39, 695–697 (2003). [CrossRef]

9. A.V. Oppenheim and R. W. Schafer, Discrete-Time Signal Processing, Prentice Hall (1989)

10. B.P. Bogert, M.J.R. Healy, and J.W. Tukey, “The Quefrency Analysis of Time Series for Echoes: Cepstrum, Pseudoautocovariance, Cross-Cepstrum, and Saphe Cracking,” Proc. Symposium Time Series Analysis,M. Rosenblatt, Ed., John Wiley and Sons, New York, 209–243 (1963).

11. T. Barwicz, M.A. Popovic, P.T. Rakich, M.R. Watts, H.A. Haus, E.P. Ippen, and H.I. Smith, “Microring-resonator-based add-drop filters in SiN: fabrication and analysis,” Opt. Express 12, 1437–1442 (2004). [CrossRef] [PubMed]

Fast recursive algorithm for broadband APFs using complex cepstrums

Abstract

1. Introduction

2. Group delay of APFs, minimum-phase denominator polynomial and DGD

3. Relationship between minimum-phase response and cepstral coefficients

4. Summary of algorithm

5. Simulation

6. Conclusion

References

Cited By

Figures (3)

Equations (23)

Optics Express