Nesrine Changuel
email: nesrine.changuel@alcatel-lucent.com

Bessem Sayadi
email: bessem.sayadi@alcatel-lucent.com

Michel Kieffer
email: kieffer@lss.supelec.fr

Inter-frame dependent rate and distortion models for statistical multiplexing of video programs

Keywords: R-D models, H264/AVC coding, Statistical multiplexing, video broadcasting

Statistical multiplexing (SM) is a useful technique for transmitting multiple video streams over a channel while satisfying multiple simultaneous bandwidth and quality related constraints. Design of efficient SM algorithms is based on Rate and Distortion (R-D) information extracted from successive encoded frames. Good R-D models usually take into account inter-frame dependencies specified in terms of impact of encoding parameters of a given frame over the next one. A highly accurate low complexity rate and distortion model is proposed in this paper. A minimum number of R-D measurement points per stream is shown to be sufficient to predict the effect of interframe dependencies. Design of efficient SM algorithms satisfying the above mentioned performance constraints on the basis of the proposed R-D model is described.

I. INTRODUCTION

Dynamically allocating transmission bandwidth between multiple video streams is one of the major challenges of broadcast systems such as DVB or S-DMB. Statistical multiplexing [START_REF] He | Linear rate control and optimum statistical multiplexing for H.264 video broadcast[END_REF] based on joint control of multiple video coders offers a good quality smoothing solution in this situation. Additionally, the attempt is made to maximize the number of video streams sharing the same transmission channel in the face multiple simultaneous constraints. These are related to bandwidth usage, minimum required video quality, smooth quality transitions over time and balanced quality between streams [START_REF] Changuel | Predictive control for efficient statistical multiplexing of digital video programs[END_REF]. Rate and Distortion (R-D) characteristics of the transmitted streams are known to offer a useful tool for solving the resulting constrained optimization problem.

The R-D characteristics for a video stream can be modeled either by considering each frame as an independent unit or by taking into account the inter-frame dependencies. The latter is used for performing motion compensation between a reference and a predicted frame. Several parametric models representing R-D behavior of video coders for independent frames have been proposed. For example [START_REF] Xie | Sequence-based rate control for constant quality video[END_REF] introduces a quadratic model used to solve the rate control problem for constant quality video. Nevertheless, this model does not take into account the control parameters of the video encoder. In [START_REF] Milani | A rate control algorithm for the H.264 encoder[END_REF], a ρ-domain model is proposed, where ρ indicates the number of null coefficients in a DCT encoded/quantified image block. This model offers a reduced complexity solution due to its linear structure but it does not consider the inter-frame dependencies.

The quality of motion-compensated inter-Predicted (P) or Bidirectionally-predicted (B) frames is strongly impacted by the quality of the corresponding reference frame(s). Thus dependent R-D models are more complex than if independence between successive frames is considered. Moreover, the number of R-D measurement points required for a good quality fit is higher for such models. The dependent R-D model presented in [START_REF] Lin | Bit-rate control using piecewise approximated rate-distortion characteristics[END_REF] incorporates the influence of current and previous quantization parameters as well as the effect of texture and motion information. However, a large number of R-D measurements is needed for a good parametrical fit -12 for the distortion model and 16 for the rate model. The inter-frame dependency is also considered in [START_REF] Minghui | Rate-distortion optimized bit allocation for error resilient video transcoding[END_REF] where R-D models for an I-Frame and a P-Frame are presented. Nevertheless, these models do not take into account the quantization parameter (QP) of the considered video coder.

The dependent R-D models proposed in [START_REF] Lin | Bit-rate control using piecewise approximated rate-distortion characteristics[END_REF] and [START_REF] Minghui | Rate-distortion optimized bit allocation for error resilient video transcoding[END_REF] are based on experimental analysis. An analysis of inter-frame dependency based on the statistical properties of DCT coefficients is proposed by [START_REF] Wang | Theoretic analysis of inter frame dependency in video coding[END_REF]. The pdf of DCT coefficients follows a Laplace distribution and it is experimentally shown to be a function of the energy of the residual predicted frames , which is directly impacted by the reconstruction quality of the reference frame.

This paper proposes a new dependent R-D Model that offers a good prediction performance and a moderate computational complexity i.e., requiring much less R-D measurement points per picture than in [START_REF] Lin | Bit-rate control using piecewise approximated rate-distortion characteristics[END_REF]. IT focuses on the dependency between P frames and their reference I or P frames which are assumed to be unique in the video sequence. It analysis the impact of quantization parameters of the reference frame on the rate and distortion of the predicted frames.

Section II introduces the statistical multiplexing problem. Section III introduces a R-D model considering independent successive frames. It is shown to be sufficiently accurate for R-D modeling of I-Frames. Then, in section IV a dependent R-D model is introduced. Multiple regimes constituting this model are explained in section V in terms of an independent regime switching auto-regressive model with variable quantization parameters. The theoretical R-D curves obtained with this model are shown in section V to be very satisfying with a dependency behavior similar to that obtained in the experimental R-D curves of H.264 encoded frame. The concluding section VI presents some future perspectives on R-D models based on very few parameters estimation.

II. PROBLEM STATEMENT Statistical multiplexing is an efficient technique for maximizing the number of transmitted programs in broadcast systems. Figure 1 presents a broadcast system in which N video programs are encoded in parallel, the compressed bitstreams are multiplexed and transmitted over a communication channel at a constant rate R c . The video programs are encoded using a video coder, e.g., H264/AVC [START_REF]Advanced video coding for generic audiovisual services[END_REF]. In this work, only two types of frames are considered: Intra (I) and Predictive (P) frames. Apart from bandwidth constraints, the control of each encoder involved in the statistical multiplexing has to ensure that:

• Each decoded program offers at least a minimum specified quality (quality constraint),

• programs are decoded with more or less the same quality (fairness constraint),

• for each program, the quality of the reconstructed video varies smoothly (smoothness constraint). Satisfying all constraints simultaneously may sometimes be difficult.Due to the non-stationary content of each program. Variations may be due, e.g to scene changes or to high activity within a scene. In [START_REF] Changuel | Predictive control for efficient statistical multiplexing of digital video programs[END_REF], a predictive control scheme multiplexing H.264/AVC encoded video programs has been presented, taking all previously mentioned constraints into account. The control may be performed at a frame level (a distinct QP is selected for each frame), or considering several frames simultaneously (all sharing the same QP). If it is done at a frame level, for the j-th frame, [START_REF] Changuel | Predictive control for efficient statistical multiplexing of digital video programs[END_REF] introduces a control window containing the previous, current, and w -2 future frames. The quantization parameters Q j = (Q 1j , Q 2j , . . . , Q N j ) are evaluated for the j-th GOP by solving the following constrained optimization problem for the rate and quality control:

Q j , . . . , Q j+w-2 = arg min Qj ,...,Qj+w-2 N i=1 D ij (Q ij ) (1) subject to                N i=1 R ij (Q ij ) R c D ij (Q ij ) D max |P ij (Q ij ) -P j | ∆P p |P ij (Q ij ) -P ij-1 ( Q ij-1 )| ∆P s , |P ij+k (Q ij+k ) -P ij+k-1 (Q ij+k-1 )| ∆P s with k = 1 . . . w -2 and i = 1, . . . , N . (2) R ij (Q ij ) and D ij (Q ij )
denote the rate and distortion for the j-th frame in the i-th program when encoded with the quantization parameter Q ij . D max is the maximum tolerated distortion to satisfy the minimum quality constraint. The quality of the j-th frame in the i-th program is evaluated by the Peak Signal-to-Noise Ratio (PSNR)

P ij (Q ij ) = 10 log 10 255 2 D ij (Q ij ) (3)
A controlled variation ∆P S of the PSNR between successive frames of a given program can be allowed in view of the smoothness constraint. The fairness constraint of the PSNR between successive frames of a given program, whereas for the fairness constraint, the difference of PSNR between programs is bounded by ∆P p . Q i(j-1) denotes the choice of the quantization parameter for the j -1-th frame in the i-th program.

Performance of the control depends on the models for D ij (Q ij ) and R ij (Q ij ). frame j as Intra coded, if the future frames (j + 1, . . . , j + w -2) are considered independently, their corresponding R-D models will be tuned by the quantization parameters applied to each frame (Q j+1 , . . . , Q j+w-2 ). Nevertheless, for a given program, the R-D models of future frames are impacted by the choice of the quantization parameter at frame j. Therefore, R-D models should involve the quantization parameter used for the current frame and at least the one used in the corresponding reference frame. In the sequel, models for D ij and R ij are proposed as function of Q ij and the QP at the reference frame, usually

Q ij-1 .
In what follows, the program index i is omitted.

III. INDEPENDENT RATE AND DISTORTION MODEL

aim of this section is to provide a independent R-D model for I frames. Experimental R-D curves for P frames show that the previous model may not be efficiently extended for these frames. The R-D curves of a P frame depend significantly on the QP of the corresponding reference frame.

A. R-D model for I frames

This section introduces a R-D model suitable for I frames, which are encoded independently of any reference. Exponential R-D models presented in [START_REF] Changuel | Predictive control for efficient statistical multiplexing of digital video programs[END_REF] are considered

R j (Q j ) = a R exp(b R Q j ), (4)
and proposed model is tested on different frames in different video sequences. Figure 2 shows the rate and distortion functions for the first intra frame in Foreman.cif. Due to the nonprecision of the obtained models with two encoding trials, we represent it with 3 and 7 measurements while Q j vary. One sees that 3 measurements are enough to efficiently fit the R-D models.

D j (Q j ) = a D exp(b D Q j ) (5

B. Experimental R-D curves for P frames

To determine whether the models (4) and (5) may be extended to P frames, several R-D curves have been plotted as a function of Q j-1 , the QP of the reference frame, considering several values Q j , the QP of the P frame. The reference frame is an I frame here.

Figure 3 shows the logarithm of the rate as a function of Q j-1 , parametrized with Q j , for the first P-frame in Foreman.cif sequence. Figure 4 shows the logarithm of the distortion for the same frame. The R-D curves depend significantly on Q j-1 . The parameters of the models (4) and (5) would thus depend on Q j-1 .

Dependent R-D models have thus to be introduced, as done in [START_REF] Lin | Bit-rate control using piecewise approximated rate-distortion characteristics[END_REF] and [START_REF] Minghui | Rate-distortion optimized bit allocation for error resilient video transcoding[END_REF].

IV. DEPENDENT RATE AND DISTORTION MODEL

In order to propose dependent R-D models, several regimes are first identified on the experimental curves of Figure 3

j-1 Log(Rate) Q j =20 Q j =25 Q j =30 Q j =35 Q j =40
Q j =20 Q j =25 Q j =30 Q j =35 Q j =40
Fig. 4. logarithm of the distortion as a function of Q j-1 , parametrized with Q j in Foreman sequence

A. Rate regimes

All rate curves of Figure 3 decrease when Q j-1 decreases. When the quality of the reference frame increases, the motion compensation is usually more efficient, and the texture (motion compensation residuals) contains less energy, and requires thus less bits to be represented. Besides this tendency, three regimes may be identified for the rate curves.

When Q j-1 Q j , the rate remains approximatively constant. As long as the quality of the reference frame remains good, it allows to get a texture containing mainly elements which cannot be efficiently predicted from the reference frame, whatever its quality. The rate is thus mainly determined by the amount of elements which are difficult to predict in the P frame.

When Q j-1 is much larger than Q j (Q j-1 Q j + K),
the rate remains approximatively constant. In this regime, the reference frame has been coarsely quantized, and is thus of poor quality. The motion compensation is not efficient and the amount of information in the texture is not significantly affected by a poor or very poor reference frame. The rate is thus essentially determined by Q j .

The intermediate regime, when Q j-1 ∈ [Q j , Q j + K], corresponds to the situation described at the beginning of this section. Increasing the quality of the reference frame reduces the amount of information of the texture, and thus the rate.

B. Distortion regimes

All distortion curves of Figure 4, decreases when Q j-1 decreases, as for the rate curves. The interpretation is the same: having a better reference increases the efficiency of the motion-compensation, and reduces the amount of information in the texture. Even if the texture is coarsely quantized, the reconstructed frame may be of good quality thanks to the good quality of the reference frame. Now two regimes may be observed.

When Q j-1 Q j , the distortion remains approximatively constant, since it is essentially determined by Q j .

When Q j-1 Q j , the distortion decreases when Q j-1 decreases. In this regime, the motion-compensation residuals of frame j are coarsely quantized compared to frame j -1. When the motion-compensation is performed, many blocks are well predicted, and for those blocks, the motion-compensation residuals are negligible when compared to the quantization step corresponding to Q j . Many blocks have a distortion that is thus determined by Q j-1 .

The fact that the R-D curves do not converge to the same point when Q j-1 is small is due to the fact that independently of the value of Q j-1 , some block are not well predicted during motion-compensation. The corresponding texture for those blocks are thus of quite high energy. The R-D characteristics for those blocks are actually determined by Q j , and have a behavior similar to I frames.

C. Proposed dependent model

Based on the analysis of Sections IV-A and IV-B, the following piecewise linear models are proposed for

log (R j (Q j , Q j-1 )) and log (D j (Q j , Q j-1 )) log (R j (Q j , Q j-1 )) =    α r j + β r j Q j if Q j-1 Q j α r j + β r j Q j-1 if Q j Q j-1 Q j + K α r j + β r j (Q j + K) if Q j-1 Q j + K (6)
Where R j is the distortion at the predicted frame. and

log (D j (Q j , Q j-1 )) = α d j + β d j Q j-1 if Q j-1 Q j α d j + β d j Q j else (7)
Where D j is the distortion at the predicted frame. β r j and β d j are constants, and α r j and α d j correspond to the information linked to the current QP Q j . This relationship is similar to that shown in (4) and (5) for the rate and the distortion. In the logarithmic scale the relationship between α d j , α r j and Q j is so represented by a linear function. This induces to an R-D model involving a total of seven parameters.

V. ALTERNATIVE DEPENDENT MODEL

The main difficulty with the preceding models (6) and ( 7) is that they require several coding trials of the P frame and of its reference frame with several distinct values of the pair (Q j-1 , Q j ) to allow their parameters to be identified. The transform DCT coefficients can be classified into DC and AC coefficients. The DC coefficient is the mean value of the image block and carries most of the energy in the image block. The AC coefficients carry energy depending on the amount of detail in the image block. Most of the energy is compacted in the DC coefficient and a few ac coefficients. In this section, an alternative approach is investigated, considering a model in the transform domain of the dependency between the DC coefficients of blocks of the P frame and of its reference frame.

The rate due to signaling and motion vectors is not considered here. Several authors agrees on the Gaussian model of the DC coefficient in the transform domain [START_REF] Reininger | Distributions of the two-dimensional dct coefficients for images[END_REF], [START_REF] Smoot | Study of dct coefficient distributions[END_REF]. For AC coefficients, there is less consensus, e.g., [START_REF] Reininger | Distributions of the two-dimensional dct coefficients for images[END_REF] and [START_REF] Smoot | Study of dct coefficient distributions[END_REF] consider a Laplace distribution, whereas [START_REF] Kamaci | Frame bit allocation for the h.264/avc video coder via cauchy-density-based rate and distortion models[END_REF] use a Cauchy distribution.

A. Independent regime-switching autoregressive model

Here, an independent regime-switching autoregressive (IRS-AR) model

Y j = a Xj Y j-1 + b Xj U j , (8)
is considered to represent the DC value Y j of a given block of frame j and Y is the DC value of the corresponding collocated block in the reference frame j -1 (possibly after motion compensation). In (8), X j is a sequence of independent and identically distributed (iid) binary-valued random variables with Pr (X j = 0) = 1 -ρ and Pr (X j = 1) = ρ, U j is a sequence of iid zero-mean and unit variance Gaussian variables. The sequences X j and U j are assumed to be independent. This model allows to take into account the fact that most of the time (when X j = 0), the DC value of collocated blocks is quite similar, in which case a 0 is close to one, and that at some time instants (when X j = 1), there is few correlation, in which case a 1 is close to zero (a 1 = 0 in what follows).

It is assumed that the variance of Y j does not depend on j. This imposes some constraints on the parameters a 0 , b 0 , and b 1 . If X j = 0, one gets

σ 2 y = a 2 0 σ 2 y + b 2 0 = b 2 0 1 -a 2 0 (9)
and if X j = 1, one obtains

σ 2 y = b 2 1 . (10)
Combining ( 9) and (10), one gets

b 2 1 = b 2 0 1 -a 2 0 . (11)
The IRS-AR model is thus characterized by three parameters a 0 , b 0 , and ρ, which are assumed to be known in what follows. This model is only a very coarse approximation of the dependency between DC coefficients of successive frames. Nevertheless, at low bitrate, most of rate is due to the low frequency coefficients. The rate and distortion performance of scalar quantizers applied to Y j-1 and Y j is studied in the following section.

B. Rate and Distortion performance

Assume that Y j-1 is simply quantized with a scalar uniform midtread quantizer q j-1 with step size ∆ j-1 (intra coding). For Y j a predictive coding is performed with

Y j = a Xj q j-1 (Y j-1 ) (12)
as prediction for Y j (the value of X j is assumed to be known).

The prediction residual

E Xj j = Y j -Y j = a Xj (Y j-1 -q j-1 (Y j-1 )) + b Xj U j (13)
is then quantized with a stepsize ∆ j . Depending on X j , (13) may become

E 0 j = a 0 (Y j-1 -q j-1 (Y j-1 )) + b 0 U j when X j = 0 and E 1 j = b 1 U j when X j = 1.
The aim of the remainder of this section is to provide rate and distortion curves for the model (8) as a function of ∆ j-1 and ∆ j . In the H.264/AVC standard [START_REF] Itu | Advanced video coding for generic audiovisual services[END_REF], [START_REF] Richardson | H.264/MPEG-4 part 10: Transform & qantization[END_REF], characteristics of the quantizers depend on a quantization parameter Q. The relation between the quantization stepsize ∆ and Q may be approximated as

∆ (Q) = 2 Q-4 6 /P F, (14)
where P F is a constant which value depends of the subband, see [START_REF] Ma | Rate-distortion analysis for h.264/avc video coding and its application to rate control[END_REF].

The distortion for Y j may be written as

D Xj j (∆ j-1 , ∆ j ) = +∞ -∞ (x -q j (x)) 2 f E X j j (x) dx, (15)
where f E X j j (x) is the probability density function of E Xj j . The rate required to represent the quantized Y j is evaluated as the entropy of the output of the quantizer fed with

E Xj j R Xj j (∆ j-1 , ∆ j ) = - +∞ k=-∞ P k (∆ j-1 , ∆ j ) log(P k (∆ j-1 , ∆ j ))
(16) where

P k (∆ j-1 , ∆ j ) = (k+ 1 2 )∆j (k-1 2 )∆j f E X j j (x) dx. (17)
When X j = 1, E 1 j is zero-mean Gaussian with variance b 2 1 = σ 2 y and does not depend on ∆ j-1 ,

f E 1 j (x) = 1 2πσ 2 y exp - x 2 2σ 2 y . (18)
is thus known. High-rate approximations for [START_REF] Sayood | Introduction to Data Compression[END_REF] and (16) are easily obtained, see [START_REF] Sayood | Introduction to Data Compression[END_REF]. At medium to low rates (large values of ∆ j compared to b 1 such high-rate approximation becomes coarse, but [START_REF] Sayood | Introduction to Data Compression[END_REF] and (16) contain only few significant terms. When X j = 0, f E 0 j (x) is the convolution of the pdfs of a 0 (Y j-1 -q j-1 (Y j-1 )) and of b 0 U j . One may show that

f E 0 j (x) = 1 8πσ 2 y +∞ k=-∞ exp - (x + k∆ j-1 a 0 ) 2 2σ 2 y G(x, ∆ j-1 , k), (19) where
G(x, ∆ j-1 , k) = erf 2a 0 x + ∆ j-1 1 -2k 1 -a 2 0 2 √ 2σ y 1 -a 2 0 -erf 2a 0 x -∆ j-1 1 + 2k 1 -a 2 0 2 √ 2σ y 1 -a 2 0 . (20)
Now, since f E 0 j (x) and f E 1 j (x) are known, one may evaluate numerically D Xj j (∆ j-1 , ∆ j ) and R Xj j (∆ j-1 , ∆ j ) using (15) and (16). The expectation of the rate and distortion with respect to X j is then

D j (∆ j-1 , ∆ j ) = (1 -ρ) D 0 j (∆ j-1 , ∆ j ) + ρD 1 j (∆ j ) (21) R j (∆ j-1 , ∆ j ) = (1 -ρ) R 0 j (∆ j-1 , ∆ j ) + ρR 1 j (∆ j ), (22)
since D 1 j and R 1 j do not depend on ∆ j-1 .

C. Experimental results

Figures [START_REF] Lin | Bit-rate control using piecewise approximated rate-distortion characteristics[END_REF] and (6) represent log D j (Q j-1 , Q j ) and log R j (Q j-1 , Q j ), with a 0 = 0.99, b 2 0 = 3, and ρ = 0.1, meaning that 10% of DC coefficients may not be efficiently predicted from the previous frame. The conversion between ∆ and Q is done using 14. The chosen quantization parameters are Q j-1 ∈ [15, . . . , 50] for the reference Y j-1 and Q j ∈ {20, 25, 30, 35, 40} for Y j .

Theoretical R-D curves shown in Figures [START_REF] Lin | Bit-rate control using piecewise approximated rate-distortion characteristics[END_REF] and (6) have the same behavior as the R-D curves obtained using H.264/AVC video coder in Figures [START_REF] Xie | Sequence-based rate control for constant quality video[END_REF] and (4). The same regimes are observed. However, the two sets of rate curves do not have the same values. A constant element is missing in the theoretical curves to have the same values as in the experimental curves. This is due to the rate of the motion compensation residues which is not considered in the theoretical model. This later has an almost constant evolution while varying the quantization parameter, which perfectly explains the missing part in the total rate in the theoretical curves. Theoretical and experimental distortion curves are very close. In this case, all distortion elements are considered and the proposed model is able to accurately reproduce the distortion values from the H.264/AVC video coder.

Q j =20 Q j =25 Q j =30 Q j =35 Q j =40
Q j =40 Q j =35 Q j =30 Q j =25 Q j =20
Fig. 6. Distortion dependency between current and reference DCT coefficient Thus, the model (8), defined by three parameters is a viable alternative to the much more complicated models (6) and (7) involving seven parameters.

The estimation of the parameters of the model (8) may be done with several coding trials for different values of the pair (Q j-1 , Q j ), as would be done for the models (6) and [START_REF] Wang | Theoretic analysis of inter frame dependency in video coding[END_REF]. Nevertheless, since the parameters of the model proposed in this section have a physical meaning, one expects to be able to perform their identification without performing several coding trials.

VI. CONCLUSION

In this work, we propose a new rate and distortion models for the statistical multiplexing system using H.264/AVC video encoders. These models take into account the interframe dependencies used for performing motion compensation between a reference and a predicted frame even I and P frames or P and P frames. In this paper, the dependencies specified in terms of impact of encoding parameters of a given frame over the next one is proved experimentally and interpreted theoretically by considering the quantization of a independent regime-switching autoregressive source model with varying quantization parameters. Theoretical R-D curves obtained for this source model show a behavior which is very similar to the curves obtained using H.264/AVC video. The proposed theoretical models is based on three parameters a 0 , b 0 , and ρ. Further work, will be dedicated to the estimation of these parameters based on several coding trials for different values of the pair (Q j-1 , Q j ).

Fig. 1 .

 1 Fig. 1. Structure of a statistical multiplexer

  )where j denotes the index of the current frame and Q j is the QP used to encoded it. At least two encoding trials per frame are necessary to estimate a D , b D , a R and b R . The

Fig. 2 .

 2 Fig. 2. Rate and distortion model fitting when 3 and 7 experimental values are considered (Foreman)

  the refrence frame: Q

Fig. 3 .

 3 Fig.3. logarithm of the rate as a function of Q j-1 , parametrized with Q j in Foreman sequence

1

  Quantization parameter of the reference DCT: Q j-1 Log(Rate)

Fig. 5 .6

 5 Fig. 5. Rate dependency between current and reference DCT coefficient

ACKNOWLEDGMENT

The authors are deeply grateful to V.Kumar for his helpful comments and continual encouragement.


Céline Elie-Caille
email: celine.caille@clipproteomic.fr

Jean-Yves Rauch

Alain Rouleau

Wilfrid Boireau

Preparation of flat gold terraces for protein chip developments

A simply method to prepare flat gold terraces on mica for AFM biomolecular characterization is described. The procedure includes preheating of the substrate, a metal deposition followed by an annealing step, all of these steps lead at elevated temperature (300 to 420°C). This approach allowed to prepare large flat gold terraces (200 to 500 nm), that constituted ideal substrates for visualization and characterization of self-assembly monolayer of biomolecules at the nanoscale. We illustrated this potential of characterization with the reconstitution of a protein monolayer.

Introduction

Repartition and orientation of immobilized proteins on bio-devices are very important to assure their high performance. There have been few methods to evaluate surface coverage, molecules repartition and even tridimensional structure and orientation of immobilized proteins. One of them is time-of-flight secondary ion mass spectrometry (TOF-SIMS) that is able to analyze upper surface of one layer of molecules 1,2 . Atomic force microscopy is one of the methods of choice to achieve the protein surface coverage, the biofilm homogeneity with a molecular resolution until protein structure, without any "denaturating" sample preparation or labeling 3 .

Gold surfaces provide convenient supports for a number of reasons. Gold is rather chemically inert to biomacromolecules; but, on the other hand, it is accessible to chemisorption via the formation of stable gold-thiol bonds. Based on that, a protein having a unique cystein, i.e. a unique sulfhydryl residue at the periphery of the 3D structure, will naturally self-assemble on gold substrate. Bearing this sulfhydryl moiety, cytochrome b5 is able to perform direct selfassembly onto gold surface by chemisorption process leading to a packed and oriented self assembled monolayer (SAM) 2 .

The classical way to prepare thin gold films is vacuum evaporation of gold onto glass, silicon or mica, at low pressure. This deposition method induces metal films presenting globular gold grains measuring 25 to 30 nm in diameter 4 . Such rough surfaces present limitations in resolution when gold immobilized biomolecules have to be imaged and characterized.

Atomically flat gold substrates seem ideal to better characterize surface bound biomolecules using scanning microscopy methods. Mica substrate is often preferred as it provides a clean and highly flat surface. When gold is deposited on such an atomically flat surface, the metal film can be stripped from mica and used on this "mica-like face" 5 . The other possibility is to thermally anneal the gold covered mica since this method allows improving the flatness of the gold film. Nevertheless the first "stripping" method often needs the deposition of another layer on gold, in order to propose a new "support" to the gold stripped film. This method is moreover delicate to realize due to the stripping step of metal film in the case of a mechanical pealing. Stripping can be also obtained through a chemical treatment using THF. This last possibility presents also limitations especially in terms of surface contamination by mica flasks 6 .

We propose here a new method to prepare gold films presenting large flat terraces on mica.

Our method is rapid, highly controlled in terms of pressure and temperature and limitation of surface contamination, since the whole process happens under a non-ruptured vacuum (substrate preheating, metal deposition and annealing). We present the relevancy of such a flat gold surface in the characterization of a protein monolayer specifically self-assembled on it.

Material and methods

Materials

Pure 99,99% gold metal and mica sheets (G250-1, from Agar Scientific) were used. Engineered cytochrome b5 was prepared as described in the reconstitution procedure. Octylglucopyranoside (OG) and dithiothreitol (DTT) were purchased from Sigma-Aldrich.

Metal deposition and annealing run

Gold thin films were deposited on mica substrates to prepare flat terraces, by electron beam evaporation in an Alliance Concept EVA450 vacuum system at a pressure of 4* x 10 -6 mbar. Indeed, we wanted to obtain very thin grain size, in order to mimic the flat state of the mica's surface while realizing flat gold terraces. Before gold deposition run, the substrates were first heated at a stabilized temperature of 300°C for half an hour and then etched with argon ion beam powered by a hallow cathode of 3cm diameter at a pressure of 3,6x10 -4 mbar. Positive argon ions were neutralized by flow rates of negative species before impact of the substrates. An Argon flow rate was 5 sccm in cathode and neutralizer, and the pressure in the chamber was obtained by a high-speed turbo molecular pump Leybold of 1200 l.min -1 coupled with a primary rotative pump. Finally, evaporation was obtained by heating the pure 99,99% gold metal in a liner of graphite. In a first time, deposition was made on the shadow between liner and substrates until the right and regular speed rate control by an Inficon XTC2 quartz crystal was reached.

Gold was deposited either in one or two steps. For sample I, deposition of gold film has been made in two steps: first, 150 nm was deposited at a rate of 5 nm.sec -1 which will be reduced to 0.05 nm.sec -1 to achieve the thickness of 200 nm. For samples II, III and IV, the gold films were obtained in one step at low gold deposition rate.

For all samples, heating, etching and gold layer depositions were all realized at 300°C. Samples I and II stayed in the same chamber for a run of 30 min and 8h of annealing at 300°C. Samples III, IV and V were annealed in another vacuum system at 2.5*x10 -6 mbar for 3h, 20h or 32h at 420°C respectively. Sample VI was obtained after a classical gold deposition of 150 nm (without heating) followed by an annealing step of 20h at 420°C.

Reconstitution of cytochrome b5 monolayer

Engineered cytochrome b5 was derived from human microsomal cytochrome b5 by genetic engineering resulting in the substitution of (i) the 26 C-terminal amino-acid residues by the -NGHHHH-COOH sequence and (ii) the serine 23 in Hb5(His) 4 by cystein as previously described 7 .

Modified

AFM characterization

The AFM used was a Nanoscope III (Veeco, Santa Barbara, CA). Imaging was performed in contact mode using NPS-oxide sharpened silicon nitride probes (Veeco). For the feedback controls, typical values of set-point for imaging were between 0.5 to 1.5 V, depending on scan size.

Image J software

We used Image J software to determine the area of gold terraces obtained on different substrates.

Results & discussion

We chose deposition of gold using electron beam evaporation method, since it offers possibility to work under elevated temperature until 300°C before, during and after the gold deposition. This According to the same Thornton diagram, it also means that with sputtering conditions the gold layer grew in columnar structure. The evaporation way with low deposition rate is the best one to form gold flat terraces. In the same field, the evaporation of gold with high speed was not chosen, because the roughness of the coating increased with the thickness and the flat terraces disappeared. The only way we found to better mimic the flatness of the mica substrate and obtain flat gold terraces, was evaporation beam with a temperature of 300°C and a pressure less than 5*x10 -5 mbar.

It appears that long bake-outs of the mica helped to outgas contamination which prevents direct stacking of Au atoms onto the mica surface 9 . We decided to use etching with neutral argon ion gun, in order to clean the surface without high particles energies and to avoid disturbance of the mica substrate.

Morphology and size of terraces

We have evaporated gold on pre-heated mica substrate and have tested different conditions of metal deposition and temperature. Different gold deposition rate, thickness and temperature conditions were experimented (Table 1). Indeed, the rate of metal deposition influences the roughness of the gold film.

A classical evaporation without heating induces a rough gold film consisting of grains with a diameter around 30 nm 4 . While decreasing the gold deposition rate the surface topography gets smoother 5 .

Thermal annealing of the evaporated films dramatically improves the flatness of the gold substrates, producing atomically flat terraces 10,11 .

For sample I, the gold thickness is 200 nm and only a short annealing time at 300°C (30 min) is applied. Small terraces with irregular shape were obtained (Fig1A,B). The sample II was obtained after 50 nm gold deposition followed by a long annealing time at 300°C (8h). In this case, terraces present round shaped morphology (Fig1C,D). In the case of sample III and IV, 150 nm of gold were deposited but annealing time was 3h or 20h at 420°C. In this last case, we obtained better results.

Indeed, after 3h annealing, terraces are larger than sample I and present an homogeneous shape, (Fig1E,F). After 20h annealing, the surface still presents these individual terraces, but several terraces seemed have fused together to give very large ones (Fig1G,H).

Concerning the size, we managed to prepare terraces from 22 500 nm 2 for sample I (Fig2A,D), to 55 000 nm 2 for sample III (Fig2B,D) and until more than 100 000 nm 2 (Fig2C,D) for sample IV. Thus, the tendency for the terraces to grow in size with annealing time appeared clearly (Fig2D).

Nevertheless, there is a maximum heating duration beyond which the gold surface is completely disorganized like "melt" on/with mica (sample V and Fig3A). AFM images do not present gold terraces anymore and the surface topography is a sort of intermediary between classical rough gold film and highly flat mica. Nogues and Wanunu (2004) 12 also reported this phenomena, saying that overexposure (long time) of gold films to high temperatures had to be avoided. Indeed, when mica was heated to temperatures above 500 °C, Derose et al 9 remarked that the surface deteriorated becoming quite rough. It appears that the mica disintegrated at these temperatures causing poor epitaxy.

We also made two observations: first, there is a minimum gold thickness allowing preparation of relatively large terraces. Indeed, when only a deposition of 50 nm of gold is realized at 300°C (sample II), even followed by an annealing time of 8h at 300°C, small round and "donuts" shape terraces were obtained (Fig1C,D). Second, if the gold deposition is realized without heating the substrate, but followed by an annealing time at 420°C during several hours (sample VI), we did not succeed in gold terrace realization. The gold substrate presents in this case a rough surface presenting gold grains of around 30 nm in diameter (fig3B), that corresponds to the roughness of "classical" gold films.

Our work present then an easy way to prepare flat gold surfaces under highly controlled conditions of pressure and temperature. Moreover, this approach, using equipments in a "clean room", presents the advantage to protect the gold interface against contamination by polluting species in air. This last point is of deep importance when the substrate has to be efficiently and homogeneously functionalized by proteins for a biosensor development.

Terraces flatness

Higher resolution images revealed nicely individual gold terraces. These terraces are smooth on 100 (for the smallest) to 200 nm in length or more and present a maximum height variation of only around 5Å (Fig4). The morphology exhibited by such gold surfaces is therefore ideal for clear visualization, distinction and characterization of immobilized molecules and proteins.

Visualization and counting of immobilized protein

We developed flat surfaces in order to better characterize the protein monolayer immobilized at a gold surface. Indeed, classical evaporated gold films present a high roughness, with a surface composed of small grains in the same size range than proteins. The topography of such a rough surface prevents a clear visualization and characterization of a protein layer grafted on it. In our aim to characterize the reconstituted biomolecular film at the interface of sensor surface, the development of flat surface could allow to count individual molecules and determine their surface coverage, eventually their tridimensional structure and orientation. As a model, we worked with a protein containing a single cystein residue allowing its binding to gold. In Fig5 we show clearly the difference between before (Fig5A) and after (Fig5B) incubation of proteins on the gold terraces. Binding of proteins induced surface changes and appearance of lot of motifs on the gold surface. As a comparison, a classical evaporated gold film, without annealing, is represented before (Fig5C) and after (Fig5D) protein adsorption on the surface. Proteins are undistinguishable showing the difficulty to characterize the protein layer on this kind of rough surface. In opposite, on surface presenting terraces, proteins are visible, distinguishable and it is possible to control the density, homogeneity, to count proteins and determine their surface coverage. We counted around 1200 proteins per µm Images also revealed that the protein layer is homogeneous on gold surface, with proteins regularly distributed on the substrate. This AFM investigation revealed that neither particular concentration of molecules on sides of terraces nor aggregation of proteins appeared on the surface.

Conclusion

A new rapid, reproducible and highly controlled procedure for the preparation of clean flat gold terraces was presented. We demonstrated the usefulness of such flat surfaces for reconstitution of a protein monolayer and its easily characterization by AFM in terms of surface coverage, density and repartition of immobilized proteins. In order to use such flat gold surfaces in a biosensor, we are currently working on gold terraces formation on glass. With collaborators working on theoretical modeling of surface behavior, we are studying the ability of these flat surfaces to generate plasmons.

This last development would open the way of these innovative atomically gold substrate investigation through their plasmon interrogation, giving informations on kinetics and mass loading. Thus, these flat gold surfaces could represent a really interesting flat and "active" substrate for biochip developments.

Figure legends

Fig1: Gold terraces on mica obtained with different conditions of metal deposition rate, thickness and temperature. For each sample, a large scan (around 1µm) and one zoom (around 500 nm) are presented. (A,B) Sample I obtained with 200 nm gold deposited in two-steps: 150 nm at 5 nm/s then 50 nm at 0.05 nm/s, followed by 30 min 300°C annealing. (C,D) Sample II obtained with 50 nm gold at 0,1 nm/s, followed by 8h 300°C

annealing. Sample III (E,F) obtained with 150 nm gold deposited at 0,1nm/s followed by 3h annealing at 420°C.

Sample IV (G,H) obtained with 150 nm gold deposited at 0,1nm/s followed by 20h at 420°C. Scale bar: 100 nm.

  temperature of 300°C corresponds to a transition structure of densely packed fibrous grains in Thornton diagram 8 enabling to keep the state of mica's surface and to obtain the flat gold terraces. On this diagram, it is also possible to understand that evaporation with lower pressure and little energies of particles -i.e. less than 1eV-corresponds to a better solution than sputtering with a pressure of 5*x10 -3 mbar of argon and particles energies from 5 to 10 eV. According to Thornton diagram, it is important to choose evaporation conditions with lower pressure and a temperature T/Tm in the range of 0,2 to 0,3. T corresponds to the temperature of the condensation vapor, i.e. the temperature of the substrate, and Tm is the gold liquid temperature, i.e. 1063°C.

Fig2:

  Fig2: Gold terrace areas of (A) sample I, (B) sample III and (C) sample IV. (D) Area of terraces function of annealing time. Scale bar: 100 nm.

Fig3:

  Fig3: Effect of heating during gold deposition and annealing time. (A) Gold surface after an annealing time of 32h at 420°C (preheating of substrate at 300°C, 150 nm gold thickness at 300°C). (B) Gold deposition without heating the substrate (150 nm gold thickness followed by annealing at 420°C 20h). Scale bar: 100 nm.

Fig4:Fig5:

  Fig4: Flatness of individual gold terraces. (A) Gold terraces prepared with 50 nm gold evaporated at 0,1nm/s, followed by 8h 300°C annealing (sample II). (B) Cross-section along the white line in A. (C) Gold terraces prepared with 150 nm gold evaporated at 0,1nm/s followed by 3h annealing (sample III). (D) Cross-section along the white line in B. (E) Gold terraces prepared with 150 nm gold at 0,1nm/s followed by 20h at 420°C (sample IV). (F) Cross-section along the white line in E. Scale bar: 50 nm.

Table 1 :

 1 Various conditions of gold deposition and annealing runs.

	Sample	Heating parameters	Etching at 300°C	Gold thickness and	Annealing
	number		(min)	deposition rate at 300°C	parameters
	I			150 nm at 5 nmsec -1	30 min at 300°C
				50 nm at 0.05 nmsec -1
	II	30 min at 300°C	1min at 3.6.10 -4 mbar	50 nm at 0.1 nmsec -1	8h at 300°C
	III				3h at 420°C
	IV				20h at 420°C
	V			150 nm at 0.1 nmsec -1	32h at 420°C
	VI	/	/		20h at 420°C

  human cytochrome b5 bears a unique cystein, i.e. a unique sulfhydryl residue, and is

	processed in saline phosphate buffer (PBS, 100 mM at pH 7.4 with NaCl 50 mM). Protein is reduced
	by a reducing agent in excess (dithiothreitol, DTT) (20/1 by moles) during 10 min. Reduced

cytochrome b5 (20 µl at 1 µM) is able to react directly by chemisorption onto gold substrate. During 15 min incubation, the reconstitution of a packed and oriented protein monolayer by direct selfassembly is obtained. A washing step by buffer first (20 µl 3 times) and secondly by a detergent (Octylglucopyranoside (OG), 20 µl at 40 mM) was realized after protein immobilization in order to remove not covalently bond proteins.

  2 , meaning 1.2 10 9 molecules/mm 2 thus 2 fmol/mm 2 . Aoyagi S. et al (2008) 2 mentioned that it is possible to perform SAMs with various density in a window of 1 to 200 fmol/mm 2 . Our AFM results then fit in this range.

Acknowledgments

We thank the centrale de technologie Mimento from FEMTO-ST institute in Besançon to have allowed realization of these different gold substrates.


Romain Couillet
email: romain.couillet@supelec.fr

Mérouane Debbah
email: merouane.debbah@supelec.fr

Jack W Silverstein

Asymptotic Capacity of Multi-User MIMO Communications

This paper introduces two new formulas to derive explicit capacity expressions of a class of communication schemes, which include single-cell multi-user MIMO and singleuser MIMO with multi-cell interference. The extension of a classical theorem from Silverstein allows us to assume a channel Kronecker model between the base stations and the cellular terminals, provided that they all embed a large number of antennas. As an introductory example, we study the singleuser MIMO setting with multi-cell interference, in the downlink. We provide new asymptotic capacity formulas when single-user decoding of the incoming data or MMSE decoding are used. Simulations are shown to corroborate the theoretical claims, even when the number of transmit/receive antennas is not very large.

I. INTRODUCTION

In the last decade, while mobile networks were expected to run out of power and frequency resources, Foschini [START_REF] Foschini | On Limits of Wireless Communications in a Fading Environment when Using Multiple Antennas[END_REF] and Telatar [START_REF] Telatar | Capacity of multi-antenna Gaussian channels[END_REF] introduced the notion of MIMO (multiple input multiple output) systems and predicted a growth of capacity performance of min(n R , n T ) times the single antenna capacity for an n T -antenna transmitter and an n R -antenna receiver. However, this tremendous multiplexing gain can only be achieved for large SINR (signal-to-interference plus noise ratio) and without signal correlation at the channel ends. In case of correlation due to antenna spacing or to poorly scattering environments, the theoretically achievable gains are still not completely available. Moreover, in present multi-cell wireless mobile networks, neither base stations nor users cooperate; this leaves the device manufacturers with the dilemma of increasing the signal processing capabilities of the transmit/receive units to result into non significant throughput gains when adjacent cells interfere one another. Besides, due to limited computational constraints, suboptimal linear techniques such as MMSE (minimum mean square error) decoding are used at the receiver [START_REF] Paulraj | Introduction to Space-Time Wireless Communications[END_REF], in place of optimal single-user decoders.

In this work, we derive the per-antenna channel capacity of MMSE receivers against optimal single-user decoders in multi-cell networks, when the number of antennas at the transmitters and receivers is large. The capacity here is defined as the supremum of the achievable rates between a base station and a specific user (in uplink or in downlink) interfered by other cells. We model all transmission channels by the well-spread Kronecker model [START_REF] Chuah | Capacity Scaling in MIMO Wireless Systems under Correlated Fading[END_REF]. Few major contributions propose to study the capacity performance of point-to-point communications with interference. In [START_REF] Shamai | Information-Theoretic Considerations for Symmetric, Cellular, Multiple-Access Fading Channels -Part I[END_REF], the authors carry out the performance analysis of TDMA-based networks with intercell interference. In [START_REF] Zaidel | Multicell Uplink Spectral Efficiency of Coded DS-CDMA With Random Signatures[END_REF], a random matrix approach is used to study large CDMA-based networks with inter-cell interference. In the MIMO context, [START_REF] Moustakas | MIMO Capacity Through Correlated Channels in the Presence of Correlated Interferers and Noise: A (Not So) Large N Analysis[END_REF] provides an analytic solution to our problem in two-cell networks, using replica methods [START_REF] Guo | Multiuser Detection and Statistical Physics[END_REF]. These methods are however tedious since they require heavy combinatorial calculus. We propose in the following a more direct approach, based on analytical tools of random matrix theory. Note also that [START_REF] Peacock | Eigenvalue distributions of sums and products of large random matrices via incremental matrix expansions[END_REF] provides a deterministic capacity expression in our context, when the receive covariance matrices R k are assumed to share the same eigenspace 1 .

Although this specific work is dedicated to the study of point-to-point MIMO systems with multi-cell interference, the potential applications to the mathematical results we introduce cover a larger class of problems, in which channel capacities express as the log determinant of a sum of Gram matrices X i X H i , where X i is a large matrix modelled as Kronecker. For example, beside the uplink/downlink multi-cell singleuser MIMO, our results are applicable to single-cell multiuser MIMO communications in the uplink, evaluation of the capacity region of multiple access channels and broadcast channels [START_REF] Couillet | A deterministic equivalent approach for the capacity analysis of multi-user MIMO channels[END_REF].

The remainder of this work is structured as follows: in Section II, we provide mathematical preliminaries and we introduce a new important theorem, for which we provide a sketch of proof. In Section III, we introduce the system model. In Section IV, the point-to-point capacity of the channel between a base station and a user, interfered by other cells, is derived when single-user decoding or MMSE decoding are performed at the receiver. In Section V, we provide simulation results of the previously derived theoretical formulas. Finally, in Section VI, we give our conclusions.

Notation: In the following, boldface lower-case symbols represent vectors, capital boldface characters denote matrices (I N is the N × N identity matrix). X ij denotes the (i, j) entry of X. The Hermitian transpose is denoted (•) H . The operators tr X, |X| and X represent the trace, determinant and spectral norm of matrix X, respectively. The symbol x n a.s.

-→ x denotes almost sure convergence of x n to x. The notation F Y stands for the empirical distribution of the eigenvalues of the Hermitian matrix Y.

II. MATHEMATICAL PRELIMINARIES

Part of this work is dedicated to the introduction of a novel theorem, from which the multi-cell downlink and uplink capacities will be given compact expressions. This theorem generalizes Silverstein and Bai's formula [START_REF] Silverstein | On the empirical distribution of eigenvalues of a class of large dimensional random matrices[END_REF] to multiple random matrices with separable variance profiles, i.e. following the Kronecker model, and unfolds as follows, Theorem 1: Let K, N ∈ N be some positive integers. Let

B N = K k=1 R 1 2 k X k T k X H k R 1 2 k (1)
be an N × N matrix with the following hypothesis for all k ∈ {1, . . . , K},

1) T k is n k × n k Hermitian nonnegative definite, n k ∈ N * , qith unit normalized trace, 2) R 1 2 k is the N × N Hermitian nonnegative definite square root of the nonnegative definite matrix R k with unit normalized trace, 3) X k is N × n k with i.i.d. complex Gaussian entries with variance 1/n k . For k ∈ {1, . . . , K}, let c k = n k /N . Also denote, for z ∈ C \ R + , m N (z) = 1 N (B N -zI N ) -1 .
Then, as all n k and N grow large (while

K is fixed), with ratio c k m N (z) -m (0) N (z) a.s. -→ 0 (2)
where

m (0) N (z) = 1 N tr K k=1 τ k dF T k (τ k ) 1 + τ k c k e k (z) R k -zI N -1 (3)
and the set of functions {e i (z)}, i ∈ {1, . . . , K}, form the unique solution to the K equations

e i (z) = 1 N tr R i K k=1 τ k dF T k (τ k ) 1 + τ k c k e k (z) R k -zI N -1 (4) such that sgn(ℑ[e i (z)]) = sgn(ℑ[z]).
The function m N (z) is the Stieltjes transform of the random variable with cumulative distribution function F BN . The complete proof of a more general expression of this theorem is given in an extended version of the present article [START_REF] Couillet | A deterministic equivalent approach for the capacity analysis of multi-user MIMO channels[END_REF].

Remark 1: This theorem allows us to derive Stieltjes transforms of large matrices independently of the realization of the X k matrices. In wireless communications, this allows one to characterize the performances of a multi-user or multicell communication based only on the transmit and receive correlations R k and T k . This further helps to estimate channel capacity thanks to the Shannon transform, Theorem 2: Let B N be a random Hermitian matrix as defined in Theorem 1 with the additional assumption that, for

all N , n k , tr(T k ) = n k , tr(R k ) = N , and let x > 0. Then, for large N , n k , V(x) -V (0) (x) a.s.
-→ 0, where

V(x) = log 1 + b x dF BN (b) (5)
and

V (0) (x) = 1 N log det I N + 1 x K k=1 R k τ k dF T k (τ k ) 1 + c k e k (-x)τ k + K k=1 1 c k log (1 + c k e k (-x)τ k ) dF T k (τ k ) + x • m (0) N (-x) -1 (6)
A proof of this result (with less stringent hypothesis) is provided in [START_REF] Couillet | A deterministic equivalent approach for the capacity analysis of multi-user MIMO channels[END_REF]. From both results, one can derive the capacity formulation of a large range of multi-user/multi-cell network models. In the following, we cast attention on single-cell networks with multi-cell interference.

III. SYSTEM MODEL

In this work we derive capacity expressions of wireless channels between a multi-antenna transmitter and a multiantenna receiver, the latter of which is interfered by several multi-antenna transmitters. This scheme is well-suited to multi-cell wireless networks with orthogonal intra-cell and interfering inter-cell transmissions, both in downlink and in uplink. The following scenarios encompass in particular

• multi-cell uplink: the base station of a cell indexed by i ∈ {1, . . . , K} receives data from one user in this cell 2and is interfered by K -1 users transmitting on the same physical resource from remote cells indexed by j ∈ {1, . . . , K}, j = i.

• multi-cell downlink: the user being allocated a given time/frequency resource in a cell indexed by i ∈ {1, . . . , K} receives data from its dedicated base-station and is interfered by K -1 base stations in neighboring cells indexed by j ∈ {1, . . . , K}, j = i. In the following, in order not to confuse both scenarios, only the downlink scheme is considered. However, one must keep in mind that the provided results can easily be adapted to the uplink case.

Consider a wireless mobile network with K ≥ 1 cells indexed from 1 to K, controlled by non-physically connected base stations. On a particular time/frequency resource, each base station serves only one user; therefore the base station and the user of cell j will also be indexed by j. Without loss of generality, we focus our attention on user 1, equipped with n R ≫ K antennas and hereafter referred to as the user or the receiver. Every base station j ∈ {1, . . . , K} is equipped with n Tj ≫ K antennas. We additionally denote c j = n Tj /n R .

Denote s j ∈ C nT j , E[s j s H j ] = I nT j , the signal transmitted by user j, y ∈ C nR the signal received by the base station and n ∼ CN(0, σ 2 I nR ) the noise vector received at the base station. The fading MIMO channel between base station j and the user is denoted H j ∈ C nR×nT j . Moreover we assume that H j has a separable variance profile, i.e. can be decomposed as

H j = R 1 2 j X j T 1 2 j ( 7
)
with R j ∈ C nR×nR the (Hermitian) correlation matrix at the receiver with respect to the channel H j , T j ∈ C nT j ×nT j the correlation matrix at transmitter j and X j ∈ C nR×nT j a random matrix with Gaussian independent entries of variance 1/n Tj .

Remark 2: Note that in this model, and contrary to what is often assumed, R j , the correlation matrix at the receiver, explicitly depends on j. In the uplink scenario, this assumption is of particular relevance in the sense that base stations are usually placed in areas clear of scatterers. In these circumstances, the solid angle from which the signals from user j originate influences the signal correlation at the receive antenna array. Hence the dependence of the receive correlation matrices on j. Note moreover that, in this model, the transmit power assumption E[s j s H j ] = I nT j is not restrictive in the sense that the transmit power correlation of user j can be included into the matrix T j . However, the Kronecker model has two major drawbacks: (i) the inner matrix X j implicitly assumes a high density of scatterers3 in the communication link and (ii) the correlations on both sides must be inter-independent and independent of the realizations of X j , which is inaccurate to some extent.

With the assumptions above, the communication model unfolds

y = H 1 s 1 + K j=2 H j s j + n (8)
where s 1 is the useful signal (from base station 1) and s j , j ≥ 2, constitute interfering signals.

IV. MULTI-CELL MIMO CAPACITY

A. Optimal Decoding

If the receiving user considers the signals from the K -1 interfering transmitters as pure noise and knows the exact value of the SNR (signal-to-noise ratio) σ -2 , then base station 1 can transmit with arbitrarily low decoding error at a perreceive antenna rate C SU (σ 2 ) given by

C SU (σ 2 ) = 1 n R log 2 |I nR + 1 σ 2 K j=1 H j H H j | - 1 n R log 2 |I nR + 1 σ 2 K j=2 H j H H j | (9)
Assume that n R and the n Ti , i ∈ {1, . . . , K}, are large compared to K and such that no eigenvalue of R i or T i is too large. As in Theorem 1, we define the function m i,(0) as the approximated Stieltjes transforms of K j=i H j H H j , i ∈ {1, 2},

m i,(0) (z) = 1 n R tr K k=i t k dF T k (t k ) 1 + t k c k e i k (z) R j -zI nR -1 (10)
where, for all i ∈ {1, . . . , K}, e i j (z) is solution of the fixedpoint equation

e i j (z) = 1 n R tr R j K k=1 t k dF T k (t k ) 1 + t k c k e i k (z) R k -zI -1 (11)
From Theorem 2, we then have approximately

C SU (σ 2 ) = 1 n R log det I + 1 σ 2 K k=1 R k τ k dF T k (τ k ) 1 + c k e 1 k (-σ 2 )τ k - 1 n R log det I + 1 σ 2 K k=2 R k τ k dF T k (τ k ) 1 + c k e 2 k (-σ 2 )τ k + K k=1 1 c k log 1 + c k e 1 k (-σ 2 )τ k dF T k (τ k ) - K k=2 1 c k log 1 + c k e 2 k (-σ 2 )τ k dF T k (τ k ) + σ 2 • [m 1,(0) (-σ 2 ) -m 2,(0) (-σ 2 )] (12)
However, this capacity expression assumes no specific power allocation at the transmitter. If the transmit covariance matrices T j are replaced by T 1 2 j P j T 1 2 j with P j the signal covariance matrix at transmitter j, we remark that (9) can be rewritten as

C SU (σ 2 ) = 1 n R log |I+ 1 σ 2 A -1 2 R 1 2 1 X 1 T 1 2 1 P 1 T 1 2 1 X H 1 R 1 2 1 A -1 2 | (13) with A = I nR + 1 σ 2 j>1 H j P j H H j .
From [START_REF] Tulino | Impact of antenna correlation on the capacity of multiantenna channels[END_REF], the optimal power allocation for base station 1 consists in aligning the eigenvectors of P 1 to those of T 1 and in choosing the eigenvalues p 1 , . . . , p nT 1 as

   p i = 0 (1/α i ) -1 ≤ 1 nT 1 nT 1 l=1 (1 -α l ) p i = 1-αi 1 n T 1 Pn T 1 l=1 (1-α l ) otherwise (14) for α -1 i equal to 1 + 1 σ 2 h H i A -1 2 I N + 1 σ 2 A -1 2 H -i P 1 H H -i A -1 2 -1 A -1 2 h i (15) h i is the i th column of H 1 , and
H -i is H 1 with column i removed.
Denoting

h i = R 1 2
1 x i , we have x i centered Gaussian with covariance T 1ii /n T1 I nR and independent of

R 1 2 1 A -1 2 I N + 1 σ 2 A -1 2 H -i P 1 H H -i A -1 2 -1 A -1 2 R 1 2
1 .

Therefore, asymptotically on n R , from lemmas in [START_REF] Silverstein | On the empirical distribution of eigenvalues of a class of large dimensional random matrices[END_REF],

α i = 1 + T 1ii σ 2 n R tr R 1 A + 1 σ 2 H -i P 1 H H -i -1 -1 (16) = 1 + T 1ii σ 2 n R tr R 1 A + 1 σ 2 H 1 P 1 H H 1 -1 -1 (17) =   1 + T 1ii σ 2 n R tr R 1   I + 1 σ 2 K j=1 H j P j H H j   -1    -1 (18)
= 1 + T 1ii e 1 1 (-σ 2 ) -1 (19)
This leads to the power allocation

     p i = 0, T 1ii e 1 1 (-σ 2 ) ≤ 1 nT 1 nT 1 l=1 1 -1 + T 1 ll e 1 1 (-σ 2 ) -1 p i = 1-(1+T1 ii e 1 1 (-σ 2 )) -1 1 n T 1 Pn T 1 l=1 1-(1+T1 ll e 1 1 (-σ 2 ))
-1 otherwise (20) In most cases however, the diagonal entries of T 1 all equal 1, and therefore the optimal power allocation policy is uniform array power allocation.

B. MMSE Decoder

Achieving C SU requires non-linear processing at the receiver, such as MMSE successive interference cancellation. A suboptimal linear technique, the MMSE decoder, is often used instead. The communication model in this case reads

y =   k j=1 H j H H j + σ 2 I nR   -1 H H 1   k j=1 H j s j + n   (21)
and each entry of y will be processed individually. This technique makes it possible to transmit data reliably at any rate inferior to the per-antenna MMSE capacity C MMSE ,

C MMSE (σ 2 ) = 1 n R nT 1 i=1 log 2 (1 + γ i ) (22)
where, denoting h j ∈ C nT j the j th column of H 1 and R 1 2

1 x j = h j , the SINR γ i expresses as

γ i = h H i K j=1 H j H H j + σ 2 I nR -1 h i 1 -h H i K j=1 H j H H j + σ 2 I nR -1 h i (23) = h H i   K j=1 H j H H j -h i h H i + σ 2 I nR   -1 h i (24) = x H i R 1 2 i   K j=1 H j H H j -h i h H i + σ 2 I nR   -1 R 1 2 i x i (25)
where Equation (24) comes from a direct application of the matrix inversion lemma. With these notations, x i has i.i.d. complex Gaussian entries with variance T 1ii /n Ti and the inner matrix of the right-hand side of (25) is independent of x i (since the entries of H 1 H H 1h i h H i are independent of the entries h i ). Applying Lemma 3.1 in [START_REF] Silverstein | On the empirical distribution of eigenvalues of a class of large dimensional random matrices[END_REF], for n Ti large, approximately

γ i = T 1ii n T1 tr R 1   K j=1 H j H H j -h i h H i + σ 2 I nR   -1 (26)
From Lemma 2.1 in [START_REF] Bai | On the signal-to-interference-ratio of CDMA systems in wireless communications[END_REF], the rank 1 perturbation (-h i h H i ) does not affect asymptotically the trace in (26). And therefore, approximately,

γ i = T 1ii n T1 tr R 1   K j=1 H j H H j + σ 2 I nR   -1 (27)
Noting that e 1 (z) in Section IV-A corresponds to the normalized trace in Equation ( 27) (this is shown precisely in the proof of Theorem 1 [START_REF] Couillet | A deterministic equivalent approach for the capacity analysis of multi-user MIMO channels[END_REF]), we finally have the compact expression for C MMSE ,

C MMSE (σ 2 ) = 1 n R nT 1 i=1 log 2 1 + 1 c 1 T 1ii e 1 (-σ 2 ) (28)
In practice, when no power allocation strategy is applied, T 1ii = P the average power per transmit symbol, and the capacity becomes

C MMSE = c 1 • log 2 (1 + P c1 e 1 (-σ 2 )).

V. SIMULATION AND RESULTS

In the following, we apply the results ( 12) and (28) to the downlink of a two-cell network. The capacity analyzed here is the achievable rate on the link between base station 1 and the user, the latter of which is interfered by base station 2. The relative power of the signal received from user 2 is on average Γ times that of user 1. Both base stations 1 and 2 are equipped with linear arrays of n T antennas and the user with a linear array of n R antennas. The correlation matrices T i at the transmission and R i at the reception, i ∈ {1, 2}, are modeled thanks to a generalization of Jake's model including solid angles of transmit/receive power, i.e. for instance,

T i ab = θ (i) max θ (i) min exp 2π • i • d Ti ab λ cos(θ) dθ (29)
with d Ti ab the distances between antennas indexed by a, b ∈ {1, . . . , n Ti } for transmitter i, (θ (i) min , θ (i) max ) the angles over which useful power (i.e. power that will be received by user 1) is transmitted, and λ the wavelength.

In Figure 1, we took n R = 16, Γ = 0.25 and we consider optimal single-user decoding at the receiver. For every realization of T i , R i , 1000 channel realizations are processed to produce the simulated ergodic capacity and compared to the theoretical capacity (28). Those capacities are then averaged over 100 realizations of T i , R i , varying in the random choice of θ In Figure 2, with the same assumptions as previously, we apply MMSE decoding at the base station. Here, a slight difference is observed in the high SNR regime between theory and practice. This was somehow expected, since the large n R approximations in Silverstein's lemmas [START_REF] Silverstein | On the empirical distribution of eigenvalues of a class of large dimensional random matrices[END_REF] are very loose for σ 2 close to R -in the sense of the Euclidean distance. To cope with this gap, many more antennas must be used. We also observe a significant difference in performance between optimum and linear MMSE decoders, especially in the high SNR region. Therefore, in wireless networks, when interfering cells are treated as Gaussian correlated noise at the cell-edge, i.e. where the interference is maximum, the MMSE decoder provides tremendous performance loss.

VI. CONCLUSION

In this paper, we introduced an important theorem relating the Stieltjes transform of a class of large matrices to a deterministic approximate. Based on this formula, we provided compact expressions for the optimal and MMSE-decoder capacities of point-to-point MIMO with inter-cell interference, for random channel matrices with separable variance profile, both in downlink and in uplink. The simulations show perfect match with the theoretical formulas in the low-SNR region, even if fewer antennas are used at the transmitters and receivers. As for the high SNR region, a large number of antennas must be used to reach an accurate match between theory and Monte-Carlo simulations.

  max with constraint θ (i) max -θ (i) min = π/2, while d Ti ab = 10λ|a -b| at the transmitters, d R ab = 2λ|a -b| at the receiver. The SNR ranges from -5 dB to 30 dB, and n T ∈ {8, 16}. We observe here that Monte-Carlo simulations perfectly match the capacity obtained from Equation[START_REF] Tulino | Impact of antenna correlation on the capacity of multiantenna channels[END_REF].

Fig. 1 .

 1 Fig. 1. Capacity of point-to-point MIMO in two-cell uplink, optimal singleuser decoding, n R = 16, n T ∈ {8, 16}, Γ = 25%.

  VII. ACKNOWLEDGMENT Silverstein's work is supported by the U.S. Army Research Office under Grant W911NF-05-1-0244. Debbah's work is partially supported by the European Commission in the framework of the FP7 Network of Excellence in Wireless Communications NEWCOM++.

Fig. 2 .

 2 Fig. 2. Capacity of point-to-point MIMO in two-cell uplink, with MMSE decoder, n R = 16, n T ∈ {8, 16}, Γ = 25%.

this assumption is too strong for our current study since it boils down to supposing that all entities have the same geometrical antenna pattern and that transmit/receive energy comes from the same solid angle for all users.

this user is allocated a given time/frequency resource, which is orthogonal to time/frequency resources of the other users in the cell; e.g. the multi-access protocol is OFDMA.

the number and distance between scatterers must be of the same order as the number and distance between the transmit and receive antennas.


Allou Samé

Latifa Oukhellou

Etienne Côme

Patrice Aknin
email: aknin@inrets.fr

Mixture-model-based signal denoising

Keywords: Denoising, Asymmetrical noise, Regression, Gaussian mixture model, EM algorithm, GEM algorithm

This paper proposes a new signal denoising methodology for dealing with asymmetrical noises. The adopted strategy is based on a regression model where the noise is supposed to be additive and distributed following a mixture of Gaussian densities. The parameters estimation is performed using a Generalized EM (GEM) algorithm. Experimental studies on simulated and real signals in the context of a diagnosis application in the railway domain reveal that the proposed approach performs better than the least-squares and wavelets methods.

Introduction

Denoising or noise reduction is a subject of research in both signal and image processing. The aim is to remove the additive noise from the signal without distorting it. Various linear and non linear methods have been proposed to solve this problem. They lead to different solutions depending on the noise model and the signal properties which are closely linked to applications.

One of the most popular and recent denoising approaches is based on wavelet analysis [START_REF] Mallat | A theory for multiresolution signal decomposition: the wavelet representation[END_REF][START_REF] Mallat | A wavelet tour of signal processing[END_REF]. It consists in applying thresholding algorithms to the wavelet coefficients in order to extract the most significant of them. This method allows an unknown function to be reconstructed from noisy data, but its theoretical framework involves a normal noise density and requires an optimal choice of the threshold. For this purpose, several thresholding rules have been proposed and successfully applied (Dohono et al. (1996), [START_REF] Ogden | Data dependent wavelet thresholding in nonparametric regression with change-point applications[END_REF].

Many statistical approaches like least-squares method have also been proposed. In this paper, we explore a strategy for denoising based on mixture models which are widely used in clustering (Dempster et al. 1977, McLachlan and[START_REF] Mclachlan | The EM Algorithm and Extentions[END_REF]. The basic idea of the proposed approach is to fit a regression model on the signal with the assumption that the error term (or additive noise) follows a mixture of normal distributions. The parameters of such a model are then estimated by means of a Generalized Expectation Maximization (GEM) algorithm.

It should be noticed that this approach is different from the wavelet method, which is dedicated to reconstruct an unknown signal in a white Gaussian noise. Here, we assume a parametric regression model of the signal and a mixture of normal distributions for the noise. In many situations where the normal assumption cannot be satisfied, the normal mixtures seem to be suitable because they can represent a wide variety of noise densities including both symmetrical and non-symmetrical cases [START_REF] Marron | Exact mean integrated squared error[END_REF].

The paper is organized as follows. At first, we present the regression model and introduce some preliminary notations. Then, we detail the parameter estimation procedure via a GEM algorithm. Section four gives the performances of the approach illustrated on simulated examples and section five is concerned with an application in the railway domain.

Additive noise and mixture model

A signal is supposed to be represented by an independent sample ((x 1 , y 1 ), . . . , (x n , y n )) where variables x and y, defined on IR, represent the dependent variable (for example the time) and the independent variable (the signal at time x), respectively. Polynomial regression functions have been chosen, for their simplicity and their coherence with many real signals. An observed signal is then modeled by

y = a T x + ε, (1)
where a = (a 0 , a 1 , a 2 , . . . , a m ) T , x = (1, x, x 2 , . . . , x m ) T , m is the polynomial order and ε is an additive noise whose distribution does not depend on x.

Usually, ε is assumed to be normally distributed with zero mean and variance σ 2 , and the estimation of the model parameters proceeds by solving a classical least-squares problem, but in certain practical problems the presence of other sources of noise results in an asymmetrical non-normal noise which makes the least-squares method unsuitable.

To solve this problem, the noise ε is assumed to be centered and distributed according to a mixture of K normal densities

f (ε) = K k=1 π k N (ε; µ k , σ 2 k ) (2)
where N (•; µ k , σ 2 k ) is the normal distribution density with mean µ k and variance σ 2 k , the π k ′ s are the proportions of the mixture. The proportions π k and the means µ k satisfy the constraints

K k=1 π k = 1 and K k=1 π k µ k = 0. The component density N (•; µ k , σ 2
k ) can thus be regarded as one source of noise. The well known flexibility of mixture distributions allows the regression model to handle both symmetrical and non-symmetrical noises [START_REF] Bartolucci | The use of mixtures for dealing with nonnormal regression errors[END_REF].

The next section shows how the parameters of the proposed model can be estimated.

Parameters estimation via the GEM algorithm

The model parameters are estimated by maximizing the conditional likelihood or equivalently the conditional log-likelihood defined by

L(Ψ) = n i=1 log p(y i |x i ; Ψ), (3)
where Ψ = (a, π 1 , . . . , π K , µ 1 , . . . , µ K , σ 2 1 , . . . , σ 2 K ) is the parameter vector to be estimated. It can be easily verified that, conditionally on x, the variable y has the mixture distribution density

p(y|x; Ψ) = K k=1 π k N (y ; a T x + µ k , σ 2 k ).
The (observed) log-likelihood is then written

L(Ψ) = n i=1 log K k=1 π k N y i ; a T x i + µ k , σ 2 k , (4)
where

x i = (1, x i , x 2 i , . . . , x m i ).
This equation clearly highlights that the proposed model is simply a constrained Gaussian mixture of regressions [START_REF] Aitkin | Mixture models, outliers and the EM algorithm[END_REF]Wilson 1980, DeSarbo and[START_REF] Desarbo | A maximum likelihood methodology for clusterwise linear regression[END_REF] in which the coefficient vector a is common for each component.

The Expectation-Maximization (EM) algorithm [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF]) is used to perform the maximization of the log-likelihood (4) which can not be solved directly. Let us recall that the EM algorithm requires the specification of complete data whose log-likelihood can be maximized more easily than the observed log-likelihood (4). In this situation, the complete data are obtained by adding to each observation (x i , y i ) its component membership variable z i ∈ {1, . . . , K}. Using the binary coding of z i , z ik = 1 if z i = k and z ik = 0 otherwise, the complete data log-likelihood is written as

L C (Ψ) = n i=1 K k=1 z ik log π k N y i ; a T x i + µ k , σ 2 k , (5)
Given an initial parameter vector Ψ (0) , the EM algorithm consists in alternating the two following steps until convergence.

E step : Expectation

This step consists in evaluating the expectation of the complete data loglikelihood conditionally on the observed data and the current parameter vector Ψ (q) :

Q(Ψ; Ψ (q) ) = E L C (Ψ)|x 1 , . . . , x n , y 1 , . . . , z n ; Ψ (q) = n i=1 K k=1 t (q) ik log π k N y i ; a T x i + µ k , σ 2 k (6)
where L C (Ψ) is the complete data log-likelihood and

t (q) ik = E(z ik |x i , y i ; Ψ (q) ) = P (z ik = 1|x i , y i ; Ψ (q) ) = π (q) k N (y i ; a (q) T x i + µ (q) k , σ (q) k 2 ) K ℓ=1 π (q) ℓ N (y i ; a (q) T x i + µ (q) ℓ , σ (q) ℓ 2 ) (7)
is the posterior probability that y i originates from the kth mixture component, given x i . Thus, the E step simply requires the computation of the posterior probabilities t (q)

ik .

M step : Maximization

This step consists in computing the parameter vector Ψ (q+1) that maximizes with respect to Ψ, the quantity Q(Ψ; Ψ (q) ) which, in our situation, can be written

Q(Ψ; Ψ (q) ) = k,i t (q) ik log π k - 1 2 k,i t (q) ik log(2π) + log σ 2 k + (y i -(a T x i + µ k )) 2 σ 2 k . ( 8
)
As in the classical Gaussian mixture situation, it can be verified that the proportions maximizing Q are given by

π (q+1) k = n i=1 t (q) ik n . (9)
The simultaneous maximization with respect to the parameters a, (µ k ) and (σ 2 k ) cannot be performed analytically. But it can be shown that, for fixed parameters (σ 2 k ), the optimum parameter vector [a T , µ 1 , . . . , µ K ] is given by

[a T , µ 1 , . . . , µ K ] T = K k=1 H T k H k -1 K k=1 H T k Y k , (10)
with

H k = 1 σ k • diag t (q) 1k , . . . , t (q) nk • X k Y k = 1 σ k • diag t (q) 1k , . . . , t (q) nk • Y,
where diag(e 1 , . . . , e n ) is the diagonal matrix with diagonal elements e 1 , . . . , e n , and

X k =         x T x T n 0 . . . 0 1 0 . . . 0         and Y =         y 1 . . . y i . . . y n         .
Notice that X k is a n × (K + m + 1) matrix with all elements equal 0 except those in the columns 1, . . . , m + 1 and k + m + 1.

On the other hand, for fixed parameters a and (µ k ), it can be shown that the parameter σ 2 k maximizing Q is given by

σ 2 k = n i=1 t (q) ik y i -(a T x i + µ k ) 2 n i=1 t (q) ik . ( 11
)
Thus, we propose to compute the parameter vector Ψ (q+1) as follows:

(a) compute a (q+1) , µ The strategy described by (a) and (b), which fulfil Q(Ψ (q+1) ; Ψ (q) ) ≥ Q(Ψ (q) ; Ψ (q) ), is commonly used when the maximization cannot be directly performed and also known as a conditional-M step [START_REF] Mclachlan | The EM Algorithm and Extentions[END_REF]. It provides a Generalized EM algorithm (GEM) [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF] whose convergence properties are the same as those of the EM algorithm. Usually, the GEM algorithm is iterated for a number of random initial parameters Ψ (0) and only the parameter that provides the largest log-likelihood value is adopted.

Experiments on simulated signals

This section is designed to evaluate the performance of the proposed GEM algorithm as regards the quality of estimation, by using simulated signals.

For this purpose, results obtained with GEM are compared with those of the and wavelets methods. Let us remind that least-squares and wavelet denoising methods [START_REF] Donoho | Density estimation by wavelet thresholding[END_REF] are based on the assumption that the noise has a zero mean normal distribution.

Simulation protocol

Each simulated signal is generated according to a polynomial regression function g(x) (the denoised signal that is to be recovered) and an asymmetrically distributed noise. For all simulated signals, we consider a sample size n = 600 and the value x i = i for the dependant variable (i = 1, ..., n). In the following, we describe the setting of the polynomial coefficients, the type of asymmetrical noise distribution, the criteria used to measure the estimation quality, and the performed experiments.

Polynomial coefficients

focus here only on second degree polynomial functions as the denoised signal, in accordance with the real measured signals described in section 5. Since no significant influence of the regression coefficients on the quality of estimation has been observed, only one polynomial function with coefficients vector a = (1.2; 3 • 10 -5 ; -5 • 10 -5 ) is considered in all the described simulations below, i.e., g(x) = 1.2 + 3 • 10 -5 x -5 • 10 -5 x 2 . Relatively small values (-5•10 -5 and 3•10 -5 ) of the polynomial coefficients have been chosen to get simulated signals whose amplitude is very close to the real signals amplitude.

Asymmetrical noise distributions

Two different types of asymmetric distributions have been adopted for the noise generation:

• a two-components Gaussian mixture distribution with density

f (ε; π 1 , π 2 , µ 1 , µ 2 , σ 2 1 , σ 2 2 ) = π 1 N (ε; µ 1 , σ 2 1 ) + π 2 N (ε; µ 2 , σ 2 2 ) (12)
with mean 2 k=1 π k µ k = 0;

• a three-parameters Weibull distribution

f (ε; κ, λ, α) = κ λ ε -α λ κ-1 exp - ε -α λ κ ]α;+∞[ (ε), (13)
with α = -λΓ(1+ 1 κ ), Γ being the gamma function. This particular choice of the parameter α ensures a zero mean distribution.

Figures 1 and2 show examples of signals simulated according to the adopted polynomial regression and the two noise distributions (Gaussian mixture and Weibull).

Evaluation criteria

The quality of estimation provided by GEM, least-squares and wavelets is jointly evaluated by the following criteria measuring the quality of the estimated (denoised) signal and the quality of the noise estimation, respectively:

• The L 2 distance between the true regression function g(x i ) with parameter a and the estimated regression function g(x i ) with parameter a, defined by

D 2 (g, g) = 1 n n i=1 [a T x i -a T x i ] 2 ; ( 14
KL(f, f ) = IR f (x) log f (x) f (x) dx. ( 15
)
The Kullback-Leibler divergence is computed using Monte Carlo integration method since its evaluation is not analytically straightforward in the case of Gaussian mixtures.

Performed experiments

Four sets of experiments have been performed for this simulation study.

• The first set of experiments studies the quality of the estimates when the ratio σ 2 1 /σ 2 2 of the Gaussian mixture variances varies in the set {1/5, 1/4, 1/3, 1/2, 1}, for equal proportions π 1 = π 2 = 0.5, variance σ 2 1 = 2, and ∆ = |µ 1µ 2 | = 4.

• The second set studies the quality of the estimates when the ratio π 1 /π 2 of the Gaussian mixture proportions varies in the set {1/3, 1/2, 1, 2, 3}, for equal variances σ 2 1 = σ 2 2 = 4, and ∆ = 4.

• The third set studies the quality of the estimates when the distance ∆ between the means of the Gaussian mixture components varies in the set {0, 5, 10, 16, 21}, for proportions π 1 = 0.7, π 2 = 0.3 and variances σ 2 1 = 2, σ 2 2 = 10.

• The fourth set of experiments studies the quality of the estimates when the parameter λ of the Weibull distribution varies in the set {5, 10, 15, 20, 25, 30}, for a parameter κ = 1.4.

Results

Table 1 represents, for the four sets of experiments, the D 2 distance and the KL divergence obtained with the methods GEM (applied to K = 2 and K = 3 components), least-squares and wavelets. Notice that each value of D 2 and KL is an average over 30 different simulated samples (Monte-Carlo simulations). As it can be seen, GEM applied with K = 2 or K = 3 performs globally better than its competitors, for any modification of noise distribution parameters σ 2 1 /σ 2 2 , π 1 /π 2 , ∆ or α, in terms of the D 2 distance. Not surprisingly, the KL divergence for GEM is distinctly better for any situation; this behavior illustrates that the mixture model is more suitable for dealing with various situation of asymmetrical noise distribution. When the distance ∆ or the Weibull parameter λ increases, it can be observed that the D 2 distance for all methods increases. In practice, this phenomenon, even more pronounced for least-squares and wavelets, is due to an increase of the noise in addition to the noise asymmetry. It can also be observed that the number of components K = 3 is more adapted for Weibull distribution.

Application to defect diagnosis in the railway domain

In the railway domain, signal denoising is an important preprocessing task for defect diagnosis. Most of the automatic diagnosis systems are indeed based on the analysis of measurement signals which are often corrupted with noise.

Our application concerns the diagnosis of the track/vehicle transmission system in the French high speed lines [START_REF] Aknin | Track circuit diagnosis by automatic analysis of inspection car measurements[END_REF]. The identification of two sources of noise led us to consider that the noise is distributed according to a mixture of two normal densities, rather than a single normal density.

The first task consists in segmenting the complete signal into local arches with a piecewise approach, and the second task is focused on the denoising procedure on each arch. Figure 3 shows an example of a measurement signal. The evaluation of the proposed approach performances on these experimental signals is not possible because the original signals and noise distributions are not available. For this reason, we choose to test the proposed denoising algorithm on signals whose original denoised signals are given by an electrical model, which allows to achieve realistic simulations of the system [START_REF] Aknin | Track circuit diagnosis by automatic analysis of inspection car measurements[END_REF]. We use two different denoised signals called DS1 and DS2, composed of n = 778 obervations (see figure 4). The same noise distributions (Gaussian mixture and Weibull) as for the simulation study are added to DS1 and DS2. The quality of denoised signals reconstruction is evaluated using the D 2 distance. Unlike the cases studied in section 4, the signals considered in this part are generated with non polynomial regression Table 1: D 2 distance and KL divergence in relation to σ 2 1 /σ 2 2 , π 1 /π 2 , ∆ and parameter λ of the Weibull distribution, obtained with GEM (applied with K = 2 and K = 3), least-squares and wavelets. functions but we try to approximate them with second degree polynomial functions.

σ 2 1 σ 2 2 D 2 KL GEM GEM LS Wav. GEM GEM LS Wav. (K=2) (K=3) (K=2) (K=3) 1/5 0.
Results obtained for DS1 and DS2 using the same noise distribution parameters as in the simulation study are shown in tables 2 and 3. It will be observed that, once again, our proposed GEM algorithm performs better than least-squares and wavelets. The relative behavior of the compared methods are almost the same as for the simulation study.

In rare situations where variances have small values, wavelets denoising has better performances than GEM and least-squares. This remark can be illustrated, for instance, in the case of the denoised signal DS2, by setting σ 2 1 = 5 • 10 -4 , π 1 = π 2 = 0.5, ∆ = 0.05 and choosing values of σ 2 2 corresponding to different ratios σ 2 1 /σ 2 2 (see Table 4). The good performances of wavelets in this case, can be attributed to the non polynomial nature of the denoised signal DS2.

Conclusion

An original methodology for signal denoising was proposed in this paper. This methodology is based on a polynomial regression model where the noise is supposed to be additive and distributed following a Gaussian mixture distribution. A Generalized EM (GEM) algorithm is developed to estimate both the noise distribution and the polynomial regression coefficients. The experimental study on simulated signals reveals that the GEM algorithm provides reasonably good results for relatively low degree of the regression function, compared to the standard least-squares and wavelets methods. When the noise is distributed according to certain asymmetric distributions like Weibull, the experimental study shows that the number of component K = 3 is more adapted than K = 2. This point suggests to adopt a suitable strategy for choosing the optimal number of components of the mixture as regards of the likelihood criterion and the complexity of the model. For this purpose the adaptation of Bayesian Information Criteria (BIC) [START_REF] Schwarz | Estimating the number of components in a finite mixture model[END_REF]) could be an interesting direction of work. Final results on physical signals in the context of a real application in the railway domain show that the GEM algorithm can represent an efficient mean for denoising real signals in presence of asymmetrical noise.

  11), where the parameters a (q+1) and µ (q+1) k have been inserted in the right hand side.

Figure 1 :

 1 Figure 1: (left) Example of a signal simulated according to a second degree polynomial function and a Gaussian mixture distribution, and its true denoised signal; (right) corresponding noise mixture distribution and its component densities

Figure 2 :

 2 Figure 2: (left) Example of a signal simulated according to a second degree polynomial function and a Weibull distribution, and its true denoised signal; (right) corresponding noise Weibull distribution

Figure 3 :Figure 4 :

 34 Figure 3: Experimental signal recorded on a track/vehicle transmission system

Figure 5 :

 5 Figure 5: Denoised signal corresponding to the experimental signal of figure 3.

Table 2 :

 2 D 2 distance in relation to σ 2 1 /σ 2 2 , π 1 /π 2 , ∆ and parameter λ of the Weibull distribution; GEM, least-squares and wavelets are applied to noised signals based on DS1.

	σ 2 1 σ 2 2	GEM	D 2 GEM	LS	Wav.	π 1 π 2	GEM	D2 GEM	LS	Wav.
		(K=2) (K=3)				(K=2)	(K=3)
	1/5	0.025	0.026	0.040	0.137	1/3	0.024	0.025	0.032	0.127
	1/4	0.019	0.019	0.031	0.115	1/2	0.025	0.028	0.031	0.144
	1/3	0.019	0.019	0.026	0.091	1	0.027	0.028	0.031	0.165
	1/2	0.023	0.024	0.027	0.104	2	0.029	0.031	0.034	0.136
	1	0.014	0.015	0.016	0.067	3	0.039	0.039	0.043	0.133
	∆		D 2			λ		D2
		GEM	GEM	LS	Wav.		GEM	GEM	LS	Wav.
		(K=2) (K=3)				(K=2)	(K=3)
	0	0.012	0.012	0.015	0.055	4	0.012	0.018	0.029	0.112
	5	0.018	0.017	0.031	0.098	6	0.046	0.039	0.069	0.239
	10	0.030	0.030	0.054	0.175	8	0.069	0.063	0.118	0.422
	16	0.031	0.031	0.091	0.318	10	0.122	0.093	0.151	0.601
	21	0.043	0.043	0.125	0.512	12	0.153	0.126	0.207	0.793

Table 3 :

 3 D 2 distance in relation to σ 2 1 /σ 2 2 , π 1 /π 2 ,∆ and parameter λ of the Weibull distribution; GEM, least-squares and wavelets are applied to noised signals based on DS2.

	σ 2 1 σ 2 2	GEM	D 2 GEM	LS	Wav.	π 1 π 2	GEM	D2 GEM	LS	Wav.
		(K=2) (K=3)				(K=2)	(K=3)
	1/5	0.021	0.021	0.028	0.109	1/3	0.027	0.028	0.029	0.118
	1/4	0.022	0.022	0.035	0.159	1/2	0.029	0.030	0.036	0.166
	1/3	0.021	0.024	0.027	0.121	1	0.038	0.045	0.042	0.150
	1/2	0.019	0.019	0.022	0.090	2	0.023	0.023	0.025	0.150
	1	0.015	0.016	0.016	0.086	3	0.031	0.032	0.039	0.171
	∆		D 2			λ		D2
		GEM	GEM	LS	Wav.		GEM	GEM	LS	Wav.
		(K=2) (K=3)				(K=2)	(K=3)
	0	0.014	0.015	0.014	0.083	4	0.022	0.018	0.031	0.115
	5	0.017	0.018	0.024	0.086	6	0.040	0.032	0.065	0.285
	10	0.026	0.026	0.049	0.207	8	0.053	0.047	0.076	0.458
	16	0.025	0.025	0.064	0.339	10	0.109	0.091	0.159	0.645
	21	0.044	0.044	0.133	0.591	12	0.171	0.141	0.206	1.043

Table 4 :

 4 Results obtained for DS2+noise in the case of small values of the mixture component variances : σ 2 1 = 5.10 -4 .

	σ 2 1 σ 2 2			D 2 (×10 -4 )
		GEM	GEM	LS Wav.
		(K=2) (K=3)
	1/5	0.56		0.57	0.57 0.39
	1/4	0.56		0.56	0.56 0.30
	1/3	0.56		0.56	0.56 0.30
	1/2	0.54		0.54	0.54 0.22
	1	0.54		0.54	0.54 0.16
	3
	2
	y (Ampere)
	1
	0	5000	10000	15000	20000	25000
				x (Samples)

0 . . . 0 1 0 . . . 0 . . . . . . . . . . . . . . . . . . . . . . . . x T i 0 . . . 0 1 0 . . . 0 . . . . . . . . . . . . . . . . . . . . . . . .