Machine learning for optical modulation format recognition

Machine Learning for Optical Modulation Format Recognition Ricard Boada Farràs Bachelor Thesis S136096 May 2014 DTU Fotonik, Technical University of Denmark, Kgs. Lyngby, Denmark Supervised by: Idelfonso Tafur Monroy Robert Borkowski 2 3 Machine Learning for Optical Modulation Format Recognition Ricard Boada Farràs Bachelor Thesis S136096 May 2014 DTU Fotonik, Technical University of Denmark, Kgs. Lyngby, Denmark Supervised by: Idelfonso Tafur Monroy Robert Borkowski 4 Name of the thesis Machine Learning for Optical Modulation Format Recognition Author(s): Ricard Boada Farràs Supervisor(s): Idelfonso Tafur Monroy (DTU Fotonik) Robert Borkowski (DTU Fotonik) This report is a part of the requirements to achieve the Bachelor of Science (BSc) in Telecommunications at Technical University of Denmark. The report represents 15 ECTS points. Department of Photonics Engineering Technical University of Denmark Oersted Plads, Building 343 DK-2800 Kgs. Lyngby Denmark www.fotonik.dtu.dk Tel: (+45) 45 25 63 52 Fax: (+45) 45 93 65 81 E-mail: [email protected] 5 Abstract Next generation of optical networks will necessitate the use of advanced functionalities in the transceivers. In particular, modulation format recognition will enable dynamically switched energy efficient networks that can adapt to traffic demand and optimize resource utilization. This thesis presents the first comparison of clustering algorithms applied to modulation format recognition in Stokes space. 2-, 4-, 8-PSK and 8-, 12-, 16-QAM modulation formats are recognized and the requirements in terms of OSNR performance, accuracy and complexity are analyzed. A novel technique is also proposed for modulation format recognition in Stokes space, based on maximum likelihood between the received signal and the characteristics of the modulations targeted, showing OSNR requirements lower than those in the previously published literature. 6 7 Acknowledgements First of all, I would like to express my gratitude to Professor Idelfonso Tafur Monroy for giving me the opportunity to join the Metro-Access & Short Range Communications Group at DTU Fotonik and helping me with my stay in Denmark. Thanks to him and also to my co-supervisor, PhD Robert Borkowski, for all their support during the thesis. I feel highly indebt for all the feedback and advice they gave me on my daily work. I would also like to acknowledge Joan M. Gené Bernaus, from Polytechnic University of Catalonia, for encouraging me to do my thesis in DTU and helping me before and during the whole project. I wish to thank my department colleagues and especially my office workmates, Los Boludos del 358, for all the good moments since I arrived in Denmark. Further thanks to my classmates, Els Competitius, from Barcelona, who have always been there during these four years. Finally, none of this would have been possible without the patience and belief of my parents, my sister, my family and my friends. Thank you so much. 8 9 Table of Contents Abstract ............................................................................................................................. 5 Acknowledgements .......................................................................................................... 7 Table of Contents ............................................................................................................. 9 Acronyms ....................................................................................................................... 11 List of Figures ................................................................................................................. 12 List of Tables .................................................................................................................. 14 1 Introduction .................................................................................................. 15 1.1 Motivation ............................................................................................. 15 1.2 Problem statement ................................................................................. 17 1.3 State of the Art ...................................................................................... 18 1.4 Methodology ......................................................................................... 20 1.5 Contributions ........................................................................................ 21 1.6 Thesis outline ........................................................................................ 21 2 Polarization demultiplexing ......................................................................... 23 2.1 Polarization demultiplexing in Stokes space ........................................ 24 2.2 Representation of Stokes parameters in the Poincaré sphere ................ 26 3 Machine Learning ........................................................................................ 29 3.1 Clustering algorithms ............................................................................ 30 3.1.1 K-Means ............................................................................................ 31 3.1.2 Expectation Maximization (EM) ....................................................... 33 3.1.3 Density Based Spatial Clustering of Application with Noise (DBSCAN) ………………………………………………………………………36 3.1.4 Ordering Points to Identify the Clustering Structure (OPTICS) ....... 39 3.1.5 Spectral Clustering ............................................................................ 41 3.1.6 Gravitational Clustering .................................................................... 43 16 The ability of deciding settings in each node reduces the importance of the control plane, used to communicate the connection parameters between different devices. This is normally done in a parallel link to the data one, normally slower, that makes the whole communication more difficult to operate since the information carried in this control plane must be known by the receiver in advance to properly demodulate the data. Cognitive optical networks (CON) aim to provide a solution for the control of these heterogeneous optical networks, by simultaneously supporting different switching paradigms and protocols [4]. A cognitive network has three main building blocks [5]:  Network monitoring elements (NME), which inform about the current status of the network;  Software adaptable elements (SAE), which adapt their configuration according to current conditions of the network;  Cognitive processes, which learn and make decisions according to current and past observations on the system. In this thesis we will focus on the two main SAEs: a reconfigurable transmitter, capable of signal generation using arbitrary modulation formats, and an autonomous adaptable receiver, able to ensure that signals transmitted are optimally acquired and demodulated without prior knowledge on the signal parameters and modulation format of the incoming transfer of data. There are several different scenarios in which cognitive networks enhance the performance of the system. For example, they make possible to mix protocols and physical layer interfaces in the same network, they provide the capability to use a large selection of different QoS requirements by identifying bottlenecks or changing priorities in the nodes, or also might be used for security purposes by analyzing patterns of operation that could risk the safety of the system. Introduction of cognition into the optical layer will enable energy efficiency. By autonomously adapting and identifying the transmission parameters it will be possible to optimize the resource utilization of the system, leading to a more power efficient solution than that offered by traditional networks. Three main frameworks have analyzed and pushed forward the idea of cognition in application to fiber-optic communication systems: COGNITION [6] by Zervas et al., CON [7] by Wei et al. and CHRON [8]. 17 Cognitive optical networks (CON) Performance monitoring Cognitive processes Software adaptability Figure 1.2. Interaction between different elements in a cognitive network. 1.2 Problem statement With the advent of highly reconfigurable networks capable of adapting the modulation format of the conveyed signal to the conditions of the link, seems very important to have a mechanism to autonomously detect the characteristics of the signal at the receiver side without any extra information, i.e. reducing the importance of a control plane to send control commands. In particular, it is very important to detect correctly the modulation format used in the communication in order to be able to demodulate it without any prior information. It also helps for Modulation Format Opaque algorithms, like equalization, which work more effectively when the modulation format is known in advance. This functionality is known as Modulation Format Recognition (MFR) [9]. MFR has been extensively explored for wireless applications in the field called ‘cognitive radio’ [10] with applications ranging from access networks to military networks [11], where the goal is to optimize the spectrum utilization by using the most optimal modulation format given the channel conditions. However it has not been much researched in optical networks until very recently due to conventionally static nature of these networks. One can think of modulation format recognition in several contexts. An immediate application is Optical Burst Switching (OBS) for heterogeneous networks where the connections are turned up and down very fast. The need for rapid connection setups motivates the employment of MFR as the use of the control plane and human 18 operation would noticeably slow down the operation of the whole system. It also helps to decrease the overhead produced when many connections are initialized. Modulation format recognition can be performed in several different ways on optical communications as stated in the state of the art, from applying a simple iterative Machine Learning algorithm to the constellation, to using statistical analysis or employing power distributions. Figure 1.3. Example of two modulation formats and their corresponding Poincaré sphere. 1.3 State of the Art Researchers are increasingly studying new methods to be applied on modulation format recognition for optical communications. So far, the state of the art includes six different techniques: a. Modulation format identification based on constellation diagram measurements [12]. This approach is tested with a burst phase modulated Radio-over-Fiber communication and consists on using the information provided by the constellation to recognize the format used. The two criteria employed are the number of levels in the constellation histogram to classify between PSK and QAM modulations and the number of clusters formed, identified using K-Means. This information is compared with a knowledge database containing characteristic data of the recognizable modulations. One of the main drawbacks of this method is its need for entire receiver-side processing (polarization demultiplexing, equalization and carrier phase recovery) before MFR, as illustrated in Figure 1.3 (a). b. Artificial Neural Networks (ANN) applied on MFR [13]. This technique uses a feature based classification technique named Artificial Neural Networks with features extracted from an asynchronous amplitude 19 histogram (AAH) to identify the modulation format. It works with direct detection and it suffers from the need for prior training. c. Modulation format recognition with high-order statistics (HOS) in Stokes space [14]. This MFR is performed in Stokes space. Firstly, a main classification between 16QAM and the other possible formats (OOK, BPSK, and QPSK) is done using a variational learning algorithm. The discrimination between these three latter modulations consists on evaluating signal cumulants on the 2D projection of the Stokes space constellation and, by considering a Gaussian Mixture Model (GMM), clustering this projection with a hierarchical clustering algorithm. Figure 1.3 (b) shows the position within the DSP flow where the MFR is applied. d. Stokes space-based optical MFR with a machine learning technique [15]. This approach considers the Stokes space representation of the signal, where each modulation format shows a different number of clusters. A clustering method is used to determine this number of clusters and hence the modulation format of the incoming signal. The algorithm chosen is a variation of Expectation Maximization, a commonly used Machine Learning technique. This method allows for MFR at an early stage in the receiver, as displayed in Figure 1.3 (c); does not require training and provides valuable information for modulation format dependent methods, like equalization. e. Power distribution based MFR for digital coherent receivers [16]. Modulation formats are recognized by identifying special features from the received signal power distribution. Although this technique produces very good results comparing to the FEC limit, the features extracted from the power distribution depend on the OSNR, which requires monitoring, and thus, the knowledge of the modulation in advance. The position in the DSP flow of a coherent receiver is shown in Figure 1.4 (d). f. Modulation format identification using physical layer characteristics [17]. The technique tested distinguishes between 15 possible singlepolarized modulations by extracting information from the complex electric field magnitude histogram. Differentiation between different formats is done in separate stages by comparing measurements taken directly from the magnitude of the electric field. This procedure is tested with WDM signals with 40000 samples. Figure 1.4 (e) shows the coherent receiver DSP flow of this approach. 20 In this thesis we aim to continue the work started by Borkowski et al. [15]. Different clustering algorithms will be tested in Stokes space to assess their performance in recognizing six different polarization division multiplexed (PDM) modulation types. Moreover, a new algorithm for modulation format recognition will be presented and compared in equal terms with the other methods analyzed. I/Q imbalance compensation Low pass filter CD compensation Timing recovery Polarization demultiplexing Equalization Carrier phase recovery Sampled signal Constellation MFR I/Q imbalance compensation Low pass filter CD compensation Timing recovery Polarization demultiplexing HOS MFR Sampled signal I/Q imbalance compensation Low pass filter CD compensation Timing recovery Sampled signal Stokes space MFR I/Q imbalance compensation Low pass filter CD compensation Timing recovery Polarization demultiplexing Sampled signal Power distribution MFR I/Q imbalance compensation Low pass filter CD compensation Timing recovery Polarization demultiplexing Sampled signal Physical Layer MFR (a) (b) (c) (d) (e) Figure 1.4. Coherent receivers’ DSP flows of state-of-the-art techniques on optical modulation format recognition. 1.4 Methodology The methodology used in this project consists firstly of a deep background research on the topic of modulation format recognition in the optical domain, by reviewing all the state-of-the-art methods used in the literature (in journals, papers, reference books, …). A study about the main features of cognitive networks and polarization demultiplexing techniques has also been done. The most part of the time has been devoted to learn about different clustering techniques and to the analysis of their suitability for modulation format recognition 21 purposes. The implementation of the algorithms has been done using the software MATLAB. Finally, the different techniques have been compared entirely by simulation in three categories analyzing their performance, reliability and complexity. 1.5 Contributions After carefully reviewing the state-of-the-art, with this thesis we are presenting the first comparison of clustering algorithms in application to modulation format recognition in fiber-optic networks so far in the literature. Six different techniques (KMeans, Expectation Maximization, DBSCAN, OPTICS, Spectral Clustering and Gravitational Clustering) have been contrasted in terms of their OSNR performance, accuracy and complexity to assess their feasibility in a real working environment. Furthermore, we are introducing a novel algorithm for identifying the modulation format, solely from the input data set, without assuming any particular model for the data. Its performance has been compared with the other algorithms in the same three categories yielding some promising results to be further examined. 1.6 Thesis outline The structure of this document is organized as follows. In Chapter 2, a brief introduction to polarization demultiplexing and a more detailed explanation about Stokes space are presented. Some ideas regarding the representation of states of polarization in the Poincaré sphere are also given. Chapter 3 introduces Machine Learning by giving an overall view of all available techniques and their most popular applications. Focusing more on clustering, six different algorithms are described: K-Means, Expectation Maximization, DBSCAN, OPTICS, Spectral Clustering and Gravitational Clustering. Finally, two methods to analyze the clustering outcomes are mentioned. Chapter 4 is devoted to describe the novel proposed algorithm in detail, specially focusing in the characteristics of the signals used to identify them and mentioning the advantages introduced. In Chapter 5, after stating the simulation setup used to contrast the methods, the results obtained are assessed in terms of OSNR performance, accuracy and complexity. The presented algorithm is also compared in these three categories. Finally, in chapter 6, we summarize the main conclusions of the work done and suggest some lines on future work. 22 23 2 Polarization demultiplexing Polarization multiplexing is a type of optical transmissions that uses two orthogonal polarization states, as shown in Figure 2.1, to transmit two different optical signals at the same wavelength. It is a widely used method to double the spectral efficiency. Both polarization states can be excited independently and, ideally, would propagate separately. However, the state of polarization (SOP) rotates along the fiber and this causes the signal to arrive misaligned to the receiver, with both states mixed. This is mainly a consequence of birefringence, which is the optical property by which materials show a different refractive index depending on the polarization and direction of light. Other fiber impairments, like polarization mode dispersion (PMD) – polarization states travel at different speeds –, differential group delay (DGD) – difference in propagation time between X and Y polarization states – or polarization dependent losses (PDL) – difference in losses between both polarization states – are also responsible for the SOP rotation. In this scenario, a demultiplexing process is required in order to recover the original transmitted optical signal and, therefore, the information sent. At the receiver, the incoming signal is split into two orthogonal optical waves by the polarization beam splitter (PBS), with a system of coordinates different from that of the transmitter. Thus, each of the two orthogonal optical signals at the receiver is actually a linear combination of the transmitted optical waves. EY EX Figure 2.1. Orthogonal polarization states 24 To perform polarization demultiplexing, several different techniques exist, such as the commonly used constant modulus algorithm (CMA). Due to the misalignment in polarization states introduced by the fiber, simple modulation schemes can be observed as multi-level signals at the receiver side. To avoid this, CMA introduces the idea of enforcing constant modulus to the received data in order to properly demultiplex the signal, by mapping resulting symbols with expected symbols. This is true only for typical modulation formats like phase-shift keying (PSK), where the modulus is constant for all symbols. A variation of this algorithm, radius directed equalization with multiple modulus algorithm (RDE-MMA), considers several levels to overcome this problem, allowing to deal with more advanced modulation formats. Polarization demultiplexing in Stokes space is another method for the same purpose, with the advantage it does not require knowledge of the modulation format and it is faster in convergence. Therefore, working in Stokes space is very suitable for Modulation Format Recognition purposes. 2.1 Polarization demultiplexing in Stokes space Any state of polarization of a polarized lightwave can be formally described through the Stokes parameters, and represented unequivocally with a different vector orientation in Stokes space. Given an arbitrary plane defined by 𝑥 and 𝑦, normal to the propagation axis of the fiber, 𝑧, any electromagnetic field can be described easily as shown in (2.1), where 𝑎𝑥 and 𝑎𝑦 are the amplitudes, and 𝜙𝑥 and 𝜙𝑦 are the phases of the two orthogonal components. where 𝐸 󰇍  =𝑒𝑥𝑥+𝑒𝑦𝑦 𝑒𝑥=𝑎𝑥𝑒𝑗(𝜔𝑡+𝜙𝑥) 𝑒𝑦=𝑎𝑦𝑒𝑗(𝜔𝑡+𝜙𝑦) (2.1) Jones calculus is a simple formalism introduced by R. Clark Jones (1941) for describing polarized light in vector format. Moreover, linear optical elements can also be represented by Jones matrices. The previous equation, (2.1), can be expressed in terms of Jones vectors as shown in (2.2), where the first term 1 √2 is used for normalization. 𝐸=1 √2(𝑒𝑥 𝑒𝑦)= 1 √2 (𝑎𝑥𝑒𝑗(𝜔𝑡+𝜙𝑥) 𝑎𝑦𝑒𝑗(𝜔𝑡+𝜙𝑦)) (2.2) 25 The effect of an assumed ideal optical fiber can be modelled by means of a unitary Jones matrix independent of the optical frequency as described in (2.3). The unitary condition states that the determinant of the matrix is 1, and this makes any orthogonal wave to preserve its orthogonality through the fiber. Thus, to recover the initial alignment of the signal it is necessary to compute the inverse matrix 𝐻−1. 𝐻=(𝑎 𝑏 −𝑏∗𝑎∗) (2.3) For example, launching the polarization state 𝐽=(1,0)𝑇 causes to have 𝐽′= (𝑎,−𝑏∗) at the end of the fiber. Applying 𝐻−1 to this signal yields the initial Jones vector transmitted, since |𝑎|2+|𝑏|2=det(𝐻)=1. The derivation is shown mathematically in (2.4). We successfully recover the initially transmitted data. 𝐽𝑟𝑒𝑐=𝐻−1[ 𝐻 𝐽𝑡𝑥 ]=(𝑎∗−𝑏 𝑏∗𝑎)[(𝑎 𝑏 −𝑏∗𝑎∗)(10)]= =(|𝑎|2+|𝑏|2 𝑎𝑏∗−𝑏∗𝑎)=(10) (2.4) We have seen that only by identifying the polarization states of the signal at the receiver it is possible to construct the inverse matrix 𝐻−1 and thus the original wave. The equations in (2.5) define the Stokes parameters in terms of received data in Jones vector form. Stokes parameters were first introduced by Sir George Gabriel Stokes in 1852, and allow to represent any state of polarized light (even unpolarized), unlike Jones vector which only allows polarized light. They can be derived from the transmitted intensities when the beam is passed through different degree polarizers [18]. 𝑺=(𝑆0 𝑆1 𝑆2 𝑆3)=12 ( 𝑒𝑥𝑒𝑥∗+𝑒𝑦𝑒𝑦∗ 𝑒𝑥𝑒𝑥∗−𝑒𝑦𝑒𝑦∗ 𝑒𝑥∗𝑒𝑦+𝑒𝑥𝑒𝑦∗ −𝑗𝑒𝑥∗𝑒𝑦+𝑗𝑒𝑥𝑒𝑦∗ ) =12 ( 𝑎𝑥2+𝑎𝑦2 𝑎𝑥2−𝑎𝑦2 2𝑎𝑥𝑎𝑦cos(𝜙𝑦−𝜙𝑥) 2𝑎𝑥𝑎𝑦sin (𝜙𝑦−𝜙𝑥) ) (2.5) In (2.5) Stokes parameters are shown without normalizing, which is generally done with respect to 𝑆0. We can also note that the parameters are independent of the carrier frequency 𝜔 due to they operate in interpolarization differences. For the same reason, they are also independent of the polarization mixing and the phase offset. Each parameter conveys different information:  𝑆0 represents the total power of the light beam  𝑆1 represents 0º linear polarized light  𝑆2 represents 45º linear polarized light  𝑆3 represents circularly polarized light 32 𝝁𝑘=1 |𝑀𝑘|∑𝒙(𝑖) 𝑖∈𝑀𝑘 (3.2) c. Check for convergence. The cost function (3.3) is computed as the average distance from each point to its assigned cluster centroid. When the difference between current and last iteration values of this function goes below a specified threshold, the algorithm is said to converge. Another possible definition for convergence is determining if none of the points has changed its label (cluster) from previous assignments. 𝐽=1 𝑁∑∑‖𝒙(𝑗)−𝝁𝑘‖ 𝑗∈𝑀𝑘 𝐾 𝑘=1 (3.3) As stated above, one of the main advantages of K-Means is its simplicity, measured in number of operations as 𝑜(𝑁), where N is the number of training samples. This is why it is typically used for initializing more complex algorithms, such as Expectation Maximization, or for clustering an already transformed data space as with Spectral Clustering. Figure 3.2. Example of K-Means performance on a 2D space with two clusters. However, K-Means also has some drawbacks. One of them is caused by the random initialization, which can make the algorithm fall into a local optimum of the cost function instead of the global minimum. To circumvent this problem a commonly used solution is to run several initializations and choose the clustering produced by the one yielding a smaller cost function. This clearly increases the computational -2 0 2 0 2 4-2-4 µ1 µ2 33 complexity but increases the reliability. Another popular problem presented by KMeans is the need to previously specify the number of clusters, which in our case is the desired outcome. This can be solved by running the algorithm with every number of clusters targeted, each one corresponding to a different modulation format, and selecting the one that provides better results in terms of some figure of merit like the cost function or any other clustering indices such as Silhouette, which is further explained in Section 3.2. Finally, K-Means only allow a certain shape for clusters. For instance, if we had two ring-shaped clusters, one inside the other, it would be impossible to separate them using this algorithm. Selection of initial partition with #clusters = K START Cluster assignment Recompute centroids Convergence achieved? No Yes END Figure 3.3. K-Means flow chart. 3.1.2 Expectation Maximization (EM) The Expectation Maximization algorithm (Dempster et al, 1977) can be seen as a maximum-likelihood estimator with some latent variables, i.e. variables whose value is unknown [23]. In the most usual case it is applied to estimate parameters from data that can be fit to a simple linear superposition of Gaussian components, also called Gaussian Mixture Model (GMM). The goal in this case is to find the Gaussian parameters (means, covariance matrices and mixing coefficients) that define the model, as shown in Figure 3.4. 34 In the following paragraphs we assume we have a data set {𝒙1,…,𝒙𝑁} consisting of N observations of a random D-dimensional Euclidean variable x. We also define a Kdimensional binary random variable z where 𝑧𝑖𝜖{0,1}, indicating to which cluster the current sample belongs. Using the previous defined notation a GMM can be described as in (3.4), where 𝜋𝑘 symbolize the mixing coefficients, as described in (3.5). 𝑝(𝑥)=∑𝜋𝑘𝒩(𝑥 | 𝜇𝑘,Σ𝑘) 𝐾 𝑘=1 (3.4) Figure 3.4. Example of Expectation Maximization applied with a Gaussian Mixture Model formed by two clusters. The colors of the points represent the choice if a hard clustering assignment was done. The EM algorithm is based on two main steps: an expectation step (E), where given the actual estimated values of the parameters defining the model the probability of each point belonging to each cluster is computed, and a maximization step (M), where new values for the model parameters are computed given these probabilities. Both steps are repeated until convergence. There are some important concepts to understand EM:  Mixing coefficient (πk): They define the weight of each Gaussian in the total mixture. It is normally computed as the ratio of points for each cluster. 𝜋𝑘=𝑝(𝑧𝑘𝑖=1) ∀𝑘=1:𝐾 (3.5)  Responsibility (γk): It is the probability of each point belonging to each cluster. 𝛾𝑘=𝑝𝑟𝑜𝑏(𝑧𝑘=1|𝒙) ∀𝑘=1:𝐾 (3.6) Schematically the procedure for a GMM is the following (depicted in Figure 3.5): 1. Initialization of the parameters: means µk, covariances Σk and mixing coefficients πk. Also compute the value of the log likelihood. 2. E step: Evaluate the responsibilities using current parameters values, (3.7). -2 0 2 0 2 4-2-4 µ1 µ2 σ1 σ2 35 𝛾(𝑧𝑛𝑘)=𝜋𝑘𝒩(𝑥𝑛|𝝁𝑘,𝚺𝑘) ∑𝜋𝑗𝒩(𝑥𝑛|𝝁𝑗,𝚺𝑗) 𝐾 𝑗=1 (3.7) 3. M step: Estimate parameters using current responsibilities, (3.8). 𝝁𝑘𝑛𝑒𝑤=1 𝑁𝑘∑𝛾(𝑧𝑛𝑘)𝒙𝑛 𝑁 𝑛=1 𝚺𝑘𝑛𝑒𝑤=1 𝑁𝑘∑𝛾(𝑧𝑛𝑘)(𝒙𝑛−𝝁𝑘𝑛𝑒𝑤)(𝒙𝑛−𝝁𝑘𝑛𝑒𝑤)𝑇 𝑁 𝑛=1 𝜋𝑘𝑛𝑒𝑤=𝑁𝑘 𝑁 where 𝑁𝑘=∑𝛾(𝑧𝑛𝑘) 𝑁 𝑛=1 (3.8) 4. Evaluate the log likelihood and check for convergence. If the algorithm has not converged, steps 2, 3 and 4 are repeated until convergence is achieved, (3.9). ln𝑝(𝑿|𝝁,𝚺,𝝅)=∑ln{ ∑𝜋𝑘𝒩(𝑥𝑛|𝝁𝒌,𝚺𝒌) 𝐾 𝑘=1 } 𝑁 𝑛=1 (3.9) The initialization step is fundamental in order to find the best choice of parameters, those that fit best the data set. Depending on the initialization the algorithm can get trapped in a local optimum far from the optimum global solution and, therefore, result in a bad clustering of the input data. For this reason, K-Means is commonly used to initialize the model. The convergence decision also determines the outcome of the algorithm. It consists in checking if the algorithm has converged by fixing a certain threshold for the log likelihood function and checking it iteratively. The threshold setting is thus a very important variable in order to ensure the correct operation of the algorithm. EM is often compared to K-Means as they both are very similar. However, while K-Means performs a hard assignment of each point to a certain cluster, EM provides a set of probabilities (called responsibilities) for each point to belong to each cluster. This procedure is normally called soft assignment. In [14] there is the derivation to express K-Means as a limit of EM. 36 Initialization of parameters (µk, Σk, πk) START Expectation step γ(znk) Maximization step (µknew, Σknew, πknew) Convergence achieved? No Yes END Log likelihood evaluation ln p(X|µ, Σ, π) Figure 3.5. EM flow chart. 3.1.3 Density Based Spatial Clustering of Application with Noise (DBSCAN) Some of the usual problems when dealing with clustering algorithms are: i) it is often difficult to determine the input parameters to be used for a specific database; ii) they are normally very computationally costly; and iii) they are restricted to only some cluster shapes, as shown in Figure 3.7. The DBSCAN algorithm (Ester et al, 1996) is a density-based technique that offers a solution to these three main problems [25]. The algorithm consists in identifying the clusters by referring to the density of the database elements, with minimal knowledge of the domain by the user. Moreover, it differentiates points of noise and reliable points to be included in a cluster. To understand the procedure it is necessary to introduce some definitions before:  ε-neighborhood of a point: It comprises all the points within a distance below than ε from the point.  Directly density-reachable: One point 𝑞 is directly density-reachable from another point 𝑝 if 𝑞 belongs to the ε-neighborhood of 𝑝 and there are more than 37 𝑀𝑖𝑛𝑃𝑡𝑠 points in this neighborhood. Then 𝑝 is called a core point and 𝑞 is a border point. This can be seen in Figure 3.6 (a).  Density-reachable: One point 𝑞 is density-reachable from another point 𝑝 if there is a succession of points 𝑝1,…,𝑝𝑛, such that 𝑝1=𝑝 and 𝑝𝑛=𝑞 and 𝑝𝑖+1 is directly density-reachable from 𝑝𝑖. It is depicted in Figure 3.6 (b).  Density-connected: A point 𝑞 is density-connected to a point 𝑝 if there is a point 𝑜 such that from it 𝑝 and 𝑞 are both density-reachable given ε and 𝑀𝑖𝑛𝑃𝑡𝑠. The density-connected property can be observed in Figure 3.6 (c). q p q p p q o (a) (b) (c) Figure 3.6. Definitions necessary for density-based algorithms with respect to MinPts=3. (a) p is a core point, q is a border point. q is directly density reachable from p; (b) q is density reachable from p, but p is not density reachable from q; (c) p and q are density-connected to each other by o. With these definitions it is possible to define a cluster as the aggregation of points such that for each pair of points in the cluster they are density-connected with respect to ε and 𝑀𝑖𝑛𝑃𝑡𝑠. If one point 𝑝 belongs to a cluster and another point 𝑞 is densityreachable from 𝑝 given ε and 𝑀𝑖𝑛𝑃𝑡𝑠, then 𝑞 also belongs to the same cluster. Noise points are identified as the set of points in the dataset that do not belong to any of the clusters. -2 0 2 0 2 4-2-4 Figure 3.7. Data set recognizable using a density-based method, non-detectable by other techniques like K-Means. 38 The algorithm DBSCAN is briefed as follows (also in Figure 3.8): 1. Arbitrary select a training sample 𝒑 of the data set and retrieve all points density-reachable from it given ε and 𝑀𝑖𝑛𝑃𝑡𝑠. 2. If 𝑝 is a core point (more than 𝑀𝑖𝑛𝑃𝑡𝑠 points in its neighborhood), then expand a new cluster: a. Obtain the neighbors in the neighborhood and assign them to the same cluster. b. For each neighbor still not classified or marked as noise, check for its neighborhood and repeat the same process in a. and b. Each of these points will be density reachable from 𝑝. Otherwise, no points are density-reachable from 𝑝 and the algorithm selects the next training sample of the data set. 3. DBSCAN finishes when it reaches the last sample in the data set. All points have to be visited. Arbitrary select a starting sample, p START Expand Cluster ClusterId ClusterId = 1 Unclassified point? Last point of the data set? Retrieve next point in the data set Increase ClusterId ClusterId+=1 END Core point? Mark point as Noise Yes Yes Yes No No No Figure 3.8. Flow chart of DBSCAN. The value of ε and 𝑀𝑖𝑛𝑃𝑡𝑠 can be determined dynamically depending on the signal, as explained in [15, section 4.2], by obtaining, for each point the distance to the k nearest neighbor and looking for the elbow of this function, which corresponds to ε. 39 𝑀𝑖𝑛𝑃𝑡𝑠 may be computed by taking the average number of neighbors each point has at distance ≤𝜀. 3.1.4 Ordering Points to Identify the Clustering Structure (OPTICS) OPTICS (Ankerst et al, 1999) is another density-based clustering method, a more generalized view of the DBSCAN algorithm [26]. It introduces the idea of, instead of hard-assigning clusters, producing a special order of the database containing all the valuable information to extract a certain clustering structure. In OPTICS the same definitions given before for DBSCAN and shown in Figure 3.6 are applied together with two new ones, core distance (CDist) and reachability distance (RDist) with respect to another object, as defined in (3.10) and (3.11), respectively. On one hand, core distance (3.10) is the smallest distance ε’ between p and an object of its ε-neighborhood so that p would be considered a core object with respect to 𝜀′<𝜀. On the other hand, reachability distance (3.11) with respect to o can be thought as the minimum needed distance for point p to be directly density-reachable from o if o is a core object. Figure 3.9 shows graphically these definitions. Figure 3.9. Core distance and reachability distances with MinPts=4. RDist(p, o) exemplifies the case of RDist equal to CDist, while RDist(p’, o) represents the second explained case. 𝐶𝐷𝑖𝑠𝑡(𝑝)={𝑢𝑛𝑑𝑒𝑓𝑖𝑛𝑒𝑑,𝑖𝑓 |𝑁𝜀(𝑝)|<𝑀𝑖𝑛𝑃𝑡𝑠 𝑀𝑖𝑛𝑃𝑡𝑠−𝑑𝑖𝑠𝑡𝑎𝑛𝑐𝑒(𝑝),𝑜𝑡ℎ𝑒𝑟𝑤𝑖𝑠𝑒 (3.10) 𝑅𝐷𝑖𝑠𝑡(𝑝,𝑜)={ 𝑢𝑛𝑑𝑒𝑓𝑖𝑛𝑒𝑑,𝑖𝑓 |𝑁𝜀(𝑝)|<𝑀𝑖𝑛𝑃𝑡𝑠 max( 𝐶𝐷(𝑝),𝑑𝑖𝑠𝑡𝑎𝑛𝑐𝑒(𝑝,𝑜) ),𝑜𝑡ℎ𝑒𝑟𝑤𝑖𝑠𝑒 (3.11) The ordering that OPTICS produces consists in storing each point together with its core-distance and its reachability-distance. The procedure is similar to that of DBSCAN, as shown in Figure 3.10. 1. For each unvisited point in the data set, compute its core distance and set its reachability distance as undefined. 2. If the core distance has a defined value, fill a list with its neighbors. For each element in the list, compute the core distance and if it is not RDist(p’,o) oCDist(o) p’ p RDist(p,o) ε 40 undefined, add its neighbors to the list. Retrieve the next point in the list. Repeat this process until the list is empty. The process of adding items to the list consists in computing the reachability distance. If the point already contained a reachability distance and the new one is smaller, or none was defined, this value is updated. Arbitrary select a starting sample, p START ClusterId = 1 Unprocessed point? END Compute CDist and set RDist as undefined CDist != undefined? Update list with neighbours (also compute Rdist) Retrieve next point in the data set Seeds list empty? Retrieve next point in the seeds list All data set visited? No No Yes No Yes No Yes Yes Figure 3.10. OPTICS flow chart. With the information given by OPTICS (core distance and reachability distance for each point) it is possible to extract the associated clustering structure. The idea comes from the reachability plot, which is a representation of the reachability distance for each point in the data set where is possible to differentiate clusters by identifying the valleys. The OPTICS reference paper [26] shows this diagram and gives an algorithm to assign the points to its corresponding cluster given the ordered list and the input parameters ε’ (assuming 𝜀′≤𝜀) and 𝑀𝑖𝑛𝑃𝑡𝑠. 41 3.1.5 Spectral Clustering Spectral Clustering is one of the most recent clustering algorithms [27], [28]. It offers some advantages in front of more traditional clustering techniques, as for instance it does not require the assumption of a certain density function like EM, assuming a Gaussian mixture model, and it is not vulnerable to local minima as K-Means. The main feature of spectral methods is to perform a different mapping producing a more compacted and tight clustering structure, over which a simple clustering method such as K-Means can be applied to easily identify clusters. This change in representation is obtained by computing as many eigenvectors as clusters to be found in the data set on a certain distance-dependent matrix. (a) (b) Figure 3.11. Spectral clustering dimensionality reduction. Comparison before (a) and after (b) spectral technique, for 25 dB PDM BPSK signal. There are many variants of spectral techniques applied to clustering, different in the way they compute the eigenvectors and on which matrix they apply. The work here is done according to Ng, Jordan and Weiss [28]. The starting point is to build an undirected graph weighted with some distance measure between points. With this in mind, some definitions are needed:  Weighted adjacency matrix (W): In this matrix (NxN), the value in 𝑤𝑖𝑗 represents a certain distance measure between points i and j. The matrix is symmetric and is also called similarity or affinity matrix 𝑊=(𝑤11 ⋯ 𝑤1𝑁 ⋮ ⋱ ⋮ 𝑤𝑁1 ⋯ 𝑤𝑁𝑁) (3.12)  Degree matrix (D): It is a matrix (NxN) containing in its diagonal the degree of the corresponding vertices (3.13). The degree of a vertex is defined as the sum of all its adjacent edges weights. 48 This coefficient gives results in the range between -1 and 1, where ±1 correspond to a very tight clustering structure with high inter-cluster distance with respect to the intra-cluster, and 0 occurs when both distances are very similar and no appreciable difference can be seen between clusters. The number of clusters is decided by choosing the one that gives a higher average Silhouette coefficient for all points. In Figure 3.16 we can observe the Silhouette coefficient plot for the detection of a 20 dB PDM BPSK signal with K-Means. Two clusters are detected in the data set, which effectively corresponds to a PDM BPSK signal in Stokes space. Figure 3.16. Modulation format recognition using Silhouette coefficient. On the other hand, we also want to quantify how good a certain clustering method is performing in terms of the percentage of points assigned to the ideal cluster. For this purpose we need a tool to compare the optimum label set with the clustered one that does not account for permutations, i.e. we cannot control if the numbers of ideal clusters will identically match the clustered ones but we want to recognize them as the same group of points. For example, point a can ideally belong to cluster A and be assigned to cluster B, but in fact cluster A (before) and B (after) are the same (contain same points). One of these metrics is Adjusted Mutual Information (AMI). AMI measures the agreement between two sets of labels, ignoring permutations. It does not assume any predefined cluster structure and has the advantage that random assignments always give close-to-zero values – it is normalized against chance. To compute the AMI index between two partitions of the data set, 𝑈= 𝑈1,…,𝑈𝑅 and 𝑉=𝑉1,…,𝑉𝐶 with different number of clusters R and C, first we have to build a contingency table as defined in (3.24). 0 0,1 0,2 0,3 0,4 0,5 0,6 0,7 0,8 0,9 1 2 4 8 16 32 60 Silhouette coefficient Number of clusters Silhouette coefficient for 20 dB BPSK signal 49 𝑚𝑖𝑗=|𝑈𝑖∩𝑉𝑗| ∀𝑖=1:𝑅,∀𝑗=1:𝐶 (3.24) From information theory, the mutual information between two random variables can be defined as in (3.25) with 𝑝(𝑖) being the probability that an object is assigned into 𝑈𝑖 and 𝑝(𝑖,𝑗) the probability that an object belongs to both 𝑈𝑖 and 𝑉𝑗. 𝑀𝐼(𝑈,𝑉)=∑∑𝑝(𝑖,𝑗)·log(𝑝(𝑖,𝑗) 𝑝(𝑖)𝑝(𝑗)) |𝑉| 𝑗=1 |𝑈| 𝑖=1 (3.25) Vinh, Epps and Bailey [31] derived the formula (3.26) to find the expected mutual information between two sets of assignments. With this expected value it is possible to compute the AMI value as written in (3.27), where H(U) and H(V) represent the entropies of both partition sets. 𝐸[𝑀𝐼(𝑈,𝑉)]=∑∑ ∑ 𝑛𝑖𝑗 𝑁log(𝑁·𝑛𝑖𝑗 𝑎𝑖𝑏𝑗)𝑎𝑖!𝑏𝑗!(𝑁−𝑎𝑖)!(𝑁−𝑏𝑗)! 𝑁!𝑛𝑖𝑗!(𝑎𝑖−𝑛𝑖𝑗)!(𝑏𝑗−𝑛𝑖𝑗)!(𝑁−𝑎𝑖−𝑏𝑗+𝑛𝑖𝑗)! min (𝑎𝑖,𝑏𝑗) 𝑛𝑖𝑗=max (𝑎𝑖+𝑏𝑗−𝑁,0) 𝐶 𝑗=1 𝑅 𝑖=1 (3.26) 𝐴𝑀𝐼= 𝑀𝐼−𝐸[𝑀𝐼] max(𝐻(𝑈),𝐻(𝑉))−𝐸(𝑀𝐼) (3.27) The higher the AMI value is, the greater the percentage of points correctly assigned to their cluster is. Therefore, we will use this index to compare the algorithms in terms of reliability, when the algorithm is capable of detecting the conveyed modulation format. 50 51 4 Proposed algorithm for modulation format recognition based on Maximum Likelihood in Stokes space In this section we propose a novel algorithm to be applied on modulation format recognition for optical communications. The algorithm works in Stokes space and is especially suited for three dimensional data representable in the Poincaré sphere. 4.1 Proposed algorithm Stokes space transformed data display a different number of clouds of points (clusters) for different polarization multiplexed modulation formats. However, these clusters do not maintain a fixed position in the Poincaré sphere, but they can suffer from the rotation introduced by effects such as birefringence or other fiber impairments, as explained in Chapter 2. Specifically, the state of polarization rotates around the axe defined by S1, which represents the linear polarization, as depicted in Figure 1.2. S1 S3 S2 Figure 4.1. Rotation suffered in the Poincaré sphere around S1. 52 If we could determine the angle of rotation the received signal has with respect to a certain initial position of the centroids, then we would be able to identify the transmission parameters by comparing the incoming signal with the stored centroids positions rotated by the optimum angle computed. This is essentially what the proposed algorithm does. In order to compute the angle of rotation, we can select the value which minimizes the Euclidean distance between the input data and the rotated version of the centroids for each modulation. In this proposal, the optimization is done with respect to QPSK centroids as they are common for all the modulation formats targeted, except for BPSK, in which case the rotation is also optimal since the algorithm tends to maintain the two also existing clusters in QPSK and leave two unpopulated clusters. QPSK centroids are located at the four points of the outer surface of the sphere which intersect with the plane defined by 𝑆1=0, as shown in Figure 4.2. S0 S1 S2 S3 µ1 1 0 1 0 µ2 1 0 -1 0 µ3 1 0 0 -1 µ4 1 0 0 1 Figure 4.2. QPSK initial centroids' positions used as a reference to obtain the optimum angle. This technique requires a database storing the initial position of the clusters for the different modulation formats targeted, which is fixed and does not depend on the received signal. The algorithm proposed is summarized in the following steps, and it is schematically displayed in Figure 4.4. 1. We compute the Stokes parameters for the received data. This results in a rotated version of the clusters in the Poincaré sphere. 2. The rotation angle of the incoming signal is obtained by optimizing the Euclidean distance between the input data, 𝒙𝒊 and the QPSK rotated centroids, 𝝁𝒋(𝜃), as described in (4.1). 𝜃=argminθ{ ‖𝒙𝒊−𝝁𝒋(𝜃)‖2 } ∀𝑖=1:N,∀𝑗=1:4 (4.1) 53 where 𝝁𝑗(𝜃) is obtained as follows: 𝝁𝒋(0)=[𝑆1,𝑆2,𝑆3] 𝑪=(𝑆2,𝑆3)→𝑪(𝜃)=𝑪𝑒𝑗𝜃=(𝑆2′,𝑆3′) 𝝁𝒋(𝜃)=[𝑆1,𝑆2′,𝑆3′] The optimization is computed as an unconstrained nonlinear optimization. 3. For each modulation format targeted, the corresponding stored cluster centroids are rotated by the obtained optimum angle 𝜃 . This step could also be done by rotating the data instead, which would be more computationally expensive. 4. For each point in the received data set, we perform a hard assignment, as described in (4.2), to the closest centroid rotated by 𝜃 according to the modulation analyzed, and store the labels given to the samples in 𝒄. 𝑐(𝑖)=𝑎𝑟𝑔 𝑚𝑖𝑛𝑗‖𝒙𝑖−𝝁𝑗 (𝜃)‖ ∀𝑗=1:𝐾 (4.2) 5. Once all the assignments have been done according to the different modulation formats analyzed, we decide which clustering structure fits best to the data set using the Silhouette coefficient. This metric chooses the format which maximizes the relation between the inter-cluster distance and intra-cluster distance as described in Section 3.2. Figure 4.3. Poincaré sphere for a 20 dB QPSK signal with π/16 radians rotation. 54 START Rotation angle optimization compared to QPSK Rotate centroids by optimal angle Assignment of points into closest clusters Evaluate cluster structure using Silhouette END Figure 4.4. Proposed algorithm flow chart. 4.2 Conclusions The technique presented in this section was initially thought to classify between six different polarization multiplexed modulations (BPSK, QPSK, 8PSK, 8QAM, 12QAM and 16QAM). Nevertheless, it can be easily adapted to detect more modulations by simply updating the knowledge database with the location of centroids for the new formats. The algorithm presented benefits from the advantages the recent conceived Stokes space-based DSPs include, like polarization demultiplexing, compensation for chromatic dispersion and cross-polarization modulation, or polarization dependent losses, as explained in Chapter 2. It takes profit as well from the benefits of treating signals in Stokes space, which make signals independent of polarization mixing, carrier frequency offset or phase offset. 55 5 Simulation The simulation results obtained by each algorithm are described in this section. First we will state the common simulation setup used for all methods and the conditions under which each test has been done. Next, we will justify the OSNR range in which the analysis has been made, and finally, the results will be commented. 5.1 Simulation setup In this section, the previously defined clustering algorithms have been contrasted with the novel proposed algorithm, except for gravitational clustering, which did not yield good results as it is extremely sensitive to the parameters selection. They have been analyzed in three different categories (Table 5.1) in order to make a fair comparison between them. Firstly, they have been compared in terms of the OSNR lower bound from which they start to detect the modulation format correctly. Another comparison metric used is the reliability in terms of the Adjusted Mutual Information (AMI) index, and finally the clustering time spent by each technique to apply modulation format recognition. Parameter Unit OSNR lower bound dB AMI % Clustering time min 2000 samples per signal Table 5.1. Categories to be simulated for each modulation format. The modulation formats tested are all dual-polarized: BPSK, QPSK, 8PSK, 8QAM, 12QAM and 16QAM. The signals have been tested in a different OSNR range depending on their nature; from 5 to 25 dB for PSK signals and from 15 to 35 dB for QAM signals, in 1 dB steps. QAM signals are more advanced modulations and need a higher OSNR in order to be demodulated correctly by the receiver. All simulations have been performed with a simulated received signal of only 2000 samples. 56 In Table 5.2 the specific parameters needed for each algorithm are shown. These values have been carefully selected according to the targeted domain – points in the Poincaré sphere – and algorithm used. As shown, the algorithm introduced in Section 4 does not require any input parameter. Analyzed algorithm Parameters K-Means MaxIter = 25 Expectation Maximization MaxIter = 25 DBSCAN Eps =0.08; MinPts=25; OPTICS MinPts = 25; Eps2 = 0.02; Spectral Clustering Sigma=0.4; Proposed algorithm - Table 5.2. Parameters selection for the simulation test. All the simulations have been run in the same computer, equipping Windows 8 and MATLAB R2013b (Intel® Core™2 Duo CPU P8700 @ 2.53 GHz, with 6 GB RAM). 5.2 FEC limit Not all the modulation formats results should be analyzed in the same way, since the number of possible symbols per modulation is different, and thus the complexity to properly acquire and demodulate signals varies. In order to set a reference for all the formats we will use the forward error correction (FEC) threshold as a bound. This limit represents the maximum bit error rate (BER) allowed in received signals to be recovered, considering we have transmitted them using a 7% overhead FEC code [32]. It corresponds to a BER equal to 3.8e-3 (- 2.42 dB). The FEC limits in terms of minimum OSNR required for different modulation formats, considering a 28Gbaud signal, can be seen graphically in Figure 5.1 and are described in Table 5.3. BPSK, QPSK and 16QAM values are computed theoretically, while 8PSK, 8QAM and 12QAM are simulated. BPSK QPSK 8PSK star 8QAM 12QAM 16QAM 9.73 dB 12.04 dB 17.38 dB 17.14 dB 17.93 dB 18.75 dB Table 5.3. FEC threshold for different modulation formats (28 Gbaud, [email protected]). 57 Figure 5.1. FEC limit in terms of the minimum required OSNR for different modulations, computed with 28Gbaud at 0.1nm linewidth. 5.3 Simulation results As mentioned in Table 5.1, three different analyses concerning the same input data have been carried out. The first study, shown in Figure 5.2, regards the minimum optical signal to noise ratio needed for the algorithms to correctly detect the format of the received signal. In all cases, the value shown in the results corresponds to the first OSNR in the simulation range that succeeds in the detection of the correct number of clusters. This does not mean that all subsequent OSNR values are also correctly recognized, although a high percentage of them are. AMI index will help to reflect this event with low values. In the graphic the FEC limit corresponding to each modulation format is also displayed as a reference. As we can observe, most algorithms except OPTICS, manage to recognize PSK signals correctly in almost all the OSNR range considering the threshold imposed by the FEC limit. The reason OPTICS might fail, reminding that DBSCAN is nothing more than a more specific version of OPTICS, is that the choice of parameters is often difficult and clearly determines the outcome of the algorithm. For QAM signals the results are not as close to the FEC limit, since the number of clusters is higher and it is more difficult to separate them. Density-based methods – DBSCAN and OPTICS – fail to identify 16QAM in the OSNR range considered. According to Figure 5.2, the proposed algorithm and Spectral clustering offer the best balance among all the clustering methods analyzed, by giving OSNR lower bounds very close to the FEC limit for all modulation formats. BPSK 9,73 dB QPSK 12,04 dB 8PSK 17,38 8QAM 17,14 12QAM 17,93 16QAM 18,75 -6 -5,5 -5 -4,5 -4 -3,5 -3 -2,5 -2 -1,5 -1 2 4 6 8 10 12 14 16 18 20 log(BER) OSNR (dB) FEC limit BPSK QPSK 8PSK 8QAM 12QAM 16QAM 64 [15] R. Borkowski, D. Zibar, A. Caballero, V. Arlunno, and I. T. Monroy, “Stokes Space-Based Optical Modulation Format Recognition for Digital Coherent Receivers,” IEEE Photonics Technol. Lett., vol. 25, no. 21, pp. 2129–2132, Nov. 2013. [16] J. Liu, Z. Dong, K. P. Zhong, A. P. T. Lau, C. Lu, and Y. Lu, “Modulation Format Identification Based on Received Signal Power Distributions for Digital Coherent Receivers,” in Optical Fiber Communication Conference, 2014, p. Th4D.3. [17] E. J. Adles, M. L. Dennis, S. Member, W. R. Johnson, T. P. Mckenna, C. R. Menyuk, J. E. Sluz, R. M. Sova, M. G. Taylor, and R. A. Venkat, “Blind Optical Modulation Format Identification From Physical Layer Characteristics,” J. Light. Technol., vol. 32, no. 8, pp. 1501–1509, 2014. [18] A. Kumar and A. Ghatak, Polarization of Light with Applications in Optical Fibers, Illustrate. SPIE Press, 2011, pp. 75–135. [19] B. Szafraniec, B. Nebendahl, and T. Marshall, “Polarization demultiplexing in Stokes space.,” Opt. Express, vol. 18, no. 17, pp. 17928–39, Aug. 2010. [20] C. M. Bishop, “Pattern Recognition and Machine Learning (Information Science and Statistics),” Aug. 2006. [21] A. K. Jain, M. N. Murty, and P. J. Flynn, “Data clustering: a review,” ACM Computing Surveys, vol. 31, no. 3. pp. 264–323, 1999. [22] J. Macqueen, “Some methods for classification and analysis of multivariate observations,” Proc. Fifth Berkeley Symp. Math. Stat. Probab., vol. 1, pp. 281–297, 1967. [23] A. P. Dempster, N. M. Laird, and D. B. Rubin, “Maximum Likelihood from Incomplete Data via the EM Algorithm,” J. R. Stat. Soc. Ser. B, vol. 39, pp. 1–38, 1977. [24] C. M. Bishop, Pattern Recognition and Machine Learning, vol. 4. Springer-Verlag New York, Inc., 2006, pp. 1–738. [25] M. Ester, H. P. Kriegel, J. Sander, and X. Xu, “A Density-Based Algorithm for Discovering Clusters in Large Spatial Databases with Noise,” in Second International Conference on Knowledge Discovery and Data Mining, 1996, pp. 226–231. [26] M. Ankerst, M. M. Breunig, H.-P. Kriegel, and J. Sander, “OPTICS: ordering points to identify the clustering structure,” ACM SIGMOD Rec., vol. 28, pp. 49–60, 1999. [27] U. Luxburg, “A tutorial on spectral clustering,” Stat. Comput., vol. 17, no. 4, pp. 395–416, Aug. 2007. [28] A. Y. Ng, M. I. Jordan, and Y. Weiss, “On spectral clustering: Analysis and an algorithm,” in Advances in Neural Information Processing Systems, 2001, pp. 849–856. [29] W. E. Wright, “Gravitational clustering,” Pattern Recognition, vol. 9, no. 3. pp. 151–166, 1977. [30] P. J. Rousseeuw, “Silhouettes: A graphical aid to the interpretation and validation of cluster analysis,” Journal of Computational and Applied Mathematics, vol. 20. pp. 53–65, 1987. [31] N. X. Vinh, “Information Theoretic Measures for Clusterings Comparison : Variants , Properties , Normalization and Correction for Chance,” J. Mach. Learn. Res., vol. 11, pp. 2837–2854, 2010. [32] I.-T. S. Group, “ITU-T Rec. G.975.1 (02/2004) Forward error correction for high bit-rate DWDM submarine systems,” Standard, pp. 1–58, 2005.