Chromatic dispersion equalization FIR filter design based on discrete least-squares approximation

Dawei Wang; Hao Jiang; Guowei Liang; Qianxin Zhan; Zikang Su; Zhaohui Li; Zhaohui Li

doi:10.1364/OE.430955

1. Introduction

Chromatic dispersion (CD) is a major effect causing severe signal distortions in optical fiber communication systems. The modern coherent optical communication technology utilizes coherent receivers with digital signal processing (DSP) techniques to equalize the signal suffering from a large accumulative CD [1,2]. Many DSP methods have been proposed that operate either in the frequency domain [3,4] or in the time domain [5,6,7,8]. Due to the quasi-static nature of the currently deployed optical network, CD can be considered as a static effect that requires only static equalization. Dynamic or adaptive CD equalization is also feasible but is usually enabled by an extra CD estimation stage before the CD equalization gets applied [9]. The frequency sampling method (FSM) [3] obtains a CD equalization response by sampling the ideal CD equalization response, namely H(ω)=e^jK^(ωT)2 by following the convention in [7], with N equally spaced points and then applies it to the N-point fast Fourier transform (FFT) of the input signal. The parameter K is defined as DLλ²/(4πcT²) where D, L, λ, T stands for the dispersion parameter, the fiber length, the central wavelength, and the sampling period, respectively. The FSM has advantages in both design and equalization complexity. However, there are no definite error criteria associated with the FSM except that the CD equalization response within the entire digital band ωT∈[-π, π] is exact at those N discrete frequencies. When performed in the time domain, the CD equalization involves the digital filter design techniques associated with distinct design criteria. The criteria include the impulse-invariant method [6], the passband least-squares (LS) approximation of H(ω) [7], the minimum mean-squared error (MMSE) between the transmitted and the received symbols [8], etc. Albeit with increasingly high computational complexity for large filter orders, the time-domain CD equalization filter delivers better performance than the conventional FSM, especially for the high-order modulation formats [3]. It seems that the passband least-squares approximation [7] achieves the best results among all existing methods. Still, we believe it is useful to build connections between them. In this paper, we show that the passband least-squares approximation can be performed discretely in an equivalent form. We also present a practical design procedure based on estimating the filter taps from several passband samples, which requires only FFTs without any evaluation of functions such as the complex erf. The proposed method has no performance loss against the closed-form solution and enables non-rectangular weighting in the passband least-squares approximation. A conventional impulse-invariant CD filter is revisited, improved and a connection is made between the improved impulse-invariant and the Fourier series interpretations.

2. Filter design principles

Although positioned very close to the head of DSP chain where the signal is still oversampled, the CD equalization problem is better to be solved in the context of a multirate system [4,8]. The main reason is that all CD equalization filters will inevitably cause amplitude fluctuations of the signal’s spectrum (i.e., passband ripples) which are turned into inter-symbol interference (ISI) after signal downsampling, as illustrated in Fig. 1. The best strategy is therefore to maximize the flatness of the filter frequency response within the signal’s passband (defined by the pulse-shaping function), or alternatively, to minimize the time-domain ISI caused by the CD equalization (cf. [8]). In contrast to a full-band all-pass design, it will require fewer taps to achieve the same level of equalization and add minimum burden to the following dynamic channel equalizers for the ISI removal. Trade-offs could be taken to balance the ISI and the noise suppression incurred by the CD equalization, but they need to be justified.

Fig. 1. Principle of passband CD equalization. The ISI-free pulse shaping determines the spectrum shape in the passband, which is altered by the passband ripples of the CD equalization process. The magnitudes of various filter responses are displaced on purpose for better illustration. Only the positive frequency range (0,π) is plotted here.

Download Full Size | PDF

2.1 Impulse-invariant method

By directly inverting the all-pass CD equalization response H(ω) for ωT∈[-∞, ∞] via Fourier transform, the analog impulse response can be obtained analytically [6], i.e.,

(1)$$h(t) = \frac{1}{{2\pi }}\int_{ - \infty }^\infty {{e^{jK{{({\omega T} )}^2}}}} {e^{j\omega t}}d\omega = \sqrt {\frac{j}{{4\pi K{T^2}}}} {e^{ - j\frac{\pi }{{4K{T^2}}}{t^2}}},\textrm{ }j = \sqrt { - 1} .$$

The impulse-invariant (II) method is to take symmetrically the first N samples of h(t) around t = 0 with sampling interval T=1/f_s and construct the FIR taps from the given samples h(nT). The method has certain optimal properties in the time domain approximation. But it has also severe frequency-domain aliasing problems if the analog filter has significant frequency response above one-half the sampling frequency, which is the case for an all-pass CD response without considering the bandlimited nature of the signal. A better design strategy is hence to include the bandlimiting effect and obtain the analog impulse response via the inverse Fourier transform. With the definition $erf(x) = \frac{2}{{\sqrt \pi }}\int_0^x {{e^{ - {t^2}}}} dt$, the impulse response is given by

(2)$$\begin{aligned} {h_{\textrm{II}}}(t) &= \frac{1}{{2\pi }}\int_{ - \infty }^\infty {W{e^{jK({\omega t} )}}^{^2}{e^{j\omega t}}d\omega } \\ &\textrm{ = }\frac{1}{{4\sqrt \pi }}{e^{ - j\frac{{{t^2}}}{{4K{T^2}}}}}\left[ {erf \left( {\sqrt { - jK} \left( {T{\omega_2} + \frac{t}{{2KT}}} \right)} \right) - erf \left( {\sqrt { - jK} \left( {T{\omega_1} + \frac{t}{{2KT}}} \right)} \right)} \right], \end{aligned}$$

where W is a rectangular lowpass function taking the unit for ω∈[ω₁, ω₂] and zero otherwise. The filter taps h[n] are then obtained by sampling the improved impulse response at instances t = nT. Namely, the FIR filter taps obtained via the impulse-invariant method applied to a bandlimited CD response is written as, for -(N-1)/2 ≤ n ≤ (N-1)/2 and N is odd,

(3)$${h_{{{\rm I}{\rm I}}}}[n] = \frac{1}{{4\sqrt \pi }}{e^{ - j\frac{{{n^2}}}{{4K}}}}\left[ {erf \left( {\sqrt { - jK} \left( {T{\omega_2} + \frac{n}{{2K}}} \right)} \right) - erf \left( {\sqrt { - jK} \left( {T{\omega_1} + \frac{n}{{2K}}} \right)} \right)} \right],$$

which is seen to be the N-point double-sided Fourier series of the periodical frequency response H_p(ω) = WH(ω) with ωT∈[-π, π]. This completes the first result of this paper and the solution Eq. (3) should be compared to the D vector in [7]. By the convergence theorems of Fourier series, the truncated Fourier series represents a least-squares approximation of the frequency response H_p(ω) within the entire digital band [-π, π], i.e.,

(4)$${{\mathbf h}_{{{\rm I}{\rm I}}}} = \arg \min \frac{1}{{2\pi }}\int_{ - \pi }^\pi {{{\left|{WH(\omega ) - \sum\limits_{n ={-} {{(N - 1)} / 2}}^{{{(N - 1)} / 2}} {h[n]{e^{ - jn\omega T}}} } \right|}^2}d(\omega T)}$$

with h_II = (h_II[-(N-1)/2],…,h_II[(N-1)/2])^T and T for transpose. As long as the weighting function W does not extend beyond the Nyquist band, i.e., max{|ω₁|, |ω₂|}<f_s/2, the improved impulse-invariant filter has no aliasing problem. However, as shown later, the passband ripple is another critical factor determining the filter performance, which is exploited by the passband least-squares approximation.

2.2 Least-squares approximation

Although the impulse-invariant filter h_II can already produce a frequency response that is optimal in the least-squares sense, it does not fully fulfill the goal of minimizing the passband ripples that determine the level of ISI. By reformulating the least-squares criterion and limiting the optimization range to be the passband only, a better filter is found [7], namely,

(5)$${{\mathbf h}_{\textrm{LS}}} = \arg \min \frac{1}{{2\pi }}\int_{{\Omega _1}}^{{\Omega _2}} {{{\left|{H(\omega ) - \sum\limits_{n ={-} {{(N - 1)} / 2}}^{{{(N - 1)} / 2}} {h[n]{e^{ - jn\omega T}}} } \right|}^2}d(\omega T)} ,$$

where h_LS = (h_LS[-(N-1)/2],…,h_LS[(N-1)/2])^T, Ω₁ and Ω₂ are chosen to be the spectrum edges of the signal. It is clear that the filter given by Eq. (4) is the limit case of Eq. (5) when the frequency range approaches the whole digital band. In contrast to the filter given in Eq. (4), the approximation error of filter Eq. (5) out of the range [Ω₁, Ω₂] is not controlled and could be huge. The fact that the CD is an all-pass effect with quadratic phase response enables analytical solutions of Eq. (5), i.e., with Q given in [7, Eq. (10)],

(6)$${{\mathbf h}_{\textrm{LS}}} = {{\bf Q}^{ - 1}}{{\mathbf h}_{{{\rm I}{\rm I}}}}$$

which however entails numerical evaluation of nontrivial functions such as the error function (erf) with complex-valued arguments. Also, this type of approximation implies a rectangular weighting of the signal’s passband spectrum. However, closed-form solutions to Eq. (5) may be inaccessible for an arbitrary non-rectangular weighting that might be desirable for other reasons (we give an example in a later section). The above expression also presents a link between the impulse responses generated by the II and LS method. To improve upon the direct solution of the least-squares problem given as Eq. (5), we propose to formulate the spectral approximation in a discrete manner, i.e., constructing a linear model

(7)$${\bf C}{{\mathbf h}_{\textrm{LS}}} = {{\mathbf H}_p}$$

where H_p is a px1 vector that contains the p-point (p > N) nonzero samples of H(ω) in the passband. Explicitly, with k, p₁, p₂, M being integers,

(8)$${H_p}[k] = {e^{jK{{(2\pi k/M)}^2}}},\textrm{ } - {p_1} \le k \le {p_2}$$

where p = p₁+p₂+1, M is the FFT size and is larger than p so that -p₁/M and p₂/M are very close to the passband edges Ω₁ and Ω₂, respectively. The matrix C is a shifted partition of the full discrete Fourier transform (DFT) operator D, which is of size M X M and has matrix elements

(9)$$D[k,n] = {e^{{{ - j2\pi nk} / M}}},\textrm{ 0} \le k,n \le M - 1$$

with both k and n being integers. Note that this is simply the Fourier transform of an identity matrix. Matrix C is of size p X N, accounting for evaluating the M-point DFT at p specific frequencies of an N X 1 vector h_LS, and has elements

(10)$$C[k,n] = {e^{{{ - j2\pi (k - {p_1})\left( {n - \frac{{N - 1}}{2}} \right)} / M}}},\textrm{ }0 \le k \le p - 1,\textrm{ }0 \le n \le N - 1.$$

The FFT size M is chosen large to better match the spectrum edges and its value does not affect the filter design complexity. Whereas a large p will result in a large matrix C for a given filter order N and increases the design complexity. The least-squares estimation of the N-point FIR taps h_LS is then expressed as

(11)$${{\mathbf h}_{\textrm{LS}}} = {({{\bf C}^{\dagger} }{\bf C} + \eta {\bf I})^{ - 1}}{{\bf C}^\dagger }{{\mathbf H}_p}$$

where the symbol ${\dagger} $ is conjugate transpose, ηI is a scaled identity to avoid bad matrix inversion. This is the main result of this paper. To draw parallels with the method in [7], the Q matrix is the ${{\mathbf C}^{\dagger} }{\mathbf C}$ here and the D vector is now the ${{\mathbf C}^{\dagger} }$H_p. The proposed method has several benefits: it avoids the numerical evaluation of the complex erf functions, uses only the DFT matrix that is easy to compute, and applies to the non-rectangular weighting of the signal’s spectrum in the passband (i.e., H_p could be arbitrary). The proposed method designs the CD equalization filter digitally and requires only a set of samples of the target spectrum. Although the solution (11) appears similar to an MMSE solution, little can be said about the underlying statistics of the least-squares approximation error.

3. Numerical simulation

We conduct numerical simulations to study the performance of the proposed filter design method. As illustrated in Fig. 2, a large block of 32GBaud 16-QAM signal samples with 2-fold oversampling and root-raised cosine (RRC) pulse shaping is generated. The CD effect is introduced to the signal by directly multiplying the FFT of a large signal block with the fine-grained CD response (which contains more than 10 million samples). The CD synthesis in this manner approximates the ideal CD response with great precision and the observable bit error ratio (BER) increases shown in the later simulation results are then due to the CD equalization errors. Alternative CD synthesis methods can also be considered provided they have enough precision in approximation of H(ω). Independent and identically distributed noise samples are added to the signal for simulating a finite signal-to-noise ratio (SNR). At the receiver, a matched filtering is performed, followed by the static CD equalization implemented as an FIR filter. The BER is calculated after the signal is downsampled to the symbol rate. A roll-off factor (ROF) of 0.22 for the RRC pulse shaping is assumed throughout the simulation. The K parameter used in the simulation is about 20, corresponding to a 500 km length of fiber with a 16 ps/nm/km dispersion parameter at a carrier wavelength of 1550 nm.

Fig. 2. Simulation setup for performance evaluation of the CD equalization filter approximating the ideal response H(ω). The x(t) denotes the 16-QAM symbols and v(t) is the additive noise. The upsampling and downsampling factor denoted by L is 2.

Download Full Size | PDF

A rectangular weighting with a width equal to the signal’s bandwidth is chosen to select the passband of H(ω), i.e., from -19.52 to 19.52 GHz. Several frequency samples (denoted by p) are used to estimate the h_LS according to Eq. (11). The BER performances of various designs are shown in Fig. 3. The impulse-invariant filter h_II is improved by bandlimiting the H(ω) [cf. Eq. (3)] such that the BER finally converges for a large filter order. The performance is improved further by the least-squares filter h_LS designed either by the direct evaluation given in [7] or by performing the discrete least-squares estimation proposed here. The performance difference between the two methods is negligible. Since the bandlimited impulse-invariant filter h_II is approximating the full digital band, the same number of filter taps results in a poor passband flatness that explains the performance degradation compared to h_LS clearly seen in the figure.

Fig. 3. BER of a coherently detected 16-QAM signal at 14 dB SNR versus the order of CD equalization FIR filters designed based on the estimated h_LS, the directly evaluated h_LS, the bandlimited h_II, and the full-band h_II, respectively. The estimated h_LS is obtained by using p = 611 samples of H(Ω=ωT) with passband (-0.61π, 0.61π). That is, the spectrum edges are Ω₁=-0.61π and Ω₂=0.61π. The FFT size used for solving h_LS is M = 1000 and the factor η is 1E-11.

Download Full Size | PDF

In Fig. 4, the number of CD response samples (p number) is varied to test the optimal design condition of h_LS. It can be concluded that the filter performance starts to degrade once the number of response samples drops below the threshold that is required for the BER to converge (about 130 in this case). Above that threshold, the performance is optimal and the discrete least-squares estimation of h_LS needs not be overdetermined thanks to the ηI term. For example, it can be seen that for p = 131 and all N > p, the estimation as in Eq. (11) is in fact under-determined. The simple form of the DFT operator along with the relatively small p number shows that the design complexity of the proposed method is low.

Fig. 4. BER of a coherently detected 16-QAM signal at 14 dB SNR versus the order of CD equalization FIR filter designed based on the estimated h_LS with different p values. The rest of the simulation parameters are the same as those given in Fig. 3 except that a smaller M is used for a smaller p with the same ratio p/M=611/1000 = 0.611. For p = [305,183,131,107,95,83], M is [500,300,216,176,156,136], respectively.

Download Full Size | PDF

From the filter frequency responses shown in Fig. 5, it is clear that the design technique ensures a good level of flatness in the passband. A longer filter and more frequency response samples will always provide better flatness. The Gibbs phenomenon is more prominent for a short filter since the filter is approximating a rectangular lowpass response in magnitude. It appears that the filter has a significantly raised out-of-band response and it becomes more severe when the filter is shorter. The out-of-band behavior of the filter may raise concerns because it will boost the residual noise, after matched filtering, in the stopband in real use cases. A better strategy would be combining the matched filtering and the time-domain CD equalization and design a filter that does the two tasks jointly.

Fig. 5. Frequency response of the FIR filter designed based on the estimated h_LS for CD equalization (modulation format-independent). The simulation parameters including the M values are the same as those given in Fig. 4.

Download Full Size | PDF

We therefore apply the proposed design method to sample the joint frequency response that includes both the RRC and CD equalization effects and estimate the h_LS accordingly. The frequency response of the matched RRC filter H_RRC(ω) and the ideal CD equalization response H(ω) are multiplied together as the new target frequency response, and the passband range remains the same as (-0.61π, 0.61π) since it is only determined by the RRC shape. Denote the joint spectrum as H_tot(ω)=H_RRC(ω)H(ω), also illustrated in Fig. 2. The H_p used in Eq. (11) is then obtained by sampling H_tot uniformly within the given passband. The estimation of h_LS then follows. Directly after the coherent detection, the raw samples are filtered by h_LS, producing samples that are matched filtered and CD equalized simultaneously. As shown in Fig. 6, the required number of taps for the BER to converge is inevitably larger compared to the case of doing RRC and CD filtering separately. But we suffer no BER performance loss and gain in the filter responses, as shown in Fig. 7, which ensures a 20 dB suppression in the stopband.

Fig. 6. BER of a coherently detected 16-QAM signal at 14 dB SNR versus the order of the FIR filter designed based on the estimated h_LS for joint matched filtering and CD equalization, and for CD equalization only, respectively. The simulation parameters are the same as those given in Fig. 3.

Download Full Size | PDF

Fig. 7. Frequency response of the FIR filter designed based on the estimated h_LS for joint matched filtering and CD equalization (modulation format-independent). The η factor is 1E6 and the rest of the simulation parameters are the same as those given in Fig. 3.

Download Full Size | PDF

4. Conclusion

We have shown a discrete design method for the CD equalization FIR filter based on the passband least-squares approximation. The filter taps can be obtained from the passband samples of the ideal CD response. The method has no performance loss compared to the direct calculation of the passband LS approximation and applies to the non-rectangular weighting of the frequency response. We have also shown that, in contrast to the conventional all-pass filter, an improved impulse-invariant filter can be obtained via taking the Fourier series of the bandlimited CD response. It has no spectrum aliasing and represents the limit case of the passband LS filter.

Funding

National Key Research and Development Program of China (2018YFB1801800); National Natural Science Foundation of China (U2001601); Science and Technology Planning Project of Guangdong Province (2020B0303040001).

Disclosures

The authors declare no conflicts of interest.

Data availability

All data included in this study are available upon request by contact with the corresponding author.

References

1. M. S. Faruk and S. J. Savory, “Digital signal processing for coherent transceivers employing multilevel formats,” J. Lightwave Technol. 35(5), 1125–1141 (2017). [CrossRef]

2. A. P. T. Lau, Y. Gao, Q. Sui, D. Wang, Q. Zhuge, M. H. Morsy-Osman, M. Chagnon, X. Xu, C. Lu, and D. V. Plant, “Advanced DSP techniques enabling high spectral efficiency and flexible transmissions: Toward elastic optical networks,” IEEE Signal Process. Mag. 31(2), 82–92 (2014). [CrossRef]

3. R. Kudo, T. Kobayashi, K. Ishihara, Y. Takatori, A. Sano, and Y. Miyamoto, “Coherent optical single carrier transmission using overlap frequency domain equalization for long-haul optical systems,” J. Lightwave Technol. 27(16), 3721–3728 (2009). [CrossRef]

4. F. Pittala, I. Slim, A. Mezghani, and J. A. Nossek, “Training-aided frequency-domain channel estimation and equalization for single-carrier coherent optical transmission systems,” J. Lightwave Technol. 32(24), 4849–4863 (2014). [CrossRef]

5. G. Goldfarb and G. Li, “Chromatic dispersion compensation using digital IIR filtering with coherent detection,” IEEE Photonics Technol. Lett. 19(13), 969–971 (2007). [CrossRef]

6. S. J. Savory, “Digital filters for coherent optical receivers,” Opt. Express 16(2), 804–817 (2008). [CrossRef]

7. A. Eghbali, H. Johansson, O. Gustafsson, and S. J. Savory, “Optimal least-squares FIR digital filters for compensation of chromatic dispersion in digital coherent optical receivers,” J. Lightwave Technol. 32(8), 1449–1456 (2014). [CrossRef]

8. Y. Liu, Y. Zhang, Y. Peng, and Z. Zhang, “Equalization of chromatic dispersion using Wiener filter for coherent optical receivers,” IEEE Photonics Technol. Lett. 28(10), 1092–1095 (2016). [CrossRef]

9. D. Wang, C. Lu, A. P. T. Lau, and S. He, “Adaptive chromatic dispersion compensation for coherent communication systems using delay-tap sampling technique,” IEEE Photonics Technol. Lett. 23(14), 1016–1018 (2011). [CrossRef]

Chromatic dispersion equalization FIR filter design based on discrete least-squares approximation

Abstract

1. Introduction

2. Filter design principles

2.1 Impulse-invariant method

2.2 Least-squares approximation

3. Numerical simulation

4. Conclusion

Funding

Disclosures

Data availability

References

Data availability

Cited By

Figures (7)

Equations (11)

Optics Express