Method and apparatuses for algorithm on QAM coherent optical detection

Abstract

Blind polarization demultiplexing algorithms based on complex independent component analysis (ICA) by negentropy maximization for quadrature amplitude modulation (QAM) coherent optical systems are disclosed. The polarization demultiplexing is achieved by maximizing the signal's non-Gaussianity measured by the information theoretic quantity of negentropy. An adaptive gradient optimization algorithm and a Quasi-Newton algorithm with accelerated convergence are employed to maximize the negentropy. Certain approximate nonlinear functions can be substitutes for the negentropy which is strictly derived from the probability density function (PDF) of the received noisy QAM signal with phase noise, and this reduces the computational complexity. The numerical simulation and experimental results of polarization division multiplexing (PDM)-quadrature phase shift keying (QPSK) and PDM-16QAM reveal that the ICA demultiplexing algorithms are feasible and effective in coherent systems and the simplified ones can also achieve equivalent performance.

Description

BACKGROUND OF THE INVENTION

The coherent optical communication systems which employ polarization division multiplexing (PDM) are promising and excellent solutions for high capacity and spectral-efficient communication. The PDM scheme can simply double the transmission rate by utilizing both polarization-orthogonal tributaries at the identical wavelength as multiplexing paths in the fiber. The availability of coherent optical receivers with PDM has been validated by plenty of high capacity experiments [1-5]. Despite of its attractive benefit, there are unavoidable problems in the implementation of the PDM method. The main problematic issue is the crosstalk of the two paths due to the random polarization state variation and polarization mode dispersion (PMD). In order to successfully recover the transmitted data, polarization demultiplexing in optical domain [6-8] or electrical domain [9-11] is needed to separate the mixed signals. In a digital coherent receiver, digital signal processing (DSP) techniques can be employed for polarization demultiplexing. The blind constant modulus algorithm (CMA) and its variants [12-17] are most commonly used, but they are not specially designed for the polarization demultiplexing purpose and might cause the singularity problem [18] of “converge to the same source”, implying that the demultiplexing techniques need to be improved. Supervised adaptive least mean square (LMS) algorithm [19] and single carrier frequency domain equalization (SC-FDE) [20] methods, which require training sequence, can also be implemented. However, they are primarily for channel equalization and behave analogously as CMA in polarization demultiplexing.

Independent component analysis (ICA) was originally proposed for the blind signal separation (BSS) problems, but now it has been developed into broad applications such as voice separation, image processing, bioinformatics, etc. Though widely used in signal processing, its applications in optical communication sphere are rare. Polarization demultiplexing based on ICA method has been explored through maximum likelihood estimation [21-22], in which a gradient optimization algorithm is used with the criterion strictly matching the probability density function (PDF). An approach based on signal higher order statistics makes the ICA demultiplexing algorithm independent of modulation format [23]. There is another ICA method exploiting the magnitude boundedness of digital signal for low symbol rate case [24] and high symbol rate case [25]. The ability of separating the mixed signal components further suggests the potential application of ICA in the newest spatial-division multiplexing technology [26].

BRIEF SUMMARY OF THE INVENTION

In this application, we derive a complex-value ICA algorithm by negentropy maximization originally proposed by Tülay Adali [27-28] for the purpose of polarization demultiplexing in actual fiber optic communication systems and verify it by numerical simulation and experiment. The result shows that the demultiplexing algorithm based on ICA is effective and its computation complexity can be decreased without performance deterioration.

The demultiplexing problem will be described and the fiber PDM channel will be analyzed. The proposed ICA algorithms and their simplifications will be stated in detail. The convergence, stability and the performance in the PMD emulator will be analyzed and the experimental results will be presented.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 is an overall diagram of a coherent PDM system.

FIG. 2(a) is a PDF of QPSK signal with phase noise and OSNR of 24 dB, FIG. 2(b) is a PDF of 16QAM signal with phase noise and OSNR of 28 dB, FIG. 2(c) is a cost function for QPSK, and FIG. 2(d) is a cost function for 16QAM.

FIG. 3 shows nonlinear functions (a) |a sin h(y)|², (b), |y|⁴ (c) |y|⁶ and their PDF (d) e^{−|a sin h(y)|2}, (e) e^−|y|4, (f) e^−|y|6.

FIG. 4 shows the learning curve of the gradient algorithm.

FIGS. 5(a)-(d) show the convergence of the quasi-Newton algorithm.

FIG. 6 shows the MSE of average convergent value versus polarization rotation frequency.

FIG. 7 is a schematic diagram of the simulation system.

FIG. 8 shows the BER of the mixed and unmixed PDM-QPSK or PDM-16QAM signal versus DGD value.

FIGS. 9(a)-(f) show constellations of (a) QPSK before demultiplexing; (b) QPSK after demultiplexing by ICA; (c) QPSK after follow-up DSP; (d) 16QAM before demultiplexing; (e) 16QAM after demultiplexing; and (f) 16QAM after follow-up DSP.

FIG. 10 shows the PDM-QPSK experimental system.

FIG. 11 is an Experiment performance comparison of CMA and ICA with a PMD emulator.

DETAILED DESCRIPTION OF THE INVENTION

Problem Description and PDM Channel Model

As shown in FIG. 1, in a PDM system, independent and identically distributed complex signal a=[a_X,a_Y]^T in QAM format is transmitted over both polarization tributaries of fiber. At the receiver, these signals are polarization-split and frequency down-converted by a polarization-diverse coherent receiver in order to extract the in-phase and quadrature components in two orthogonal polarizations of the received signal.

As noted in [21], a is affected by the time-variant phase noise and frequency offset of the lasers at optical transceivers, and also affected by Erbium doped fiber amplifier (EDFA) generated amplified spontaneous emission (ASE) noise in the fiber link, which can be viewed as white Gaussian noise. So the independent components which the ICA algorithms pursue can be expressed as s_K=(a_K+n_K)e^jϕK (K=X or Y), where ϕ_K is the interference due to the signal phase and is complex Gaussian noise. Additional carrier recovery process is needed to eliminate the ill effects of ϕ_K after applying the polarization demultiplexing method.

We assume that chromatic dispersion (CD) has been completely compensated in this paper, since stationary CD can be compensated by a fixed digital equalizer using frequency-domain or time-domain truncation method [15].

PMD is the most indispensable factor to be considered in designing optical PDM communication systems. In theory, PMD is mathematically modeled as a concatenation of birefringence fiber segments with arbitrary rotations around their principal axes and stochastic differential group delay (DGD) of Maxwellian distribution. The Jones transformation matrix can be written as:

$\begin{array}{cc}H\left(\omega \right)=\prod _{i=1}^n{P}_i\left(\omega \right)S_i& \left(1\right)\end{array}$

Here, P_i(ω) is the ith section's delay matrix, and S_i is the ith scattering matrix:

$\begin{array}{cc}{P}_i\left(\omega \right)=\left[\begin{array}{cc}{e}^{j\mathrm{\omega \Delta }{\tau }_i/2}& 0\\ 0& e^{-j\omega \Delta {\tau }_i/2}\end{array}\right]& \left(2\right)\\ P_i\left(\omega \right)=\left[\begin{array}{cc}{e}^{j\omega \Delta {\tau }_i/2}& 0\\ 0& e^{-j\omega \Delta {\tau }_i/2}\end{array}\right]& \left(3\right)\end{array}$ where ω is the angular frequency, Δτ_i is the DGD value of the ith section, θ_i and φ_i, uniformly distributed in [0,2π), respectively denote the frequency-independent random rotation and phase shift of principal axes [29], so that the corresponding PMD vectors would cover the whole Poincaré sphere. In all, H(ω) in (1) is a frequency-dependent unitary matrix:

$\begin{array}{cc}H\left(\omega \right)=\left[\begin{array}{cc}{H}_1\left(\omega \right)& H_2\left(\omega \right)\\ -H_2^{*}\left(\omega \right)& H_1^{*}\left(\omega \right)\end{array}\right],\left(\omega \in \left[0,2\pi \right)\right)& \left(4\right)\end{array}$ and |H₁(ω)|²+|H₂(ω)|²=1. Due to the frequency-dependent nature of PMD, PDM transmission channel model should be accurately described as a 2×2 multiple-input multiple-output (MIMO)-finite impulse response (FIR) structure. Thus the received signal x=[x_X,x_Y]^T is a mixed and distorted version of s=[s_X,s_Y]^T, which also means that the polarization demultiplexing algorithms should be, in principle, capable of dealing with convolutional mixing. However, considering the fact that the 1^st-order PMD value of fiber link is relatively low in many practical scenarios and adaptive equalizers are usually employed for PMD compensation [15] in the following DSP processing, when CD has been well compensated previously, the mixing channel can be simplified as an instantaneous matrix in a short period of time. As was done in [21-23], ignoring the frequency-selective nature of H(ω), we express the mixing matrix as

$\begin{array}{cc}H\left(\theta ,\phi \right)=\left[\begin{array}{cc}\cos \theta & \sin e^{j\phi }\\ -\sin \theta e^{-j\omega }& \cos \theta \end{array}\right]& \left(5\right)\end{array}$ where θ and φ are the parameters to be estimated by ICA algorithm discussed below.

Based on the FIG. 1, the mixing matrix links the independent components and mixed signal, namely x=Hs. The goal of ICA demultiplexing algorithm is to seek a matrix W which is the best estimation of H⁻¹, so that the independent components can be obtained from y=Wx. Because the polarization state of fiber is time-varying and usually unstable, the separation method has to be adaptive. The applicability of ICA demultiplexing algorithm depends on its performance, robustness and computational complexity.

Complex ICA Algorithm Detection

A. The Principle and Cost Function

According to the central limit theorem, a classical result in probability theory, the statistics of a mixed signal tends to be more Gaussian than its independent components under certain conditions. As noted by A. Hyvärinen [30], “Non-Gaussian is independent.” Thus, ICA is to maximize the non-Gaussianity of y=w^Hx, where w is one of the column vectors in W. Non-Gaussianity can be measured by negentropy which is based on the information theoretic quantity of differential entropy [30]. The negentropy of a complex value y can be defined as [27]

$\begin{array}{cc}\begin{array}{c}{J}_{\mathrm{neg}}=H\left(y_{\mathrm{Gauss}}\right)-H\left(y\right)\\ =\mathrm{const}-H\left(y\right)\end{array}& \left(6\right)\end{array}$ where y_Gauss is a complex Gauss variable of the same variance as y, H is the differential entropy defined as

$\begin{array}{cc}\begin{array}{c}H\left(y\right)=E\left\{\log \left[p\left(y\right)\right]\right\}\\ =-\int p\left(y\right)\log p\left(y\right)dy\\ =-\int \int p\left(y^R,y^I\right)\log p\left(y^R,y^I\right)d{y}^{R}d{y}^I\end{array}& \left(7\right)\end{array}$ and p(y)=p(y^R,y^I) is joint PDF of complex variance y. It can be proved in information theory that a Gauss variable has the largest entropy among all random variables of equal variance, so negentropy is always positive and is zero if and only if y is Gauss. Since H(y_Gauss) is constant, maximizing the non-Gaussianity of y is equivalent to minimizing the bivariate differential entropy in (7). Therefore, when using negentropy as a criterion of non-Gaussianity, the optimal cost function isJ(w)=E{−log [p_s(w^Hx)]}  (8)where p_s is the PDF of s_K(K=X or Y). In practice, expectation operator E in (8) is usually replaced by an arithmetic average or instantaneous value.

As mentioned in Section II, s=(a+n)e^{j ϕ}, n□N(0,2σ²), then the PDF of the unmixed independent component is

$\begin{array}{cc}\begin{array}{c}{P}_s\left(s\right)=\frac{1}{M}\sum _{i=1}^M{p}_{s|a_i}\left(s|a_i\right)\\ =\frac{1}{M}\sum _{i=1}^M\frac{1}{2{\mathrm{\pi \sigma }}^2}\exp \left(-\frac{{s}^2+{a_i}^2}{2{\sigma }^2}\right)I_0\left(\frac{a_i^{*}s}{{\sigma }^2}\right)\end{array}& \left(9\right)\end{array}$ where a_i is the complex point in the constellation diagram, and I_v is the v^th-order modified Bessel function of the first kind. The detailed derivation of p_s(s) is in [21] and the derivation of a from optical signal to noise ratio (OSNR) is presented in Appendix A. The PDF and cost function images of QPSK and 16QAM are shown in FIG. 2. The OSNRs of them are 24 dB and 28 dB, respectively. Because phase noise ϕ is uniformly distributed in [0,2π), the PDF and the cost functions only rely on the modulus of constellation points. In other words, they are circular symmetric, as can be seen from (9) and FIG. 2.

B. The Gradient Optimization Algorithm

To adaptively calculate the optimal vector w which minimizes the cost function J(w), a gradient optimization method can be employed. The update rule is

$\begin{array}{cc}w\leftarrow w+\mu \frac{\partial J\left(w\right)}{\partial w^{*}}w\leftarrow \frac{w}{\sqrt{w^{H}w}}& \left(10\right)\end{array}$ where μ is a negative learning rate, the 2^nd step in the update rule is to guarantee the normality of output.

$\frac{\partial J\left(w\right)}{\partial w^{*}}$ is the complex gradient which approaches zero near convergence, and given by

$\begin{array}{cc}\frac{\partial J\left(w\right)}{\partial w^{*}}=\frac{1}{2}E\left\{x{g}^{*}\left(w^{H}x\right)\right\}& \left(11\right)\end{array}$

The explanation and derivation of complex gradient are in Appendix B. The gradient algorithm can gradually adjust the vector parameters online using the newly updated data. Like other gradient algorithms, it is fit for varying and nonstationary environments, such as fiber PDM channel with PMD interference, but the convergence rate and stability depend on the initial w and learning rate μ.

C. Acceleration Convergence by Quasi-Newton Algorithm

The gradient algorithm has the drawback that it cannot converge rapidly and accurately. Adali proposed a Quasi-Newton algorithm to accelerate its convergence based on the Lagrangian cost function [27]L(w,λ)=J(w)+λ(w^Hw−1)  (12)

The second term of (12) is the constraint condition on w, λ is the Lagrange multiplier, and J(w) is defined in (8). Newton algorithm is a 2^nd-order updating rule which converges faster than gradient algorithm. Using complex gradient and Hessian in [31], Newton update can be defined as

$\begin{array}{cc}\begin{array}{c}\Delta =-{\left(\frac{{\partial }^{2}L}{\partial {*}^{}\partial T^{}}{❘}_{w=w_n}\right)}^{-1}\frac{\partial L}{\partial {*}^{}}{❘}_{w=w_n}\\ =-L_{-1}^{}{L}_{*}^{}\\ =-{\left(J_{}+\lambda \tilde{I}\right)}^{-1}\left(J_{*}^{}+\lambda \right)\end{array}& \left(13\right)\end{array}$ where =[w₁,w₁*, w₂,w₂*]^T. With the sophisticated derivation in Appendix C, we obtain the updating rule of Quasi-Newton algorithm:w←−½E{xg*(y)}+E{g_a(y)}w+E{xx^T}E{g_b(y)}w*  (14)

Though the updating rule of Quasi-Newton algorithm is much more complicated than the gradient algorithm, it is immune to the choice of learning rate and converges faster. A practical and feasible choice is that the processor employs the Quasi-Newton algorithm in batch processing mode at the beginning of computing. Once convergence has been achieved, it may turn into a gradient optimization mode to track the variation.

D. Simplification and Approximation

The cost function and update rules of the gradient and Quasi-Newton algorithms, as mentioned above, are rigidly and exactly derived from the PDF of the independent component. However, the complexity of update rules in (10) and (14), whose complete general formulae are in Appendix B & C, almost leads to their infeasibility in real-time receivers. It is necessary to simplify the expressions.

Some approximate nonlinear functions can be substitutes for negentropy [30], because they implicitly introduce some higher-order statistics which can be viewed as measures of non-Gaussianity. A. Hyvärinen have also proved that “any sufficiently smooth even function” can be used as a cost function for ICA by either maximizing or minimizing its value [32]. It is proposed in [33] that the cost function for complex ICA can be defined asJ(w)=E{|G(w^Hx)|²}  (15)where G is a nonlinear function which has to be chosen to match the PDF of independent component. As observed from (8), the PDF matching the cost function in (15) isp_G(y)=e^−|G(y)|2  (16)

Some nonlinear functions, such as G₁(y)=a sin h(y), G₂(y)=y², G₃(y)=y³, along with their associated PDFs are shown in FIG. 3. They are circular symmetric and similar to the actual PDFs and cost functions in FIG. 2, but the cost function |y|⁶ is excessively sensitive to outliers due to its rapid rate of growth.

The gradient of the cost function in (15) is

$\begin{array}{cc}\frac{\partial J\left(w\right)}{\partial w^{*}}=\frac{1}{2}E\left\{{\mathrm{xG}}^{*}\left(y\right)g\left(y\right)\right\}& \left(17\right)\end{array}$ and the update rule of the Quasi-Newton algorithm [33] is

$\begin{array}{cc}w\leftarrow -E\left\{G^{*}\left(y\right)g\left(y\right)x\right\}+E\left\{g\left(y\right)g^{*}\left(y\right)\right\}w+E\left\{{\mathrm{xx}}^T\right\}E\left\{G^{*}\left(y\right)g^{\prime }\left(y\right)\right\}{w}^{*}& \left(18\right)\end{array}$ where g is the derivative of G and g′ is the derivative of g. The low-complexity approximate update rules are suitable for hardware implementation.

E. Calculation of Independent Components on Both Polarizations

The previously discussed algorithms are classified as “one-unit” algorithms which estimate just one of the independent components. To obtain several independent components, a conventional practice is to run one-unit algorithms several times for different weights respectively and then decorrelate them. As for the issue of polarization demultiplexing, this process can be reduced and simplified because the optimized weights w₁ and w₂ must be orthogonal due to the Jones matrix of fiber transmission in (5), namely, w₁^Hw₂=0. Therefore, once we have estimated one of them, say w₁, then w₂ can also be calculated without applying the one-unit algorithm again. To calculate a complex vector that is orthogonal to w₁, a simple way is to use the Gram-Schmidt orthogonalization algorithm:

$\begin{array}{cc}{w}_2=w_0-w_0{w}_0^H{w}_1{w}_2\leftarrow \frac{w_2}{\sqrt{w_2^H{w}_2}}& \left(19\right)\end{array}$ where w₀ is a vector that is orthogonal to the initial w₁. The orthogonalization algorithm extracts a vector of the same direction as w₁ and leaves the vector w₂ that is orthogonal to it. Additionally, it is worthwhile to note that ICA cannot make a distinction between the separated demultiplexing signals, so some special frame header information are needed to identify them.

F. Preprocessing

Before applying an ICA algorithm, it is usually profitable to preprocess the data. The preprocessing is mainly for centering and whitening x, to make x a zero-mean and uncorrelated variable, namely, E{x}=0 and E{xx^H}=I. One popular method for whitening is to use eigenvalue decomposition (EVD) of the covariance matrix E{xx^H}, which leads to very high computation-complexity. The whitening can be achieved adaptively by a gradient algorithm, in which the update rule isW←W+μ[I−Wxx^HW^H]  (20)Where W is the whitening matrix, Wx is the whited signal, and μ is the convergent rate. A rough interpretation of the rule is that [I−Wxx^HW^H] becomes zero when whitening is achieved at convergence.Numerical Simulation

A. Convergence of the Algorithms

In order to investigate the convergence of the proposed algorithms, we firstly assume that PDM channel is static, and the parameters in (5) are θ=40°, φ=60°. The OSNRs are set to 18 dB for PDM-QPSK and 24 dB for PDM-16QAM. The convergence depends on the choice of mixing matrix, but the case shown is typical. A uniformly distributed phase noise has also been attached to the symbols. The learning curves of the gradient algorithms and Quasi-Newton algorithms are shown in FIG. 4-5, respectively, which are the arithmetic means of 1000 times simulation. The nonlinear function used here and in the following part of the article is G₂(y)=y².

The curves in FIG. 4 represent the change of cost functions in (8) derived from the PDF of QPSK and 16QAM signal, which reach steady state after 200˜300 iterations. There is a tradeoff between the accuracy and convergent speed, so μ has to be chosen carefully. The curves of the same kind of signal converge in a similar fashion. They almost converge to the minima of their cost functions.

The Quasi-Newton algorithms are in batch processing mode which employs 2000 symbols of both polarizations. The convergence is much faster than the above gradient algorithms as can be seen in the horizontal ordinate in FIG. 5, but the computation complexity is even higher. So we could use the Quasi-Newton algorithms to get close to convergence, and then switch to the gradient algorithm to track it.

B. Dynamic Tracking of the Adaptive Algorithms

In order to gain further insight into the dynamic tracking behavior of the adaptive gradient algorithms, we multiply the channel Jones matrix used in Section IV.A by an endless polarization rotation matrix which is formatted as

$\begin{array}{cc}A=\left(\begin{array}{cc}\cos \omega t& \sin \omega t\\ -\sin \omega t& \cos \omega t\end{array}\right)& \left(21\right)\end{array}$ where ω is the polarization rotation angular frequency. In FIG. 6, MSE refers to the mean square error of the average convergent values relative to the minima of their corresponding cost functions. The average convergent value is also calculated via an arithmetic mean of 1000 times simulation, and the symbol rates of the QPSK and 16QAM signal are supposed to be 28 Gsym/s. Seen from FIG. 6, the four methods can track a polarization rotation frequency less than 10 Mrad/s without significant convergent stability deterioration.

C. Performance in PMD Emulator

The algorithms are tested with a PMD emulator to evaluate the demultiplexing ability in fiber link. In the simulation in FIG. 7, the optical carrier from a laser of 100 kHz line-width is split into two orthogonal polarizations by a polarization beam splitter (PBS). The two I/Q modulators are driven by 28 GSym/s two levels or four levels signal to generate PDM-QPSK or PDM-16QAM signal, the bit rates of which are therefore 107 Gbit/s and 428 Gbit/s, respectively. After that, the two polarizations tributaries are coupled again by a polarization beam coupler (PBC). An instantaneous mixing Jones matrix as in Section IV.A with a concatenated PMD emulator is employed as transmission channel. As mentioned in Section II, CD is neglected, and the PMD emulator is controlled by the mean differential group delays (DGDs) of fiber segments to simulate the effects of the 1st and higher order PMD of real fibers. The OSNRs are the same as before. At the receiver, a polarization-diverse 90° hybrid, along with a local laser of the same 100 kHz line-width are used for coherent detection.

For every DGD value, BERs (bit error ratios) are measured over 32000 symbols after ICA and follow-up DSP are applied, which includes blind adaptive equalization and phase recovery. The ICA algorithms are implemented in this way: the first 2000 symbols of each polarization are employed for the Quasi-Newton algorithms in a batch processing mode to make the demultiplexing matrix w approach to convergence, and then it switches to the gradient algorithm to track the polarization change adaptively for the rest of the symbols. The results are shown in FIG. 8 with BERs of unmixed signal. No matter which kind of ICA algorithm is applied, the one accurately derived from the signal PDF or the approximate one using nonlinear functions, the results are similar. The coincidence of the curves proves the successful polarization demultiplexing of ICA algorithms and also proves that the simplified ICA algorithms are as effective as the accurate ones even though the computation is greatly reduced.

FIG. 9 shows the constellations for the signal before/after applying ICA methods and after finishing all the DSP processing for the case of zero-DGD. The ICA algorithms transfer the signals to their intended moduli.

Experimental Results

The proposed ICA algorithms are also tested in the experimental system shown in FIG. 10. In the experiment, 40 Gbit/s PDM-QPSK signal is generated by one IQ modulator, a PBS, a PBC and an optical delay line. A 1^st-order PMD emulator and a polarization scrambler are inserted before 100 km standard single mode fiber (SSMF) to ensure that the resulting PMD has all orders. The eye diagram of the generated optical QPSK signal is in FIG. 10(a), and the OSNR of the received optical signal is set to 20 dB by adjusting the EDFA and the optical attenuator.

The DSP algorithms for experimental data are identical to those used for the simulation system except that fiber CD compensation is employed. The BERs shown in FIG. 11 are calculated over 40000 bits. The ICA and CMA have similar performance, and this result is coincident with that of [22]. The result reveals that if singularity does not occur, the ICA and CMA are comparable in polarization demultiplexing for PDM-QPSK. Even so, ICA excels CMA, because singularity cannot exist with ICA [21] and ICA is adaptive to other modulation formats.

CONCLUSION

We have derived the polarization demultiplexing algorithms and their simplifications based on ICA by negentropy maximization. It is found that they are effective for coherent detected PDM-QAM signals. It is further shown by experiment that the performance of ICA and CMA for demultiplexing PDM-QPSK is comparable, but ICA has its own advantages of immunity to singularity and modulation format independence. The ICA algorithms can also be potentially applied to eliminate the crosstalk and interferences between sub-channels in the newest spatial-division multiplexing systems which employ multi-core or few-mode fibers.

APPENDIX

A. Derivation of σ from OSNR

The power of noise n can be expressed as

$\begin{array}{cc}\begin{array}{c}E\left({n}^2\right)=\int \int {n}^2{p}_n\left(n\right)dxdy\\ =\int \int \frac{{n_x}^2+{n_y}^2}{2{\mathrm{\pi \sigma }}^2}\exp \left(-\frac{{n_x}^2+{n_y}^2}{2{\sigma }^2}\right)dxdy\\ =\int \frac{{n_x}^2}{\sqrt{2\pi }\sigma }\exp \left(-\frac{{n_x}^2}{2{\sigma }^2}\right)dx\int \frac{{n_y}^2}{\sqrt{2\pi }\sigma }\exp \left(-\frac{{n_y}^2}{2{\sigma }^2}\right)dy\\ =2{\sigma }^2\end{array}& \left(22\right)\end{array}$

Assuming the power of received signal is normalized, so received OSNR is

$\begin{array}{cc}\begin{array}{c}\mathrm{OSNR}=\frac{E\left({a}^2\right)}{E\left({n}^2\right)}\\ =\frac{1-E\left({n}^2\right)}{E\left({n}^2\right)}\\ =\frac{1-2{\sigma }^2}{2{\sigma }^2}\end{array}& \left(23\right)\end{array}$

B. Derivation of Complex Gradient the Cost Function

The complex gradient is derived in detail in [31, 34]. Here we use the results directly. Complex gradient is defined as

$\begin{array}{cc}\frac{\partial f}{\partial z}=\frac{1}{2}\left(\frac{\partial f}{\partial z^R}-j\frac{\partial f}{\partial z^I}\right);\frac{\partial f}{\partial z^{*}}=\frac{1}{2}\left(\frac{\partial f}{\partial z^R}+j\frac{\partial f}{\partial z^I}\right)& \left(24\right)\end{array}$ where z∈□, z* is the conjugate of z and f is analytic with respect to z and z* independently. So using the chain rule, we have

$\begin{array}{cc}\begin{array}{c}\frac{\partial J\left(w\right)}{\partial w_i^{*}}=\frac{1}{2}E\left\{\frac{\partial \left\{-\log \left[p_s\left(w^{H}x\right)\right]\right\}}{\partial w_i^R}+j\frac{\partial \left\{-\log \left[p_s\left(w^{H}x\right)\right]\right\}}{\partial w_i^I}\right\}\\ =\frac{1}{2}E\left\{\begin{array}{c}{g}^R\left(y\right)\frac{\partial {\left(w^{H}x\right)}^R}{\partial w_i^R}+g^I\left(y\right)\frac{\partial {\left(w^{H}x\right)}^I}{\partial w_i^R}+\\ j\left[g^R\left(y\right)\frac{\partial {\left(w^{H}x\right)}^R}{\partial w_i^I}+g^I\left(y\right)\frac{\partial {\left(w^{H}x\right)}^I}{\partial w_i^I}\right]\end{array}\right\}\\ =\frac{1}{2}E\left\{g^R\left(y\right)x_i^R+g^I\left(y\right)x_i^I+j\left[g^R\left(y\right)x_i^I-g^I\left(y\right)x_i^R\right]\right\}\\ =\frac{1}{2}E\left\{x_i{g}^{*}\left(y\right)\right\},\left(x_i=x_i^R+{\mathrm{jx}}_i^I,g\left(y\right)=g^R\left(y\right)+{\mathrm{jg}}^I\left(y\right)\right)\end{array}& \left(25\right)\end{array}$

Thus

$\begin{array}{cc}\frac{\partial J\left(w\right)}{\partial w^{*}}=\frac{1}{2}E\left\{{\mathrm{xg}}^{*}\left(w^{H}x\right)\right\}& \left(26\right)\end{array}$ where

$\begin{array}{cc}\begin{array}{c}{g}^R\left(y\right)=-\frac{\partial \log p_s\left(y\right)}{\partial y^R}\\ =\frac{-y^R}{M2{\mathrm{\pi \sigma }}^4{p}_s\left(y\right)}\sum _{i=1}^M\left\{\exp \left(-\frac{{y}^2+{a_i}^2}{2{\sigma }^2}\right)\left[\begin{array}{c}{I}_1\left(\frac{a_i^{*}y}{{\sigma }^2}\right)\frac{a_i^{*}}{y}-\\ I_0\left(\frac{a_i^{*}y}{{\sigma }^2}\right)\end{array}\right]\right\}\end{array}& \left(27\right)\end{array}$

Similarly,

$\begin{array}{cc}{g}^I\left(y\right)=\frac{-y^I}{M2{\mathrm{\pi \sigma }}^4{p}_s\left(y\right)}\sum _{i=1}^M\left\{\exp \left(-\frac{{y}^2+{a_i}^2}{2{\sigma }^2}\right)\left[I_1\left(\frac{a_i^{*}y}{{\sigma }^2}\right)\frac{a_i^{*}}{y}-I_0\left(\frac{a_i^{*}y}{{\sigma }^2}\right)\right]\right\}& \left(28\right)\end{array}$

C. Derivation of Updating Rule of Quasi-Newton Algorithm

We start the derivation from (13) and assuming Δ=_n+1−_n. We can have(_J+λĨ)_n+1=−_J*+_Jn  (29)*_J is the complex gradient [31,34], whose odd elements are defined in (11). The complex Hessian _J is given by [27],

$\begin{array}{cc}\begin{array}{c}{J}_{}=\frac{{\partial }^{2}J\left(w\right)}{\partial {*}^{}\partial T^{}}\\ =\left[\begin{array}{cccc}\frac{{\partial }^{2}J\left(w\right)}{\partial w_1^{*}\partial w_1}& \frac{{\partial }^{2}J\left(w\right)}{\partial w_1^{*}\partial w_1^{*}}& \frac{{\partial }^{2}J\left(w\right)}{\partial w_1^{*}\partial w_2}& \frac{{\partial }^{2}J\left(w\right)}{\partial w_1^{*}\partial w_2^{*}}\\ \frac{{\partial }^{2}J\left(w\right)}{\partial w_1\partial w_1}& \frac{{\partial }^{2}J\left(w\right)}{\partial w_1\partial w_1^{*}}& \frac{{\partial }^{2}J\left(w\right)}{\partial w_1\partial w_2}& \frac{{\partial }^{2}J\left(w\right)}{\partial w_1\partial w_2^{*}}\\ \frac{{\partial }^{2}J\left(w\right)}{\partial w_2^{*}\partial w_1}& \frac{{\partial }^{2}J\left(w\right)}{\partial w_2^{*}\partial w_1^{*}}& \frac{{\partial }^{2}J\left(w\right)}{\partial w_2^{*}\partial w_2}& \frac{{\partial }^{2}J\left(w\right)}{\partial w_2^{*}\partial w_2^{*}}\\ \frac{{\partial }^{2}J\left(w\right)}{\partial w_2\partial w_1}& \frac{{\partial }^{2}J\left(w\right)}{\partial w_2\partial w_1^{*}}& \frac{{\partial }^{2}J\left(w\right)}{\partial w_2\partial w}& \frac{{\partial }^{2}J\left(w\right)}{\partial w_2\partial w_2^{*}}\end{array}\right]\\ =E\left\{g_a\left(y\right)\left[\begin{array}{cccc}{x}_1{x}_1^{*}& 0& x_1{x}_2^{*}& 0\\ 0& x_1^{*}{x}_1& 0& x_1^{*}{x}_2\\ x_2{x}_1^{*}& 0& x_2{x}_2^{*}& 0\\ 0& x_2^{*}{x}_1& 0& x_2^{*}{x}_2\end{array}\right]\right\}+\\ E\left\{\left[\begin{array}{cccc}0& x_1^2{g}_b\left(y\right)& 0& x_1{x}_2{g}_b\left(y\right)\\ x_1^{*2}{g}_b^{*}\left(y\right)& 0& x_1^{*}{x}_2^{*}{g}_b^{*}\left(y\right)& 0\\ 0& x_2{x}_1{g}_b\left(y\right)& 0& x_2^2{g}_b\left(y\right)\\ x_2^{*}{x}_1^{*}{g}_b^{*}\left(y\right)& 0& x_2^{*2}{g}_b^{*}\left(y\right)& 0\end{array}\right]\right\}\end{array}& \left(30\right)\end{array}$ where

$\begin{array}{cc}\begin{array}{c}{g}_a\left(y\right)=4\left[\frac{\partial g^R\left(y\right)}{\partial y^R}+\frac{\partial g^I\left(y\right)}{\partial y^I}+j\left(\frac{\partial g^R\left(y\right)}{\partial y^I}-\frac{\partial g^I\left(y\right)}{\partial y^R}\right)\right]\\ =\frac{8g\left(y\right)}{y}+\frac{2h\left(y\right){y}^2}{{\mathrm{\pi \sigma }}^2}\end{array}\begin{array}{c}{g}_b\left(y\right)=4\left[\frac{\partial g^R\left(y\right)}{\partial y^R}-\frac{\partial g^I\left(y\right)}{\partial y^I}+j\left(\frac{\partial g^R\left(y\right)}{\partial y^I}+\frac{\partial g^I\left(y\right)}{\partial y^R}\right)\right]\\ =\frac{2h\left(y\right)\left[{\left(y^R\right)}^2-{\left(y^I\right)}^2+j2y^R{y}^I\right]}{{\mathrm{\pi \sigma }}^4}\end{array}& \left(31\right)\\ h\left(y\right)=\frac{1}{{\mathrm{Mp}}_s\left(y\right)}\sum _{i=1}^M\left\{\exp \left(-\frac{{y}^2+{a_i}^2}{2{\sigma }^2}\right)\left[\left({\frac{a_i}{\sigma y}}^2+\frac{1}{{\sigma }^2}-\frac{g\left(y\right)}{y}\right)I_0\left(\frac{a_{i}y}{{\sigma }^2}\right)-\frac{2a_i}{y}\left(\frac{1}{{y}^2}+\frac{1}{{\sigma }^2}-\frac{g\left(y\right)}{2y}\right)I_1\left(\frac{a_{i}y}{{\sigma }^2}\right)\right]\right\}& \left(32\right)\end{array}$

Assuming x has been whitened, and after removing the even rows of the matrix _J and _J*, it results in

$\begin{array}{cc}\left(H_J+\lambda I\right)w_{n+1}=-\frac{1}{2}E\left\{{\mathrm{xg}}^{*}\left(y\right)\right\}+E\left\{{\mathrm{xx}}^H{g}_a\left(y\right)\right\}{w}_n+E\left\{{\mathrm{xx}}^T{g}_b\left(y\right)\right\}{w}_n^{*}\approx -\frac{1}{2}E\left\{{\mathrm{xg}}^{*}\left(y\right)\right\}+E\left\{g_a\left(y\right)\right\}{w}_n+E\left\{{\mathrm{xx}}^T\right\}E\left\{g_b\left(y\right)\right\}{w}_n^{*}& \left(33\right)\end{array}$

At the convergent point, (_J+λĨ) will become real, as proved in [27]. Therefore, the fix-point update is

$\begin{array}{cc}w\leftarrow -\frac{1}{2}E\left\{{\mathrm{xg}}^{*}\left(y\right)\right\}+E\left\{g_a\left(y\right)\right\}w+E\left\{{\mathrm{xx}}^T\right\}E\left\{g_b\left(y\right)\right\}{w}^{*}& \left(34\right)\end{array}$

REFERENCES

• [1] X. Zhou, et al., “High Spectral Efficiency 400 Gb/s Transmission Using PDM Time-Domain Hybrid 32-64 QAM and Training-Assisted Carrier Recovery,” IEEE/OSA Journal of Lightwave Technology, vol. 31, pp. 999-1005, 2013.
• [2] J. Yu, et al., “Transmission of 200 G PDM-CSRZ-QPSK and PDM-16 QAM With a SE of 4 b/s/Hz,” IEEE/OSA Journal of Lightwave Technology, vol. 31, pp. 515-522, 2013.
• [3] P. J. Winzer, “High-Spectral-Efficiency Optical Modulation Formats,” IEEE/OSA Journal of Lightwave Technology, vol. 30, pp. 3824-3835, 2012.
• [4] M. Yiran, Y. Qi, T. Yan, C. Simin, and W. Shieh, “1-Tb/s Single-Channel Coherent Optical OFDM Transmission With Orthogonal-Band Multiplexing and Subwavelength Bandwidth Access,” IEEE/OSA Journal of Lightwave Technology, vol. 28, pp. 308-315, 2010.
• [5] Y. Ma, Q. Yang, Y. Tang, S. Chen, and W. Shieh, “1-Tb/s single-channel coherent optical OFDM transmission over 600-km SSMF fiber with subwavelength bandwidth access,” OSA Optics Express, vol. 17, pp. 9421-9427, 2009.
• [6] X. Yuan, J. Zhang, Q. Jing, Y. Zhang, and M. Zhang, “A novel scheme for automatic polarization division demultiplexing,” in Proc. Asia Communications and Photonics Conference and Exhibition (ACP), 2011, pp. 83092P.
• [7] A.-L. Yi, L.-S. Yan, J. Ye, W. Pan, B. Luo, and X. S. Yao, “A Novel Scheme for All-Optical Automatic Polarization Division Demultiplexing,” in Proc. Asia Communications and Photonics Conference and Exhibition (ACP), 2009, pp. TuAA1.
• [8] X. S. Yao, L. S. Yan, B. Zhang, A. E. Willner, and J. Jiang, “All-optic scheme for automatic polarization division demultiplexing,” OSA Optics Express, vol. 15, pp. 7407-7414,2007.
• [9] Z. Yu, X. Yi, Q. Yang, M. Luo, J. Zhang, L. Chen, et al, “Polarization demultiplexing in stokes space for coherent optical PDM-OFDM,” OSA Optics Express, vol. 21, pp. 3885-3890, 2013.
• [10] C. Gong, X. Wang, N. Cvijetic, and G. Yue, “A Novel Blind Polarization Demultiplexing Algorithm Based on Correlation Analysis,” IEEE/OSA Journal of Lightwave Technology, vol. 29, pp. 1258-1264,2011.
• [11] B. Szafraniec, B. Nebendahl, and T. Marshall, “Polarization demultiplexing in Stokes space,” OSA Optics Express, vol. 18, pp. 17928-17939, 2010.
• [12] L. Ling, T. Zhenning, Y. Weizhen, S. Oda, T. Hoshida, and J. C. Rasmussen, “Initial tap setup of constant modulus algorithm for polarization de-multiplexing in optical coherent receivers,” in Proc. Optical Fiber Communication Conference and Exposition and the National Fiber Optic Engineers Conference (OFC/NFOEC), 2009, pp. 1-3.
• [13] A. Vgenis, C. S. Petrou, C. B. Papadias, I. Roudas, and L. Raptis, “Nonsingular Constant Modulus Equalizer for PDM-QPSK Coherent Optical Receivers,” IEEE Photonics Technology Letters, vol. 22, pp. 45-47, 2010.
• [14] K. Kikuchi, “Polarization-demultiplexing algorithm in the digital coherent receiver,” in Digest of the IEEE/LEOS Summer Topical Meetings, 2008, pp. 101-102.
• [15] X. Zhou and J. Yu, “Digital signal processing for coherent optical communication,” in Proc. 18th AnnualWireless and Optical Communications Conference (WOCC), 2009, pp. 1-5.
• [16] O. Zia-Chahabi, R. Le Bidan, C. Laot, and M. Morvan, “A self-reconfiguring constant modulus algorithm for proper polarization demultiplexing in coherent optical receivers,” in Proc. Optical Fiber Communication Conference and Exposition (OFC/NFOEC), 2012 and the National Fiber Optic Engineers Conference, 2012, pp. 1-3.
• [17] I. Roudas, A. Vgenis, C. S. Petrou, D. Toumpakaris, J. Hurley, M. Sauer, et al, “Optimal Polarization Demultiplexing for Coherent Optical Communications Systems,” IEEE/OSA Journal of Lightwave Technology, vol. 28, pp. 1121-1134,2010.
• [18] S. J. Savory, “Digital Coherent Optical Receivers: Algorithms and Subsystems,” IEEE Journal of Selected Topics in Quantum Electronics, vol. 16, pp. 1-16, 2010.
• [19] Y. Han and G. Li, “Coherent optical communication using polarization multiple-input-multiple-output,” OSA Optics Express, vol. 13, pp. 7527-7534, 2005.
• [20] C. C. Do, C. Zhu, A. V. Tran, S. Chen, T. Anderson, D. Hewitt, et al., “Chromatic dispersion estimation in 40 Gb/s coherent polarization-multiplexed single carrier system using complementary Golay sequences,” in Proc. Optical Fiber Communication Conference and Exposition and the National Fiber Optic Engineers Conference(OFC/NFOEC), 2012, pp. 1-3.
• [21] Pontus. Johannisson, Henk. Wymeersch, Martin. Sjodin, A. Serdar Tan, et al, “Convergence Comparison of the CMA and ICA for Blind Polarization Demultiplexing,” IEEE/OSA Journal of Optical Communications and Networking, vol. 3, pp. 493-501, 2011.
• [22] H. Zhang, Z. Tao, L. Liu, S. Oda, T. Hoshida, and J. C. Rasmussen, “Polarization demultiplexing based on independent component analysis in optical coherent receivers,” in Proc. 34th European Conference on Optical Communication (ECOC) 2008, pp. 1-2.
• [23] X. Xiaobo, F. Yaman, Z. Xiang, and L. Guifang, “Polarization Demultiplexing by Independent Component Analysis,” IEEE Photonics Technology Letters, vol. 22, pp. 805-807, 2010.
• [24] T. Oktem, A. T. Erdogan, and A. Demir, “Adaptive Receiver Structures for Fiber Communication Systems Employing Polarization-Division Multiplexing,” IEEE/OSA Journal of Lightwave Technology, vol. 27, pp. 5394-5404, 2009.
• [25] Turgut M. Öktem, Alper T. Erdogan, and Alper Demir, “Adaptive Receiver Structures for Fiber Communication Systems Employing Polarization Division Multiplexing: High Symbol Rate Case,” IEEE/OSA Journal of Lightwave Technology, vol. 28, pp. 1536-1546, 2010.
• [26] I. B. Djordjevic, “Spatial-Domain-Based Hybrid Multidimensional Coded-Modulation Schemes Enabling Multi-Tb/s Optical Transport,” IEEE/OSA Journal of Lightwave Technology, vol. 30, pp. 2315-2328, 2012.
• [27] M. Novey and T. Adali, “Complex ICA by Negentropy Maximization,” IEEE Transactions on Neural Networks, vol. 19, pp. 596-609, 2008.
• [28] M. Novey and T. Adali, “Complex Fixed-Point ICA Algorithm for Separation of QAM Sources using Gaussian Mixture Model,” in Proc. IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), 2007, pp. 11-445-11-448.
• [29] P. K. A. Wai, C. R. Menyuk, and H. H. Chen, “Stability of solitons in randomly varying birefringent fibers,” Optics. Letters, vol. 16, pp. 1231-1233, 1991.
• [30] A. Hyvarinen and E. Oja, “Independent component analysis: algorithms and applications,” Neural Networks, vol. 13, pp. 411-430, 2000.
• [31] A. van den Bos, “Complex gradient and Hessian,” IEE Proceedings—Vision, Image and Signal Processing, vol. 141, pp. 380-383, 1994.
• [32] E. Bingham and A. Hyvarinen, “A fast fixed-point algorithm for independent component analysis of complex valued signals,” International Journal of Neural Systems, vol. 10, pp. 1-8, 2000.
• [33] T. Adali, L. Hualiang, M. Novey, and J. F. Cardoso, “Complex ICA Using Nonlinear Functions,” IEEE Transactions on Signal Processing, vol. 56, pp. 4536-4544, 2008.
• [34] D. H. Brandwood, “A complex gradient operator and its application in adaptive array theory,” IEE Proceedings H-Microwaves, Optics and Antennas, vol. 130, pp. 11-16, 1983.

Claims

1. A method for polarization demultiplexing in an optical system for maximizing a signal's non-Gaussianity using independent component analysis, wherein the non-Gaussianity is measured by negentropy based on information quantity of differential entropy, the method comprising: receiving an optical beam containing the signal;generating a local laser beam;coherently detecting the signal using the local laser beam;employing an optimal cost function in a processor;employing a quasi-Newton update algorithm in the processor until the algorithm converges; andthereafter employing a gradient optimization update algorithm to track signal variation, to minimize the optimal cost function in the processor,wherein the negentropy is characterized by the probability density function (PDF) of the signal, and the minimizing the cost function includes substituting the PDF with a nonlinear function.

2. The method of claim 1, wherein the optimal cost function is J(w)=E{−log[ps(wHx)]} or J(w)=E{|G(wHx)|2}.

3. The method of claim 2, wherein a probability density function matching the cost function is pG(y)=e−|G(y)|2.

4. The method of claim 1, wherein the gradient optimization update algorithm is

5. The method of claim 1, wherein the gradient optimization update algorithm gradually adjusts vector parameters using newly updated data.

6. The method of claim 1, wherein the quasi-Newton update algorithm is

7. The method of claim 1, wherein the Quasi-Newton update algorithm is employed by the processor in a batch processing mode.

8. The method of claim 1, further comprising running a one-unit algorithm a plurality of times for different weights respectively; and decorrelating the results of one-unit algorithm, wherein a plurality of independent components are obtained.

9. The method of claim 1, further comprising employing unique frame headers to identify separated signals.

10. The method of claim 1, further comprising preprocessing data for whitening.

11. The method of claim 10, wherein the preprocessing includes using eigenvalue decomposition.

12. The method of claim 10, wherein the preprocessing is achieved adaptively by a gradient algorithm.

13. An optical system for maximizing a signal's non-Gaussianity using independent component analysis, wherein the non-Gaussianity is measured by negentropy based on information quantity of differential entropy, the optical system comprising a transmitter and a receiver, wherein the transmitter includes a laser and transmits a signal over a fiber optic medium using the laser,wherein the receiver carries out steps comprising:receiving an optical beam containing the signal;generating a local laser beam;coherently detecting the signal using the local laser beam;employing an optimal cost function;employing on the signal a quasi-Newton update algorithm until the algorithm converges;employing after the quasi-Newton update algorithm a gradient optimization update algorithm to minimize the optimal cost function; andrunning a one-unit algorithm a plurality of times for different weights respectively and decorrelating the results, wherein a plurality of independent components are obtained.

14. The optical system of claim 13, wherein the receiver carries out further steps comprising employing unique frame headers to identify separated signals.

15. The optical system of claim 13, wherein the receiver carries out further steps comprising preprocessing data for whitening.