Member, IEEE Vincent Gouldieff

Member, IEEE Jacques Palicot

Member, IEEE Steredenn Daumont

Blind Modulation Classification for Cognitive Satellite in the Spectral Coexistence context

Keywords: Blind Signal Processing, Cognitive Satellite, Coexistence Scenario, Automatic Modulation Classification, M th -Power nonlinear Transformation

The scarcity of the usable satellite spectrum has led to a need for advanced communication techniques. Among them, multi-satellite networks have been promoted for their ability to reuse the frequency bands in a more aggressive manner. This is made possible by sharing both the spectral and spatial degrees of freedom. Meanwhile, the concept of Cognitive Satellite (CS) was deeply investigated to make the best of these novel architectures. In this context, very few blind signal processing methods were proposed to tackle the co-channel interference issues entailed by such systems. Herein, we develop a novel low-complexity algorithm for blind Automatic Modulation Classification (AMC) and parameters estimation in the strongly interfering scenario. The proposed method rely on an Analytical study of the M th -Power nonlinear Transformation (AM PT) of the co-channel mixture. Theoretical analysis and extensive simulations show the effectiveness of the proposed method in various scenarios.

I. INTRODUCTION

N EXT generation of Satellite Communications (SatComs) systems mainly aims at enhancing the global network throughput while ensuring Quality of Service (QoS) [START_REF] Vidal | Next generation High Throughput Satellite system[END_REF]. The advent of multi-satellite and hybrid networks using multibeam technologies has led to numerous spectrum sharing techniques [START_REF] Alegre-Godoy | Offered Capacity Optimization mechanisms for Multi-beam Satellite Systems[END_REF]. To push the global capacity optimization into a corner, the concept of Cognitive Radio [START_REF] Palicot | Radio Engineering: from Software Radio to Cognitive Radio[END_REF] has been considered for SatComs to make the best use of the available diversities (e.g., spectral, temporal, spatial diversities). The paradigm of Cognitive Satellite [START_REF] Hoyhtya | Application of Cognitive Radio techniques to Satellite Communication[END_REF] was born. Thereby, strong co-channel interference may occur between the primary link and the secondary links at both the uplink and the downlink. For instance, such a scenario is encountered when two or more multibeam satellites share the same geographic region [START_REF] Sharma | Cognitive Beamhopping for Spectral Coexistence of Multibeam Satellites[END_REF] or when a cellular network is built on a dual-satellite system [START_REF] Christopoulos | User Scheduling for Coordinated Dual-Satellite Systems with Linear Precoding[END_REF].

In the cognitive context, several features of the interfering signals may be dynamically chosen so that the spectral efficiency is maximized [START_REF] Soysal | Channel Estimation and Adaptive M-QAM in Cognitive Radio Links[END_REF]. In order to cut the overhead down, these parameters are not transmitted and it becomes mandatory to estimate them at the receiver [START_REF] Palicot | A New Concept for Wireless Reconfigurable Receivers[END_REF]. Classical methods of blind signal processing applied on the entailed co-channel mixture rapidly show their limits. Their inability to correctly infer the Signal-of-Interest (SoI) is all the more true when V. Gouldieff and J. Palicot are with the Signal, Communication and Embedded Electronics (SCEE) Research Team, CentraleSupelec/IETR, Avenue de la Boulaie -CS 47601, 35573 Cesson-Sevigne Cedex, France. E-mail : vincent.gouldieff@supelec.fr, jacques.palicot@supelec.fr.

V. Gouldieff and S. Daumont are with Zodiac Data Systems SAS, 5 Esplanade Anton Phillips, 14460 Colombelles, France. Formerly, she was with the SCEE Research Team. E-mail : steredenn.daumont@zodiacaerospace.com. the composite power of the interferers grows. Actually, while the SINR (Signal to Interference-plus-Noise Ratio) of the SoI weakens, the available statistical properties of the interferers are not exploited as useful information for the inference of the unknown features of the co-channel signals.

Consequently, classical 1 methods for Automatic Modulation Classification (AMC) applied on co-channel mixtures of signals perform poorly. This fact has recently motivated the scientific community to investigate the Multiuser2 AMC (MAMC) issue. In the blind scenario, MAMC has long been seen as a real challenge which may explain why the Multiuser case is so sparse in the AMC literature (cf. Section III).

Thereby, very few MAMC methods were proposed, whether it was considered in the single-node or in the distributed multinode contexts 3 . The latter would be preferred in interconnected networks provided for the purpose and in the military domain, where the data exchange between the nodes is not an issue. However, in the context of CR, the savings in terms of overhead would be mitigated by the consequent signaling between the nodes, needed to fusion the collected information. This fact has motivated the study of MAMC techniques in the terrestrial cognitive single-node context. To the best of our knowledge, Spooner was the first to propose an algorithm for MAMC in the single-node scenario, in 1996 [START_REF] Spooner | Classification of Co-channel Communication Signals using Cyclic Cumulants[END_REF]. Since then, it seems that no efficient and practically feasible algorithm has been proposed in the literature to face the blind single-node MAMC issue, which is precisely our goal.

Throughout this article, we actually are interested in the blind MAMC of asynchronously mixed co-channel signals received by a single antenna. More precisely, we are prone to consider the (almost-)totally spectral overlapping scenario. This case study is of real interest for three main reasons. First, the (almost-)totally overlapping scenario would appear in the context of Cognitive Satellite where time and frequency are slotted for an easier management of the resource. Secondly, the asynchronous assumptions can not be relaxed since perfect time/frequency alignment between the overlapping signals is generally not reachable. Finally, the considered scenario is a major issue in the field of MAMC: it actually seems that no proposed signal processing method is currently able to face it.

The MAMC algorithm we propose in this article relies on an Analytical study of the M th -Power Nonlinear Transformation (AM PT) [START_REF] Gouldieff | Blind Digital Modulation Classification based on Mth-Power Nonlinear Transformation[END_REF]. As further discussed, this transform is chosen for its ability to blindly recognize the constellations of the cochannel signals at a small computational expense and with a good robustness against the unknown parameters. The whole method is classically structured around a two step architecture compounded of a preprocessing unit and a recognition algorithm, as depicted in Fig. 1. The paper is organized as follows: In Section II, we describe the considered scenarios and formulate the working hypothesis. In Section III a brief review of (M)AMC methods is unfolded, while in Sections IV to VII, we develop the underlying theory of AM PT for AMC and MAMC. In Section VIII, the proposed AM PT-based algorithms in the co-channel context are detailed. Finally, the performance of the proposed methods is studied via extensive simulations throughout Section IX.

II. PROBLEM STATEMENT Thereby, we describe the scenarios entailing asynchronous co-channel interference of satellite signals, and we illustrate them with specific case-studies. Then, we derive the associated mathematical model and assumptions in the blind context.

A. System model and interference scenarios

Small satellite dishes (such as the one used in VSATs) are increasingly considered at the end users for their cost and space requirement. Due to their non-directivity, they willingly allow multi-satellite communications, but they also inherently induce unwanted interference at both the forward and the reverse links [START_REF] Weerackody | Sensitivity of Satellite Interference to Locations of Small-Aperture Terminals[END_REF]. Thereby, we consider Geostationary (GEO) bent-pipe satellites used for Fixed Satellite Services (FSS).

1) Scenario A -Interference at the end user (Forward link):

We may consider the scenario induced by the transmission of multiple gateways through closely located satellites. Interference occurs when their footprints and frequency plan overlap at the downlink (Fig. 2). It equivalently happens with a multibeam satellite when the frequency reuse policy is aggressive [START_REF] Chatzinotas | Joint Precoding with Flexible Power Constraints in Multibeam Satellite Systems[END_REF]. 3) Case-studies for Scenarios A and B in the literature: Scenario A captures the case of coexistent multibeam GEO satellites in Ka band, as shown in [START_REF] Sharma | Cognitive Beamhopping for Spectral Coexistence of Multibeam Satellites[END_REF] and [START_REF] Christopoulos | User Scheduling for Coordinated Dual-Satellite Systems with Linear Precoding[END_REF]. In these articles, it is proposed to mitigate co-channel interference via coordination between the gateways. This assumption is in practice very restrictive and then we relax it in the present work.

Scenario B accurately models the interference produced by cognitive SatComs at the reverse link, as shown in scenarios D and E of the CoRaSat approach [START_REF] Liolis | Cognitive Radio Scenarios for Satellite Communications: the CoRaSat Approach[END_REF]. Thereby, some portions of the C and Ku bands are used on a primary basis by FSS. Non-directive terminals may produce significant interference due to their opportunistic link with a secondary satellite.

B. Mathematical model of the involved co-channel mixture

We derive the expression of the mixture at the VSAT for the forward link (Scenario A). The derivation for Scenario B is similar and leads to an analogous model with the same mathematical structure. First, the uplink signal at gateway u is expressed in baseband as

x u (t) = k∈Z s u (k) • h u (t -kT u ), (1)
where s u is the i.i.d. symbol stream carried by constellation C u ∈ C, h u the pulse-shaping function and T u the symbol time.

The physical link between each gateway u and the VSAT receiver (including the uplink, the bent-pipe satellite and the downlink) is classically modeled for fixed satellite systems as a delayed single-path channel h p u (t) such that

h p u (t) = a u • e i(2πfut+φu) • δ(t -τ u ), (2)
where a u is the amplitude, f u the carrier frequency at the downlink, φ u the initial phase and τ u the delay. You may refer to Section II.C for the detailed assumptions related to the aforementioned parameters.

Denoting

x(t) = [x 1 , • • • , x U ] and h p u (t) = [h p 1 , • • • , h p U ]
T with U the number of signals sharing the same band, we get at the receiver the following co-channel mixture of signals:

x(t) = x(t) • h p u (t) + ν(t), (3)
where ν(t) stands for Additive White Gaussian Noise (AWGN) with unknown variance σ 2 ν .

At the receiver, signal x(t) is sampled (to output x(n)) and blindly preprocessed (to output y(n)), as described in Fig. 1. The latter task includes a translation to pseudo base-bandi.e. the resultant Power Spectral Density (PSD) is located at a frequency f ry which is close to zero -and a normalization step -as the one proposed in [START_REF] Gouldieff | Blind Digital Modulation Classification based on Mth-Power Nonlinear Transformation[END_REF] -so that the noise-free mixture has unit power. A low-pass pre-filter is also considered for a better rejection of adjacent channels and out-of-band noise.

The discrete time samples y(n) obtained at the output of the preprocessing unit are then equivalently given by

y(n) = U u=1 y u (n) + ν(n), (4)
where, y u (n) corresponds to the u th signal component:

y u (n) = a u •e i(2πfrun+φu) k∈Z s u (k)h u (nT e -kT u -τ u ), (5)
where T e denotes the sampling period and f ru the residual Carrier Frequency Offset (CFO) for signal u. Thereby, we classically consider that h u are Root Raised Cosine (RRC) filters with unknown Roll-Off factor β u .

C. Assumptions on the parameters in the realistic context

The satellite channels, in addition to the absence of synchronization at the gateways and satellites, inherently produce impairments between the downlink signals at the receiver.

1) Effects induced by the Local Oscillators (LO):

The LO have an impact on carrier frequency and symbol time. In fact, denoting u the summed LO imprecisions of link u, the received carrier frequency at the VSAT is f u = (1 + u ) • f c , with f c the target downlink carrier frequency. In the same way, we have T u = (1 + u ) • T , with T the target symbol rhythm.

Considering a standard value of around 0.1 ppm for u [START_REF] Comtech | CDM-625 Advanced Satellite Modem[END_REF], a bandwidth B of the order of the MHz and the use of C, Ku or Ka band lead to the following mathematical properties:

• For all u, T u is considered to all be almost equal T .

• For all u, the CFO f ru is small compared to T such that the "almost-overlapping" scenario holds. Without loss of generality (wlog), we assume in this work that the f ru are all different, even slightly, with probability one. A graphical representation of these effects is shown in Fig. 4 for a mixture of 3 signals. The CFOs have been voluntarily overstated for a better readability. The dashed line represent the shape of the low-pass pre-filter in the preprocessing unit. Last, P u stands for the power of signal u and P y for the summed power in band B, so that P y -P ν = 1. At the receiver, we assume that a proper algorithm, such as the one in [START_REF] Mazo | Jitter Comparison of Tones Generated by Squaring and by Fourth-Power Circuits[END_REF], estimates T . The residual CFOs are assumed as unknown: they are estimated by the proposed M PT method.

2) Other assumptions on signal y u : In (5), we assume that φ u and τ u are uniformly distributed, respectively in [0; 2π[ and [0; T [, and unknown. Note that unlike several AMC methods, we do not assume that T e = T (i.e. we do not perform any time synchronization at the receiver, which would be impractical in the interfering context). Instead, we assume that y(n) is properly oversampled by a factor ρ = T /T e which respects Shannon's sampling theorem, i.e. T -1 e ≥ B.

3) Assumptions related to s u and C: Symbol sets s u are independent and identically distributed, and respectively carried by constellations C u ∈ C. The C u may all be different. All the C u are zero mean, unit variance and unknown. Set C, however, is assumed to be known. As depicted in Fig. 5, we consider PSK, QAM and "hybrid" constellations for the sake of a sufficient diversity. Moreover, these constellations are typically used in SatComs (see e.g., IESS-310 and ETSI EN 301 210).

i -i 1 -1 i -i 1 -1 i -i 1 -1 i -i 1 -1 i -i 1 -1 i -i 1 -1 BPSK QPSK 8PSK 8AMPM R8QAM 16QAM
Fig. 5. Representation of the considered constellations in the complex plan

D. Overview of the assumptions

Table I sums up the assumptions proposed in the current Section as well as the blind aspects of the proposed method. Note that "Estimated with AM PT" means that we propose in this article an algorithm to blindly estimate the considered parameter. " * " means that no particular assumption is taken.

III. AMC AND MAMC: A PRACTICAL REVIEW

In this Section, we unfold a brief review of the general methods used in the literature for AMC and MAMC.

A. A brief review of AMC methods in the single-signal case

For the sake of a better understanding, it seems crucial to quickly draw a hierarchical inventory of classically used AMC algorithms in the single-signal case before outlining the formerly proposed MAMC methods.

Classically, the first step is to demodulate -e.g., normalize, filter and synchronize in time and/or frequency -the received signal. Then, a Likelihood-Based (LB) or a Feature-Based (FB) recognizer performs the constellation classification [START_REF] Dobre | Survey of Automatic Modulation Classification Techniques: Classical Approaches and New Trends[END_REF]. While LB algorithms are optimal in the likelihood sense -they use exact or approximated Likelihood Ratio Tests (LRT) as a decision criterion -FB methods are generally preferred for their robustness to model mismatch and their close-to-optimal performance at a lower computational complexity [START_REF] Hazza | An Overview of Feature-based Methods for Digital Modulation Classification[END_REF].

In the blind scenario, the performance of the whole AMC method relies on both the quality of the preprocessing stepwhich may be intricate in the blind context -and on the ability of the features to discriminate the constellations of a given set (Fig. 1). Several features have been proposed in the literature to classify both analog and digital modulations [START_REF] Hazza | An Overview of Feature-based Methods for Digital Modulation Classification[END_REF], such as Higher-Order Moments (HOM) and Higher-Order Cumulants (HOC) [START_REF] Swami | Hierarchical Digital Modulation Classification using Cumulants[END_REF]. However, these features are very sensitive to time and frequency synchronizations: in the blind scenario, any estimation inaccuracy would lead to a substantial loss in the performance of the whole classifier. Mainly to improve the robustness against the unknown in the blind scenario, Cyclic Moments (CM), Cyclic Cumulants (CC) [START_REF] Dobre | Cyclostationarity-Based Modulation Classification of Linear Digital Modulations in Flat Fading Channels[END_REF] and spectral nonlinearities of the signal [START_REF] Gouldieff | Blind Digital Modulation Classification based on Mth-Power Nonlinear Transformation[END_REF] were exploited as strong features for the blind classification of digital modulations.

B. Proposed MAMC methods in the co-channel context

Some of the features listed in Subsection III.A were also studied in the context of spectral overlapping communication signals in the single-node context (i.e. with a single antenna).

Spooner proposed a Minimum Distance (MD) classifier based on Cyclic Cumulants [START_REF] Spooner | Classification of Co-channel Communication Signals using Cyclic Cumulants[END_REF] [21] [START_REF] Spooner | The Cumulant Theory of Cyclostationary Time-Series, Part II: Development and Applications[END_REF] as a solution to the MAMC issue. Thereby, the features used for the classification are chosen by puncturing the two-dimensional n th -order Cyclic Temporal Cumulant Function (CTCF) at some specific frequencies and cyclic-frequencies for different values of n.

The method aims at being blind and only requires the knowledge of a few parameters. It mostly uses the signal bandwidth and the pulse-shaping filter characteristics as discriminants: thus, it may not be directly applicable in the present case study.

Zaerin et al. proposed a second approach to tackle the single-antenna MAMC issue [START_REF] Zaerin | Multiuser Modulation Classification based on Cumulants in AWGN channel[END_REF]. It relies on a MD classification based on some ratios of classical Higher Order Cumulants (HOCs). While this approach shows a very low complexity, its performance would quickly degrade in the blind scenario with realistic impairments: actually, it is assumed in [START_REF] Zaerin | Multiuser Modulation Classification based on Cumulants in AWGN channel[END_REF] a perfect knowledge of all the parameters involved in the mixture. Moreover, the method performs well only if the residual f ru , delays τ u and phases φ u are null which does not hold in practice. This method does not meet our constraints either.

C. AM PT: a natural choice for efficient MAMC

We finally justify the use of Analytical M th -Power nonlinear Transformation for MAMC in the blind co-channel context.

Note that the assumptions related to timing and frequency imply that we are operating in one of the most challenging MAMC case study. In fact, if the T u were different, CM/CC based methods would work successfully, as in [START_REF] Spooner | Classification of Co-channel Communication Signals using Cyclic Cumulants[END_REF]. In the case of partially overlapping signals -i.e. when the residual CFOs are of the order the symbol rhythm -filter-bank may be use to separate the interferers before the AMC step [START_REF] Zaerin | Multiuser Modulation Classification based on Cumulants in AWGN channel[END_REF]. However, to the best of our knowledge, the realistic scenario proposed in Section II has never been specifically carried out.

1) A strong relation with Cyclostationarity-based methods:

Cyclic Moments (CM) and Cyclic Cumulants (CC) seem the most promising in the blind multiuser context for their robustness to Carrier Frequency Offsets (CFO), Phase Offsets and Time Offsets [START_REF] Hazza | An Overview of Feature-based Methods for Digital Modulation Classification[END_REF] that inherently exist in the mixture. As a consequence, no intricate demodulation step is needed and the AMC is almost directly performed on the received signal.

The Cyclic Moment of y is especially defined as

R α y (τ ) n,q = ∞ -∞ n i=1 y ( * )i (t + τ i )e -i2παt dt (6)
where n is the order and ( * ) i represents the q conjugations. In the literature, n th -order/q-conjugated CMs with q = n/2 are classically used to generate features which are for instance robust to phase noise [START_REF] Dobre | Cyclostationarity-Based Robust Algorithms for QAM Signal Identification[END_REF]. With such a transform, a strong spectral line appears at each multiple of the cyclic-frequency.

In the presence of multiple co-channel signals, it is mandatory to choose q = n/2 so that the features will not overlap in the cyclic spectral space. As further defined in Section IV, the M PT of y may be seen as a derivative of the CM of y with n = M , q = 0 and τ = 0.

2) M PT-based AMC in the literature: The use of M PT for "classical" AMC was first proposed in [START_REF] Reichert | Automatic Classification of Communication Signals using Higher Order Statistics[END_REF] and recently applied in the Compressed Sensing context [START_REF] Lim | Automatic Modulation Recognition for Spectrum Sensing using Nonuniform Compressive Samples[END_REF]. An improved version was also studied in [START_REF] Gouldieff | Blind Digital Modulation Classification based on Mth-Power Nonlinear Transformation[END_REF] to distinguish between PSK, QAM and "hybrid" constellations. Thereby, it was proposed to base the recognition on the theoretical power of the generated spectral lines -through a method called AM PT -and not only on their cardinality. In fact, the fusion of these features for several M proved to provide stronger statistics for the classification. In the present work, we generalize the AM PTbased AMC to co-channel mixtures: the further theoretical developments lead us to a practical blind MAMC algorithm.

3) A moderate complexity: Since MAMC shows a practical interest in embedded Smart Modems, the complexity matters. M PT-based MAMC relies on the computation of a few PSDs, which complexity is of the order of the FFT. The complexity of the other processing blocs is mainly due to numerical methods on scalars (Section XIII) and they are in practice negligible.

Compared to the existing literature, the proposed method is slightly less complex than the Cyclic-Cumulants-based MAMC in [START_REF] Spooner | Classification of Co-channel Communication Signals using Cyclic Cumulants[END_REF]. However, AM PT-based method is more complex than [START_REF] Zaerin | Multiuser Modulation Classification based on Cumulants in AWGN channel[END_REF] (relying on Classical Cumulants which are low-complexity) but as seen the blind performance results (Section IX), the trade-off is clearly in favor of our method.

IV. BASICS ON M th -POWER NONLINEAR TRANSFORM In the scope of this article, we define the M th -Power Transform of y, denoted as M PT y (f ), as the M th -root of the Power Spectral Density (PSD) of y M (t). Then we have

M PT y (f ) F y M (t) (f ) 2 M , (7)
where F(•)(f ) stands for the classical Fourier Transform (FT). As we will further see, the M th -root is not necessarily required, but it nonetheless allows the power of the extracted features to be comparable when several values of M are used. Then, the M th -root is preferred for a better fusion of the features in the Minimum Distance (MD) sense. However, in the case of a hierarchical tree-based classification, the M th -root is not needed since the features are compared with precomputed thresholds instead of with one another.

Also remark that the discrete-time model has to be considered in practice: the M PT of y(n) -where n is the time index -is practically performed thanks to the classical Discrete-Time Fourier Transform (DTFT). Then, as exposed in [START_REF] Gouldieff | Blind Digital Modulation Classification based on Mth-Power Nonlinear Transformation[END_REF], the M PT is computed from discrete-time signal y(n) according to

M PT y (f ) = N e -2 M • Ne-1 n=0 y M (n) • e -2iπnf /Ne 2 M , (8)
where we recall that N e stands for the number of samples.

In the next Section, we derive the analytical expression of the power and position of the spectral lines occurring in the M PT function in the single-signal scenario. These results are then extended to the co-channel interference context -i.e. for U > 1 (almost)-totally overlapping signals -in Section VI.

V. SPECTRAL LINES IN THE SINGLE-SIGNAL CONTEXT

In this section, we derive the theoretical power of the spectral lines generated by the M PT functions [START_REF] Soysal | Channel Estimation and Adaptive M-QAM in Cognitive Radio Links[END_REF] with a single-signal (i.e. when U = 1). These powers -jointly taken with their respective spectral positions -represent the features that are further exploited to classify the incoming constellation.

A. Theoretical power of the spectral lines

We assume that y(t) is a continuous noise-free infinitelength signal, whose discretized version is given by (4). For the sake of coherence with Section VI, we keep index u in the following notation so that y u (t) = y(t).

The M PT function is not affected by residual phases and delays as proposed in [START_REF] Lim | Automatic Modulation Recognition for Spectrum Sensing using Nonuniform Compressive Samples[END_REF]. For conciseness, these parameters are both set to 0 without any loss of generality. It is obvious that residual CFO f ru shifts the spectral line by M •f ru , while amplitude a u multiplies its power by a 2 u . For ease of notation, we set f ru to 0 and a u to 1 in what directly follows.

First, we develop the general literal expression of continuous time signal y u (t) raised to power M (please refer to Appendix A for further justification of these results). Starting with

y M u (t) = k∈Z s u (k) • h u (t -k • T ) M , (9)
and following the work in [START_REF] Mazo | Jitter Comparison of Tones Generated by Squaring and by Fourth-Power Circuits[END_REF], y M u (t) is fully written as

y M u (t) = i∈I M M i y u,i (t), (10)
where, in [START_REF] Spooner | Classification of Co-channel Communication Signals using Cyclic Cumulants[END_REF], the following notation was adopted:

IM = i ∈ Z M , M j=1 ij = M, i k = 0 =⇒ ∀k ≥ k, i k = 0 , (11) M i M ! M j=1 (i j !) , (12)
y u,i (t) = k∈K i M M j=1 s ij u (k j ) • h ij u (t -k j • T ) , (13)
and where, in (13), we have

K i M = k ∈ Z M , kj =0 if ij =0 i<j =⇒ ki<kj else . (14)
In [START_REF] Spooner | Classification of Co-channel Communication Signals using Cyclic Cumulants[END_REF], some terms in the right-hand side sum are such that E [y u,i (t)] = 0. The presence of a spectral line at each frequency which is a multiple of the symbol rhythm is made obvious according to the Poisson Summation Formula, since

E [y u,i (t)] = 1 T k∈K 0i M µ i u,k m∈RZ e i•mt • H i u,k (m) , (15)
where we denoted

K 0i M = k ∈ K i M s.t. k 1 = 0 , (16)
µ i u,k = E   M j=1 s ij u (k j )   , (17)
H i u,k (f ) = F   M j=1 h ij u (t -k j • T )   . (18)
From ( 7), [START_REF] Spooner | Classification of Co-channel Communication Signals using Cyclic Cumulants[END_REF] and [START_REF] Comtech | CDM-625 Advanced Satellite Modem[END_REF], we deduce the analytical power of the spectral line generated by M PT yu around null frequency. We denote this power λ u,M , whose analytical expression is

λ u,M = 1 T i∈I M M i • k∈K 0i M µ i u,k • H i u,k (0) 2 M (19)
As we can see, several terms participate to the theoretical power of this spectral line. In fact, for a given M and for a transmitted constellation C u , some specific couples of vectors (i, k) are such that neither µ i u,k nor H i u,k (0) is null. Theoretical expression of µ i u,k can be obtained from the i th -order moment µ (i) u of the emitted constellation C u under the i.i.d. hypothesis for the source symbols. We actually have

µ i u,k = M j=1 E s ij u (k j ) = M j=1 µ (ij ) u = µ i u . (20)
Some values of µ

(i) u
are tabulated in Table II for the constellations considered thereby. For odd values of M , the moments are all null and they have been willingly removed.

Meanwhile, obtaining an analytical expression for H i u,k (0) is quite tedious for the pulse-shaping functions classically considered -such as the root-raised cosine filter used thereby. In practice, for a given radio standard, we would numerically compute H i u,k (0) and store it for different sets of parameters. As an example and in order to illustrate the behaviour of ( 19), we develop in the next two subsections the simplified expression of λ u,M for M = 2 and M = 4.

(i) u FOR Cu ∈ C Mod. → BPSK QPSK 8PSK 8AMPM R8QAM 16QAM µ (2) u 1 0 0 -1 5

B. Simplified expression of λ u,M for M = 2

As previously proposed, y 2 u (t) easily splits into two sums:

y 2 u (t) = y u,20 (t) + 2 • y u,11 (t), (21)
where

y u,20 (t) = k1 s 2 u (k 1 ) • h 2 u (t -k 1 T ) y u,11 (t) = k1<k2 s u (k 1 )s u (k 2 )•h u (t -k 1 T ) h u (t -k 2 T )
Term y u,20 (t) may show a spectral line if E s 2 u = 0, which is true for non-π/2-rotation-invariant constellations such as BPSK or R8QAM. However, with i.i.d. symbols, E[s u (k 1 )s u (k 2 )] = 0 whatever C u and no spectral line occurs.

If we consider that h u is a square-root raised cosine filter and whatever the roll-off β u , we eventually get H 20 u,0 (0) = T . Then, in this context, the analytical expression of λ u,2 is

λ u,2 = 1 T • µ (2) u • T + 0 • k>0 H 11 u,0k (0) = µ (2) u . (22)
C. Simplified expression of λ u,M for M = 4

Likewise, the theoretical expression of λ u,4 is drawn in the same context. Now, several terms in the development of y 4 u contribute to the power of the spectral line. Herein, we get Among the right-hand side terms in [START_REF] Dobre | Cyclostationarity-Based Robust Algorithms for QAM Signal Identification[END_REF], only y u,4000 and y u,2200 may show a non-null expectation and contribute to the power of the spectral line. Meanwhile, terms y u,3100 , y u,2110 and y u,1111 show null expectations: they have no effect on the power of the spectral line whatever the constellation.

The development of theoretical λ u,4 according to ( 19) is a bit more tedious than for M = 2 but follows the exact same principle. We eventually get after some simplifications

λ u,4 = 1 T µ (4) u H 4000 u,0 (0)+6 µ (2) u 2 k>0 H 2200 u,0k00 (0) 1 2
.

No simple analytical formulas were found neither for H 4000 u,0 (0) nor for H 2200 u,0k00 (0) in the situation where h u is a square-root raised cosine (SRRC) filter. Anyway, these values are numerically computable and they can be tabulated to further reduce the complexity of the classification algorithm.

VI. SPECTRAL LINES IN THE CO-CHANNEL CONTEXT

We now consider the co-channel interference scenario: to this end, we adopt the assumptions proposed in Section II. To develop the analytical expression of the power and position of the spectral lines, we first consider that the mixture is composed of U continuous-time and free of noise signals.

A. Preliminary considerations

By expanding the M th -power of the received mixture of U signals according to the multinomial theorem, we get

y M (t) = m∈M M U M m U u=1 y mu u (t), (25)
where m = (m 1 , • • •, m U ), where the multinomial coefficients are classically developed as in section V as

M m M ! U u=1 m u ! , (26)
and where M M U denotes the following set:

M M U = m = (m 1 , • • • , m U ) ∈ N U , U u=1 m u = M . (27)
B. Cardinality of the set of spectral lines for a given M Some right-hand side terms of the sum in ( 25) may produce a spectral line. In what follows, we derive the expression of the subsets of M M U prone to generate spectral lines, mainly depending on the set of constellations {C u }.

Let first remark that all zero-mean π-rotational symmetric constellations don't produce any spectral line if at least one of the m u is odd. Thus, for usual constellations such as the one considered in the scope of this article, the sum in ( 25) is performed on the following reduced set whatever

{C u } u∈{1,•••,U } : M M,2 U = m ∈ M M U ∩ (2N) U , (28)
which can be formulated using the notation of ( 27) as

M M,2 U = 2 • M M/2 U . ( 29
)
The number of terms in the sum in [START_REF] Lim | Automatic Modulation Recognition for Spectrum Sensing using Nonuniform Compressive Samples[END_REF] participating to the theoretical power of the spectral lines that may appear in the M PT of y is the cardinality of M M/2 U . This set is known, in combinatorics, as the set of all the M/2-multicombinations [START_REF] Ryser | Combinatorial Mathematics[END_REF]. The cardinality of this set for U ≥ 1 is expressed as

Card M M/2 U = U M/2 , (30)
where (( • )) -reads "U multichoose M/2" -is analytically given in terms of classical combinations and for U ≥ 1 as

U M/2 = U + M/2 -1 M/2 . ( 31
)
We finally get the number of expected spectral lines around the null frequency for a given {M,

U } ∈ N * 2 , Card M M,2 U = (U + M/2 -1)! (U -1)! • (M/2)! . ( 32
)
Some values of Card M M,2 U for small values of M and U are gathered in Table III: Note that it is still possible to refine the set given by ( 28) when considering the order of the rotational symmetry of the emitted constellations. Actually, as this order grows, some standard moments of the constellation become null and some spectral lines do no longer appear. Then, if we assume that C u is 2π/X u -rotational symmetric -with X u the highest possible -then the sum in ( 25) is performed on reduced set

M M,X U = m ∈ M M U ∩ XN . ( 33
)
The analytical formula of Card M M,X U whatever X, U and M is rather inconvenient and proves not to be crucial. Instead, we assume that the number of spectral lines is always given by (32) whatever the emitted constellations. However, some spectral lines may appear with a null power.

C. Analytical power of the spectral lines

Following the previous development, it is possible to derive the theoretical power of each term in the sum in [START_REF] Lim | Automatic Modulation Recognition for Spectrum Sensing using Nonuniform Compressive Samples[END_REF]. In fact, for a given vector m, we may get (see Appendix B.)

λ m = M ! U u=1 m u ! • U u=1 a 2 u • λ u,mu mu 2 2 M , (34)
where, as defined previously, λ u,mu stands for the theoretical power of the spectral line generated by y mu u if it was alone.

Formula (34) is of major interest since it shows that the theoretical powers of the spectral lines in the interfering case can simply be derived from the ones in the single-signal case.

D. Position of the spectral lines

Let us now consider the spectral shift entailed by the residual carrier frequencies f ru on the spectral lines. While replacing y u in (25) by its expression, literally given by ( 5), it is straightforward that the spectral line generated for a given m simply appears at frequency f m , such that

f m = U u=1 m u • f ru . ( 35
)
Similarly to the way we normalize the power of the spectral lines -as the M th -root of the PSD -we define the normalized spectrum scale f such that f = f M . Then, the normalized position of the spectral line indexed by vector m is

f m 1 M U u=1 m u • f ru . ( 36
)
E. Theoretical features for a mixture of two signals

According to the previous development, we gathered in Table IV the expected number, powers and positions of the spectral lines for a few values of M when U = 2.

TABLE IV PROPERTIES OF THE FEATURES FOR

A MIXTURE OF 2 SIGNALS M ↓ Peaks Norm. f m Theoretical Power λm 2 2 f r1 , f r2 a 2 1 λ 1,2 , a 2 2 λ 2,2 4 3 f r1 , f r1 +f r2 2 , f r2 a 2 1 λ 1,4 , a 1 a 2 6λ 1,2 λ 2,2 , a 2 2 λ 2,4 6 4 f r1 , 2f r1 +f r2 3 , f r1 +2f r2 3 , f r2 a 2 1 λ 1,6 , 3 15a 4 1 a 2 2 λ 2 1,4 λ 2,2 3 15a 2 1 a 4 2 λ 1,2 λ 2 2,4 , a 2 2 λ 2,6
In order to verify the proposed analytical development, we jointly represent, Fig. 6, the theoretical λ m computed according to Table IV and the M PT functions of normalized y(n) for a few M . Thereby, a mixture of a two signals was considered, respectively with constellations BPSK and R8QAM, and residual CFOs of -20 Hz and +40 Hz. The Signal-to-Interferer Ratio (SIR) of the strongest signal is 3.5 dB. The noise-free mixture was up-sampled by a factor ρ = 4 and N e was set to 4 • 10 5 samples (which is high enough to verify the "infinite-length" observation hypothesis). As expected, the theoretical position f m and power λ m of each spectral line in the M PT matches the simulated ones.

F. Towards the noisy discrete-time finite-length model

Up to now, we considered the theoretical references obtained for continuous-time, infinite-length, noise-free signals.

With noisy mixtures and small observation lengths, the power of the spectral lines is no longer deterministic and the spectral noise becomes quite substantial. In Section VII, we derive the Probability Density Functions (PDFs) for both the spectral lines and the background noise. This development aims at:

• Analytically derive the references in this realistic context.

• Get the theoretical probabilities of good detection of the spectral lines and false alarm rate in the blind scenario. • Compute the theoretical Correct Classification Rate (CCR) in various AMC and MAMC scenarios.

VII. PROBABILISTIC CONSIDERATIONS

In this Section, we consider the effects of a noisy discretetime finite-length model for the received signal, namely:

• The spectral background noise power impact on the power of the spectral lines. • The statistical behaviour of the power of both the background noise and the spectral lines.

A. Motivation and preliminaries

Mazo rightly pointed out the first effect in [START_REF] Mazo | Jitter Comparison of Tones Generated by Squaring and by Fourth-Power Circuits[END_REF]: performing the M PT computation on a finite-length signal inherently leads to the integration of the background noise in a bandwidth of size fe

Ne around f m . The mean background noise quickly increases with both M and U and decreases with the composite SNR, which motivates to study its impact on measured λ m . Mazo proposed a general framework to compute the mean background noise power, based on the expression of the PSD as the Fourier Transform of the autocorrelation function of y. While this method is easy to apply for M = 2 and U = 1 -which was done by Reichert in [START_REF] Reichert | Automatic Classification of Communication Signals using Higher Order Statistics[END_REF] for Gaussian filtersit becomes infeasible for higher values of M and U due to the exponentially growing number of cross terms appearing in the autocorrelation function of y. Besides, it does not give any information on the distributions of the spectral power, which is needed to unfold the analytical expression of both the good detection and false alarm rates in the blind scenario.

Instead, we propose to draw the PDFs of the spectral power by considering that the mean background noise power around the spectral lines in (M PT y ) M (f ), denoted by η M , is easily estimable. Following [START_REF] Reichert | Automatic Classification of Communication Signals using Higher Order Statistics[END_REF], we may consider the distributions under two different hypothesis, namely:

• H 0 -M PT y (f ) belongs to the background noise • H 1 -M PT y (f ) is a spectral line.

B. Distribution of the power of the background noise (H 0 )

It can be shown that the probabilistic behaviour of (M PT y ) M (f ) under hypothesis H 0 follows a centered chisquared distribution with two degrees of freedom -see e.g., [START_REF] Reichert | Automatic Classification of Communication Signals using Higher Order Statistics[END_REF] [27] [START_REF] Jenkins | Spectral Analysis and Its Applications[END_REF]. Then, the PDF for each frequency bin around the null frequency and under hypothesis H 0 is derived as

p M PT M y (f ) = x|H 0 = 1 η M • e -x η M , (37)
where we recall that η M stands for the mean background noise power surrounding the spectral lines. η M is given by

η M = M PT M y (f ) f ∈B M , (38)
where B M is such that for all f in B M , M PT y (f ) shows no spectral line and is almost identically distributed. This second assertion is true when f is close to the null frequency. Note that "< • >" classically stands for the mean value.

C. Distribution of the power of the spectral lines (H 1 )

Deriving the PDF for the spectral lines is also needed to get the analytical expression of the references. It also allows to get the theoretical probability of good detection of the spectral lines and the Correct Classification Rate (CCR) of the (M)AMC method with a Minimum Distance classifier.

The power of (M PT y )

M/2 (f m ) proves to follow a Rice distribution [START_REF] Reichert | Automatic Classification of Communication Signals using Higher Order Statistics[END_REF]. Equivalently, (M PT y ) M (f m ) follows a noncentral Chi-squared distribution with 2 degrees of freedom:

p M PT M y (f m ) = x = 1 η M •e - x+λ M m η M •I 0 2 xλ M m η M , (39)
where I 0 (•) stands for the modified Bessel function of the first kind with order zero. From (39), we get (see Appendix C.)

p [M PTy(fm) = x] = M x M -1 ηM •e - x M +λ M m η M •I0 2 x M λ M m ηM . (40
λ ref m = M √ η M •e - λ M m η M •Γ 1+ 1 M • 1 F 1 1+ 1 M , 1, λ M m η M , (41)
where Γ stands for Euler's Gamma function and p F q (a, b, x) the p, q-order generalized hypergeometric function.

We may remark that if λ M m is large compared to η M , then the following approximation holds:

1 F 1 1+ 1 M , 1, λ M m η M e λ M m η M Γ 1 + 1 M • λ M m η M 1 M . (42)
Then, λ ref m λ m . In other words, the PDF given by ( 40) is almost centred around λ m . In that case, it is not mandatory to use (41) to refine the references and using the expression given by (34) gives sufficient precision for the references.

E. Blind detection of the Spectral Lines

In the blind context, the first step is to detect the spectral lines in the computed (M PT y ) M functions. The probability of detecting a false spectral line at a given frequency f with threshold ξ M is derived according to (37) as

P F,s,M = +∞ ξ M 1 η M • e -x η M dx = e - ξ M η M . (43)
The detection of the spectral lines can be handled on a reduced part of the spectrum: this set is centered around 0 and its bounds mainly depend on the maximum expected value of the CFOs. If we assume that the maximum normalized CFO for all the users is ±f max , then the peak detection is performed on n b = 2M N e f max frequency bins. Moreover, the power of the background noise bins proves to be identically distributed in

[-f max ; f max ]\{f m } if f max << T -1
. Then, we get the following false alarm rate in the considered bandwidth:

P F,M 1 -1 -e - ξ M η M n b . (44)
Fig. 7 shows the theoretical probability of false detection on the whole spectrum-of-interest (44) and the simulated one for M = {2, 4}. Here, n b was set to 50M , U to 2, ρ to 4, N e to 4 • 10 4 and SN R to 10 dB. These probabilities are represented as a function of the normalized threshold ξ M /η M . Fig. 7. Theoretical and simulated false alarm rates P F,2 and P F,4

From (44), we get the expression of ξ M for a given P F,M :

ξ M -η M • ln 1 -(1 -P F,M ) 1 n b . (45)
Generally -mainly to get the least possible outliers and have a simple post-processing algorithm -we choose ξ M such that P F,M is as small as possible while correctly detecting most of the spectral lines. For a given ξ M , the probability of good detection of the spectral line occurring at frequency f m is derived using the Cumulative Distribution Function (CDF) of the non-central Chi-squared distribution. We may get:

P D,m = Q 1   2λ M m η M , 2ξ M η M   (46)
where Q 1 (•, •) stands for the first-order Marcum Q-function.

VIII. BLIND MAMC BASED ON AM PT

Throughout this section, we propose to blindly solve the inverse problem of estimating parameters U and {f ru , a u , C u } u∈{1,•••,U } , given functions M PT y for some M.

A. General scheme of the proposed blind MAMC method

We first draw the block diagram of the whole method with the references to specific equations or algorithms in Fig. 8.

B. Preprocessing Unit

As previously discussed, the preprocessing unit transforms the incoming signal x(n) into a normalized, pseudo-centred and pre-filtered signal y(n). For instance, this step can be performed according to the method proposed in [START_REF] Gouldieff | Blind Digital Modulation Classification based on Mth-Power Nonlinear Transformation[END_REF].

The preprocessing unit also estimates some parameters of the mixture if needed, such as the oversampling factor ρ.

C. Detection of the spectral lines in the M PTs

The decision starts with the detection of the spectral lines in the computed M PTs. As stated in Section VII, proper threshold ξ M has to be chosen to find a compromise between the good detection rate and the false-alarm rate. We proposewith Algorithm 1. -to detect the spectral lines with a chosen false-alarm rate P F,M . The spectral lines are then sorted by decreasing Peak-to-Background Noise Ratio (PBNR). Thereby, we develop the method to compute the features λ M for some M . This task is not an issue when the number of users U and the residual carrier frequencies f ru are known. In this situation, we simply compute the powers given by [START_REF] Palicot | A New Concept for Wireless Reconfigurable Receivers[END_REF] for each f m (35) such that m is in M M,2 U (28). In the blind scenario, U and the residual CFOs f ru are not known. Then, we have to estimate them to correctly run the computation of the features.

To proceed, we propose an algorithm which iteratively builds a basis { f ru } and estimate U based on the set of detected spectral lines. We first work on the spectral lines that show high PBNRs to make the construction robust to outliers and to forgotten spectral line.

On overview of the method is proposed in Algorithm 2. We denote by n SL the number of spectral lines detected with Algorithm 1. The n th spectral line is represented by its frequency f n and the order M n of the M PT. F n represents the frequency basis obtained with the n first spectral lines.

U ← U + 1 5: F n-1 ← F n-1 ∪ f n 6:
Find F n by orthogonalizing F n-1

7:

end if 8: end for At the end of the algorithm, we get the estimated number of users U and a basis F n SL -also denoted as ( f 1 , • • • , f U )that describes the positions of the detected spectral lines. The computation of λ m is straightforward using [START_REF] Palicot | A New Concept for Wireless Reconfigurable Receivers[END_REF] and (35).

E. Computation of the references λ ref m

Analytical expressions of the references may be obtained according to [START_REF] Swami | Hierarchical Digital Modulation Classification using Cumulants[END_REF] and (34), and refined according to (41) if needed. Note that for small values of M and U and for medium to high SNRs, computing the references using (34) gives sufficient precision and greatly simplifies the algorithm in charge of estimating {a u } and {C u }.

F. Finding {a u } and {C u } by solving an optimization problem

In the blind scenario, sets {a u } and {C u }, for all u in {1, • • •, U } are unknown. Our goal is to estimate these parameters in the Minimum Distance (MD) sense.

1) Problem statement:

We assume that the λ m have been correctly computed for all m ∈ M M,2 U . Then, in the MD sense, we look for set { a u , C u } such that

{ a u , C u } = arg min {au,Cu} s.t. a 2 1 +•••+a 2 U =1 (C1,••• ,C U )∈C U M m∈M M,2 U λ m -λ ref m 2 , (47)
where we recall that λ ref m stands for the references, computed according to (34) and refined with (41) if assumed necessary.

2) Splitting the optimization problem: Solving (47) is performed by splitting the optimization problem: for each possible set {C u } ∈ C U , we first find { a u } such that

{ a u } = arg min {au},||a|| 2 2 =1 M m∈M M U ( λ m -λ ref m ) 2 (48)
Then, the set of constellations { C u } is the one which minimizes (47) according to the respective estimated sets of amplitudes { a u }.

G. Specific case study: U = 2 with the 2/4PT classifier Among the co-channel scenarios, the situation where two signals coexist is encountered in several contexts such as in dual-satellite systems. This scenario was studied by Zaerin et al. in [START_REF] Zaerin | Multiuser Modulation Classification based on Cumulants in AWGN channel[END_REF] with simplifying assumptions. In Section IX, we mainly study the performance of our method in specific context and we confront it with Zaerin's method.

Then, with two co-channel signals inferred by a 2PT/4PT classifier and under the approximation λ ref m λ m , the following multivariate polynomial P has to be minimized

P(a 1 , a 2 ) = λ 20 -a 2 1 λ 1,2 2 + λ 02 -a 2 2 λ 2,2 2 + λ 40 -a 2 1 λ 1,4 2 + λ 22 -a 1 a 2 6λ 1,2 λ 2,2 2 + λ 04 -a 2 2 λ 2,4 2 (49)
This problem is easily solved by neglecting the effects of the cross term -i.e. the fourth right hand side term in (49) -on the solution { a 1 ; a 2 } which minimizes (49). This last approximation leads to an analytical formula for { a 1 ; a 2 }.

IX. SIMULATION RESULTS

In this penultimate Section, we show the effectiveness of the proposed method throughout numerical simulations in various scenarios. To this end, this Section is split into three main subsections. The first part deals with the good detection rate of the number of signals and their respective CFOs, which is a crucial step in the blind scenario. Then, we expose some results in the dual-signals blind context through the Correct Classification Rate obtained with the proposed 2PT-4PT classifier. Last, the performance of the proposed method are compared with the one obtained with Zaerin's method [START_REF] Zaerin | Multiuser Modulation Classification based on Cumulants in AWGN channel[END_REF].

Otherwise stated, we consider that all the assumptions proposed in Section II hold. All the mixed sources are pulseshaped by Square-Root Raised Cosine (SRRC) filters with roll-off 0.25. We consider that the adjacent channels and the Gaussian noise that lay outside band B are perfectly rejected by a proper bandpass filter. The oversampling factor ρ is arbitrarily set to 4. The baud-rate R is set to 1e 6 symbols per second for all the sources. The residual CFOs are randomly and uniformly chosen in band [-1; +1] kHz. Moreover, we use the classical definitions of the composite SNR: SNR dB = 10 • log 10 ( U u=1 Pu Pν

), where P u stands for the power of signal u and P ν for the power of the noise.

A. Detection of the number of signals and the residual CFOs

We first consider the problem of detecting the number of signals U in the considered bandwidth as well as the respective residual CFOs of these signals. The good detection of these parameters is in fact mandatory for the classification algorithm to properly work. We have to mention that when U is correctly estimated, then the estimated residual CFOs are also correct.

Fig. 9 shows the Correct Detection Rate (CDR) of the number of users versus the composite SNR when the users transmit BPSK constellations with same powers. To this end, we use both 2PT and 4PT and we jointly run Algorithm 1 for the detection of the spectral lines, and Algorithm 2 for the estimation of U and {f ru }. Thereby, 2e 3 symbols were used for the detection, and the ratio ξ M /η M was set to 10 so that the overall false alarm rate -i.e. in both 2PT/4PT -is 0.5%. First, Fig. 10 shows the Correct Classification Rate (CCR) arbitrarily for signal 1 when the Signal-to-Interferer Ratio (SIR) of both signals is equal to 0 dB. This scenario can be seen as a "worst case" in MAMC, as stated in [START_REF] Zaerin | Multiuser Modulation Classification based on Cumulants in AWGN channel[END_REF].

Composite SNR (dB)

Composite SNR (dB) In Table V, we also expose the confusion matrix of the classifier in the same context as in Fig. 10 with SNR = 10 dB. Note that "0%" was replaced by " • " for a better readability. As expected, the worst CCR is for both QPSK and 16QAM constellations since they are the closest in the 2PT/4PT feature-space. Meanwhile, the CCR for other types of constellation is satisfying even at low SNRs. Last, Fig. 11 shows the mean CCR for both signals in the average realistic scenario, i.e. in the blind context and with a SIR randomly chosen in set [0; 6] dB (signal 1 is assumed to be the strongest wlog). Other parameters are randomly chosen as previously stated. The mean performance is slightly better for BPSK, 8AMPM and R8QAM constellations, while an expected degradation of the CCR for QPSK and 16QAM is encountered. If considering a channel coded transmission with usual code rate, the CCR obtained at the minimum SNR allowing proper multiuser demodulation is truly acceptable. C. Comparison with Zaerin's MAMC method [START_REF] Zaerin | Multiuser Modulation Classification based on Cumulants in AWGN channel[END_REF] In this penultimate part, we propose to compare the MAMC method proposed in [START_REF] Zaerin | Multiuser Modulation Classification based on Cumulants in AWGN channel[END_REF] with our 2PT/4PT-based classifier.

1) Preamble: Note that the method in [START_REF] Zaerin | Multiuser Modulation Classification based on Cumulants in AWGN channel[END_REF] is the closest to the one we propose thereby. However, it is not designed for the blind situation and several strong assumptions are taken in the co-channel case. Especially, it is assumed that the signal is perfectly sampled at the symbol rhythm, which doesn't hold in practice. Moreover, the method is designed for rectangular pulse-shaping, for null CFO and with a perfect synchronisation between the interfering signals. Any offset, no matter how small, entails consequent degradations in the CCR. Instead, the method we propose is designed to handle these offsets and no degradation of the CCR is noticed.

To remain fair, we propose to compare both methods in the non-blind situation with no phase offset and no delay between the signals. However, we assume non-null residual CFOs when simulating both methods. We assume perfect synchronization for [START_REF] Zaerin | Multiuser Modulation Classification based on Cumulants in AWGN channel[END_REF] while the over-sampling factor is set to 4 for the proposed AM PT-based method. Rectangular pulse-shaping is assumed for [START_REF] Zaerin | Multiuser Modulation Classification based on Cumulants in AWGN channel[END_REF] while we use SRRC filters with roll-off 0.25 for the proposed AM PT-based method.

2) Simulation results: Fig. 12 shows the Correct Classification Rate (CCR) of both methods for a mixture of two signals when C = {BPSK, QPSK, 16QAM} and N = 8000.

Composite SNR (dB) The proposed AM PT-based method clearly outperforms the method in [START_REF] Zaerin | Multiuser Modulation Classification based on Cumulants in AWGN channel[END_REF], even in the context of perfectly synchronized signals. Moreover, note that the performance of the MAMC method in [START_REF] Zaerin | Multiuser Modulation Classification based on Cumulants in AWGN channel[END_REF] would significantly drop with delays between the users or when the amplitudes of the signals are different. This method shows two main advantages. First, it makes no assumption regarding time and frequency synchronization at the receiver. Thus, no intricate demodulation scheme is needed and the constellation recognition is performed in a blind manner almost directly on the received signal. Secondly, its complexity is of the order of the Fourier Transform, which makes the whole algorithm feasible in cognitive receivers.

We showed that AM PT-based MD-sense classification was able de distinguish between several constellations even if they are very close to each other in the feature sense. The theoretical mathematical background was proposed and it can be easily extended to other type of signals as long as they show non-linearities when raised to some power M .

Further work should be carried out regarding the use of other non-linear transformations that may be better adapted to the set of constellations we want to classify. Also, the study of advanced MAMC methods based on AM PT should be regarded, especially in the situation where the physical channel is not a simple impaired AWGN channel. This scenario would be relevant for cognitive satellite terrestrial radio (CSTR), land mobile SatComs or for terrestrial cellular networks.

Fig. 1 .

 1 Fig. 1. The common Feature-Based AMC architecture

Fig. 2 .Fig. 3 .

 23 Fig. 2. Forward link interference at the VSAT due to overlapping footprints

2 𝐵Fig. 4 .

 24 Fig. 4. Representation of the PSDs of the co-channel signals

y 4 u

 4 = yu,4000 + 16yu,3100 + 6yu,2200 + 36yu,2110 + 24yu,1111.[START_REF] Dobre | Cyclostationarity-Based Robust Algorithms for QAM Signal Identification[END_REF]

Theoretical 6 m, 2 Theoretical 6 m, 4 Theoretical 6 m, 6 Fig. 6 .

 2466 Fig. 6. Theoretical power λm and simulated M PTy(f )

  ) D. References in the discrete-time finite-length context Based on the above development, the references λ ref m used to compare the constellations in the Minimum Distance (MD) sense can be analytically developed has the expected value of the distribution of the power of the spectral line (40), yields

Fig. 8 .

 8 Fig. 8. The proposed AM PT-based blind MAMC algorithm

Algorithm 2

 2 Estimation of U and { f ru } by decreasing PBNR 1: Initialize F 0 = ∅ and U = 0 2: for n = 1 : n SL do 3: if f n is not a linear combination of F n-1 (35) then 4:

Fig. 9 .

 9 Fig. 9. Correct Detection Rate of the number U of overlapping BPSK signals as a function of the composite SNR, for U ∈ {0, 1, 2, 3}

Fig. 10 .

 10 Fig. 10. Correct Classification Rate (CCR) for user 1 in the context of two overlapping signals -SIR = 0 dB -Blind scenario -N = 1e 4

Fig. 11 .

 11 Fig. 11. Correct Classification Rate (CCR) for both signals in the context of two overlapping signals -SIR ∈ [0; 6] dB -Blind scenario -N = 1e 4

Fig. 12 .

 12 Fig.12. Mean Correct Classification Rate for AM PT-based method and[START_REF] Zaerin | Multiuser Modulation Classification based on Cumulants in AWGN channel[END_REF]

  X. CONCLUSION We have proposed a new method for the blind modulation classification of co-channel overlapping signals in the cognitive satellite context. It is based on an Analytical study of the M th -Power nonlinear Transformation (AM PT) of the mixture.

TABLE II THEORETICAL

 II i th -ORDER MOMENTS µ

TABLE V CONFUSION

 V MATRIX WITH SN R = 10dB FOR THE PROPOSED 5-CLASS PROBLEM -SIR = 0 DB -BLIND SCENARIO -N = 1E 4

				Emitted Constellation
			BPSK QPSK 8AMPM R8QAM 16QAM
		BPSK	100%	.	.	.	.
	Classifier Output	QPSK 8AMPM R8QAM 16QAM	. . . .	92.8% 1.8% . 5.4%	0.6% 97.1% . 2.3%	. . 100% .	17.4% 4.5% . 78.1%

By "classical", we include all the AMC methods that are designed for a single signal corrupted by noise. Thereby, it is assumed that no interferer is spectrally overlapping with this signal whatsoever.

The terminology "Multiuser AMC" (MAMC) was first proposed in 2012 byZaerin et al.,[START_REF] Zaerin | Multiuser Modulation Classification based on Cumulants in AWGN channel[END_REF]. MAMC refers to the AMC methods that are designed to face partial or total spectral overlapping of two or more signals.

"Single-node" and "multi-node" refer to the number of terrestrial stations that are assumed to receive and process the mixture. In the satellite context, a "node" is composed of a single antenna since no spatial diversity is available.

APPENDIX A JUSTIFICATION OF SECTION V. Denote x u (k) = s u (k) • h u (t -k • T ) for better reading. We first apply the multinomial theorem to expand the power of the (infinite) sum as an infinite sum of infinite products.

In the sum, only M terms i j may be non-zero at the same time since the i j are all positive integers. Then, we denote:

• M the cardinality of the non-zero terms (0

Thus, we can rewrite (50) as:

It is then possible to invert the sums by re-indexing terms i independently from the k j . This operation leads to:

Following the fact that 0! = 1 and x 0 u = 1, the set in the first sum can be performed in the condensed generic set I M [START_REF] Gouldieff | Blind Digital Modulation Classification based on Mth-Power Nonlinear Transformation[END_REF] and the second sum in generic set K i M [START_REF] Liolis | Cognitive Radio Scenarios for Satellite Communications: the CoRaSat Approach[END_REF]. M can simply be replaced by M in both the right-hand side products.

Then, we take the expectation of y u,i (t). Remarking that that h ij u (t -k j T ) is deterministic for a given {i, k, t} yields:

The Poisson Summation Formula (PSF) is given by

where G(f ) stands for the Fourier Transform of g(t).

Remarking that (53) can be written as

where

makes it possible to apply the PSF to g(t -kT ), yields [START_REF] Comtech | CDM-625 Advanced Satellite Modem[END_REF].

APPENDIX B JUSTIFICATION OF SECTION VI.

Lets unfold the expression of λ m from ( 7) and [START_REF] Lim | Automatic Modulation Recognition for Spectrum Sensing using Nonuniform Compressive Samples[END_REF]:

where * stands for the classical convolution -and considering the properties of the Fourier Transform, (57) may be written as:

At frequency f m , the terms in the sum in (58) are all null if m = m and if the frequencies

Then, the convolution is a simple product at frequency f m . Yields:

Then, (34) is simply obtained by distributing the absolute value in the product and by remarking that

APPENDIX C JUSTIFICATION OF SECTION VII.

The Probability Density Functions (PDF) of the PSD of y M under H 0 and H 1 were developed in the literature related to spectral analysis [START_REF] Jenkins | Spectral Analysis and Its Applications[END_REF]. First, (40) can be obtained from (39) with a simple variable change. Then, (41) is derived as:

Replacing I 0 by its expression and inverting both sum and integral operators yields:

Replacing the integral by its expression may give:

Denoting by • (k) the rising factorial and noting that Γ 1+k+

and that k! = 1 (k) , we finally get expression (41).


Olga Kupriianova
email: olga.kupriianova@lip6.fr

Christoph Lauter
email: christoph.lauter@lip6.fr

Metalibm: a Mathematical Functions Code Generator

Keywords: code generation, elementary functions, mathematical libraries

There are several different libraries with code for mathematical functions such as exp, log, sin, cos, etc. They provide only one implementation for each function. As there is a link between accuracy and performance, that approach is not optimal. Sometimes there is a need to rewrite a function's implementation with the respect to a particular specification.

In this paper we present a code generator for parametrized implementations of mathematical functions. We discuss the benefits of code generation for mathematical libraries and present how to implement mathematical functions. We also explain how the mathematical functions are usually implemented and generalize this idea for the case of arbitrary function with implementation parameters. Our code generator produces C code for parametrized functions within a known scheme: range reduction (domain splitting), polynomial approximation and reconstruction. This approach can be expanded to generate code for black-box functions, e.g. defined only by differential equations.

Introduction

Each time we evaluate mathematical functions in some programming language a corresponding function from a mathematical library (libm) is called. There are several examples of existing libms: glibc libm [START_REF] Loosemore | The GNU C Library Reference Manual for version 2.19[END_REF], crlibm by ENS-Lyon [START_REF] Daramy-Loirat | A library of correctly rounded elementary functions in double-precision[END_REF], libmcr by Sun1 , libultim by IBM2 , etc. They differ not only by the developer company or supported language but also by final accuracy. The common fact for all the versions is that they provide only one manually coded implementation of each supported function and precision. As the codes for elementary functions are used in various applications, the implementations are done for the widest possible domain, for maximum possible accuracy, etc. However, such generalization complicates the algorithms and leads to poor performance: for some particular tasks there is no need to compute a precise result. For example, most physical measurements have only small number of digits after the decimal point, so computation of the result with 53 mantissa bits (about 15 decimal digits) is a waste of time. Another example is a small domain known beforehand, which means that that the result is a finite number (no overflows/underflows are possible). In this case, handling the special values like NaNs (Not-A-Number [START_REF]IEEE Computer Society, IEEE Standard for Floating-Point Arithmetic[END_REF]) or infinities in the beginning of the implementation can be skipped, hence the implementation gets faster.

There is also a link between speed and accuracy, so when we process big amounts of data and the needed accuracy is only about several bits, there is no reason to compute a precise result [START_REF] Muller | Elementary Functions: Algorithms and Implementation[END_REF].

So, there is a growing need to provide several different versions (flavors) for each libms function. Manual implementation of different function flavors is almost impossible due to the quantity of all the possible parameters and coding time. Thus, we propose to write a code generator that produces parametrized implementations. The name of this prototype is Metalibm3 . Besides producing different versions of standard libm functions, Metalibm generates implementations for composite functions as well.

The paper is organized as following: in Section 2 we explain in general the generation of "black-box" parametrized functions, Section 3 explains the workflow of the generator, in Section 4 we show how to detect a known type of function, and finally Section 5 shows the importance of the Metalibm, future work on the project and its possible application.

A Black-box Function Generator

The Metalibm code is a collection of scripts written in Sollya4 , a software tool for safe floating-point (FP) development with plenty of rigorous numerical algorithms [START_REF] Chevillard | Sollya: An Environment for the Development of Numerical Codes[END_REF].

Among the existing versions of mathematical libraries we are mostly interested in improvement of the current glibc libm. It runs on all *nix-powered machines (from supercomputers to mobile phones), so tends to be the most used one. Thus, Metalibm generates implementations in C. For each input set of parameters Metalibm generates the C code within the same scheme. On Fig. 1 there is an illustration of Metalibm routine. Metalibm computes and stores the constants with the needed accuracy, and as we use table-driven methods [START_REF] Muller | Elementary Functions: Algorithms and Implementation[END_REF], it computes and stores the tables. Then Metalibm generates the code for range reduction (domain splitting), polynomial approximation on each of the small domains, and the final reconstruction procedure. The purpose of these procedures is explained later in the paper.

As all the computations for the code generation are done in Sollya, a function to be implemented can be considered as black-box. On the step of range reduction we need to evaluate the function in some values. Sollya provides elementary functions and theirs combinations. The generation of code for some "exotic functions" like functions purely defined by differential equations (e.g. Dickman's) gets possible as soon as we have a corresponding Sollya implementation, even if it is bound to Sollya only dynamically [START_REF] Chevillard | Users manual for the Sollya tool[END_REF].

Precision of all the constants and accuracy of the interim computations are specially selected in order to obtain the final result of the specified accuracy (target accuracy parameter). Besides the generated code, Metalibm verifies the final results accuracy and generates a Gappa proof [START_REF] De Dinechin | Certifying the Floating-Point Implementation of an Elementary Function Using Gappa[END_REF].

A simple example of the parametrization file is provided in Listing 1.1. The sense of all the parameters will be explained later with the implementation details of Metalibm.

Function Generation Workflow

In order to implement a mathematical function on a given domain we have to care first about special cases (NaNs and infinities). So, the first or even precomputing step is always filtering the special inputs, handling exceptions, too large or too small inputs unless the domain is so small that these special cases cannot occur. When designing algorithms for mathematical functions evaluation we usually start with writing the inputs and output in a form of floating-point (FP) numbers 2 E m and trying to separate all the factors into two groups: a power of two n = 2 E that represents the results exponential part and one non-integer number with values in a small range to represent significand. This technique of emphasizing a non-integer part with a small range (results significand) is usually called range reduction [START_REF] Muller | Handbook of Floating-Point Arithmetic[END_REF]. The next step is building an approximation. There are different approximation techniques, but we only consider polynomial ones. The larger the range of argument is, the higher polynomial degrees are required. This means that the computation time will be high and the error analysis becomes more difficult with the growing number of FP operations. This is the main reason why we need argument reduction. In Metalibm, we use a parameter maxDegree to limit the maximal degree of constructed approximations. These polynomials are computed with Remez algorithm that is implemented in Sollya [START_REF] Chevillard | Users manual for the Sollya tool[END_REF], [START_REF] Chevillard | Évaluation efficace de fonctions numériques[END_REF].

All the transformations on the first step fully depend on mathematical properties of the implemented function, i.e. e a+b = e a e b [START_REF] Tang | Table-driven implementation of the exponential function in IEEE floating-point arithmetic[END_REF]. It can happen that it is impossible to reduce the range using only mathematical properties. In this situation the required domain for the implementation is divided into subdomains. On each of the subdomains the argument range is small and the degree of the approximating polynomial gets lower. In this case, one has to build approximations for each of the subdomains and then the reconstruction step gets more complex. In Metalibm we bound minimal size of subdomain by a parameter minWidth.

Once the generation process is launched, Metalibm checks whether it is possible to build a polynomial for the specified function, domain and other parameters. If it is not possible, it tries to reduce argument for the set of known functions (exponential, logarithm, periodic). Then there is symmetry detection and expression decomposition for composite functions. When it is still not possible to build a polynomial approximation for the reduced domain, domain splitting is performed. Metalibm adapts computational precision [START_REF] Ch | Arrondi correct de fonctions mathématiques[END_REF] in order to get the needed accuracy, so in some cases it uses double-double or triple-double arithmetic [START_REF] Schewchuk | Adaptive Precision Floating-Point Arithmetic and Fast Robust Geometric Predicates[END_REF], [START_REF] De Dinechin | Fast correct rounding of elementary functions in double precision using double-extended arithmetic[END_REF].

Some Technical Details

As it was already mentioned, Metalibm tries to detect some known properties of the function to perform range reduction. For example, let us have a look on how Metalibm tests hypothesis that the function is exponential and finds its base.

In order to detect the exponential function and to perform the appropriate argument reduction we accept the hypothesis that the function has a form of f (x) = b x on the implementation domain I with the unknown base b. It means that we can determine the base from the following:

b = exp ln f (x) x = const ∀x ∈ I.
The base b can be computed in Sollya without any information about the function f (x). Sollya will provide the needed value for the function. This is why we were talking about "black-box" functions: we do not know what code are we generating, but we have a mean to evaluate some function values.

As we have accepted the function type as f (x) = b x , it means that for another argument x 1 from the implementation domain the function value is f (x 1 ) = b x1 . Of course the base b stays the same if the hypothesis was correct and we know that we can perform the exponential argument reduction now. If

the value ε = b x f (x) -1 I ∞
is sufficiently small we accept the initial hypothesis. The detection of other types of functions is done in the analogous way.

Conclusions

The Metalibm code generator is still under development, but it already produces code for basic functions. The next goals are an optimized domain splitting procedure, producing vectorizable implementations and addition of range reduction for other functions.

As it was told, the generated code tends to improve and expand current glibc libm. So, we will try to provide generated code on the whole code generator to the GNU community for integration with the glibc libm.

Fig. 1 .

 1 Fig. 1. Metalibm scheme

  (x); // we want to get the code for exp(x) 2 dom = [-70, 70]; // on domain [-70, 70] 3 target = 2 -42 ; // the final

http://www.math.utah.edu/cgi-bin/man2html.cgi?/usr/local/man/man3/libmcr.3

http://www.math.utah.edu/cgi-bin/man2html.cgi?/usr/local/man/man3/libultim.3

http://lipforge.ens-lyon.fr/www/metalibm/

http://sollya.gforge.inria.fr/