K Hassan

I Dayoub

W Hamouda
email: hamouda@ece.concordia.ca

M Berbineau

Automatic Modulation Recognition Using Wavelet Transform and Neural Networks in Wireless Systems

Modulation type is one of the most important characteristics used in signal waveform identification. In this paper, an algorithm for automatic digital modulation recognition is proposed. The proposed algorithm is verified using higher-order statistical moments (HOM) of continuous wavelet transform (CWT) as a features set. A multilayer feed-forward neural network trained with resilient backpropagation learning algorithm is proposed as a classifier. The purpose is to discriminate among different M-ary shift keying modulation schemes and the modulation order without any priori signal information. Pre-processing and features subset selection using principal component analysis is used to reduce the network complexity and to improve the classifier's performance. The proposed algorithm is evaluated through confusion matrix and false recognition probability. The proposed classifier is shown to be capable of recognizing the modulation scheme with high accuracy over wide signal-to-noise ratio (SNR) range over both additive white Gaussian noise (AWGN) and different fading channels.

Introduction

Blind signal interception applications have a great importance in the domain of wireless communications. Developing more effective automatic digital modulation recognition (ADMR) algorithms is an essential step in the interception process. These algorithms yield to an automatic classifier of the different waveforms and modulation schemes used in telecommunication systems (2G/3G and 4G).

In particular, ADMR has gained a great attention in military applications, such as communication intelligence (COMINT), electronic support measures (ESM), spectrum surveillance, threat evaluation, and interference identification. Also recent and rapid developments in software-defined radio (SDR) have given ADMR more importance in civil applications, since the flexibility of SDR is based on perfect recognition of the modulation scheme of the desired signal.

Modulation classifiers are generally divided into two categories. The first category is based on decision-theoretic approach while the second on pattern recognition [START_REF] Dobre | Survey of automatic modulation classification techniques: classical approaches and new trends[END_REF]. The decision-theoretic approach is a probabilistic solution based on a priori knowledge of probability functions and certain hypotheses [START_REF] Wei | Maximum-likelihood classification for digital amplitude-phase modulations[END_REF][START_REF] Dobre | Likelihood-based algorithms for linear digital modulation classification in fading CHAN-NELS[END_REF]. On the other hand, the pattern recognition approach is based on extracting some basic characteristics of the signal called features [START_REF] Ebrahimzadeh | Intelligent digital signaltype identification[END_REF][START_REF] Zhao | Modulation recognition of SDR receivers based on WNN[END_REF][START_REF] Wong | Automatic digital modulation recognition using artificial neural network and genetic algorithm[END_REF][START_REF] Nandi | Algorithms for automatic modulation recognition of communication signals[END_REF][START_REF] Swami | Hierarchical digital modulation classification using cumulants[END_REF][START_REF] Dobre | Robust QAM modulation classification algorithm using cyclic cumulants[END_REF][START_REF] Mobasseri | Digital modulation classification using constellation shape[END_REF][START_REF] Prakasam | Modulation identification algorithm for adaptive demodulator in software defined radios using wavelet transform[END_REF][START_REF] Maliatsos | Interclass and intraclass modulation recognition using the wavelet transform[END_REF]. This approach is generally divided into two subsystems: the features extraction subsystem and the classifier subsystem [START_REF] Wong | Automatic digital modulation recognition using artificial neural network and genetic algorithm[END_REF]. However, the second approach is more robust and easier to implement if the proper features set is chosen.

In the past, much work has been conducted on modulation identification. The identification techniques, which had been employed to extract the signal features necessary for digital modulation recognition, include spectral-based feature set [START_REF] Nandi | Algorithms for automatic modulation recognition of communication signals[END_REF], higher order cumulants (HOC) [START_REF] Swami | Hierarchical digital modulation classification using cumulants[END_REF][START_REF] Dobre | Robust QAM modulation classification algorithm using cyclic cumulants[END_REF], constellation shape [START_REF] Mobasseri | Digital modulation classification using constellation shape[END_REF], and wavelets transforms [START_REF] Prakasam | Modulation identification algorithm for adaptive demodulator in software defined radios using wavelet transform[END_REF][START_REF] Maliatsos | Interclass and intraclass modulation recognition using the wavelet transform[END_REF]. With their efficient performance in pattern recognition problems (e.g., modulation classification), many studies have proposed the application of artificial neural networks (ANNs) as classifiers [START_REF] Ebrahimzadeh | Intelligent digital signaltype identification[END_REF][START_REF] Zhao | Modulation recognition of SDR receivers based on WNN[END_REF][START_REF] Wong | Automatic digital modulation recognition using artificial neural network and genetic algorithm[END_REF][START_REF] Nandi | Algorithms for automatic modulation recognition of communication signals[END_REF].

In [START_REF] Hong | Identification of digital modulation types using the wavelet transform[END_REF], Hong and Ho studied the use of wavelet transform to distinguish among QAM, PSK, and FSK signals.

In their work, they have used a wavelet transform to extract 2 EURASIP Journal on Advances in Signal Processing the transient characteristics in a digital modulated signal. It has been shown that when the signal-to-noise ratio (SNR) is greater than 5 dB, the percentage of correct identification is about 97%.

In [START_REF] Wong | Automatic digital modulation recognition using artificial neural network and genetic algorithm[END_REF], Wong and Nandi have proposed a method for ADMR using artificial neural networks and genetic algorithms. In their study, they have presented the use of resilient backpropagation (RPROP) as a training algorithm for multi-layer perception (MLP) recogniser. The genetic algorithm is used in [START_REF] Wong | Automatic digital modulation recognition using artificial neural network and genetic algorithm[END_REF] to select the best feature subset from the combined statistical and spectral features set. This method requires carrier frequency estimation, channel estimation, and perfect phase recovery process.

Using the statistical moments of the probability density function (PDF) of the phase, the authors in [START_REF] Lichun | Comments on "Signal classification using statistical moments[END_REF] have investigated the problem of modulation recognition in PSKbased systems. It is shown that the nth moment (n even) of the signal's phase is a monotonically increasing function of the modulation order. On the basis of this property, the study in [START_REF] Lichun | Comments on "Signal classification using statistical moments[END_REF] formulates a general hypothesis testing problem to develop a decision rule and to derive an analytical expression for the probability of misclassification. Similarly, El-Mahdy and Namazi [START_REF] El-Mahdy | Classification of multiple M-ary frequency-shift keying signals over a Rayleigh fading channel[END_REF] developed and analyzed different classifiers for M-ary frequency shift keying (M-FSK) signals over a frequency nonselective Rayleigh fading channel. The classifier in [START_REF] El-Mahdy | Classification of multiple M-ary frequency-shift keying signals over a Rayleigh fading channel[END_REF] employs an approximation of the likelihood function of the frequency-modulated signals for both synchronous and asynchronous waveforms. Employing adaptive techniques, Liedtke [START_REF] Liedtke | Adaptive procedure for automatic modulation recognition[END_REF] proposed an adaptive procedure for automatic modulation recognition of radio signals with a priori unknown parameters. The results of modulation recognition are important in the context of radio monitoring or electronic support measurements. A digital modulation classification method based on discrete wavelet transform and ANNs was presented in [START_REF] Wu | Automatic digital modulation recognition using wavelet transform and neural networks[END_REF]. In this paper, an error backpropagation learning with momentum is used to speed up the training process and improve the convergence of the ANN. This method was developed in [START_REF] Wu | Automatic digital modulation recognition based on ART2A-DWNN[END_REF] by combining adaptive resonance theory 2A (ART2A) with discrete wavelet neural network. It was shown through simulations that high recognition capability can be achieved for modulated signals corrupted with Gaussian noise at 8 dB SNR. Three different automatic modulation recognition algorithms have been investigated and compared in [START_REF] Dayoub | A blind modulation type detector for DPRS standard[END_REF]. The first is based on the observation of the amplitude histograms, the second on the continuous wavelet transform and the third on the maximum likelihood for the joint probability densities of phases and amplitudes.

In [START_REF] Pedzisz | Automatic modulation recognition of MPSK signals using constellation rotation and its 4th order cumulant[END_REF], Pedzisz and Mansour derived and analyzed a new pattern recognition approach for automatic modulation recognition of M-PSK signals in broadband Gaussian noise. This method is based on constellation rotation of the received symbols and fourth-order cumulants of the in-phase distribution of the desired signal. In [START_REF] Li | Automatic digital modulation recognition based on euclidean distance in hyperspace[END_REF], the recognition vector of the decision-theoretic approach and that of the cumulantbased classification are combined to compose a higher dimension hyperspace to get the benefits of both methods. The composed vector is applied to a radial basis function (RBF) neural network, yielding to more reasonable reference points. The method proposed in [START_REF] Li | Automatic digital modulation recognition based on euclidean distance in hyperspace[END_REF] was shown to cover large number of modulation schemes in AWGN channels even under low SNR. In [START_REF] Tadaion | Signal activity detection of phase-shift keying signals[END_REF], Tadaion et al. have derived a generalized likelihood ratio test (GLRT), where they have suggested a computationally efficient implementation thereof. Using discrete wavelet decompositions and adaptive network-based fuzzy inference system, a comparative study of implementation of feature extraction and classification algorithms was presented in [START_REF] Avci | An expert discrete wavelet adaptive network based fuzzy inference system for digital modulation recognition[END_REF].

Also in [START_REF] Su | Real-time modulation classification based on maximum likelihood[END_REF], Su et al. described a likelihood testbased modulation classification method for identifying the modulation scheme of SDR in real-time without pilot transmission. Unlike prior works, the study in [START_REF] Su | Real-time modulation classification based on maximum likelihood[END_REF] converts an unknown signal symbol to an address of a look-up table where it loads the precalculated values of the test functions for the likelihood ratio test to produce the estimated modulation scheme in real-time.

In this paper we focus on the continuous wavelet transform (CWT) to extract the classification features. One of the reasons for this choice is due to the capability of the transform to precisely introduce the properties of the signal in time and frequency [START_REF] Addison | Illustrated Wavelet Handbook[END_REF]. The extracted features are higher order statistical moments (HOM) of the continuous wavelet transform. Our proposed classifier is a multi-layer feed-forward neural network trained using the resilient backpropagation learning algorithm (RPROP). Principal component analysis-(PCA-) based features selection is used to select the best subset from the combined HOM features subsets. This classifier has the capability of recognizing the M-ary amplitude shift keying (M-ASK), M-ary frequency shift keying (M-FSK), minimum shift keying (MSK), Mary phase shift keying (M-PSK), and M-ary quadratic amplitude modulation (M-QAM) signals and the order of the identified modulation. The performance of the proposed algorithm is examined based on the confusion matrix and false recognition probability (FRP). The AWGN channel is considered when developing the mathematical model and through most of the results. Some additional simulations are carried to examine the performance of our algorithm over several fading channel models to assess the performance of our algorithm in a more realistic channel.

The remainder of the paper is organized as follows. Section 2 defines the mathematical model of the proposed problem and presents CWT calculations of different considered digitally modulated signals. Section 3 describes the process of feature extraction using the continuous wavelet transform. Section 4 focuses on features set pre-processing and subset selection, besides the structure of the artificial neural network and the learning algorithm. The results, algorithm performance analysis, and a comparative study with some existing recognition algorithms are presented in Section 5. Conclusions and perspectives of the research work are presented in Section 6.

Mathematical Model

In this study, the properties of the continuous wavelet transform are used to extract the necessary features for modulation recognition. The main reason for this choice is due to the capability of this transform to locate, in time and frequency, the instantaneous characteristics of a signal. More simply, the wavelet transform has the special feature of multiresolution analysis (MRA). In the same manner as Fourier transform can be defined as being a projection on the basis of complex exponentials, the wavelet transform is introduced as projection of the signal on the basis of scaled and time-shifted versions of the original wavelet (so-called mother wavelet) in order to study its local characteristics [START_REF] Addison | Illustrated Wavelet Handbook[END_REF]. The importance of wavelet analysis is its scale-time view of a signal which is different from the time-frequency view and leads to MRA.

The continuous wavelet transform of a received signal s(t) is defined as [START_REF] Addison | Illustrated Wavelet Handbook[END_REF] CWT(a, τ)

= +∞ -∞ s(t)ψ * a,τ (t)dt, ( 1 )
where a > 0 is the scale variable, τ ∈ R is the translation variable, and * denotes complex conjugate. This defines the so-called CWT, where CWT(a, τ) define the wavelet transform coefficients. The Haar wavelet is chosen as the mother wavelet where it is given by [START_REF] Addison | Illustrated Wavelet Handbook[END_REF]

ψ(t) = ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ 1, if 0 ≤ t < T 2 , -1, if T 2 ≤ t < T, 0, otherwise. (2)
The main purpose of the mother wavelet is to provide a source function to generate ψ a,τ (t), which are simply the translated and scaled versions of the mother wavelet, known as baby wavelets, as follows [START_REF] Addison | Illustrated Wavelet Handbook[END_REF]:

ψ a,τ (t) = 1 √ a ψ t -τ a . ( 3
)
Let the received waveform r(t), 0 ≤ t ≤ T s be described as

r(t) = channel [s(t)]. (4)
where T s is the symbol duration and channel is the channel function which includes the channel effect on the signal. For additive white Gaussian noise (AWGN) channel, the received waveform is described as

r(t) = s(t) + n(t), ( 5 )
where n(t) is a complex additive white Gaussian noise. The signal s(t) can be presented as [START_REF] Hong | Identification of digital modulation types using the wavelet transform[END_REF] s(t) = s(t)e j(2π fct+θc) , [START_REF] Wong | Automatic digital modulation recognition using artificial neural network and genetic algorithm[END_REF] where f c is the carrier frequency, θ c is the carrier initial phase, and s(t) is the baseband complex envelope of the signal s(t), defined by

s(t) = √ s N i=1
C i e j(wit+ϕi) g Ts (t -iT s ), [START_REF] Nandi | Algorithms for automatic modulation recognition of communication signals[END_REF] with N being the number of observed symbols, g Ts (t) is the pulse shaping function of duration T s , s is the average signal power, and C i = A i + jB i is the complex amplitude.

In our work we will focus on different M-ary shift keying modulated signals digitalized in RF or IF stages (the carrier frequency is unknown) with respect to SDR principles. That is, it is essential to know that the recognition is done without any priori signal information.

Presenting and calculating the wavelet transform of digitally modulated signals using different modulation schemes will clarify the role of wavelet analysis in feature extraction procedure. The wavelet analysis concept will be studied using only one family of wavelets (Haar wavelet). All the results and figures of CWT presented in this section are obtained using the Haar wavelet. Nevertheless, in our simulations we will extend our results to other families including Daubechies, Morlet, Meyer, Symlet, and Coiflet.

By extending the work of Hong and Ho [START_REF] Hong | Identification of digital modulation types using the wavelet transform[END_REF], from (1)-( 3), [START_REF] Wong | Automatic digital modulation recognition using artificial neural network and genetic algorithm[END_REF], and [START_REF] Nandi | Algorithms for automatic modulation recognition of communication signals[END_REF], the magnitude of continuous wavelet transform is given by

|CWT(a, τ)| = 4S i √ s √ a(w c + w i ) sin 2 (w c + w i ) aT s 4 , (8)
where

S i = |C i | = A 2 i + B 2 i is the amplitude of the ith symbol.
The normalized signal is defined as follows:

s(t) = s(t) |s(t)| = s(t)e j(wct+θc) . (9)
In what follows, the continuous wavelet transform of the normalized signal will be taken into consideration. Knowing that the amplitude of the normalized signal is constant and from [START_REF] Swami | Hierarchical digital modulation classification using cumulants[END_REF], it is clear that the signal normalization will only affect the wavelet transform of nonconstant envelope modulations (i.e., ASK and QAM), and will not affect wavelet transform of constant envelope ones (i.e., FSK, MSK, and PSK). Note that there will be distinct peaks in the wavelet transform of the signal and that of the normalized one resulting from phase changes at the times where the Haar wavelet covers a symbol change. In what follows, we consider the magnitude of the wavelet transforms for different modulation schemes.

Given the complex envelope of QAM signal

s QAM (t) = N i=1 A i + jB i g Ts (t -iT s ), (10)
where (A i , B i ) are the assigned QAM symbols, the corresponding wavelet transform is given by

CWT QAM (a, τ) = 4S i √ aw c sin 2 w c aT s 4 . ( 11
)
It is clear from [START_REF] Prakasam | Modulation identification algorithm for adaptive demodulator in software defined radios using wavelet transform[END_REF] that for a certain scale value, the |CWT| is a multi-step function. Considering the normalized QAM signal:

s QAM (t) = N i=1 e j arctan(Bi/Ai) g Ts (t -iT s ). ( 12
)
The |CWT| is constant since the signal loses its amplitude information. Figure 1 shows the multi-step CWT magnitude of 64-QAM signal and the constant CWT magnitude of normalized 64-QAM signal (as a function of n the translation sampling index).

Let us consider the complex envelope of ASK signal

s ASK (t) = N i=1 A i g Ts (t -iT s ), (13)
where 8), the wavelet transform of ASK signal is given by

A i ∈ {2m -1 -M, m = 1, 2, . . . , M}. From (
|CWT ASK (a, τ)| = 4A i √ aw c sin 2 w c aT s 4 . ( 14
)
It is clear from ( 14) that for a certain scale, the |CWT| of ASK signal is a multi-step function since the amplitude A i is a variable. As for the normalized ASK signals

s ASK (t) = N i=1 sign(A i )g Ts (t -iT s ), (15)
and its corresponding |CWT| is constant. Figure 2 shows CWT magnitude of both 16-ASK signal and its normalized version.

When considering the complex envelope of PSK signals

s PSK (t) = S N i=1 e jϕi g Ts (t -iT s ), (16)
where ϕ i ∈ {(2π/M)(m -1), m = 1, 2, . . . , M}, the wavelet transform is given by

|CWT PSK (a, τ)| = 4 √ S √ aw c sin 2 w c aT s 4 . ( 17
)
It is clear from [START_REF] Wu | Automatic digital modulation recognition using wavelet transform and neural networks[END_REF] that for a certain scale value, the |CWT|of PSK signals is almost a constant function. Given the normalized signal

s PSK (t) = N i=1 e jϕi g Ts (t -iT s ), (18)
the |CWT|is shown to be constant. Also, normalization will not affect wavelet transform of PSK signals since it is a constant envelope signal. Figure 3 shows the constant CWT magnitudes of 16-PSK signal and its normalized version.

For FSK, the complex envelope is defined by:

s FSK (t) = S N i=1 e j(wit+ϕi) g Ts (t -iT s ), (19)
where w i ∈ {w 1 , w 2 , . . . , w M } and ϕ i is the initial phase. From [START_REF] Dayoub | A blind modulation type detector for DPRS standard[END_REF], the wavelet transform of FSK signal is given by and the |CWT|of FSK signal is a multi-step function with w i being a variable. Also, the FSK normalized signal is given by

|CWT FSK (a, τ)| = 4 √ S √ a(w c + w i ) sin 2 (w c + w i ) aT s 4 , (20)
s FSK (t) = N i=1 e j(wit+ϕi) g Ts (t -iT s ). ( 21
)
One can show that |CWT| of the normalized FSK is a multistep function. This is clear from Figure 4, where we show the CWT magnitudes for 16-FSK and its normalized version. Finally, we consider MSK as a special case of continuous phase-frequency shift keying (CPFSK) with modulation index 0.5. The CWT magnitude of MSK signal is expected to be a two-step function similar to 2-FSK signal (Figure 5). (iii) The statistical properties including the mean, the variance, and higher order moments (HOM) of wavelet transforms are different from modulation scheme to another. These statistical properties also differ depending on the order of modulation, since the frequency, amplitude, and other signal properties may change depending on the modulation order.

Features Extraction

(iv) There are distinct peaks in wavelet transforms of different modulated signals and their normalized ones when the Haar wavelet covers a symbol change. Note that the median filtering helps in removing these peaks which will affect |CWT|statistical properties.

According to the above observations, we propose a feature extraction procedure as follows. The CWT can extract features from a digitally modulated signal. These features can be collected by examining the statistical properties of wavelet transforms of both the signal and its normalized one. Since median filtering affects the statistical properties, these properties will be calculated with and without applying filtering. Based on our simulations, we noted that moments of order higher than five will not improve the overall performance of our algorithm. Therefore, in what follows, we consider moments of order up to five to calculate the HOM of wavelet transforms.

Figure 6 shows the processing chain of features extraction. As shown, the digitalized received signal is first normalized then the CWT of the received signal and the normalized one are obtained where the first subset of features will be the HOM (up to 5). A median filter is then applied to cut off the peaks in the corresponding wavelet transforms. Finally the HOM of these two filtered transforms will form the other features subset. This large number of features may contain redundant information about the signal. However, these features will surely have the necessary information to distinguish between different modulations. In order to select a smaller number of features a subset selection algorithm is proposed.

Classifier

The considered ADMR approach is divided into two subsystems: the features extraction subsystem and the classifier subsystem as shown in Figure 7. The ADMR problem

(after features extraction) can be considered as a data classification problem. When the proper features are extracted, one can choose any good algorithm for classification, that is, the classification process is independent from the features extraction process. Some works use the thresholds and decisions trees to classify modulation schemes [START_REF] Prakasam | Modulation identification algorithm for adaptive demodulator in software defined radios using wavelet transform[END_REF][START_REF] Hong | Identification of digital modulation types using the wavelet transform[END_REF], and others employ ANNs to achieve that [START_REF] Ebrahimzadeh | Intelligent digital signaltype identification[END_REF][START_REF] Zhao | Modulation recognition of SDR receivers based on WNN[END_REF][START_REF] Wong | Automatic digital modulation recognition using artificial neural network and genetic algorithm[END_REF][START_REF] Nandi | Algorithms for automatic modulation recognition of communication signals[END_REF].

ANNs were widely employed in the last decades, and they are among the best solutions for pattern recognition and data classification problems. ANNs were proven to increase the recognition performance of modulated signals. For instance the authors in [START_REF] Nandi | Algorithms for automatic modulation recognition of communication signals[END_REF] introduced two algorithms for analog and digital modulations recognition based on the spectral features of the modulated signal. It was shown that the first decision-theoretic algorithm has a poorer performance than the second ANN-based one. In this study, the proposed classifier is a multi-layer feed-forward neural network.

Artificial Neural Network.

ANN is an emulation of biological neural system. ANN is configured through a learning process for a specific application, such as pattern recognition. ANNs with their remarkable ability to derive meaning from complicated or imprecise data can be used to extract patterns that are too complex to be noticed by other computer techniques.

ANN usually consists of several layers. Each layer is composed of several neurons or nodes. The connections among each node and the other nodes are characterized by weights. The output of each node is the output of a transfer function which its input is the summed weighted activity of all node connections. Each ANN has at least one hidden layer besides the input and the output layers. There are two known architectures of ANNs: the feed-forward neural networks and the feedback ones. There are several popular feed-forward neural network architectures such as multi-layer perceptrons (MLPs), radial basis function (RBF) networks, and selforganizing maps (SOMs). We had chosen MLP feed-forward networks in our work because of their simplicity and effective implementations; also they are extensively used in pattern recognition and data classification problems.

Artificial Neural Network Size.

The network size includes the number of hidden layers and the number of nodes in each hidden layer. The network size is an important parameter that affects the generalization capability of ANN. Of course, the network size depends on the complexity of the underlying scenario where it is directly related to network training speed and recognition precision. In this paper the network size has been chosen through intensive simulations.

An improvement can be carried out to our work by using an algorithm that can automatically optimize the neural network size by balancing the minimum size and the good performance, since it is harder to manually search the optimal size. There are several techniques that help to approach the optimal size; some of them starts with huge network size and try to prune it toward the optimal size [START_REF] Castellano | An iterative pruning algorithm for feedforward neural networks[END_REF], others start with small network size and try to increase it toward the optimal size [START_REF] Fahlman | The cascade-correlation learning architecture[END_REF], and some works combine both the pruning and the growing algorithms [START_REF] Lin | Neural Fuzzy Systems: A Neuro-Fuzzy Synergism to Intelligent Systems[END_REF].

Cascade-correlation algorithm (CCA) attempts to automatically choose the optimal network size [START_REF] Fahlman | The cascade-correlation learning architecture[END_REF]. Instead of just adjusting the weights in a network of fixed topology, CCA begins with a minimal network, and then automatically adds new hidden nodes one by one, creating a multi-layer structure. For each new hidden node, CCA attempts to maximize the magnitude of the correlation between the new node's output and the residual error signal which CCA is trying to eliminate.

Features Subset Selection.

The large number of extracted features causes that some among them share the same information content. This will lead to a dimensionality problem. The obvious solution is the features selection, that is, reducing the dimension by selecting some features and discarding the rest. A features space with a smaller dimension will allow more accurate classification (regardless the classifier) due to data organization and projecting data to another space in which the discrimination is more obvious. The output of the features selection process is the input of the feed-forward neural network. Then, features selection also affects the neural network convergence and allows speeding its learning process and reducing its size. Among several possible features selection algorithms, we will investigate principal component analysis (PCA) and linear discriminate analysis (LDA).

PCA constructs a low-dimensional representation of the data (extracted features) that describes as much of the variance in that data as possible. PCA is mathematically defined as an orthogonal linear transformation that transforms the data to a new space such that the greatest variance by any projection of the data comes to lie on the first dimension (called the first principal component), the second greatest variance on the second dimension, and so on [START_REF] Kwak | Principal component analysis based on L1-norm maximization[END_REF]. This moves as much of the variance as possible into the first few dimensions. The values in the remaining dimensions, therefore, tend to be highly correlated and may be dropped with minimal loss of information. PCA is the simplest of the true eigenvector-based multivariate analyses.

Let us suppose that X is the input data (extracted features). PCA attempts to find the linear transformation W which maximizes W T COV (X-X) W, where COV (X-X) is the covariance matrix of the zero-mean data. It can be shown that W is formed of the first d principal eigenvectors (i.e., principal components) corresponding to the greatest d eigenvalues of the covariance matrix. The selected features are given by

P = W * X -X . ( 22
)
LDA is a supervised technique that attempts to maximize the linear separability between data points (features) belonging to different classes (targeted modulation schemes) [START_REF] Fukunaga | Introduction to Statistical Pattern Recognition[END_REF]. It does so by taking into consideration the scatter between-classes besides the scatter within-classes, that is, finds a linear transform so that the between-classes variance is maximized, and the within-classes variance is minimized. The within-classes scatter S w and the between-classes scatter S b are defined as

S w = c∈C p c xc∈c x c -μ c x c -μ c * , S b = c∈C μ c -μ μ c -μ * , ( 23
)
where C is the set of possible classes (modulation schemes), p c is the prior of class c ∈ C, x c is a data point of class c, μ c is the mean of class c and μ represents the mean of all classes. LDA attempts to find the linear transformation W which maximizes the so-called Fisher criterion:

J(W) = W * S b W W * S w W . ( 24
)
LDA seeks to find directions along which the classes are best separated. On the other side, PCA is based on the data covariance which characterizes the scatter of the entire data.

Although one might think that LDA should always outperform PCA (since it deals directly with class separation), empirical evidence suggests otherwise [START_REF] Martinez | PCA versus LDA[END_REF]. For instance, LDA will fail when the discriminatory information is not in the mean but rather in the variance of the data.

Here, a modulation recognition performance comparison shows that LDA slightly outperforms PCA in the poor recognition region, and the performance of the two algorithms rapidly converges as the SNR goes high. Anyway, we will use PCA due to its simplicity and direct implementation.

Training Algorithm. The classification process basically consists of two phases: training phase and testing phase.

A training set is used in supervised training to present the proper network behavior, where each input to the network is introduced with its corresponding correct target. As the inputs are applied to the network, the network outputs are compared to the targets. The learning rule is then used to adjust weights and biases of the network in order to move the network outputs closer to the targets until the network convergence. The training algorithm is mostly defined by the learning rule, that is, the weights update in each training epoch. There are a number of efficient training algorithms for ANNs. Among the most famous is the backpropagation algorithm (BP). An alternative is BP with momentum and learning rate to speed up the training. The weight values are updated by a simple gradient descent algorithm

Δw i j (t) = -ε δE δw i j (t) + μΔw i j (t -1). ( 25
)
The learning rate, ε, scales the derivative, and it has a great influence on training speed. The higher learning rate is, the faster convergence is but with possible oscillation. On the other hand, a small learning value means too many steps are needed to achieve convergence. A variant of BP with adaptive learning rate can be used. The learning rate is adaptively modified according to the observed behavior of EURASIP Journal on Advances in Signal Processing the error function. A BP algorithm employs the momentum parameter, μ, to scale the influence of the previous step on the current. The momentum parameter is believed to render the training process more stable and to accelerate convergence in shallow regions of the error function. However, as practical experience has shown, this is not always true. It turns out in fact, that the optimal value of the momentum parameter is equally problem-dependent as the learning rate.

In this paper, we consider the resilient backpropagation algorithm (RPROP) [START_REF] Riedmiller | Direct adaptive method for faster backpropagation learning: the RPROP algorithm[END_REF]. Basically, RPROP performs a direct adaptation of the weight update based on local gradient information. Only the sign of the partial derivative is used to perform both learning and adaptation. In doing so, the size of the partial derivative does not influence the weight update.

The adaptive update-value Δ i j for RPROP algorithm was introduced as the only factor that determines the size of the weight update. Δ i j evolves during the learning process based on the local behavior of the error function E, according to the following learning rule:

Δ i j (t) = ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ η + * Δ i j (t -1), if δE δw i j (t -1) δE δw i j (t) > 0, η - * Δ i j (t -1), if δE δw i j (t -1) δE δw i j (t) < 0, Δ i j (t -1), else. ( 26
)
where 0 < η -< 1 < η + . The direct adaptation works as follows. Whenever the partial derivative of the corresponding weight changes its sign, which implies that the last update was too large, and the algorithm jumped over a local minimum, the update-value is decreased by the factor η -. If the derivative retains its sign, the update-value is slightly increased (η + ) in order to accelerate convergence in shallow regions.

Once the update-value for each weight is updated, the actual weight update follows a very simple rule as shown in the following equations:

Δw i j (t) = ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ -Δ i j (t), if δE δw i j (t) > 0, +Δ i j (t), if δE δw i j (t) < 0, 0, otherwise ,
w i j (t + 1) = w i j (t) + Δw i j (t). ( 27
)
If the partial derivative is positive (i.e., increasing error), the weight is decreased by its update-value. If the derivative is negative, the update-value is added.

To summarize, the basic principle of RPROP is the direct adaptation of the weight update-value. In contrast to learning rate-based algorithms, RPROP modifies the size of the weight update directly based on resilient updatevalues. As a result, the adaptation effort is not blurred by unforeseeable gradient behavior. Due to the clarity and simplicity of the learning rule, there is only a slight expense in computation compared with ordinary backpropagation. Besides fast convergence, one of the main advantages of RPROP lies in the fact that no choice of parameters and initial values is needed at all to obtain optimal or at least nearly optimal convergence times [START_REF] Riedmiller | Direct adaptive method for faster backpropagation learning: the RPROP algorithm[END_REF]. Also, RPROP is known by its high performance on pattern recognition problems.

After pre-processing and features subset selection, the training process is triggered. The initiated feed-forward neural network is trained using RPROP algorithm. Finally, the test phase is launched and the performance is evaluated through confusion matrix and false recognition probability. Some authors try to explain their results through receiver operating characteristic (ROC) which is more suitable for decision-theoretic approaches where thresholds normally classify modulation schemes.

Results and Discussion

The proposed algorithm was verified and validated for various orders of digital modulation types including ASK, PSK, MSK, FSK, and QAM. Table 1 shows the parameters used for simulations. Testing signals of 100 symbols are used as input messages for different values of SNR and channel effects (AWGN channel is used unless otherwise mentioned).

The wavelet transforms were calculated, and the median filter was applied to extract the features set. Then, preprocessing and features subset selection of 100 realizations of each modulation type/order is performed as a preparation of ANN training. The performance of the classifier was examined for 300 realizations of each modulation type/order, and the results are presented using the confusion matrix and false recognition probability (FRP).

The problem of modulation recognition will be investigated with three scenarios: (i) inter-class recognition (identify the type of modulation only), (ii) intra-class recognition (identify the order of known type of modulation), and (iii) full-class recognition (identify the type and order of the modulation at the same time), as shown in Figure 8.

Performance over AWGN Channel.

The proposed classifier has shown an excellent performance over AWGN channel even at low SNR. Table 2 shows that full-class recognition of modulation schemes (16-QAM, 32-QAM, 64-QAM, 2-PSK, 8-PSK, 4-ASK, 8-ASK, 4-FSK, 8-FSK, and MSK) is achieved with high percentage when the SNR is not lower than 4 dB. Repeating the previous simulations for lower SNR values shows that the full-class recognition gives the lowest percentage for PSK signals.

Simulation results in Table 3 show that when the SNR is not lower than 3 dB, the percentage of correct interclass recognition of ASK, FSK, MSK, PSK, and QAM modulations (case I) is higher than 99%. For lower SNR values, our results show that the inter-class recognition gives the lowest percentage for PSK and FSK, but the inter-class modulation recognition will remain robust for lower SNR values for QAM and ASK signals. We note that, reducing the modulation pool used in simulations to QAM, ASK, and FSK (case II) shows a high percentage of correct interclass modulation recognition for lower SNR value (-2 dB), as shown in Table 4.

Our results show that the intra-class recognition of modulation order using the proposed classifier gives different results depending on the modulation type. For instance, our simulations show that this recognition will be better for ASK and QAM signals than other modulation types, where a high percentage of correct modulation recognition is evident. This property can help in building an adaptive modulation system that assures high quality of service.

Tables 5 and6 show the percentage of correct intraclass modulation recognition at very low SNR for QAM and ASK modulations, respectively. Also Tables 7 and8 show the percentage of correct intra-class modulation recognition for FSK at SNR = 2 dB and PSK at SNR = 4 dB, respectively. The above results demonstrate that our algorithm can achieve high percentage with low SNR for non-constant envelop signals, while it can still achieve the same performance but with higher SNR for constant envelope signals.

Figure 9 shows the performance of the FRP for several recognition cases, where each graph represents FRP when the SNR is not lower that certain value. A minimum SNR for which the FRP is less than 1%, SNR min has been considered in these results. Accordingly, the SNR min for inter-class recognition (Case I) is 3 dB, for inter-class recognition (Case II) is -2 dB, for intra-class PSK recognition is 4 dB, and for intra-class FSK recognition is 2 dB. Generally one can notice that the performance depends on the studied scenario, and it will drop down rapidly for SNRs less than SNR min . This also justifies the SNR values used in producing the results in Tables 2345678and the corresponding high percentage of recognition observed since these SNRs represent the SNR min for each case.

Algorithm Parameters Optimization.

We note that the scaling factor of the CWT has a great effect on the final performance of the classifier. Through extensive simulations, the optimum scaling factor was found to be 10 samples.

Extensive simulations show that the optimal ANN structure to be used for this algorithm is a two hidden layers network (excluding the input and the output layer), where the first layer consists of 10 nodes and the second of 15 nodes.

Let us examine the effect of the number of received symbols, N s , on the algorithm performance. The results of this investigation are shown in Figure 10, where the FRP for several recognition cases is shown at a prescribed N s . Similar to the definition of SNR min , we define N min as the minimum N s value for which FRP is less than 1%. We found that N min for inter-class recognition (Case I) is 100 symbols, for full-class recognition is 100 symbols, for intraclass FSK recognition is 75 symbols, and for intra-class QAM recognition is 50 symbols.

Generally one can notice that the performance depends on the studied scenario, and it will drop down rapidly for number of symbols less than N min .

Figure 11 shows a performance comparison between the two features selection algorithms PCA and LDA. The FRP for inter-class modulation recognition (case II) was examined versus SNR when using each selection algorithm. It is clear that LDA slightly outperforms PCA in the poor recognition region (when SNR < SNR min ). But the two algorithms have the same performance when SNR > SNR min , that is, when the recognition algorithm is well performing. However, in our work we have preferred PCA due to its simplicity and direct implementation. The confusion matrix shows a high percentage of correct full-class modulation recognition when SNR is not lower than 4 dB. The confusion matrix shows a high percentage of correct inter-class modulation recognition (case I) when SNR is not lower than 3 dB. The confusion matrix shows a high percentage of correct inter-class modulation recognition (case II) when SNR is not lower than -2 dB. The confusion matrix shows a high percentage of correct QAM intra-class recognition when SNR is not lower than -6 dB. The confusion matrix shows a high percentage of correct ASK intra-class recognition when SNR is not lower than -4 dB.

So far our results are based on Haar wavelet. Now we examine the proposed algorithm using different wavelet families seeking the optimal wavelet filter to be used. In particular, we provide in Table 9 the total recognition The confusion matrix shows a high percentage of correct PSK intra-class recognition when SNR is not lower than 4 dB.

percentage using several wavelet filters in the case of full-class recognition for SNR = 1 dB. Using Haar wavelet, our previous results show that the SNR min for full-class recognition is 4 dB. That is the reason why the algorithm performance has been investigated at SNR = 1 dB. The poor performance of the algorithm when using Haar wavelet is obvious in comparison to other wavelet families. However, the Haar wavelet, compared to other wavelets, enjoys the simplicity and the easiness of its mathematical modeling. Table 9 shows that the best performance will be found when using Meyer, Morlet, and Biorthgonal 3.5 wavelets. Note that the choice of the best wavelet filter depends on the algorithm implementation and computational complexity of the CWT.

Performance over Fading Channels.

Most of the existing works in the literature had examined their methods over AWGN channel. Here, we also developed our mathematical model and tested our algorithm over this channel. It is clear that it be will more realistic to examine the proposed algorithm performance over fading channels.

The performance of our algorithm has been evaluated in the case of full-class recognition when the SNR is not lower than 4 dB over different fading channels. The examined channel models were derived for the standards and specifications: GSM/EDGE channel models (3GPP TS 45.005 V7.9.0 (2007-2)) [START_REF]3GPP Specification Series[END_REF], COST 207 channel models [START_REF]Digital Land Mobile Radio Communications-COST 207[END_REF], and ITU-R 3G channel models (ITU-R M.1225 (1997-2)) [START_REF]Guidelines for Evaluation of Radio Transmission Technologies for International Mobile Telecommunications-2000[END_REF]. Simulations results in Table 10 show a high modulation recognition percentage over fading channels. This confirms that our algorithm has a robust performance regardless of the channel model used. Finally in Figure 12, we present the FRP as a function of the mobile speed (x) in a GSM system considering the channel model of rural area (RAx, 6 taps) as introduced in the 3GPP GSM/EDGE channel model. It is clear from these results that the algorithm performance is better over fading channels than static ones. These results confirm the special capabilities of wavelet analysis relative to conventional analysis.

Performance Comparison.

The comparison among different modulation recognition algorithms is not straightforward. This is mainly because of the fact that there are no available standard digital modulation databases. Hence, different works have applied their algorithms to cases of their own choosing [START_REF] Mobasseri | Digital modulation classification using constellation shape[END_REF].

Also, the different studies considered different modulation pools and different simulation configurations which will result in different and incomparable performances. Some algorithms need a priori information of the signal, for example, carrier frequency [START_REF] Ebrahimzadeh | Intelligent digital signaltype identification[END_REF][START_REF] Wong | Automatic digital modulation recognition using artificial neural network and genetic algorithm[END_REF], frequency offset [START_REF] Dobre | Robust QAM modulation classification algorithm using cyclic cumulants[END_REF], and channel information [START_REF] Dobre | Likelihood-based algorithms for linear digital modulation classification in fading CHAN-NELS[END_REF].

The authors in [START_REF] Ebrahimzadeh | Intelligent digital signaltype identification[END_REF] employed HOC and HOM of the baseband recovered complex envelope. A feed-forward neural network trained with RPROP algorithm was used as a classifier. The modulation schemes M-PSK, M-ASK, and M-QAM were examined. The performance of their classifier is higher than 93% for SNR > 4 dB and 98% when SNR > 8 dB. In [START_REF] Prakasam | Modulation identification algorithm for adaptive demodulator in software defined radios using wavelet transform[END_REF], the developed algorithm was verified using wavelet transform and histogram computations to identify M-PSK, M-QAM, M-FSK and GMSK. When SNR is above 5 dB, the probability of detection of this algorithm is more than 97.8%. In [START_REF] Pedzisz | Automatic modulation recognition of MPSK signals using constellation rotation and its 4th order cumulant[END_REF], the perfect classification between M-PSK signals can be obtained when SNR > 8 dB (considering 256 symbols) and SNR > 10 dB (considering 100 symbols).

A comparison between our results and the abovementioned ones will show the high performance of our algorithm. The recognition probability of M-ASK, M-PSK, M-QAM, M-FSK, and MSK is higher than 99% when SNR is not lower than 4 dB. M-PSK signals classification percentage is higher than 99% when SNR is not lower than 4 dB (considering 100 symbols). Initially our work was based on that introduced in [START_REF] Hong | Identification of digital modulation types using the wavelet transform[END_REF]. Simulations on the same modulation pool examined in [START_REF] Hong | Identification of digital modulation types using the wavelet transform[END_REF] (16-QAM, 4-PSK, and 4-FSK) show that our recognition percentage is higher than 99% when SNR is not lower than -3 dB. This outperforms the algorithm proposed in [START_REF] Hong | Identification of digital modulation types using the wavelet transform[END_REF], where the percentage of correct recognition was about 97% when SNR is greater than 5 dB. Not very much modulation recognition studies investigated the robustness of their methods over fading channels. The authors in [START_REF] Dobre | Likelihood-based algorithms for linear digital modulation classification in fading CHAN-NELS[END_REF][START_REF] El-Mahdy | Classification of multiple M-ary frequency-shift keying signals over a Rayleigh fading channel[END_REF] considered the Rayleigh fading channel model.

It is essential to focus on the fact that in our algorithm the different modulated signals are digitalized in RF or IF stages (the carrier frequency is unknown) with respect to SDR principles. The recognition is done without any priori signal information, and our algorithm shows robustness over fading channels.

Conclusion

We presented a wavelet-based algorithm for automatic modulation recognition. The proposed algorithm is capable of recognizing different modulation schemes with high accuracy at low SNRs. Our classifier has high full-class modulation recognition performance when the SNR is not lower than 4 dB. We found that the percentage of correct inter-class recognition for MSK, FSK, ASK, PSK, and QAM is high when the SNR is not lower than 3 dB. Also, the percentage of the correct intra-class recognition of modulation order was found to be high for low SNRs, and the minimum value of SNR for which the high percentage of intra-class recognition is still reachable depends on the modulation type, and could reach a very low value (-6 dB for intra-class QAM recognition). In addition, we have shown that our algorithm offers excellent performance over both AWGN and several fading channel models.

Figure 1 :Figure 2 :

 12 Figure 1: Multi-step wavelet transform of QAM64 signal and constant wavelet transform of its normalized version.

Figure 3 :

 3 Figure 3: Constant wavelet transform of PSK16 signal and its normalized version.

5 |CWT|Figure 4 :

 54 Figure 4: Multi-step wavelet transform of FSK16 signal and its normalized version.

Figure 5 :

 5 Figure 5: Multi-step wavelet transform of MSK signal and its normalized version.

Figure 6 :Figure 7 :

 67 Figure 6: The processing chain of different features subsets extraction.

Figure 8 :

 8 Figure 8: Modulation recognition scenarios including inter-class, intra-class, and full-class recognition.

Figure 9 :

 9 Figure 9: False recognition probability versus SNR. (a) Inter-class recognition (Case I). (b) Inter-class recognition (Case II). (c) Intraclass PSK recognition. (d) Intra-class FSK recognition.

  recognition Full-class recognition Intra-class FSK recognition Intra-class QAM recognition

Figure 10 :

 10 Figure 10: False recognition probability versus number of symbols (inter-class recognition, full-class recognition, intra-class QAM recognition, and Intra-class FSK recognition).

Figure 11 :

 11 Figure 11: Performance comparison between PCA and LDA features selection algorithm (False Recognition Probability versus SNR for Inter-class recognition case II).

FRPFigure 12 :

 12 Figure 12: False recognition probability versus mobile speed in typical case for rural area for full class recognition and SNR = 0 dB.

Table 1 :

 1 Modulation parameters.

	Parameter	Value
	Sampling frequency, Fs	1.5 MHZ
	Carrier frequency, Fc	150 kHZ
	Symbol rate, Rs	12500 Symbol/s
	No. symbols, Ns	100 Symbols
	Simulation parameters of digital modulation used in training, validation,
	and evaluation of the proposed algorithm.

Table 2 :

 2 Confusion matrix at SNR = 4 dB.

		QAM16	QAM32	QAM64	PSK2	PSK8	ASK4	ASK8	FSK4	FSK8	MSK
	QAM16	99.2			0.1		0.7
	QAM32		99.3		0.2		0.5
	QAM64		0.2	99.6			0.2
	PSK2				99.4	0.6
	PSK8				0.7	99.3
	ASK4						99.5	0.5
	ASK8						0.6	99.4
	FSK4								99.4	0.4	0.2
	FSK8								0.6	99.1	0.3
	MSK								0.2	0.1	99.7

Table 3 :

 3 Confusion Matrix at SNR = 3 dB.

		QAM	PSK	ASK	FSK	MSK
	QAM	99.7	0.2	0.1
	PSK	0.4	99.2	0.2	0.1	0.1
	ASK	0.3	0.1	99.4	0.1	0.1
	FSK		0.1		99.6	0.3
	MSK		0.1		0.3	99.6

Table 4 :

 4 Confusion matrix at SNR =-2 dB.

		QAM	ASK	FSK
	QAM	99.5	0.5
	ASK	1	99
	FSK	0.1	0.1	99.8

Table 5 :

 5 Confusion Matrix at SNR =-6 dB.

		QAM16	QAM32	QAM64
	QAM16	99.5	0.4	0.1
	QAM32	0.2	99.3	0.5
	QAM64	0.1	0.1	99.8

Table 6 :

 6 Confusion Matrix at SNR =-4 dB.

		ASK2	ASK4	ASK8
	ASK2	99.4	0.4	0.2
	ASK4	0.1	99.1	0.8
	ASK8	0.1	0.6	99.3

Table 7 :

 7 Confusion Matrix at SNR = 2 dB.

		FSK2	FSK4	FSK8
	FSK2	99.7	0.2	0.1
	FSK4	0.1	99.2	0.7
	FSK8	0.1	0.8	99.1

The confusion matrix shows a high percentage of correct FSK intra-class recognition when SNR is not lower than 2 dB.

Table 8 :

 8 Confusion Matrix at SNR = 4 dB.

		PSK2	PSK4	PSK8
	PSK2	99.9	0.1
	PSK4		99.8	0.2
	PSK8		0.5	99.5

Table 9 :

 9 Algorithm performance and wavelet family.

	Wavelet filter	Full-class recognition for SNR = 1
	Haar	80.2%
	Daubechies 2	91.1%
	Daubechies 3	93.7%
	Daubechies 5	92.2%
	Daubechies 8	93.3%
	Biorthogonal 1.3	92.0%
	Biorthogonal 2.2	92.4%
	Biorthogonal 2.8	96.2%
	Biorthogonal 3.5	98.1%
	Biorthogonal 6.8	92.7%
	Coiflet 1	90.7%
	Coiflet 2	92.7%
	Coiflet 3	93.7%
	Coiflet 4	88.7%
	Coiflet 5	94.7%
	Symlet 2	90.7%
	Symlet 3	93.7%
	Symlet 5	93.7%
	Symlet 7	92.7%
	meyer	98.7%
	Morlet	98.6%

Algorithm performance of full-class recognition for SNR = 1 dB using different wavelet families.

Table 10 :

 10 Confusion Matrix at SNR = 4 dB.

	Channel model	Full-class recognition
	GSM/EDGE channel models
	Typical case for rural area (RA100), 6 taps	99.7%
	Typical case for hilly terrain (HT80), 12 taps	99.1%
	Typical case for hilly terrain (HT80), 6 taps	98.3%
	Typical case for urban area (TU60), 12 taps	99.2%
	Typical case for urban area (TU60), 6 taps	97.7%
	Typical case for very small cells (TIx)	99.2%
	COST 207 channel models
	Rural Area (RA100), 6 taps	99.2%
	Typical Urban (TU60), 12 taps	99.7%
	Bad Urban (BU60), 12 taps	97.9%
	Hilly Terrain (HT80), 12 taps	99.3%
	ITU-R 3G channel models
	Indoor office (IA5)	98.4%
	Outdoor to indoor and pedestrian (PA10)	99.7%
	Vehicular -high antenna (VA100)	97.7%
	Satellite, LOS (SA100LOS)	99.3%
	Satellite, NLOS (SA100NLOS)	99.1%

Algorithm performance of full-class recognition for SNR = 4 dB over different fading channels.


Raphaël Chenouard

Laurent Granvilliers
email: laurent.granvilliers@univ-nantes.frraphael.chenouard@irccyn.ec-nantes.fr

Ricardo Soto

High-Level Modeling of Component-based CSPs

Most of modern constraint modeling languages combine rich constraint languages with mathematical notations to tackle combinatorial optimization problems. Our purpose is to introduce new componentoriented language constructs to manipulate hierarchical problems, for instance for modeling engineering system architectures with conditional sub-problems. To this end, an object-oriented modeling language is associated with a powerful constraint language. It offers the possibility of defining conditional components to be activated at solving time, declaring polymorphic components whose concrete types have to be determined, and overriding model elements. We illustrate the benefits of this new approach in the modeling process of a difficult embodiment design problem having several architectural alternatives.

Introduction

Constraint programming (CP) is a generic framework allowing users to state various kinds of constraint satisfaction problems (CSPs). Several specialized paradigms emerged to tackle specific concepts, for instance conditional CSPs [START_REF] Gelle | Solving Mixed and Conditional Constraint Satisfaction Problems[END_REF] and composite CSPs [START_REF] Sabin | Configuration as Composite Constraint Satisfaction[END_REF].

A CSP is conditional when its sub-problems may not participate in the solutions. The activation of sub-problems is subjected to constraints to be verified during the solving process. A composite CSP can be seen as a hierarchy of sub-problems. It captures the inherent structure of a component-based system. Combining both features, composite and conditional CSPs make it possible to formulate models of complex systems and to study some variations in their architectures. This class of problems finds many applications in the configuration and design areas [START_REF] O'sullivan | Constraint-Aided Conceptual Design[END_REF][START_REF] Gelle | Solving Mixed and Conditional Constraint Satisfaction Problems[END_REF][START_REF] Junker | Handbook of Constraint Programming, chapter Configuration[END_REF].

Modeling such problems is not trivial, as it is necessary to handle various concerns e.g., the component hierarchy, the activation of sub-problems, and even the definition of families of components sharing some properties. In this paper, we present a new approach to elegantly represent component-based CSPs. We extend the s-COMMA modeling language [START_REF] Soto | The Design of COMMA: An Extensible Framework for Mapping Constrained Objects to Native Solver Models[END_REF] by adding a new set of language constructs. We firstly introduce the concept of conditional (optional) component and we provide a simple construct to employ them in constraint models.

Depending on constraint-based conditions, it is possible to dynamically modify the problem architectures by activating new components. Second, components can be polymorphic. It is possible to define families of components, whose types can be refined at solving time. Third, as usual in object-oriented languages, model elements can be overridden in sub-types to facilitate the implementation of sub-classes. Finally, the language also makes it possible to handle components off the shelf or catalogs by means of compatibility constraints, which are sets of tuples to be assigned to a given set of variables.

Modeling complex real-world CSPs is a really hard task and s-COMMA is the result of several years of experience tackling that kind of problems. Now, let us illustrate such new modeling features by presenting a real-world CSP from Dassault Aviation about the design of an air conditioning system (ACS) for aircrafts.

Example

The ACS in an aircraft roughly mixes hot air from the turbojet (main air) and cold air from atmosphere (ram air) to inject air at right temperature and pressure in the cabin (see [START_REF] Chenouard | Solving an Air Conditioning System Problem in an Embodiment Design Context Using Constraint Satisfaction Techniques[END_REF] for more details). Two topologies connecting several components are depicted in Fig. 1. On the left, only one heat exchanger is used to transfer the calorific energy between the main air and the ram air. On the right, a second heat exchanger is used to pre-cool the main air. Several other topological changes may also be considered in a more complete analysis of this complex system. These architectures are hierarchical, since every component can be refined as a complex system. Moreover, the choice of topology may depend on conditions. In this problem, the goals may be to simultaneously study:

the embodiment design of the ACS, which in particular requires the solving of nonlinear constraints over continuous variables related to the geometry or the physical behavior of the system, -the feasibility of each architectural design alternative. Thus, we need a high-level modeling language being able to represent at least components subject to constraints, discrete and continuous variables, and conditional components.

Outline

The remaining of this paper is organized as follows. s-COMMA and the new modeling features are presented in Section 2. The architecture supporting the modeling language is shortly described in Sect. 3. The related work and conclusions follow.

Modeling

The core of the s-COMMA modeling language is briefly described now. A model is a collection of classes to be linked through usual composition and inheritance relations. The main class corresponds to the model entry, as follows.

[ 2 0 0 , 1 0 0 0 ] ; // " i n K o f i n p u t c o l d a i r 30 r e a l THotOut i n [ 2 0 0 , 1 0 0 0 ] ; // " i n K o f o u t p u t h o t a i r 31 r e a l TColdOut i n [ 2 0 0 , 1 0 0 0 ] ; // " i n K o f o u t p u t c o l d a i r 32 . . .
In this model, the aircraft is composed of an air conditioning system acs and a turbojet tj, being instances of other classes. A constraint block Airflows is defined to link the temperatures of air flows between the atmosphere, the turbojet, the ACS, and the cabin. A set of equality constraints is simply stated between air temperature attributes of acs and tj, and the flight conditions, namely the atmosphere air temperature TRam and the cabin air temperature TCab. These two parameters can be defined as constants in a data file for given flight conditions, as follows:

1 // Data f i l e 2
r e a l TRam := 2 2 0 . 0 5 5 ; // a t m o sp h e r e a i r t e m p e r a t u r e i n K a t 1 0 5 0 0m 3 r e a l TCab := 2 8 0 ; // e x p e c t e d c a b i n a i r t e m p e r a t u r e i n K

The ex1 attribute corresponds to the heat exchanger in the first topology, and to the first one in the second topology. That means that every topology at least requires one exchanger. The presence of one or two exchangers will be more precisely discussed in Sect. 2.2.

For the sake of clarity, only a small extract of the full model is presented above. The other components are defined in the same way. In the rest of the section we focus on the new component-oriented language constructs.

Catalogs

In most configuration and design problems, it is necessary to reuse some components from the shelf, given as catalogs stating tuples of values for their main characteristics. For instance, a heat exchanger is made of exchange surfaces to be chosen from a set of known shapes. This can simply be modeled by means of a Boolean formula that restricts the possible values of several variables such as rh, bh, and betah. This is a natural representation of catalogs, but it is not suitable to deal with a higher number of constraints, as the model becomes confuse. To avoid this, it is possible to compact the model by using a compatibility constraint.

Conditional components

In various CP problems, it is necessary to model elements whose activation depends on the satisfaction of some conditions. For instance, consider the heat exchanger whose properties vary depending on the temperature of the input hot air (THotIn). This is commonly modeled by adding a set of implications that activate the involved constraints. The model becomes more complicated if we need to model the activation of a new object depending on a condition. For instance, in the second ACS topology, an additional heat exchanger is required when the input hot temperature exceeds 600 K. One manner to state that is to include the second heat exchanger within the class and to activate the corresponding constraints. Additionally, it is necessary to deactivate the object attributes (when the condition is not satisfied) in order to avoid useless splitting operations during the search. This can roughly be done by assigning to the variables its minimum domain value. Unfortunately, the resultant model is too verbose and hard to understand. To make the definition of such a conditional formulation more concise and understandable, we introduce a new conditional statement. This new construct considers a name, a condition, and a sequence of elements to activate if the condition is satisfied. For instance, the constraints defined within the AirflowsOne constraint block are only activated if the OneExchanger conditional block evaluates to true. The integration of conditional objects is performed in the same way. This feature is more powerful since the set of variables and constraints embedded in the object can be activated at once, as for instance in the second ACS topology. The second heat exchanger is activated only if the THotIn attribute exceeds 600 K. It is also possible to include additional constraints acting over the conditional object, for example to establish the new links between air flows, and to compute the new heat efficiency. Note that the conditional block can also be located in a higher class (considering the composition's hierarchy), specifying that they apply on a given object. For instance, now the conditional block is stated within the Aircraft class, but it remains acting over the acs object. The main difference here is to modify a single object rather than the class itself. That allows one to manipulate instances of a same class with few modifications.

Polymorphic components

Abstract or non final classes in object-oriented programming languages allow one to uniformly manipulate instances of sub-classes. From a modeling viewpoint, it is common to specify sub-types corresponding to families of systems differing in some aspects. Our goal is to give users a mean for declaring some objects of a given type and to let the solving engine dynamically refine this type, i.e., to allow non deterministic choices of component types. We just define the new polymorphic keyword to be used in the declaration of a variable. As a consequence, every polymorphic object of type T occurring in a solution may have type T or a sub-type of T. In this case, we impose that the full type list from the inheritance hierarchy has to be explored. This feature is illustrated by considering three main families of heat exchangers depending on directions of the cold and hot air flows: co-current flows, counter-current flows, and cross-current flows as depicted in Fig. 2. Each type from this family has its own way to compute the exchange of heat energy between the two air flows. The family of exchangers can be naturally represented by an inheritance graph as shown in Fig. 2. It directly follows the following classes. Now let the (first) heat exchanger from the ACS to be polymorphic. Furthermore, suppose that the second heat exchanger must be of type CRCHExchanger and that, in this case, the first exchanger must be of the same type. The following piece of model results. co n s tr a i n t TypeExchangerOne { // t y p e r e s t r i c t i o n 14 typeOf ( ex1 ) = typeOf ( ex2 ) ; // f o r t h e f i r s t e x c h a n g e r 15 } 16 . . . 17 } Introducing polymorphic instances leads to the need for manipulating instance types. In the example, the type of the first exchanger is restricted using the typeOf keyword (lines 13 to 15). That may be compared with runtime type checking mechanisms in computer programming languages, such as instanceof in Java or type in Python. The main difference is the declarative nature of the primitives in our language, which are just constraints on types.

Overriding elements

The purpose of overriding mechanisms is to allow a sub-type from a hierarchy to provide a specific definition of an element -variable, constraint block, conditional block -defined in some ascending type. To this end, it suffices to declare an element using a naming correspondence.

In the ACS, the areas of the exchange surfaces are defined in the HeatExchanger base type. They are valid for all the given sub-types, except for the cold exchange surface of CRCHExchanger heat exchangers. In this case, it suffices to override the ColdExchangeSurfaceArea constraint, as follows. 3 Architecture s-COMMA has been implemented on a model-driven solver-independent platform. Such an architecture allows the automatic reformulation of a s-COMMA model into different executable solver models. A model-driven engineering approach has been implemented to perform the necessary transformations among source, intermediate, and target models. Although, those transformations are not detailed here for space reasons, an extended presentation about them can be found in [START_REF] Soto | The Design of COMMA: An Extensible Framework for Mapping Constrained Objects to Native Solver Models[END_REF][START_REF] Chenouard | Model-Driven Constraint Programming[END_REF][START_REF] Chenouard | Rewriting Constraint Models with Metamodels[END_REF]. A general scheme of the reformulation process is depicted in Fig. 3. An s-COMMA model is the input of the system. This model is mapped to an intermediate model in order to reformulate the language constructs not supported at the solver level. Those transformations are standard and are based on the previous implementations of s-COMMA. For instance, objects are flattened i.e., every variable and constraint are just incorporated into the intermediate model. Compatibility constraints are mapped to the corresponding set of Boolean expressions. Conditional blocks are transformed into a set of implications, and polymorphic elements are stated as single conditional blocks. These conditional blocks are then treated as logical formulas. Once the reformulation is finished, a solver model is generated from the intermediate one, which is launched to obtain the results. Finally, the solution set is analyzed according to a traceability model to display the solutions in a correct object-oriented format.

Related Work

Zinc [START_REF] Marriott | The Design of the Zinc Modelling Language[END_REF], MiniZinc [START_REF] Nethercote | MiniZinc: Towards A Standard CP Modelling Language[END_REF] (successors of OPL [START_REF] Van Hentenryck | The OPL Optimization Programming Language[END_REF]), and Essence [START_REF] Frisch | Essence : A constraint language for specifying combinatorial problems[END_REF] are the state-ofthe-art constraint modeling languages. They provide a rich semantics for modeling different kinds of problems. In particular, Zinc offers the possibility of tackling specific applications domains by means of user-defined predicates. However, a main problem is that no object-oriented features are provided. Hence, it is difficult to naturally capture the structure of component-based problems.

Compatibility constraints can undoubtedly be represented by these languages as Boolean formulas, and conditional blocks can be simulated by using Boolean variables and/or implications as illustrated in Sect. 2. But, the result seems to be a hack of the model rather than a natural problem formulation. Finally, there is no mechanism to manage families of components through polymorphism.

COB [START_REF] Jayaraman | Modeling Engineering Structures with Constrained Objects[END_REF] is an object-oriented constraint modeling language that can capture hierarchical problem structures. It provides a common language to post constraint models, which can be enriched with Prolog-like predicates. But, once again it lacks of compatibility constraints, conditional statements, and polymorphism. Additionally, less relevant to this paper, but not less important, the underlying architecture of COB is solver-dependent, being not possible to launch a model in different solving engines.

Conclusion

In this paper, we have presented new language constructs for modeling componentbased CSPs. These constructs allow one to define more concise and understandable models, that naturally represents complex hierarchical problems. For instance, compatibility constraints avoids the definition of confuse Boolean expressions. Conditional blocks are a mean to express variables, constraints, and even objects whose relevance on the model depends on a given condition. Polymorphism is suitable for modeling families of components sharing some properties. Such constructs have been integrated on the s-COMMA language with the spirit of designing a full component-oriented language, an important CP challenge as stated in [START_REF] Junker | Handbook of Constraint Programming, chapter Configuration[END_REF].

Fig. 1 .

 1 Fig. 1. Two possible topologies of an air conditioning system in an aircraft.

  a r i a b l e s r e l a t e d t o t h e e x c h a n g e s u r f a c e s 4 r e a l bh i n [ 0 , 0 . 1 ] ; 5 r e a l r h i n [ 0 , 0 . 1 ] ; 6 r e a l b e t a h i n [ 1 0 0 , 1 . 0 e 4 ] ; 7 . . . 8 co n s tr a i n t S u r f a c e s C a t a l o g { 9 ( r h =7.7089 e-4 and bh =6.35 e-3 and b e t a h =1.204 e +3) o r 10 ( r h =3.61315 e-3 and bh =1.905 e-3 and b e t a h =2.496 e +2) o r 11 ( r h =6.6167 e-4 and bh =1.051 e-2 and b e t a h =1.368 e +3) o r 12 ( r h =6.6992 e-4 and bh =9.525 e-3 and b e t a h= 1 . 2 5 e +3) ; 13 }

1 //

 1 Within t h e ACS c l a s s 2 co n s tr a i n t A i r f l o w s O n e { 3 ( THotIn <= 6 0 0 ) -> ex1 . TColdOut = n o z z l e . TIn ; 4 ( THotIn <= 6 0 0 ) -> THotIn = c o m p r e s s o r . TIn ; 5 } 6 co n s tr a i n t H e a t E f f i c i e n c y O n e { 7 ( THotIn <= 6 0 0 ) -> e p s i l o n = ex1 . e p s i l o n ; 8 } 9 }

1 //

 1 Within t h e ACS c l a s s 2 HeatExchanger ex2 ; // t h e s e c o n d e x c h a n g e r 3 4 co n s tr a i n t A i r f l o wsTwo { 5 ( THotIn > 6 0 0 ) -> ex1 . TColdOut = ex2 . TColdIn ; 6 ( THotIn > 6 0 0 ) -> ex2 . TColdOut = n o z z l e . TIn ; 7 ( THotIn > 6 0 0 ) -> THotIn = ex2 . THotIn ; 8 ( THotIn > 6 0 0 ) -> ex2 . THotOut = c o m p r e s s o r . TIn ; 9 } 10 co n s tr a i n t H e a t E f f i c i e n c y T w o { 11 ( THotIn > 6 0 0 ) -> e p s i l o n = ( ex1 . e p s i l o n + ex2 . e p s i l o n ) / 2 ; 12 } 13 }

1

  co n s tr a i n t d e a c t i v a t e A t t r i b u t e s { 2 n o t ( THotIn > 6 0 0 ) -> ( ex2 . THotIn = minDom( ex2 . THotIn ) ) ; 3 n o t ( THotIn > 6 0 0 ) -> ( ex2 . TColdIn = minDom( ex2 . TColdIn ) ) ; 4 n o t ( THotIn > 6 0 0 ) -> ( ex2 . THotOut = minDom( ex2 . THotOut ) ) ; 5 n o t ( THotIn > 6 0 0 ) -> . . . 6 }

1 //

 1 Within t h e ACS c l a s s 2 cond OneExchanger ( THotIn <= 6 0 0 ) { co n s tr a i n t A i r f l o w s O n e { 4 ex1 . TColdOut = n o z z l e . TIn ; 5 THotIn = c o m p r e s s o r . TIn ; 6 } 7 co n s tr a i n t H e a t E f f i c i e n c y O n e { 8 e p s i l o n = ex1 . e p s i l o n ; 9 } 10 }

1 c l a s s A i r c r a f t { 2

 2

Fig. 2 .

 2 Fig. 2. A simple inheritance hierarchy for heat exchangers.

1 c l a s s

  HeatExchanger { . . . } 2 c l a s s COCHExchanger extends HeatExchanger { . . . } 3 c l a s s CTCHExchanger extends HeatExchanger { . . . } 4 c l a s s CRCHExchanger extends HeatExchanger { . . . }

1 //

 1 Within t h e ACS c l a s s 2 polymorphic HeatExchanger ex1 ; // p o l y m o r p h i c i n s t a n c e 3 4 co n s tr a i n t A i r f l o w s { // s h a r e d c o n s t r a i n t between 5 ex1 . THotOut = t u r b i n e . TIn ; // t h e two t o p o l o g i e s 6 ex1 . TColdIn = d i f f u s e r . TOut ; 7 . . . 8 } 9 10 cond TwoExchangers ( THotIn > 6 0 0 ) { 11 CRCHExchanger HeatExchanger ex2 ; // t h e s e c o n d e x c h a n g e r 12 13

Fig. 3 .

 3 Fig. 3. Reformulation process.

  THotIn

	1	// Model f i l e
	2 main c l a s s A i r C r a f t {
	3	ACS a c s ;			// t h e ACS
	4	TurboJet t j ;			// t h e t u r b o j e t
	5
	6	co n s tr a i n t A i r f l o w s {
	7	t j . TIn = TRam ;		// t j i n p u t = a t m o sp h e r e
	8	t j . TOut = a c s . THotIn ; // t j o u t p u t = ACS i n p u t
	9	a c s . THotOut = TCab ;	// ACS o u t p u t = c a b i n
	10	a c s . TColdIn = TRam ;	// ACS i n p u t = a t m o sp h e r e
	11	. . .
	12	}
	13 }
	14
	15	c l a s s ACS {
	16	HeatExchanger ex1 ; // t h e ( f i r s t ) h e a t e x c h a n g e r
	17	Tu r b i n e t u r b i n e ;	// t h e t u r b i n e
	18
	19	r e a l THotIn	i n [ 2 0 0 , 1 0 0 0 ] ; // t e m p e r a t u r e i n K o f i n p u t h o t a i r
	20	r e a l TColdIn i n [ 2 0 0 , 1 0 0 0 ] ; //	"	i n K o f i n p u t c o l d a i r
	21	r e a l THotOut i n [ 2 0 0 , 1 0 0 0 ] ; //	"	i n K o f o u t p u t h o t a i r
	22	r e a l TColdOut i n [ 2 0 0 , 1 0 0 0 ] ; //	"	i n K o f o u t p u t c o l d a i r
	23	r e a l e p s i l o n i n [ 0 , 1 ] ;	// h e a t t r a n s f e r e f f i c i e n c y
	24	. . .
	25 }
	26
	27	c l a s s HeatExchanger {
	28	r e a l

i n [ 2 0 0 , 1 0 0 0 ] ; // t e m p e r a t u r e i n K o f i n p u t h o t a i r 29 r e a l TColdIn i n

We believe that an important research direction is about solving strategies for component-based CSPs, interesting studies are pointed out in [6,5,14,10]. The main problem is that the set of required techniques may not be present in any existing solver, such as consistency techniques for mixed variables, interval methods for continuous problems, various search strategies, decomposition algorithms, and so on. Developing a new solver is a really hard task, and the solution may be to make several tools cooperate inside a CP platform.