power consumption and small footprint of short-reach optical networks, the intensity modulation and direct detection (IM/DD) scheme is an attractive solution [2]. To keep pace with the required capacity, multilevel pulse amplitude modulation (PAM) has been utilized, and a steady increase of symbol rate has been observed, from 10 Gbaud (e.g., IEEE 802.3ae [3]) to 25 Gbaud (e.g., IEEE 802.3bs [4]) and 50 Gbaud (e.g., IEEE 802.3cu [5]), with higher than 100 Gbaud being studied for beyond 400 Gb/s operation [6].

One of the challenges with high symbol rate is the limited bandwidth of the deployed devices, such as digital-to-analog and analog-to-digital converters (DACs/ADCs), modulators and photodetectors (PDs). It induces inter-symbol interference (ISI) and degrades receiver sensitivity. Channel equalization via digital signal processing (DSP) can effectively mitigate the penalty of ISI. Feedforward equalizers (FFE) are commonly utilized linear equalizers that correct both postcursor and precursor ISI. They can be used alone [7], [8], or along with decision feedback equalizers (DFEs) and/or maximum likelihood sequence estimation (MLSE) [9]- [14]. The tap coefficients can be adjusted with least-mean squares (LMS) algorithm, where stochastic gradient descent (SGD) is conducted to minimize the mean-squared error (MSE) between training symbols and equalized symbols. LMS algorithm has attracted a lot of attention, mainly due to its computational simplicity. However, its convergence rate depends on the autocorrelation of the equalizer input signal, which is strongly influenced by the channel spectral characteristics. The convergence rate is slow when the channel has severe distortion [15], [16], which is often the case with bandlimited IM/DD systems. A large number of training symbols are thus required [9], [11], resulting in high transmission overhead, especially for burst data transmission.

Recursive least squares (RLS) algorithm is an alternative algorithm to adjust tap coefficients. The cost function is the weighted sum of squared error calculated using all the past received symbols. For the tap coefficients update, a matrix is used in place of the scalar step size in the conventional LMS algorithm, which greatly accelerates convergence [15], [16]. The convergence rate can be improved by one order of magnitude [13]. Nevertheless, the calculation of this matrix also gives rise to higher computational complexity. Another solution is to apply variable step size in tap coefficients adjustment. Rather than being constant, the step size is initialized to a larger value to speed up initial convergence, and then reduced to smaller ones as the equalizer output error decreases. Several algorithms have been proposed, such as a two-step procedure in [17], and varying step size at each iteration with values determined by the number of sign changes of an estimated gradient of MSE [18] or the squared instantaneous error [19]. The learning rate decay techniques developed for the SGD including exponential decay and staircase decay can also be applied to adjust step size, achieving similar performance with the aforementioned algorithms [8].

In this paper, the LMS-based equalization is accelerated by exploiting a novel training sequence. Unlike the traditional training sequence that is made up of independent and identically distributed (i.i.d.) samples, correlated samples generated from a first-order Markov chain (MC) are utilized. It alters the autocorrelation of the equalizer input signal and enables rapid convergence. The performance of the proposed technique is compared with the conventional i.i.d. sequence, where they are both employed for training. Faster convergence is observed in terms of tap coefficient adaptation and MSE. And lower pre-forward-error-correction (pre-FEC) bit error rate (BER) can be attained at a given training sequence length. An experimental demonstration of a 43 Gbaud four-level PAM (PAM-4) signal is conducted, where more than 70% training sequence length reduction can be achieved for 5-km standard single mode fiber (SSMF) transmission, and more than 48% reduction can be achieved for 10-km SSMF transmission.

For an N -tap equalizer adapted with LMS algorithm, the equalizer input signal y(k) is often arranged into vectors of length N :

Here N is assumed to be an odd number. Let the equalizer tap coefficients at instant k be w(k), and then the equalizer output signal is

where the superscript T denotes transpose. The tap coefficients are updated according to

where µ is the step size and e(k) is the equalizer error calculated using e (k)= x(k) -z(k), with x(k) being the training signal.

The optimal tap coefficients that achieve the minimum MSE (MMSE) is

where A = E{y(k)y(k) T } is an N × N autocorrelation matrix, and b = E{y(k)x(k)} is an N × 1 cross-correlation vector.

Since A is a symmetric matrix, it can be expressed as A = U ΛU T , where U is an orthonormal matrix whose column vectors are the eigenvectors of A, and Λ is a diagonal matrix with elements equal to the eigenvalues of A. Define a transformed tap coefficient error vector v(k) as

Then it can be proven that [15]- [17]:

As iteration increases, v(k) decreases exponentially, and decay rate of the i-th element is determined by the i-th eigenvalue λ i and the step size. To have E{v(k)} converge, the following condition has to be met

Consequently,

with λ max being the maximum eigenvalue. When there is a large difference between λ max and λ min (the minimum eigenvalue), having µ satisfy (8) slows down the convergence of the elements with small eigenvalues. Therefore, eigenvalue spread, defined as χ = λ max /λ min , determines the convergence rate of LMS algorithm. χ is generally difficult to calculate, but it is bounded by the dynamic range of the power spectrum

The equality holds as N →∞.

Both RLS algorithm and LMS adaptation with variable step size focus on amending (3), either by replacing µ with an estimate of A -1 or by having µ take different values at different stages. The proposed method, on the other hand, aims at changing the eigenvalue spread of A. As (9) shows, smaller χ can be obtained if the dynamic range of P (ω) reduces. Furthermore,

, with H(ω) and X(ω) being the channel response and the spectrum of the training signal. While H(ω) is fixed, it is still possible to modify P (ω) by changing X(ω). Since most of the bandlimited IM/DD systems behave as low-pass filters, the training sequence should have more power distributed at higher frequencies to "flatten" P (ω). A first-order MC is employed to generate such a sequence and the transition probability is

where a m and a n are symbols from the signal alphabet, β is a hyperparameter to tune the spectrum of x(k), and Z m is a normalization factor to ensure that the probabilities sum to 1. When β =0 , it becomes the conventional i.i.d. process. When β>0, symbols with larger Euclidean distances from the current symbol are transmitted at the next instant with higher probabilities. Fig. 1 shows the transition matrices for PAM-4 sequences with β =0 .03 and β =0 .07. In both cases, all the symbols have the lowest probability to remain unchanged at the next instant. Thus, the diagonal element achieves the minimum in each row. For the symbols with negative amplitudes, the most probable symbol at the next instant is 3, while it becomes -3for the symbols with positive amplitudes. Although the exponential term e β(a n -a m ) 2 is the same for the transitions from a m to a n and from a n to a m , Z m are not necessarily equal as the symbol at the current instant varies, and thus the transition matrix is not symmetric. As β increases, the symbol "furthest" from the current symbol gains higher probability to appear in the MC sequence at the following instant. Conversely, the current symbol itself has lower probability to be retransmitted. Fig. 2(a)-(d) plot the power spectral densities (PSDs) of the MC sequences when β are 0, 0.01, 0.03 and 0.07, respectively. The correlation between symbols enhances the power at high frequencies, and such spectral shaping becomes stronger with larger β.

To illustrate the effect of the proposed MC sequence in the LMS-based equalization, simulations are conducted where H(ω) is assumed to be fourth-order Bessel filters or Gaussian filters with varying bandwidths. Additive white Gaussian noise (AWGN) is added after the channel. A 15-tap T-spaced FFE is employed to compensate the ISI. Table I lists χ of i.i.d. sequence and MC sequence under various H(ω). 3-dB bandwidth of H(ω) is normalized by the Nyquist bandwidth BW Nyq (half the symbol rate). The element at the i-th row and j-th column of A is estimated by the Monte Carlo estimation:

where L is the length of the training sequence and is fixed at 10 5 in the simulations. Each kind of sequence is generated independently for 1000 times and χ reported in Table I is the mean value of these 1000 sequences. For both kinds of filters, χ of the i.i.d. sequence increases with the decrease of the 3-dB bandwidth as a result of stronger bandwidth limitation. For the MC sequence, χ depends on both β and the 3-dB bandwidth. For the MC sequence with β =0 .01, the relationship between χ and the 3-dB bandwidth is similar to that of the i.i.d. sequence. On the other hand, for the MC sequence with β =0 .07, χ increases as the 3-dB bandwidth becomes larger. In all the cases considered in Table I, there exists at least one MC sequence whose χ is smaller than that of the i.i.d. sequence and the reduction becomes greater when the bandwidth limitation becomes more severe.

All the tap coefficients of the 15-tap FFE is initialized to 0, except for the center tap, which is initialized to 1. Fig. 3(a) and (b) respectively plot the value of the center tap and the MSE versus the number of iterations when H(ω) is a fourth-order Bessel filter with 3-dB bandwidth of 0.4 BW Nyq . The training sequence is an i.i.d. sequence or an MC sequence with β =0 .07. It can be clearly seen that the smaller χ of the MC sequence results in faster convergence. It takes the center tap fewer iterations to converge and the MSE decreases at a faster speed. Furthermore, as shown in Fig. 4, there is negligible difference between the converged tap coefficients of these two kinds of training sequences. Thus, the tap coefficients trained by the MC sequence can be applied to channel equalization when i.i.d. data are transmitted.

It is worth mentioning that by exploiting the MC, the probability distribution of the symbols is also shaped. In the case of PAM-4 sequence, -3 and 3 have higher probabilities than -1 and 1. However, the probabilistic shaping itself does not contribute to a faster convergence since the spectrum remains unchanged. For the example considered in Table I, if the MC sequence is randomly shuffled before transmission, which produces a sequence with a white spectrum and the same probability distribution with the original MC sequence, similar χ to that of the i.i.d. sequence is observed. Consequently, the shuffled MC sequence has a similar convergence rate with the i.i.d. sequence and cannot accelerate the equalizer training as the unshuffled MC sequence does. III. EXPERIMENTAL SETUP Fig. 5(a) shows the experimental setup. A 43 Gbaud PAM-4 signal is generated offline and pulse shaped by a root-raisedcosine (RRC) filter with a roll-off factor of 0.1. The digital signal is converted to an analog waveform by an arbitrary waveform generator (AWG, Keysight M8195A) at 64 GSa/s with an analog bandwidth of 25 GHz. The generated signal is amplified by a 35-GHz driver amplifier and modulated onto a 1550-nm optical carrier via a 40-GHz Mach-Zehnder modulator (MZM) biased at the quadrature point. The modulated optical signal is transmitted over a 5-km or 10-km SSMF. A variable optical attenuator (VOA) adjusts the received optical power (ROP) before a 11-GHz photodetector (PD) with integrated amplifiers. The detected signal is captured for offline DSP by a digital storage scope (DSO, Keysight DSOZ254A) with an analog bandwidth of 25 GHz at 80 GSa/s. The frequency response of the transmission system with a 5-km SSMF is shown in Fig. 5(b). The 6-dB bandwidth is about 10 GHz.

The transmitted sequence consists of an i.i.d. training sequence or an MC training sequence followed by a data sequence with 1.5×10 5 i.i.d. PAM-4 samples with uniform distribution. As mentioned in Section II.C, the probability distribution of the PAM-4 symbols is not uniform in the MC sequence, resulting in higher average power. For fair comparison, the MC sequence is normalized to have the same average power with the i.i.d. sequence before transmission. In the receiver-side DSP, the received signal is first resampled to 2 samples per symbol, and then equalized by a T/2-spaced FFE for 5-km SSMF transmission or a T/2-spaced FFE + T-spaced DFE for 10-km SSMF transmission. Throughout the experiments, the center-tap initialization that has been used in the simulation is adopted for the FFE. For the DFE, the first tap is initialized to -1 and the other taps are initialized to 0. BER calculation is conducted after the channel equalization using only the data sequence.

Hyperparameter β is first fixed at 0.07 and the FFE has 61 taps. Fig. 6 shows the PSDs of the received training sequences. After transmitting over the link, the spectrum of the MC sequence has a smaller dynamic range due to the spectral shaping at the transmitter. Fig. 7 plots the MSE versus the number of iterations. Again, the MC training sequence has a faster speed of convergence, achieving lower MSE after the same number of iterations. Such improvement in MSE can be translated into improvement of BER. Fig. 8 8(g) and (h), respectively. When the length of the training sequence is increased to 4000, the residual ISI for both i.i.d. sequence and MC sequence becomes smaller and the noise dominates the performance. Still, the pre-FEC BER of the MC sequence is lower at all the measured ROPs, with slightly higher performance gain at higher ROP.

Table II summarize the reductions in training sequence length of the MC sequence compared with the i.i.d. sequence at various FEC thresholds. At a fixed ROP, the improvement of adopting the MC training sequence generally becomes larger as the considered FEC threshold becomes lower. For instance, at an ROP of -2.7 dBm, the length of the training sequence can be reduced by 72.63% at the FEC threshold of 1.25×10 -2 , while it increases

to 93.20% at the FEC threshold of 3.15×10 -3 . For all the cases listed in Table II, more than 70% reduction can be achieved with the MC training sequence.

To further verify that the change in the probability distribution does not lead to improvement, the shuffled MC sequence described in Section II.C is also utilized for training. The pre-FEC BER versus the length of the training sequence is plotted in Fig. 8(a) and (b), where it is labelled as "SHUFF". The shuffled MC sequence achieves similar BER with the i.i.d. PAM-4 sequence under a given training sequence length at an ROP of -2.7 dBm, and slightly worse BER at a lower ROP of -6.7 dBm.

The performance of the LMS-based equalization also depends on the number of FFE taps. If the training sequence is adequately long, the BER can be improved by increasing the number of FFE taps. Fig. 9 shows the pre-FEC BER versus the number of FFE taps in such a scenario, where there are 2×10 4 samples in the training sequence. Fig. 9(a) and (b) are the results for i.i.d. sequence and MC sequence, respectively. The BERs of a 61-tap FFE that have been shown before are illustrated with empty circles. The relationships between BER and the number of FFE taps are similar at ROPs of -2.7 dBm, -4.7 dBm and -6.7 dBm, where the BER decreases as the number of FFE taps becomes larger. However, this relationship would change if the length of the training sequence is limited. Fig. 10(a) and (b) depict the pre-FEC BER versus the number of FFE taps at an ROP of -2.7 dBm, when the training sequences have 150 and 800 symbols, respectively. The BER first decreases and then increases as the number of FFE taps increases from 21 to 85. This is because although increasing the number of FFE taps can improve the capability to mitigate channel distortion, it would slow the rate of convergence at the same time. No matter which FFE is selected, the MC sequence offers a better choice for training since it achieves lower BER in all the cases.

Another factor that affects the performance of the MC sequence is the hyperparameter β. The pre-FEC BERs at different β are shown in Fig. 11, in which the ROP is -2.7 dBm and the number of FFE taps is 61. Training sequences with 150 and 400 symbols are considered. In both cases, lower BER compared with the i.i.d. sequence can already be observed at β =0 .01. And it gets further reduced as β becomes larger until 0.07. After 0.07, the BER gets worse with higher β.Thisis mainly because the eigenvalue spread, on the other hand, will be increased if the spectral shaping is too strong. It can also be seen from Table I, where the MC sequences with β of 0.03 and 0.07 have higher eigenvalue spreads than that of the i.i.d. sequence when the 3-dB bandwidth of the channel is 0.6 BW Nyq .T h e performance discrepancy between MC sequences with different β gets smaller as the length of the training sequence increases. In real bandlimited IM/DD systems, β can be selected according to the electrical bandwidth of the deployed components. A larger β should be chosen as the bandwidth limitation becomes stronger.

Fig. 12 depicts the PSDs of the received signals after 10-km SSMF transmission. The spectral null at around 19.6 GHz is due to the power fading induced by chromatic dispersion. β is 0.07 in the MC sequence and its spectrum again has a smaller dynamic range. The channel equalization is fulfilled with a 77-tap FFE + 12-tap DFE. To mitigate the error propagation of the DFE, joint symbol decision proposed in [14] is adopted. When the Euclidean distance between the output of the equalizer and the output of a decision device is larger than 0.5, the symbol-wise decision is regarded as unreliable and 6 following symbols are utilized for joint decision. To calculate the equalizer output for the following symbols, all possible combinations of PAM-4 symbols are used in the DFE part. The combination that yields the smallest sum of squared distance between the equalizer outputs and the assumed symbol decisions is selected. Fig. 13 The pre-FEC BER versus the length of the training sequence at ROPs of -1.9 dBm and -4.9 dBm are plotted in Fig. 14 III. At the FEC threshold of 1.25×10 -2 , the reductions are generally larger at higher ROPs. 54.02% fewer training symbols are required when the ROP is -4.9 dBm, while it becomes 72.19% at an ROP of -1.9 dBm. Fig. 15 plots the pre-FEC BER versus β when the ROP is -1.9 dBm and the training sequence has 1000 symbols. The The sequence is generated with a first-order MC that introduces correlation between samples and enhances the power at high frequencies. Compared with the traditional i.i.d. sequence with a white spectrum, the MC sequence can reach the convergence with fewer iterations and achieve lower pre-FEC BER with a fixed training sequence length. Experimental results show that to transmit a 43 Gbaud PAM-4 signal over a system of 5-km SSMF and 6-dB bandwidth about 10 GHz, the length of the training sequence can be reduced by more than 70%. When a 10-km SSMF is utilized and the signal suffers from severe power fading induced by chromatic dispersion, the training sequence length can be reduced by more than 48%.