N Demni

M Zani

LARGE DEVIATIONS FOR STATISTICS OF JACOBI PROCESS

Keywords: Jacobi Process, Subordinated Jacobi process, Large deviations, Maximum likelihood AMS Classification: primary 60G05, 60F10, 33C45; secondary: 62F12

This paper is aimed to derive large deviations for statistics of Jacobi process already conjectured by M. Zani in her thesis. To proceed, we write in a simpler way the Jacobi semi-group density. Being given by a bilinear sum involving Jacobi polynomials, it differs from Hermite and Laguerre cases by the quadratic form of its eigenvalues. Our attempt relies on subordinating the process using a suitable random time-change. This will give an analogue of Mehler formula whence we can recover the desired expression by inverting some Laplace transforms. Once we did, an adaptation of Zani's result ([24]) in the non steep case will provide the required large deviations principle.

Introduction

The Jacobi process is a Markov process on [-1, 1] given by the following infinitesimal generator:

L = (1 -x 2 ) ∂ 2 ∂ 2 x + (px + q) ∂ ∂x , x ∈ [-1, 1]
for some real p, q, defined up to the first time when it hits the boundary. It appears as an interest rate model in finance (see [START_REF] Delbaen | An interest rate model with upper and lower bound[END_REF]) and in genetics ( [START_REF] Ethier | Markov Processes : Characterization and Convergence[END_REF]). One of the important features is that it belongs to the class of diffusions associated to some families of orthogonal polynomials, i.e. the infinitesmal generator admits an orthogonal polynomials basis as eigenfunctions ( [START_REF] Bakry | Characterization of Markov Semi-groups on R Associated to Some Families of Orthogonal Polynomials[END_REF]) such as Hermite, Laguerre and Jacobi polynomials . More precisely, if P α,β n denotes the Jacobi polynomial with parameters α, β > -1 defined by :

P α,β n (x) = (α + 1) n n! 2 F 1 -n, n + α + β + 1, α + 1; 1 -x 2 , x ∈ [-1, 1],
(see [START_REF] Lebedev | Special Functions And Their Applications[END_REF] for the definition of 2 F 1 ) then we can see that :

L P α,β n = -n(n + α + β + 1)P α,β n for p = -(β + α + 2) and q = β -α. The semi group density of the process first appeared in [START_REF] Karlin | Classical diffusion processes and total positivity[END_REF] then in [START_REF] Wong | The construction of a class of stationnary Markov[END_REF] where the author solved the forward Kolmogorov or Fokker-Planck equation ∂ 2 y [B(y)p] -∂ y [A(y)p] = ∂ t p, p = p t (x, y), where B, A are polynomials of degree 2, 1 respectively, and gave the principal solution (p 0 (x, y) = δ x (y)) using the classical Sturm-Liouville theory. This gives rise to a class of stationary Markov processes satisfying : [START_REF] Applebaum | Quantum Independent Increment Processes I : From Classical Probability to Quantum Stochastic Calculus[END_REF] lim t→∞ p t (x, y) = W (y) =

x 2

x 1 W (x)p t (x, y)dx

where W is the density function solution of the corresponding Pearson equation ( [START_REF] Wong | The construction of a class of stationnary Markov[END_REF]).

In our case, p t has the discrete spectral decomposition :

(2)

p t (x, y) =   n≥0 (R n ) -1 e -λnt P α,β n (x)P α,β n (y)   W (y), x, y ∈ [-1, 1]
where

λ n = n(n + α + β + 1), W (y) = (1 -y) α (1 + y) β 2 α+β+1 B(α + 1, β + 1)
with B denoting the Beta function and3 ([2], p. 99) :

R n = ||P α,β n || 2 L 2 ([-1,1],W (y)dy) = Γ(α + β + 2) 2n + α + β + 1 (α + 1) n (β + 1) n Γ(α + β + n + 1)n! .
Interested in total positivity, Karlin and McGregor showed that this kernel is positive for α, β > -1 ([14]). Few years later, Gasper ([13]) showed that, for α, β ≥ -1/2, this bilinear sum is the transition kernel of a diffusion and a solution of the heat equation governed by a Jacobi operator, generalizing a previous result of Bochner for ultraspherical polynomials ( [START_REF] Bochner | Sturm-Liouville and heat equations whose eigenfunctions are ultraspherical polynomials or associated Bessel functions[END_REF]). It is worth noting that λ n has a quadratic form while in the Hermite (Brownian) and Laguerre (squared Bessel) cases λ n = n. Hence, we will try to subordinate the Jacobi process by the mean of a random time-change in order to get a Mehler type formula. What is quite interesting is that the subordinated Jacobi process semi-group, say q t (x, y), is the Laplace transform of p 2/t (x, y). Thus, we deduce an expression for p t (x, y) by inverting some Laplace transforms already computed by Biane, Pitman and Yor (see [START_REF] Biane | Probability Laws Related To The Jacobi Theta and Riemann Zeta Functions, and Brownian Excursions[END_REF], [START_REF] Pitman | Infinitely Divisible Laws Associated With Hyperbolic Functions[END_REF]). This expression, more handable than (2), will allow us to compute the normalized cumulant generating function, and then to derive a LDP for the maximum likelihood estimate (MLE) for p in the ultraspherical case, i.e. q = 0 (β = α), a fact conjectured by Zani in her thesis ( [START_REF] Zani | Grandes déviations pour des fonctionnelles issues de la statistique des processus[END_REF]). Then, using a skew product representation of the Jacobi process involving squared Bessel processes, we construct a family {ν t } t of estimators for the index ν of the squared Bessel process based on a Jacobi trajectory observed till time t. This satisfies a LDP with the same rate function derived for the MLE based on a squared Bessel trajectory.

1.1. Inverse Gaussian subordinator. By an inverse Gaussian subordinator (see [START_REF] Applebaum | Quantum Independent Increment Processes I : From Classical Probability to Quantum Stochastic Calculus[END_REF]), we mean the process of the first hitting times of a Brownian motion with drift B µ s := B s + µs, µ > 0, namely,

T µ,δ t = inf{s > 0; B µ s = δt}, t, δ > 0.
Using martingale methods, we can show that for each t > 0, u ≥ 0,

E(e -uT µ,δ t ) = e -tδ( √ 2u+µ 2 -µ)
whence the density f t of T µ,δ t writes ( [START_REF] Applebaum | Quantum Independent Increment Processes I : From Classical Probability to Quantum Stochastic Calculus[END_REF]) :

(3)

f t (s) = δt √ 2π e δtµ s -3/2 exp - 1 2 ( t 2 δ 2 s + µ 2 s) 1 {s>0} .
1.2. The subordinated Jacobi Process. Let us consider a Jacobi process (X t ) t≥0 .

Then the semi-group of the subordinated Jacobi process (X T µ,δ t ) t≥0 is given by:

q t (x, y) = ∞ 0 p s (x, y)f t (s)ds = W (y) n≥0 (R n ) -1 ∞ 0 e -λns f t (s)ds P α,β n (x)P α,β n (y) = W (y) n≥0 (R n ) -1 E(e -λnT µ,δ t )P α,β n (x)P α,β n (y).
Writing

λ n = (n + γ) 2 -γ 2 where γ = α + β + 1 2
, and substituting δ = 1/ √ 2, µ = √ 2γ for α + β > -1 in the expression of f t , one gets :

E(e -λnT µ,δ t ) = e -nt so that q t (x, y) = W (y) n≥0 (R n ) -1 e -nt P α,β n (x)P α,β n (y).

The last sum has been already computed ( [START_REF] Andrews | Special functions[END_REF], p. 385) :

∞ n=0 (R n ) -1 P α,β n (x)P α,β n (y)r n = 1 -r (1 + r) a m,n≥0 a 2 m+n a+1 2 m+n (α + 1) m (β + 1) n u m v n m!n! = 1 -r (1 + r) a F 4 ( a 2 , a + 1 2 , α + 1, β + 1; u, v) (4)
where |r| < 1, a = α + β + 2, F 4 is the Appell function ( [START_REF] Exton | Multiple Hypergeometric Functions And Applications[END_REF]) and

u = (1 -x)(1 -y)r (1 + r) 2 v = (1 + x)(1 + y)r (1 + r) 2 .
The integral representation of F 4 (see [START_REF] Exton | Multiple Hypergeometric Functions And Applications[END_REF], p 51) yields :

q t (x, y) = W (y) Γ(a) 1 -r (1 + r) a ∞ 0 s a-1 e -s 0 F 1 (α + 1; u 4 s 2 ) 0 F 1 (β + 1; v 4 s 2 )ds.
Now, from a property of the function 0 F 1 (see [START_REF] Magnus | Formulas And Theorems for the Special Functions of Mathematical Physics[END_REF], p 214)

0 F 1 (c; w(1 -r)/2) 0 F 1 (d; w(1 + r)/2) = n≥0 P α,β n (r) (c) n (d) n w n , α = c -1, β = d -1,
one gets :

q t (x, y) = W (y) Γ(a) 1 -r (1 + r) a ∞ 0 s a-1 e -s n≥0 P α,β n (z) (α + 1) n (β + 1) n A n s 2n ds
where we set

z = x + y 1 + xy , A = (1 + xy)r 2(1 + r) 2
. Applying Fubini's Theorem gives :

q t (x, y) = W (y) 1 -r (1 + r) a n≥0 (a) 2n (α + 1) n (β + 1) n P α,β n (z)A n .
Letting r = e -t , then

q t (x, y) = W (y)e a-1 2 t 2 a-1 sinh(t/2) (cosh(t/2)) a n≥0 (a) 2n (α + 1) n (β + 1) n P α,β n (z) (1 + xy) 8 cosh 2 (t/2) n = W (y) tanh(t/2)e a-1 2 t 2 a-1 n≥0 (a) 2n (α + 1) n (β + 1) n P α,β n (z) (1 + xy) 8 n 1 cosh(t/2) 2n+a-1 .
Besides, from (3)

q t (x, y) = t e γt 2 √ π ∞ 0 p s (x, y) s -3/2 e -γ 2 s e -t 2 4s ds = t e γt 2 √ 2π ∞ 0 p 2/r (x, y) r -1/2 e -2γ 2
/r e -t 2 8 r dr.

Thus, noting that γ = (a -1)/2, we get :

∞ 0 p 2/r (x, y) r -1/2 e -2γ 2 /r e -t 2 8 r dr = √ 2π W (y) 2 a-1 tanh(t/2) t/2 n≥0 (a) 2n (α + 1) n (β + 1) n P α,β n (z) (1 + xy) 8 n 1 cosh(t/2) 2n+a-1 .
With regard to the integrand, one easily sees that the RHS is the Laplace transform of p 2/r (x, y)r -1/2 e -2γ 2 /r . 1.3. The Jacobi semi-group. The following results are due to Biane, Pitman and Yor (see [START_REF] Biane | Probability Laws Related To The Jacobi Theta and Riemann Zeta Functions, and Brownian Excursions[END_REF], [START_REF] Pitman | Infinitely Divisible Laws Associated With Hyperbolic Functions[END_REF]) :

∞ 0 e -t 2 8 s f C h (s) ds = 1 cosh(t/2) h , h > 0 (5) ∞ 0 e -t 2 8 s f T h (s) ds = tanh(t/2) (t/2) h , h > 0 (6)
where (C h ) and (T h ) are two families of Lévy processes with respective density functions f C h and f T h for fixed h > 0. The densities of C h and T 1 are given by ( [START_REF] Biane | Probability Laws Related To The Jacobi Theta and Riemann Zeta Functions, and Brownian Excursions[END_REF]):

f C h (s) = 2 h Γ(h) p≥0 (-1) p Γ(p + h) p! f τ (2p+h) (s) f T 1 (s) = k≥0 e -π 2 2 (k+ 1 2 ) 2 s 1 {s>0}
where τ (c) = inf{r > 0; B r = c} is the Lévy subordinator (the first hitting time of a standard Brownian motion B) with corresponding density :

f τ (2p+h) (s) = (2p + h) √ 2πs 3 exp - (2p + h) 2 2s 1 {s>0} .
Let α, β satisfy α + β + 1 > 0, thus:

p 2/r (x, y) = √ 2πrW (y)e 2γ 2 /r 2 a-1 n≥0 (a) 2n (α + 1) n (β + 1) n P α,β n (z) (1 + xy) 8 n × f T 1 ⋆ f C 2n+a-1 (r)
or equivalently (where B stands for the Beta function):

p t (x, y) = √ πW (y) 2 α+β e γ 2 t √ t n≥0 (a) 2n (α + 1) n (β + 1) n P α,β n (z) (1 + xy) 8 n f T 1 ⋆ f C 2n+α+β+1 ( 2 t ).
1.4. The ultraspherical case. This case corresponds to α = β > -1 2 and we will proceed in a slightly different way. Indeed, a = 2α + 2 and

(4) = 1 -r (1 + r) 2α+2 F 4 (α + 1, α + 3/2, α + 1, α + 1; u, v) = 1 -r (1 + r) 2α+2 1 (1 -u -v) α+3/2 2 F 1 ( 2α + 3 4 , 2α + 5 4 , α + 1; 4uv (1 -u -v) 2 )
where the last equality follows from (see [START_REF] Yu | Integrals and Series[END_REF])

F 4 (b, c, b, b; u, v) = (1 -u -v) -c 2 F 1 (c/2, (c + 1)/2, b; 4uv (1 -u -v) 2 ).
Hence,

q t (x, y) = W (y)e 2α+1 2 t 2 α+1/2 sinh(t) (cosh t -xy) α+3/2 2 F 1 ( 2α + 3 4 , 2α + 5 4 , α + 1; (1 -x 2 )(1 -y 2 ) (cosh t -xy) 2 ) = W (y)e 2α+1 2 t 2 α+1/2 sinh(t) n≥0 [(2α + 3)/4] n [(2α + 5)/4] n (α + 1) n [(1 -x 2 )(1 -y 2 )] n (cosh t -xy) 2n+α+3/2 .
Besides, for h > 0, we may write:

1 cosh t -xy h = k≥0 (h) k k! (xy) k (cosh t) k+h since xy cosh t < 1 ∀x, y ∈] -1, 1[
, ∀t ≥ 0 and where we used:

1 (1 -r) h = k≥0 (h) k k! r k h > 0, |r| < 1.
Consequently, using Gauss duplication formula,

q t (x, y) = K α W (y)e 2α+1 2 t tanh(t) n,k≥0 Γ(ν(n, k, α) + 1)(xy) k k! n! Γ(α + n + 1) (1 -x 2 )(1 -y 2 ) 4 n 1 cosh t ν(n,k,α) where ν(n, k, α) = 2n + k + α + 1/2, K α = Γ(α + 1)/[2 α+1/2 Γ(α + 3/2)].
Thus, since γ = α + 1/2 when α = β, one has:

∞ 0 p s (x, y) s -3/2 e -γ 2 s e -t 2 4s ds = √ 2πΓ(α + 1) 2 α Γ(α + 3/2) tanh(t) t W (y) n,k≥0 Γ(ν(n, k, α) + 1)(xy) k k! n! Γ(α + n + 1) (1 -x 2 )(1 -y 2 ) 4 n 1 cosh t ν(n,k,α)
or equivalently:

∞ 0 p 1/2s (x, y)e -γ 2 2s e -t 2 2 s ds √ s = √ πΓ(α + 1) 2 α Γ(α + 3/2) tanh(t) t W (y) n,k≥0 Γ(ν(n, k, α) + 1)(xy) k k! n! Γ(α + n + 1) (1 -x 2 )(1 -y 2 ) 4 n 1 cosh t ν(n,k,α)
. Using ( 5), [START_REF] Yu | Integrals and Series[END_REF], f C h et f T 1 (we take t 2 /2 instead of t 2 /8), the density is written:

p 1/2s (x, y) = √ πsΓ(α + 1) 2 α Γ(α + 3/2) W (y)e γ 2 2s
n,k≥0

Γ(ν(n, k, α) + 1) k!n!Γ(α + n + 1) xy 2 k (1 -x 2 )(1 -y 2 ) 4 n f T 1 ⋆ f C ν(n,k,α) (s). Finally (7) p t (x, y) = √ πK α e γ 2 t √ t W (y) n,k≥0 Γ(ν(n, k, α) + 1)(xy) k k!n!Γ(α + n + 1) (1 -x 2 )(1 -y 2 ) 4 n f T 1 ⋆ f C ν(n,k,α) ( 1 2t
).

2. Application to statistics for diffusions processes 2.1. Some properties of the Jacobi process. Usually in probability theory, the Jacobi process is defined on [-1, 1] as the unique strong solution of the SDE :

dY t = 1 -Y 2 t dW t + (bY t + c)dt. It is straightforward that (Y t ) t≥0 L = (X t/2
) t≥0 where X is the Jacobi process already defined in section 1 with p = 2b, q = 2c. Using the variable change y → (y + 1)/2, the equation above transforms to (t → 4t) :

dJ t = 2 J t (1 -J t )dW t + [2(c -b) + 4bJ t ] dt = 2 J t (1 -J t )dW t + [d -(d + d ′ )J t ] dt where d = 2(c -b) = q -p = 2(β + 1) and d ′ = -2(c + b) = -(p + q) = 2(α + 1)
, which is the Jacobi process of parameters (d, d ′ ) already considered by Warren and Yor ([21]). Moreover, the authors provide the following skew-product : let Z 1 , Z 2 be two independent Bessel processes of dimensions d, d ′ and starting from z, z ′ respectively. Then :

Z 2 1 (t) Z 2 1 (t) + Z 2 2 (t) t≥0 L = (J At ) t≥0 , A t := t 0 ds Z 2 1 (s) + Z 2 2 (s) , J 0 = z z + z ′ .
Using well known properties of squared Bessel processes (see [START_REF] Revuz | Continuous Martingales And Brownian Motion[END_REF]), one deduce that if d ≥ 2 (β ≥ 0) and z > 0, then J t > 0 almost surely for all t > 0. Since 1 -J is still a Jacobi process of parameters (d ′ , d), then for d ′ ≥ 2, (α ≥ 0) and z ′ > 0, J t < 1 almost surely for all t > 0. The extension of these results to the matrix Jacobi process is established in [START_REF] Doumerc | Matrix Jacobi Process[END_REF] (Theorem 3.3.2, p.36). Since 0 is a reflecting boundary for Z 1 , Z 2 when 0 < d, d ′ < 2 (-1 < α, β < 0), then both 0 and 1 are reflecting boundaries for J.

2.2. LDP in the ultraspherical case. Let us consider the following SDE corresponding to the ultraspherical Jacobi process:

(8) dY t = 1 -Y 2 t dW t + bY t dt Y 0 = y 0 ∈] -1, 1[. Let Q b
y 0 be the law of (Y t , t ≥ 0) on the canonical filtered probability space (Ω, (F t ), F) where Ω is the space of ] -1, 1[-valued functions. The parameter b is such that b ≤ -1 (or α ≥ 0), so that -1 < Y t < 1 for all t > 0. The maximum likelihood estimate of b based on the observation of a single trajectory (Y s , 0 ≤ s ≤ t) under Q b 0 (see Overbeck [START_REF] Overbeck | Estimation for continuous branching processes[END_REF] for more details) is given by

(9) bt = t 0 Y s 1 -Y 2 s dY s t 0 Y 2 s 1 -Y 2 s ds .
The main result of this section is the following Theorem.

Theorem 1. When b ≤ -1, the family { bt } t satisfies a LDP with speed t and good rate function

(10) J b (x) =    - 1 4 (x -b) 2 x + 1 if x ≤ x 0 x + 2 + (b -x) 2 + 4(x + 1) if x > x 0 > b
where x 0 is the unique solution of the equation (b -x) 2 = 4x(x + 1) = 0, x < -1.

Proof of Theorem 1 : we follow the scheme of Theorem 3.1 in [START_REF] Zani | Large deviations for squared radial Ornestein-Uhlenbeck processes[END_REF]. Set :

S t,x := t 0 Y s 1 -Y 2 s dY s -x t 0 Y 2 s 1 -Y 2 s ds,
so that for x > b (resp. x < b), P ( bt ≥ x) = P (S t,x ≥ 0) (resp. P ( bt ≤ x) = P (S t,x ≤ 0)). Therefore, to derive a large deviation principle on { bt }, we seek a LDP result for S t,x /t at 0. Let us compute the normalized cumulant generating function Λ t,x of S t,x :

(11) Λ t,x (φ) = 1 t log Q b 0 (e φSt,x ).
¿From Girsanov formula, the generalized densities are given by ( 12)

dQ b 0 dQ b ′ 0 Ft = exp (b -b ′ ) t 0 Y s 1 -Y 2 s dY s - 1 2 (b 2 -b ′ 2 ) t 0 Y 2 s 1 -Y 2 s ds .
Let

F (Y t ) = - 1 2 log(1 -Y 2 t ) . From Itô formula, F (Y t ) = t 0 Y s 1 -Y 2 s dY s + 1 2 t 0 1 + Y 2 s 1 -Y 2 s ds = t 0 Y s 1 -Y 2 s dY s + t 0 1 1 -Y 2 s ds - t 2 .
Let us denote by

D 1 (x) = {φ : (b + 1) 2 + 2φ(x + 1) ≥ 0} .
For any φ ∈ D 1 (x), we can define b(φ, x) = -1 -(b + 1) 2 + 2φ(x + 1). With the change of probability defined by [START_REF] Exton | Multiple Hypergeometric Functions And Applications[END_REF] taking b ′ = b(φ, x), the stochastic integrals simplify to (see [START_REF] Zani | Large deviations for squared radial Ornestein-Uhlenbeck processes[END_REF] p. 125 for the details):

(13) Λ t (φ, x) = 1 t log Q b(φ,x) 0 (exp({φ + b -b(φ, x))[F (Y t ) -t/2]}).
When starting from y 0 = 0, (7) reads (t → t/2) :

pt (0, y) = √ 2πK α e γ 2 t/2 √ t n≥0 Γ(2n + α + 3 2 ) 4 n n!Γ(n + α + 1) (1 -y 2 ) n+α f T 1 ⋆ f C 2n+γ (1/t),
where p = -2(α + 1) = 2b ≤ -2 and γ = -(p + 1)/2 = α + 1/2. Denote by ( 14)

D(x) = {φ ∈ D 1 (x) : G(φ, x) = b + b(φ, x) + φ < 0} .
For any φ ∈ D(x), the expectation ( 13) is finite and a simple computation gives :

Λ t (φ, x) = - φ + b -b(φ, x) 2 + 1 t log Q b(φ,x) 0 ((1 -Y 2 t ) -(φ+b-b(φ,x))/2 ) = Λ(φ, x) + 1 t log √ 2πK α(φ,x) R t (φ, x) √ t ,
where

R t (φ, x) = n≥0 Γ(2n -b(φ, x) + 1/2) 4 n n!Γ(n -b(φ, x)) B n - φ + b + b(φ, x) 2 , 1 2 e γ 2 t/2 f T 1 ⋆ f C 2n+γ ( 1 t ) , α(φ, x) = -b(φ, x) -1
and B stands for the Beta function. With regard to (1), one has for φ ∈ D(x) :

lim t→∞ Q b(φ,x) 0 ((1 -Y 2 t ) -(φ+b-b(φ,x))/2 ) = C b,φ,x 1 -1
(1 -y 2 ) -[φ+b+b(φ,x)]/2 -1 dy < ∞ by dominated convergence Theorem. Hence Λ t → Λ as t → ∞. The following lemma, which proof is postponed to the appendix, details the domain D(x) (see ( 14)) of Λ t :

Lemma 1. Denote by

φ 0 (x) = - (b + 1) 2 2(x + 1) . i)If x < (b 2 + 3)/2(b -1): then D = (-∞, φ 0 (x)). ii) If (b 2 + 3)/2(b -1) < x < -1: then D(x) = (-∞, φ 1 (x)) where φ 1 (x) is solution of G(φ, x) = 0. iii) If x > -1: then D(x) = (φ 0 (x), φ 1 (x)).
In case i) of Lemma above, Λ is steep, i.e. its gradient is infinite at the boundary of the domain (for a precise definition, see [START_REF] Dembo | Large Deviations Techniques and Applications[END_REF]). It achieves its unique minimum in φ m (x) solution of ∂Λ ∂φ (φ, x) = 0 , i.e. b(φ(x), x) = x. It is easy to see that

φ m (x) = x + 1 2 - (b + 1) 2 2(x + 1) < φ 0 (x) .
Hence, Gärtner-Ellis Theorem gives for x < b < (b 2 + 3)/2(b -1), lim

t→∞ 1 t log P ( bt ≤ x) = lim t→∞ 1 t log P (S t,x ≤ 0) = inf φ∈]∞,φ 0 (x)] Λ(φ, x) = Λ(φ m (x), x) = - 1 4 (x -b) 2 x + 1 . If b < x < (b 2 + 3)/2(b -1), notice that φ m (x) > 0 and lim t→∞ 1 t log P ( bt ≥ x) = lim t→∞ 1 t log P (S t,x ≥ 0) = inf φ∈(0,φ 0 (x)] Λ(φ, x) = Λ(φ m (x), x) = - 1 4 (x -b) 2 x + 1 .
In cases ii) and iii) of Lemma 1, Λ is not steep. Nevertheless, if the infimum of Λ is reached in

• D(x), we can follow the scheme of Gärtner-Ellis theorem for the change of probability in the infimum bound. This infimum is reached if and only if [START_REF] Lassalle | Polynômes de Jacobi généralisés[END_REF] ∂Λ ∂φ (φ 1 (x), x) > 0 , i.e. if φ m (x) < φ 1 (x).

In case x + 1 > 0, we know (see proof of Lemma 1) that φ 1 (x) < φ m (x). If x + 1 < 0, we check the sign of G(φ m (x), x). We get the following dichotomy : Let x 0 denote the unique solution of g(x)

:= 4x(x + 1) -(b -x) 2 = 0, x < -1. Since g is decreasing on ]-∞, -1] and g(b 2 +3/(2(b-1)) = (3/4)(b+1) 2 > 0 = g(x 0 ), then x 0 > (b 2 +3)/[2(b-1)]. • if (b 2 + 3)/2(b -1) < x < x 0 < -1
, the derivative ∂Λ/∂φ(φ 1 (x), x) > 0, Λ achieves its minimum on φ m (x) and lim

t→∞ 1 t log P ( bt ≥ x) = Λ(φ m (x), x) = - (x -b) 2 4(x + 1
) .

• if x 0 < x < -1 or x > -1, then ∂Λ/∂φ(φ 1 (x), x) < 0. We apply Theorem 2 of the appendix, which is due to Zani [START_REF] Zani | Large deviations for squared radial Ornestein-Uhlenbeck processes[END_REF]. Let us verify that the assumptions are satisfied, and more precisely that Λ t can take the form [START_REF] Overbeck | Estimation for continuous branching processes[END_REF]. Indeed, the only singularity

φ 1 (x) of R t comes from B (n -[φ + b + b(φ, x)]/2, 1/
2) when n = 0, and more precisely, from Γ(-[φ + b + b(φ, x)]/2) . We can write

(16) Λ t (φ, x) = Λ(φ, x) + 1 t log Γ - φ + b + b(φ, x) 2 + 1 t log √ 2πK α(φ,x) Rt (φ, x) √ t , where (17) Rt (φ
, x) = R t (φ, x)/Γ(-[φ + b + b(φ, x)]/2). Now ∀n ≥ 0, B n - φ + b + b(φ, x) 2 , 1 2 /Γ - φ + b + b(φ, x) 2
is analytic on some neighbourhood of φ 1 (x). Besides,

lim φ→φ 1 (x) ,φ<φ 1 (x) b + φ + b(φ, x) φ -φ 1 (x) = c > 0 ,
and since lim ρ→0 + ρΓ(ρ) = 1, then φ 1 (x) is a pole of order one of Γ (φ + b + b(φ, x)/2) and one writes:

1 t log Γ - φ + b + b(φ, x) 2 = - log(φ 1 (x) -φ) t + h(φ) t .
The function h is analytic on D(x) and can be extended to an analytic function on ]φ 1 (x) -ξ, φ 1 (x) + ξ[ for some positive ξ. Finally, to satisfy Assumption 1 of the appendix, we focus on Rt (φ, x)/ √ t and show that it is bounded uniformly as t → ∞. To proceed, we shall prove that this ratio is bounded from above and below away from 0. Setting A n (t) := e γ 2 t/2 f T 1 ⋆ f C 2n+γ (1/t), one has :

A n (t) √ t ≤ e γ 2 t/2 √ t k,l≥0 U k,n 1/t 0 exp - 1 2 (2n + 2k + γ) 2 s + π 2 (l + 1 2 ) 2 ( 1 t -s) ds s 3/2 = e γ 2 t/2 √ t k,l≥0 U k,n ∞ t exp - 1 2 (2n + 2k + γ) 2 s + π 2 (l + 1 2 ) 2 ( s -t ts ) ds √ s < e -2n 2 k,l≥0 U k,n e -2k 2 ∞ t exp - 1 2 (2n + 2k + γ) 2 (s -t) + π 2 (l + 1 2 ) 2 ( s -t ts ) ds √ ts = e -2n 2 k,l≥0 U k,n e -2k 2 ∞ 0 exp - 1 2 (2n + 2k + γ) 2 s + π 2 l 2 ( s t(t + s) ) ds t(t + s) with U k,n = Γ(2n + k + γ)2 2n+γ (2n + 2k + γ) k!Γ(2n + γ) .
Let Θ(x) = l∈Z e -πl 2 x = 1 + 2 l≥1 e -πl 2 x denote the Jacobi Theta function ( [START_REF] Biane | Probability Laws Related To The Jacobi Theta and Riemann Zeta Functions, and Brownian Excursions[END_REF]). Then

A n (t) √ t < e -2n 2   k≥0 U k,n e -2k 2 ∞ 0 exp - (2n + 2k + γ) 2 s 2 Θ πs 2t(t + s) ds t(t + s) + C(n, t)   where C(n, t) = 1 2 √ t k,l≥0 U k,n e -2k 2 ∞ 0 exp - (2n + 2k + γ) 2 s 2 ds √ t + s .
Recall that Θ(x) = (1/ √ x)Θ(1/x) ( [START_REF] Biane | Probability Laws Related To The Jacobi Theta and Riemann Zeta Functions, and Brownian Excursions[END_REF]), which yields :

A n (t) √ t < e -2n 2 k≥0 U k,n e -2k 2 ∞ 0 exp - (2n + 2k + γ) 2 s 2 Θ 2t(t + s) πs ds √ s + C(n) 2 √ t
where

C(n) = e -2n 2 k,l≥0 U k,n e -2k 2 ∞ 0 exp - (2n + 2k + γ) 2 s 2 ds √ s .
Since e -l 2 z < e -lz , then Θ(z) ≤ 3 for z > 1. Hence, as 2t/π ≤ 2t(t + s)/(πs), then for t large enough:

A n (t) √ t < 3e -2n 2 k≥0 U k,n e -2k 2 ∞ 0 exp - (2n + 2k + γ) 2 s 2 ds √ s + C(n) < 4C(n).
The upper bound follows since n C(n) < ∞. Besides,

Rt (φ, x) √ t > √ πΓ(1/2 -b(φ, x)) Γ(-b(φ, x))Γ{[1 -(φ + b + b(φ, x)]/2} A 0 (t) √ t = C(b, φ, x) k,l≥0 (-1) k V k ∞ 0 exp - 1 2 (2k + γ) 2 s + π 2 (l + 1 2 ) 2 ( s t(t + s) ) ds t(t + s)
where V k (t) := U k,0 e -2k(k+γ)t . One may choose t large enough independent of k such that V k (t) ≥ V k+1 (t) for all k ≥ 0. In fact, such t satisfies:

e 2(2k+γ+1)t ≥ e 2t ≥ sup k≥0 U k+1,0 U k,0 = sup k≥0 (k + γ)(2k + γ + 2) (k + 1)(2k + γ) Then: Rt √ t > C(b, φ, x)[V 0 (t) -V 1 (t)] l≥0 ∞ 0 exp - 1 2 γ 2 s + π 2 (l + 1 2 ) 2 ( s t(t + s) ) ds t(t + s) > C(b, φ, x)[γ2 γ -V 1 (t)] l≥0 ∞ 0 exp - 1 2 γ 2 s + π 2 (l + 1) 2 ( s t(t + s) ) ds t(t + s) = C(b, φ, x) 2 [γ2 γ -V 1 (t)] ∞ 0 e -γ 2 s/2 Θ πs 2t(t + s) ds t(t + s) -C(t)
where

C(t) = 1 2 √ t ∞ 0 e -γ 2 s/2 ds (t + s) < c ∞ 0 e -γ 2 s/2 ds √ s , c < 2 π .
for t large enough. Following the same scheme as for the upper bound, one gets:

Rt √ t > C(b, φ, x) 2 γ2 γ 2 π ∞ 0 e -γ 2 s/2 Θ 2t(t + s) πs ds √ s -C(t) > C(b, φ, x) 2 γ2 γ 2 π -c ∞ 0 e -γ 2 s/2 ds √ s > 0.
As a result, lim

t→∞ 1 t log P ( bt ≥ x) = Λ(φ 1 (x), x) = -(x + 2 + (b -x) 2 + 4(x + 1)),
which ends the proof of Theorem 1.

2.3. Jacobi and squared Bessel processes duality. By Itô's formula and Lévy criterion, one claims that

(Z t = Y 2 t ) t≥0 is a [0, 1]-valued Jacobi process of parameters d = 1, d ′ = -2b ≥ 2. Indeed: dZ t := d(Y 2 t ) = 2Y t dY t + Y t = 2Y t 1 -Y 2 t dW t + [(2b -1)Y 2 t + 1]dt = 2 Z t (1 -Z t )sgn(Y t )dW t + [(2b -1)Z t + 1]dt = 2 Z t (1 -Z t )dB t + [(2b -1)Z t + 1]dt.
Using the skew product previously stated, there exists R, a squared Bessel process of dimension d ′ = 2(ν + 1) = -2b and starting from R 0 = r so that:

νt := -bt -1 = log(1 -Y 2 t ) + t 2 t 0 Y 2 s 1-Y 2 s ds
is another estimator of ν based on a Jacobi trajectory observed till time t. Set t = log u, then with associated rate function :

νlog u = log[u(1 -Y 2 log u )] 2 log u 0 Y 2 s 1-Y 2 s ds = log[u(1 -Y 2 log u )] 2 u 1 Y 2 log s s(1-Y
I ν (x) = (x-ν) 2 4x if x ≥ x 1 := -(ν+2)+2 √ ν 2 +ν+1 3 1 -x + (ν -x) 2 -4x if x < x 1 .
A glance at both rate functions gives I ν (x) = J -(ν+1) (-(x + 1)) and x 0 = -(x 1 + 1).

Appendix

3.1.

A large deviations principle in a non steep case. Let {Y t } t≥0 be a family of real random variables defined on (Ω, F, P ), and denote by µ t the distribution of Y t . Suppose -∞ < m t := EY t < 0. We look for large deviations bounds for P (Y t ≥ y). Let Λ t be the n.c.g.f. of Y t : Λ t (φ) = 1 t log E(exp{φtY t }) , and denote by D t the domain of Λ t . We assume that there exists 0 < φ 1 < ∞ such that for any t sup{φ : φ ∈ D t } = φ 1 and [0, φ 1 ) ⊂ D t . We assume also that for φ ∈ D t Assumption 1. [START_REF] Overbeck | Estimation for continuous branching processes[END_REF] Λ t (φ) = Λ(φ) -α t log(φ 1 -φ) + R t (φ) t where • α > 0 • Λ is analytic on (0, φ 1 ), convex, with finite limits at endpoints, such that Λ ′ (0) < 0, Λ ′ (φ 1 ) < ∞, and Λ ′′ (φ 1 ) > 0.

• R t is analytic on (0, φ 1 ) and admits an analytic extension on a strip D γ β = (φ 1 -β, φ 1 + β) × (-γ, γ), where β and γ are independent of t.

• R t (φ) converges as t → ∞ to some R(φ) uniformly on any compact of D γ β .

Theorem 2. Under 1 For any Λ ′ (0) < y < Λ ′ (φ 1 ), ( 19) lim The rate function is continuously differentiable with a linear part. • If x + 1 < 0, ∂G ∂φ > 0 and G(•, x) is increasing. Then we see easily that G(φ 0 (x), x) < 0 iff x < (b 2 + 3)(2(b -1)), which determines cases i) and ii).

• If x + 1 > 0, φ → ∂G ∂φ is increasing hence there exists φ(x) such that ∂G ∂φ ( φ(x), x) = 0.

We compute

φ(x) = x + 1 2 - (b + 1) 2 2(x + 1) = φ m (x) .
We see that G( φ(x), x) < 0, and there exists φ 1 (x) < φ(x) such that G(φ 1 (x), x) = 0, and D(x) = (φ 0 (x), φ 1 (x)) .

t→+∞ 1 t 1 t

 11 log P (Y t ≥ y) = -sup φ∈(0,φ 1 ) {yφ -Λ(φ)} .For any y ≥ Λ ′ (φ 1 ), log P (Y t ≥ y) = -yφ 1 + Λ(φ 1 ) .

3. 2 .

 2 Proof of Lemma 1 : Note first that (b 2 + 3)(2(b -1)) < -1 if b < -1 and that the condition φ ∈ D 1 (x) ⇒ φ ≥ φ 0 (x) if x > -1 and φ ≤ φ 0 (x) if x < -1. To examine the behaviour of G, we compute ∂G ∂φ (φ, x) = 1 -x + 1(b + 1) 2 + 2φ(x + 1) .

  2 log s ) ds and {ν log u } u satisfies a LDP with speed log u and rate function J -(ν+1) (-(x + 1)).When starting at R 0 = 1, the MLE of ν based on a Bessel trajectory is given by (cf[START_REF] Zani | Large deviations for squared radial Ornestein-Uhlenbeck processes[END_REF], p. 132):

	ν1 t =	t 0	dXs Xs -2 2 t 0 ds Xs	t 0	ds Xs	=	log(X t ) t ds 2 0 Xs

(P α,β n (x)) n≥0 are normalized such that they form an orthogonal basis with respect to the probability measure W (y)dy which is not the same used in[START_REF] Andrews | Special functions[END_REF].


Mathieu Des Noes
email: mathieu.desnoes@cea.fr

Valentin Savin
email: valentin.savin@cea.fr

Jean Marc Brossier
email: jean-marc.brossier@gipsa-lab.grenoble-inp.fr

Laurent Ros
email: laurent.ros@gipsa-lab.grenoble-inp.fr

Blind Identification of the Uplink Scrambling Code Index of a WCDMA Transmission and Application to Femtocell Networks

Interference between macro and femtocells is an important issue for the development of WCDMA femtocell networks. More specifically, the uplink signal of a macro User Equipment may generate an unacceptable level of interference at the femto Base Station. To avoid this situation, interference mitigation techniques could be implemented. All the proposed techniques require the knowledge of the uplink scrambling code index of the interferer. Unfortunately, if the femto BS is in a closed access mode, there are no signalling links with the surrounding macro BSs. The femto BS has to estimate blindly this scrambling code index. An algorithm which performs a blind identification of the uplink scrambling code index of a WCDMA transmission is proposed in this article. This gives the possibility to implement interference cancelation algorithm at the femto BS.

INTRODUCTION

The development of femtocell networks is an important perspective for increasing cellular systems capacity. A femtocell is covered by a small Base Station (BS), designed for typically indoor environments (home, business) which does not require a coordinated deployment [1]. Although most of the current research activities on this topic focus on the LTE system [START_REF] Holma | LTE for UMTS: OFDMA and SC-FDMA based radio access[END_REF], current deployed systems use the UMTS-WCDMA technology [START_REF] Holma | WCDMA for UMTS[END_REF]. In this paper, an algorithm which estimates blindly the uplink scrambling code index of a User Equipment (UE) is proposed. It offers the perspective of implementing interference mitigation techniques at a femto Base Station (BS), and hence providing a solution to the problem of uplink macro to femto interference issue in WCDMA femtocell networks.

A macro BS provides the overall coverage, while femto BSs offer better indoor coverage to UEs attached to them. A femto BS can be configured to operate in open or closed access mode to visiting UEs [START_REF] Xia | Open vs. closed access femtocells in the uplink[END_REF]. In open access mode, a visiting UE is allowed to handover from the macro BS to a femto BS in order to send its data. This generates additional complexity to route the data packets and also to ensure communication security. In closed access mode, a visiting UE is not allowed to handover. This simplifies the network architecture, but this may lead to an unacceptable interference level at the femto BS. This situation occurs if a visiting UE transmits at high Fig. 1. Interference generated by a macro UE on a femto BS power, while it is located nearby the femto BS. This situation is depicted in Fig 1 . A macro UE transmits at high power because it is either at the cell edge or inside a building. Since the power control procedure concerns only the uplink with the macro BS, it will be received at the femto BS with a power much larger than the power of a femto UE. This is the well known near-far effect in CDMA systems [START_REF] Verdu | Multiuser detection[END_REF]. If the power of this interferer is too large, the femto BS may not be able to demodulate any communication with its attached UEs. This creates a "dead zone" in the network coverage. In order to avoid this situation, it is required to implement interference mitigation techniques. This subject has been deeply studied in the past two decades for CDMA systems [START_REF] Andrews | Interference cancellation for cellular systems: a contemporary overview[END_REF]. All the proposed techniques exploit the knowledge of the UE's scrambling code, which is chosen from a set of scrambling codes identified by their index. The scrambling code index of the macro UE is allocated by the macro BS and is signalled to the UE in a dedicated control channel [START_REF].0 Spreading and modulation (FDD), 3GPP[END_REF]. Unfortunately, in a closed access mode, there is no signalling link between the macro and femto BSs. Hence, a femto BS has no knowledge of the scrambling code index of an interfering macro UE. It has to estimate this index blindly. To the authors' knowledge, [START_REF] Kerr | Iterative signal processing for blind code phase acquisition of CDMA 1x signals for radio spectrum monitoring[END_REF] is the first and only article addressing this issue. The authors exploit the specificities of the framing, spreading and multiplexing procedures defined by CDMA systems to process the signal so that the resulting signal can be considered to be a linear code depending on the initial state of the scrambling code generator. This processing step can be applied to both WCDMA and CDMA 1x systems [START_REF].0 Spreading and modulation (FDD), 3GPP[END_REF] [START_REF] Gpp | -A -Physical Layer Standard for cdma[END_REF]. The received codeword is then estimated with a maxlog-MAP algorithm [START_REF] Robertson | A comparison of optimal and sub-optimal MAP decoding algorithms operating in the log domain[END_REF]. The initial state of the different spreading code generators is eventually obtained by applying a pseudoinverse of the code generator matrix. In this paper, a different decoding strategy is proposed. Exploiting the unique properties of m-sequence, it is shown that an iterative message-passing algorithm [START_REF] Yeung | An iterative algorithm and low complexity hardware architecture for fast acquisition of long PN codes in UWB systems[END_REF] can be implemented to decode the received signal after the same initial processing step as in [START_REF] Kerr | Iterative signal processing for blind code phase acquisition of CDMA 1x signals for radio spectrum monitoring[END_REF]. Then the scrambling code index can be obtained with a proper used of the "Shift and add" and decimation properties of m-sequences. In addition, in [START_REF] Kerr | Iterative signal processing for blind code phase acquisition of CDMA 1x signals for radio spectrum monitoring[END_REF], the authors focused on the CDMA 1x system, while the WCDMA system is addressed in this article. The paper is organized as follows. Section I presents the procedure used to generate the uplink WCDMA signal. The construction of the complex scrambling code is also detailed since its specific properties will be exploited by the proposed algorithm. Section II details the blind identification algorithm, which is split in 3 steps. Section III presents simulation results and Section IV concludes this paper. Notation: a sequence will be written with upper case in its BPSK representation (S(k) ∈ {-1, +1}) and with lower case in its binary representation (s(k) ∈ {0, 1}).

I. UPLINK WCDMA SIGNAL GENERATION

In this section, the modulation, spreading and multiplexing operations implemented in the uplink of the WCDMA system are first described. Then, the method standardized by 3GPP for constructing the scrambling code is detailed [START_REF].0 Spreading and modulation (FDD), 3GPP[END_REF]. It is included here for completeness, since several specific features of the scrambling will be exploited by the proposed blind estimation algorithm.

A. Spreading and multiplexing

The uplink of the WCDMA system implements a Code Division Multiple Access (CDMA) scheme. The data and control channels are mapped on a BPSK constellation and I/Q multiplexed. In the standard, these channels are referred to as the Dedicated Physical Data CHannel (DPDCH) and Dedicated Physical Control CHannel (DPCCH). The bits are first spread by a channelization sequence (C d and C c ) which ensures orthogonality between channels. The spreading factor of the DPCCH signal is equal to SF DPCCH = 256, while for DPDCH, SF DPDCH is variable and depends on the service data rate. When a DPDCH is transmitted, its channelization sequence C d is a periodic repetition of the pattern "+1, +1, -1, -1". This is due to the selection of the column number SF DPDCH /4 of the Walsh-Hadamard matrix [START_REF].0 Spreading and modulation (FDD), 3GPP[END_REF]. The channelization sequence C c used for DPCCH is the all '+1' sequence. Each spread signal is then weighted by a factor β to set its power, and I/Q multiplexed. The resulting complex signal is eventually multiplied by the scrambling code S n , where index n indicates that the scrambling code is user specific. At chip time index k, the transmitted signal can be modeled as follows:

E(k) = S n (k)(β d X DPDCH (k) + jβ c X DPCCH (k)) (1)
X DPDCH and X DPCCH models the DPDCH and DPCCH bits spread by sequences C d and C c .

In order to simplify the description, we have considered the case where only one DPDCH is transmitted. However, we note that the proposed algorithm also works when 2 DPDCHs are transmitted and I/Q multiplexed, as detailed [START_REF].0 Spreading and modulation (FDD), 3GPP[END_REF].

B. Scrambling code generation

The scrambling code s n is built from two m-sequences x n and y. Fig. 2 shows a Linear Feedback Shift Register (LFSR) sequence generator according to the Fibonacci representation [START_REF] Peterson | Introduction to spreadspectrum communications[END_REF].

D D D D D g r-1 g r-2 g 1 + + + + g 2 X(k) X(k+1) X(k+2) X(k+r) X(k+r-1) X(k+r-2)

Fig. 2. LFSR sequence with the Fibonacci representation

The two m-sequences x n and y are defined by their primitive polynomial :

g x (D) = D 25 + D 22 + 1 g y (D) = D 25 + D 24 + D 23 + D 22 + 1
These polynomials have a different definition than in the standard [START_REF].0 Spreading and modulation (FDD), 3GPP[END_REF]. This is due to the notation used in Figure 2, which is the conventional notation for LFSR sequences defined in [START_REF] Peterson | Introduction to spreadspectrum communications[END_REF] and [START_REF] Mceliece | Finite fields for computer scientists and engineers[END_REF]. The polynomials defined in the standard are obtained by exchanging coefficient g i by g r-i . The scrambling sequence is identified by an index n. Let n 23 , • • • , n 0 be the 24 bits binary representation of the index n with n 0 being the least significant bit. The registers of the sequences x n and y are reset at the beginning of each frame (every 38400 chips), with the following values :

x n (0) = n 0 , x n (1) = n 1 , . . . , x n (23) = n 23 , x n (24) = 1 y(0) = y(1) = y(2) = • • • = y(23) = y(24) = 1
The initial state of a LFSR sequence is defined by the content of its shift registers at time k = 0. As a consequence, the initial state of the sequence y at the beginning of each frame is known.

The complex scrambling code is built as follows (k = 0, 1, . . . , 38399):

S n (k) = C 1 (k) 1 + j(-1) n C 2 2⌊ k 2 ⌋ C 1 (k) = Z n (k) C 2 (k) = Z n ((k + 16777232) mod(2 25 -1)) Z n (k) = 1 -2(x n (k) ⊕ y(k)) (2)
where⌊x⌋ is the closest integer inferior or equal to x and ⊕ denotes the XOR operation. M-sequences have two interesting properties that will be exploited by the proposed algorithm [START_REF] Mceliece | Finite fields for computer scientists and engineers[END_REF]:

• "Shift and add" property : for two given delays τ 1 and τ 2 , there exists a unique τ 3 such that:

x(k + τ 1 ) ⊕ x(k + τ 2 ) = x(k + τ 3 )
• Decimation property: the decimation by a factor 2 of a m-sequence gives a shifted version of this m-sequence.

There exists a unique τ such that:

x(2k) = x(k + τ )
According to the "shift and add" property of m-sequences, it can be proved that for τ = 16777232:

x n (i + τ ) = x n (i + 4) ⊕ x n (i + 7) ⊕ x n (i + 18) y(i + τ ) = y(i + 4) ⊕ y(i + 6) ⊕ y(i + 17) (3)
This property is used to generate the sequence C 2 .

The vector (n 0 , • • • , n 23 ) is allocated by the macro BS to the macro UE in a signalling message and is unknown to the femto BS. The objective of the algorithm proposed in this paper is to perform a blind estimation of this vector.

II. BLIND IDENTIFICATION OF THE SCRAMBLING CODE

The algorithm is split in 3 steps:

• step 1 : chip level processing which results in a direct observation of a modified version of the sequence x n , denoted by x.

• step 2 : estimation of the initial shift registers state of sequence x with an iterative message-passing decoder.

• step 3 : determination of the initial state of the sequence x n with the use of a transposition matrix.

In order to ease the comprehension of the algorithm, we restrict the description to an Additive White Gaussian Noise (AWGN) channel model. The robustness to multipath channel will be discussed in section II-E and be tested by means of simulation. In addition, we also focus the description to a signal sampled at the chip rate.

A. Chip level processing

The received signal is modeled as follows (see Eq. ( 1)):

T (k) = e jθ S n (k)(β d X DPDCH (k)+jβ c X DPCCH (k))+n(k) ( 4
)
where θ is the phase rotation introduced by the channel and n(k) the additional noise modeling thermal noise as well as other sources of interference. First, a differential multiplication is performed on the received signal:

U (k) = T (k + 1)T (k) * = S n (k + 1)S n (k) * {β 2 d X DPDCH (k + 1)X DPDCH (k) +β 2 c X DPCCH (k + 1)X DPCCH (k) +jβ c β d (X DPCCH (k + 1)X DPDCH (k) -X DPDCH (k + 1)X DPCCH (k))} + I(k)
where I(k) contains all the noise cross-product terms. Let define the sequence

A(k) = C 1 (2k + 1)C 1 (2k)C 2 (2k).
The scrambling code S n (k) satisfies the following relation (see Eq. ( 2)):

S n (2k + 1)S n (2k) * = -2jA(k)

Moreover, due to the properties of the channelization sequence C d and C c , the spread signals X DPDCH (k) and X DPCCH (k) satisfy the 3 relations:

X DPDCH (2k + 1)X DPDCH (2k) = 1 X DPCCH (2k + 1)X DPCCH (2k) = 1 X DPCCH (2k + 1)X DPDCH (2k) = X DPDCH (2k + 1)X DPCCH (2k)
Exploiting these properties, if the differential signal U (k) is decimated, the negative value of its imaginary part satisfies :

V (k) = -Im (U (2k)) = 2(β 2 d + β 2 c )A(k) + w(k) (5)
The variable w(k) contains all the noise cross-product terms. From the construction of the scrambling code detailed in the previous section, it is straight forward to observe that A(k) is the BPSK representation of a binary sequence a(k) which depends on sequence x n and y (see Eq. 2):

a(k) = x(k) ⊕ ỹ(k) (6) with x(k) = x n (2k) ⊕ x n (2k + 1) ⊕ x n (2k + 16777232) ỹ(k) = y(2k) ⊕ y(2k + 1) ⊕ y(2k + 16777232)
Using the result of Eq. ( 3), we obtain:

x(k) = x n (2k) ⊕ x n (2k + 1) ⊕ x n (2k + 4) ⊕ x n (2k + 7) ⊕ x n (2k + 18) (7) ỹ(k) = y(2k) ⊕ y(2k + 1) ⊕ y(2k + 4) ⊕ y(2k + 6) ⊕ y(2k + 17)
From the "shift and add" and decimation properties, x and ỹ are shifted version of the m-sequences x n and y. The sequence a is thus a Gold sequence [START_REF] Gold | Optimal binary sequences for spread spectrum multiplexing[END_REF]. As a consequence, Eq. ( 5) tells us that V (k) is the observation of the BPSK representation of a Gold sequence with additive noise. It is known from the literature that it is possible to estimate the initial state of a LFSR sequence with a standard iterative message-passing algorithm [START_REF] Yeung | An iterative algorithm and low complexity hardware architecture for fast acquisition of long PN codes in UWB systems[END_REF][15] [START_REF] Mihaljevic | A low-complexity and highperformance algorithm for the fast correlation attack[END_REF]. The decoding strategy proposed in this paper is however different. In order to identify the initial state of the sequence x n , the proposed algorithm exploits the a priori knowledge of the initial state of the sequence y at the beginning of each frame (all '1'). At each time instant q, the receiver assumes that it is synchronized with the beginning of the frame. It is thus possible to generate the sequence ỹ by using (3). Then, the elements of vector (V (q), • • • , V (q + M -1)) are multiplied chip by chip with sequence Ỹ , the BPSK representation of the sequence ỹ:

R q (k) = V (k + q) Ỹ (k) k = 0, . . . , M -1 (8)
This operation eliminates ỹ(k) from ( 6), and hence R q (k) is a noisy observation of the BPSK representation of sequence

x.

Since x is a shifted version of sequence x, it can be decoded with an iterative message-passing algorithm which parity check matrix is matched to the sequence x (see section II-B).

The vector (R q (0), • • • , R q (M -1)) feeds this decoder, and if it fails to find a codeword, this means that either the frame synchronization assumption is not valid or a decoding failure happened. In both cases, the procedure is restarted when the next chip is received. If the decoder finds a valid codeword, the initial state of sequence x n is found by a simple matrix multiplication (see section II-C).

B. Iterative message-passing decoding

A m-sequence x satisfies the following parity check equation (g 0 = g r = 1), for all k ≥ 0:

r i=0 g r-i x(k + i) = 0
A m-sequence is thus a cyclic linear code with rate r 2 r -1 . A codeword is generated by one initial state of the shift registers. In the context of this paper, only a sequence of M variables, corresponding to M consecutive bits of the codeword, is observed. The parity check matrix of this code depends on the sequence's primitive polynomial g(D):

H =        g r • • • g 0 0 • • • • • • 0 0 g r • • • g 0 0 • • • 0 . . . . . . . . . . . . . . . . . . . . . 0 • • • 0 g r • • • g 0 0 0 • • • • • • 0 g r • • • g 0        (9)
Once the parity check matrix has been defined, it is possible to decode the received vector (R q (0), • • • , R q (M -1)) with a standard iterative message passing algorithm [START_REF] Kschischang | Factor graphs and the sum-product algorithm[END_REF] [START_REF] Wiberg | Codes and decoding on general graphs[END_REF]. In addition, as it was proposed in [START_REF] Yeung | An iterative algorithm and low complexity hardware architecture for fast acquisition of long PN codes in UWB systems[END_REF] [START_REF] Principe | Rapid acquisition of gold codes and related sequences using iterative message passing on redundant graphical models[END_REF], the use of Redundant Graphical Model (RGM) improves greatly the decoder performance. If g(D) is the sequence polynomial in GF(2), it satisfies:

g(D 2 n ) = g(D) 2 n
This property is exploited to create additional parity check equations. Polynomial g n (D) = g(D 2 n ) also generates a parity check matrix H n similar to H. Theses matrices can be concatenated to create a larger parity check matrix H RGM [START_REF] Yeung | An iterative algorithm and low complexity hardware architecture for fast acquisition of long PN codes in UWB systems[END_REF]:

H RGM =      H 0 H 1 . . . H n RGM -1      ( 10
)
where n RGM is the number of RGMs used for decoding. These RGMs increase the column weight of the parity check matrix, while keeping constant the row weight.

If the decoding is successful (all parity check equations are satisfied), the soft decision output of the decoder is converted into a binary representation with a hard decision rule. Then, according to the Fibonacci representation (Fig. 2), the first r bits of the codeword represents the content of the shift registers at initialization.

C. Determination of the initial state of the sequence x n

The decoder provides a vector of r bits representing the initial state of sequence x(k), denoted by Q x. We want to find the initial state of sequence x n (k):

Q xn = (n 0 , n 1 , • • • , n 24 )
T knowing that n 24 = 1. This task is achieved in 2 steps. The first step consists in finding the initial state of the sequence x decim (k), Q x decim , which gives x by decimation by a factor 2:

x(k) = x decim (2k) [START_REF] Arazi | Decimation of m-sequences leading to any desired phase shift[END_REF] shows that there is a fixed transposition matrix B between these two vectors and describes the method to compute this matrix:

Q x decim = BQ x
The second step eventually provides vector Q xn . It exploits the "shift and add" relation between the sequences x decim and x n (see Eq [START_REF].0 Spreading and modulation (FDD), 3GPP[END_REF]):

x decim (k) = x n (k)⊕x n (k+1)⊕x n (k+4)⊕x n (k+7)⊕x n (k+18)
Using the state transition matrix G [START_REF] Peterson | Introduction to spreadspectrum communications[END_REF], defined by :

G =          0 1 0 0 • • • 0 0 0 1 0 • • • 0 0 0 0 1 • • • 0 . . . . . . . . . 0 • • • . . . 0 0 0 0 • • • 1 1 g r-1 g r-2 g r-3 • • • g 1          (11)
we have:

Q xn = (I r + G + G 4 + G 7 + G 18 ) -1 Q x decim
where I r is the r × r identity matrix. These two steps are finally combined in a unique transposition matrix T :

Q xn = T Q x T = (I r + G + G 4 + G 7 + G 18 ) -1 B
It is important to note that matrix T shall be computed once for all and stored in memory.

D. Flow chart of the algorithm

A flow chart of the algorithm is presented in Fig. 3. The decoding operation must be implemented at the chip rate, which requires a very high data rate decoding capability.

E. Robustness to a multipath channel

The algorithm is by construction robust to multipath whose delays are larger than a chip period. This is due to the chipby-chip multiplication with sequence Y (k) in [START_REF] Kerr | Iterative signal processing for blind code phase acquisition of CDMA 1x signals for radio spectrum monitoring[END_REF]. A multipath component delayed by more than one chip will be scrambled by sequence y(k) and will be considered as an additional noise source by the decoder. In order to cope with multipath, the receiver implements a search window of length W chips. At each time q, the receiver runs the detection algorithm for time q to q + W -1. If a valid codeword is detected within this window, the detection is declared successful.

III. SIMULATION RESULTS

Let us first define the conventional synchronization hypothesis :

• H 0 : the receiver is not synchronized with the beginning of the frame.

• H 1 : the receiver is synchronized with the beginning of the frame. The performance of the algorithm are measured by the probabilities of correct detection P d , false alarm P f a and missed detection P m , defined as follows:

P d = P (Ic = 1 and Qxn = Q xn |H 1 ) P f a = P (Ic = 1|H 0 ) P m = 1 -P d
Ic is the indication function of the decoder:

Ic =
1 if all parity check equations are satisfied 0 otherwise

Qxn is the estimated initial state of sequence x n , given by the decoder output.

Performances are measured with the following simulation configurations:

• The transmitter sends a UL reference measurement channel for the 12.2 kbps service, as specified in [START_REF].0 BS Radio Transmission and Reception (FDD), 3GPP[END_REF]: • When measuring P m , the receiver is synchronized with the beginning of the frame (i.e. hypothesis H 1 is satisfied), while its is not synchronized when P f a is evaluated.

SF DPDCH =
• The decoder implements a Min-Sum (MS) messagepassing algorithm [START_REF] Wiberg | Codes and decoding on general graphs[END_REF]. The decoder stops when either all the parity check equations are satisfied or the maximum number of iteration N iter is reached.

A. AWGN channel

Simulations over 10 7 frames did not produce any false alarm. This means that P f a < 10 -6 with a good confidence level.

Fig. 4 shows the probability of missed detection as a function of the number of RGMs, for M = 4000 variables at the decoder input and N iter = 20 iterations at most. P m is measured as a function of the Signal to Noise ratio (SNR) at the input of the receiver. For an arbitrary defined target P m = 0.1, a gain of about 7 dB is obtained with 7 RGMs with respect to the case with only one graphical model. With 7 RGMs, the target P m = 0.1 is reached at a SNR around -7 dB. According to [START_REF].0 BS Radio Transmission and Reception (FDD), 3GPP[END_REF], a UE transmitting a 12.2 kbps service in an AWGN channel shall be received at its serving BS with a SNR larger than -17 dB. Hence, if it is received at the femto BS at a SNR equal to -7 dB, its power is 10 dB larger than a femto UE transmitting the same service. This macro UE is thus a strong interferer, but its scrambling code index can be detected and it can be mitigated by an appropriate interference cancelation algorithm.

Simulations not reported in this paper have also shown that there is no gain to have more than N iter = 20 iterations in the decoding algorithm. Fig. 5 shows the sensitivity of the algorithm to the number of variables M (see [START_REF] Kerr | Iterative signal processing for blind code phase acquisition of CDMA 1x signals for radio spectrum monitoring[END_REF]). There is no interest to increase the number of variables above 4000 chips. A gain of 0.5dB is achieved between M = 2000 and M = 4000. Since this gain is small, this leaves some room for a performance/complexity trade-off.

B. Multipath channel

The probability of detection is defined by the probability to detect the signal when the search window contains the beginning of the frame (see Section II-E). The performance are evaluated with 2 static multipath channels (Static A and B), having respectively 2 and 3 multipath with equal gains. The delays of the multipath are equal to {0, 2} and {0, 2, 4} chips for Static A and B channels. The search widow is set to W = 10 chips, which is larger than the channel delay spread. Fig. 6 compare the performance of the detection algorithm with AWGN, Static A and B channels. Compared to a gaussian channel, there is a degradation of about 3 dB and 8 dB at P m = 0.1 for Static A and B. This is due to the diminution of the SNR per multipath. Even if the degradation is noticeable, the algorithm is still operational when the channel contains multipath.

IV. CONCLUSION

An algorithm which performs a blind identification of the uplink scrambling code of a WCDMA transmission has been presented. It exploits the framing, spreading and multiplexing specificitions of the WCDMA standard and is thus dedicated to this system. The simulation results show that it is possible to obtain a reliable estimation at a SNR as small as -7 dB in an AWGN channel. Thus, the scrambling code of a strong interferer can be identified at the femto BS, and interference mitigation techniques can be implemented to cancel this interferer. The robustness of the algorithm to static multipath has also been validated.

[1] V. Chandrasekhar, J. Andrews, and A. Gatherer, "Femtocell networks: a survey," IEEE Communications Magazine, vol. 46, no. 9, pp. 59-67,

Fig. 3 .

 3 Fig. 3. Flow chart of the algorithm

  64, SF DPCCH = 256, β c = 1 and β d = 11/15.

Fig. 4 .Fig. 5 .

 45 Fig. 4. Sensitivity to the number of RGM -AWGN channel -N iter = 20 and M = 4000

Fig. 6 .

 6 Fig. 6. Performance with multipath channels -M = 4000


Mathieu Des Noes
email: mathieu.desnoes@cea.fr

Valentin Savin
email: valentin.savin@cea.fr

Jean Marc Brossier
email: jean-marc.brossier@gipsa-lab.grenoble-inp.fr

Laurent Ros
email: laurent.ros@gipsa-lab.grenoble-inp.fr

Blind Identification of the Scrambling Code of a Reverse Link CDMA2000 Transmission

Interference between macro and femtocells is an important issue for the development of CDMA2000 femtocell networks. More specifically, the reverse link signal of a macro User Equipment may generate an unacceptable level of interference at the femto Base Station. To avoid this situation, interference mitigation techniques could be implemented. All the proposed techniques require to know the state of the scrambling code of the interferer in the reverse link. Unfortunately, it depends on the code mask of the terminal which is unknown by the femto BS. The femto BS has to estimate blindly the state of the scrambling code. An algorithm which performs a blind identification of the scrambling code of a CDMA2000 reverse link transmission is proposed in this article. This gives the possibility to implement interference cancelation algorithm at the femto BS.

INTRODUCTION

The development of femtocell networks is an important perspective for increasing cellular systems capacity. A femtocell is covered by a small Base Station (BS), designed for typically indoor environments (home, business) which does not require a coordinated deployment [START_REF] Chandrasekhar | Femtocell networks: a survey[END_REF]. Although most of the current research activities on this topic focus on the LTE system [START_REF] Holma | LTE for UMTS: OFDMA and SC-FDMA based radio access[END_REF], current deployed systems use the UMTS-WCDMA or CDMA2000 technologies [START_REF] Holma | WCDMA for UMTS[END_REF]. In this paper, an algorithm which estimates blindly the state of the scrambling code of a terminal in the reverse link is proposed. It offers the perspective of implementing interference mitigation techniques at a femto Base Station (BS), and hence providing a solution to the problem of macro to femto interference issue in the reverse link of CDMA2000 femtocell networks. A macro BS provides the overall coverage, while femto BSs offer better indoor coverage to UEs attached to them. A femto BS can be configured to operate in open or closed access mode to visiting UEs [START_REF] Xia | Open vs. closed access femtocells in the uplink[END_REF]. In open access mode, a visiting UE is allowed to handover from the macro BS to a femto BS in order to send its data. This generates additional complexity to route the data packets and also to ensure communication security. In closed access mode, a visiting UE is not allowed to handover. This simplifies the network architecture, but this may lead to an unacceptable interference level at the femto BS. This situation occurs if a visiting UE transmits at high power, while it is located nearby the femto BS. This situation is depicted in Fig. 1. A macro UE transmits at high power because it is either at the cell edge or inside a building. Since the power control procedure concerns only the Fig. 1. Interference generated by a macro UE on a femto BS link with the macro BS, it will be received at the femto BS with a power much larger than the power of a femto UE. This is the well known near-far effect in CDMA systems [START_REF] Verdu | Multiuser detection[END_REF]. If the power of this interferer is too large, the femto BS may not be able to demodulate any communication with its attached UEs. This creates a "dead zone" in the network coverage. In order to avoid this situation, it is required to implement interference mitigation techniques. This subject has been deeply studied in the past two decades for CDMA systems [START_REF] Andrews | Interference cancellation for cellular systems: a contemporary overview[END_REF]. All the proposed techniques exploit the knowledge of the UE's scrambling code. More precisely, it is required to know the state of the UE's scrambling code generator at the beginning of each frame. The state of the scrambling code depends on its "code mask" which is used to generate a user specific delay of a long scrambling code common to all UEs. Unfortunately, in a closed access mode, there is no signalling link between the macro and femto BSs. Hence, a femto BS has no knowledge of the code mask of an interfering macro UE. It has to estimate the state of the scrambling code of a macro UE blindly. To the authors' knowledge, [START_REF] Kerr | Iterative signal processing for blind code phase acquisition of CDMA 1x signals for radio spectrum monitoring[END_REF] is the first and only article addressing this issue. The authors exploit the specificities of the framing, spreading and multiplexing procedures defined by CDMA systems to process the signal so that the resulting signal can be considered to be a linear code depending on the initial state of the scrambling code generator. This processing step can be applied to both WCDMA and CDMA2000 systems [START_REF].0 Spreading and modulation (FDD), 3GPP[END_REF] [START_REF] Gpp | -A -Physical Layer Standard for cdma[END_REF]. The received codeword is then estimated with a max-log-MAP algorithm [START_REF] Robertson | A comparison of optimal and sub-optimal MAP decoding algorithms operating in the log domain[END_REF]. The initial state of the different spreading code generators is eventually obtained by applying a pseudoinverse of the code generator matrix. Concerning the CDMA2000 system, the solution proposed in [START_REF] Kerr | Iterative signal processing for blind code phase acquisition of CDMA 1x signals for radio spectrum monitoring[END_REF] is adapted to a BPSK modulation which concerns Radio Configurations (RC) 3 to 6. However, it does not work for the 64-ary orthogonal modulation used by RC 1 and 2. These RCs are important because they are used by voice applications. This paper addresses the 64-ary orthogonal modulation case and also proposes a different decoding strategy. Exploiting the unique properties of m-sequence, it is shown that an iterative message-passing algorithm [START_REF] Yeung | An iterative algorithm and low complexity hardware architecture for fast acquisition of long PN codes in UWB systems[END_REF] can be implemented to decode the received signal after an initial processing step. Then the initial state of the scrambling code can be obtained with a proper used of the "Shift and add" and decimation properties of m-sequences. The paper is organized as follows. Section I presents the procedure used to generate the reverse link of CDMA2000 signals for RC 1 and 2. The construction of the complex scrambling code is also detailed since its specific properties will be exploited by the proposed algorithm. Section II details the blind identification algorithm, which is split in 3 steps. Section III presents simulation results and Section IV concludes this paper. Notation: a sequence will be written with upper case in its BPSK representation (S(k) ∈ {-1, +1}) and with lower case in its binary representation (s(k) ∈ {0, 1}).

I. REVERSE LINK SIGNAL GENERATION

In this section, the modulation and spreading operations implemented in the reverse link of the CDMA2000 system are first described. Then, the method standardized for constructing the scrambling code is detailed [START_REF] Gpp | -A -Physical Layer Standard for cdma[END_REF].

A. Modulation and spreading

The reverse link implements a Code Division Multiple Access (CDMA) scheme. Fig. 2 details the modulation and spreading operations for RC 1 and 2. The data are mapped on a 64-ary orthognal modulation, repeated, scrambled by the long code S n , multiplied by the complexe spreading sequence C I + jC Q and eventually filtered by the pulse shaping filter P (t).

C Q (k) C I (k) S n (k) j Tc/2 P(t) 64-ary orthogonal modulator Repetition W n (k) V(k)

Fig. 2. Modulation and spreading operations

The 64-ary orthogonal modulation consists in mapping 6 consecutive binary data d i , d i+1 , • • • , d i+5 onto a Walsh-Hadamard (WH) sequence of length 64 chips [START_REF] Gpp | -A -Physical Layer Standard for cdma[END_REF]. The WH sequence index is selected as follows:

m = d i + 2d i+1 + 4d i+2 + 8d i+3 + 16d i+4 + 32d i+5
Then each chip of the WH sequence W m is repeated 4 times and scrambled by the long code S n . The subscript n indicates that the long scrambling code is specific to user n. The scrambled signal corresponding to sequence W m is thus:

V (k) = S n (k)W m (⌊k/4⌋) (1)
⌊x⌋ is the largest integer smaller than x.

The chips of sequences S n and W m are anti-modal representation (+1,-1) of binary sequences s n and w m . Sequence V (k) is then multiplied by a complex spreading sequence C I + jC Q . These 2 sequences are built from msequences of period N = 2 15 -1 = 32767 chips [START_REF] Mceliece | Finite fields for computer scientists and engineers[END_REF]. Their characteristic polynomials are :

g I (D) = D 15 + D 13 + D 9 + D 8 + D 7 + D 5 + 1 g Q (D) = D 15 + D 12 + D 11 + D 10 + D 6 + D 5 + D 4 + D 3 + 1
The sequences are complemented with a '0' in order to create a frame of N = 32768 chips. According to the standard specifications, the initial state of each sequence is known by the receiver at the beginning of each frame. In this article, the gating procedure associated with power control is not considered. This corresponds to a gating rate equals to 1 in the standard [START_REF] Gpp | -A -Physical Layer Standard for cdma[END_REF]. A smaller gating rate (1/2 or 1/4) would require the implementation of parallel instances of the algorithm proposed in this article. The real part of the complex signal is filtered by the pulse shaping filter P (t), while the imaginary part is delayed by half a chip before pulse shaping. This implements an offset QPSK modulation. The specifications of P (t) are detailed in [START_REF] Gpp | -A -Physical Layer Standard for cdma[END_REF]. The transmitted signal can eventually be modeled as follows :

X(t) = k V (k)G k (t -kT c ) (2)
with

G k (t) = C I (k)P (t) + jC Q (k)P (t -T c /2
). T c is the chip period and is equal to 814 ns. G k (t) could be considered as a time-varying pulse shaping filter applied to sequence V (k).

B. Scrambling code generation

The scrambling code is built from a very long m-sequence a of period 2 42 chips. Fig. 3 shows the Linear Feedback Shift Register (LFSR) sequence generator according to the Galois representation [START_REF] Peterson | Introduction to spreadspectrum communications[END_REF]. Its characteristic polynomial is :

g s (D) = D 42 + D 35 + D 33 + D 31 + D 27 + D 26 + D 25 + D 22 + D 21 + D 19 + D 17 + D 16 + D 10 + D 7 + D 6 + D 5 + D 3 + D 2 + D + 1
The content of the i th shift register at time k is noted a i (k). The scrambling code is generated by the modulo-2 inner product of a 42-bit mask and the 42-bit state vector of the sequence generator:

s n (k) = r-1 i=0 a i (k)m n (i)
where ⊕ denotes the XOR operation. m n (i) is the code mask of the n th user. The vector a 0 (k), a 1 (k), • • • , a r-1 (k) is known by the receiver at some specified time since it is signalled periodically by the network. On the other hand, the code mask is unknown to the receiver and thus it cannot generate sequence s n with the schematic of Fig. 3. Due to the "Shift and add" property of m-sequences, s n is also a m-sequence with characteristic polynomial g S (D) [START_REF] Mceliece | Finite fields for computer scientists and engineers[END_REF]. It is a shifted version of sequence a: s n (k) = a(k + τ n ). τ n is specific to the n th user. It depends on the code mask and the characteristic polynomial. Using the Fibonacci representation, s n can be generated with the generator depicted in Fig. 4 [START_REF] Mceliece | Finite fields for computer scientists and engineers[END_REF]. Instead of estimating the code mask, knowing the initial state of sequence a, it is simpler to estimate the content of the shift registers of sequence s n with the Fibonacci representation. This is the objective of the algorithm proposed in this paper.

D D D g r-1 g r-2 g 1 + g 2 + + + D Code mask + S n (k) a 0 (k) a 1 (k) a r-2 (k) a r-1 (k)
D D D D D g r-1 g r-2 g 1 + + + + g 2 s n (k) s n (k+1) s n (k+2) s n (k+r) s n (k+r-1) s n (k+r-2)
Fig. 4. LFSR sequence with the Fibonacci representation

II. BLIND IDENTIFICATION OF THE SCRAMBLING CODE

The algorithm is split in 3 steps:

• step 1 : signal processing which aims at allowing the direct observation of a modified version of the sequence s n , denoted by s.

• step 2 : estimation of the initial state of sequence s with an iterative message-passing decoder.

• step 3 : determination of the initial state of the sequence s n with the use of a transposition matrix. In order to ease the comprehension of the algorithm, the description is restricted to an Additive White Gaussian Noise (AWGN) channel model. The robustness to multipath channel will be discussed in section II-E.

A. Signal processing

The received signal is down-converted to baseband and sampled at rate rate T = T c /E, where E is the oversampling factor. It is modeled as follows (see Eq. ( 2)):

R(q) = e jθ k V (k)G k (qT -kT c ) + N (q) ( 3
)
where θ is the phase rotation introduced by the channel and N (q) the additional noise modeling thermal noise as well as other sources of interference. The objective is to obtained an observation of a modified version of the sequence s n , and thus get rid off G k (t), θ and sequence W m . The pulse shaping filter defined by the standard generates inter-chip interference. As a consequence, a simple matched filter is not sufficient and an equalizer shall be implemented. Hence, the first step consists in implementing an equalizer to the time varying pulse shaping function G k (t). Let us define this equalizer as a FIR filter with L coefficients:

F k (0), • • • , F k (L -1).
It can be implemented either with a zero-forcing (ZF) or a minimum mean square error (MMSE) criterion [START_REF] Proakis | Digital communications[END_REF].

If the receiver is synchronized with the transmit signal, after equalization and downsampling, the signal is well approximated by :

Y (k) = L-1 l=0 R(kE + L/2 -l)F k (l) ≈ e jθ V (k) + N eq (k) (4)
N eq (k) is the filtered noise. If the receiver is not synchronized, the equalizer output can be regarded as a random noise sequence. This is due to the correlation properties of m-sequences c I and c Q [START_REF] Mceliece | Finite fields for computer scientists and engineers[END_REF]. The second step consists in removing the dependency upon sequence W m and phase θ. Since, the chips of sequence W m are repeated 4 times, a differential multiplication is performed within each window of 4 repeated chips. For instance, the differential multiplication using only the first two chips of the 4 chips length window will give :

Y (4k + 1)Y (4k) * = V (4k + 1)V (4k) + N dif f (k) = S n (4k + 1)S n (4k) + N dif f (k)
where N dif f (k) contains all the noise cross-product terms. In fact there are 7 possible differential products:

U i,j (k) = Y (4k + i + j)Y (4k + i) * = S n (4k + i + j)S n (4k + i) + N i,j (k) (5)
i = 0, 1, 2 and j = i + 1, • • • , 3. Let us define the sequence Si,j (k) = S n (4k + i + j)S n (4k + i), which is rewritten in GF(2) by : si,j (k) = s n (4k + i + j) ⊕ s n (4k + i). U i,j is thus a noisy observation of sequence Si,j . The objective is now to find the relationship between sequences si,j and s n .

To do so, two interesting properties of m-sequences will be exploited [START_REF] Mceliece | Finite fields for computer scientists and engineers[END_REF]:

• "Shift and add" property : for two given delays τ 1 and τ 2 , there exists a unique τ 3 such that:

x(k + τ 1 ) ⊕ x(k + τ 2 ) = x(k + τ 3 )
• Decimation property: the decimation by a factor 2 of a m-sequence gives a shifted version of this m-sequence. There exists a unique τ such that:

x(2k) = x(k + τ )
These 2 properties allow us to conclude that si,j is a shifted version of sequence s n . As a consequence, si,j has the same characteristic polynomial as s n : g s (D).

It is known from the literature that the initial state of a LFSR sequence can be estimated with a standard iterative messagepassing algorithm exploiting the sequence's characteristic polynomial [START_REF] Yeung | An iterative algorithm and low complexity hardware architecture for fast acquisition of long PN codes in UWB systems[END_REF][15] [START_REF] Mihaljevic | A low-complexity and highperformance algorithm for the fast correlation attack[END_REF]. It is thus possible to estimate the initial state of sequence si,j using its characteristic polynomial g s (D). The algorithm's principle is thus to decode the initial state of sequence si,j and to recover the state of sequence s n using a fixed state transposition matrix with sequence si,j . The successive steps of the algorithm are the followings. At each time instant q, the receiver assumes that it is synchronized with the beginning of the frame and equalizes the input signal according to Eq. ( 4). Then a differential multiplication is done, as in Eq. ( 5), to compute sequence U i,j . Since si,j is a shifted version of sequence s n , it can be decoded with an iterative message-passing algorithm which parity check matrix is matched to the characteristic polynomial of sequence s n (see section II-B). The vector (U i,j (0), • • • , U i,j (M -1)) feeds this decoder. If the decoder fails to find a codeword, this means that either the frame synchronization assumption is not valid or a missed detection happened. In both cases, the procedure is restarted when the next sample is received. If the decoder finds a valid codeword, the initial state of sequence s n is found by a simple matrix multiplication (see section II-C).

B. Iterative message-passing decoding

A m-sequence x, with characteristic polynomial g(D), satisfies the following parity check equation (g 0 = g r = 1), for all k ≥ 0:

r i=0 g r-i x(k + i) = 0
A m-sequence is thus a cyclic linear code with rate r 2 r -1 . A codeword is generated by one initial state of the shift registers. In the context of this paper, only a sequence of M variables, corresponding to M consecutive bits of the codeword, is observed. The parity check matrix of this code depends on the sequence's characteristic polynomial g(D):

H =        g r • • • g 0 0 • • • • • • 0 0 g r • • • g 0 0 • • • 0 . . . . . . . . . . . . . . . . . . . . . 0 • • • 0 g r • • • g 0 0 0 • • • • • • 0 g r • • • g 0        (6)
Once the parity check matrix has been defined, it is possible to decode the received message with a standard iterative message passing algorithm [START_REF] Kschischang | Factor graphs and the sum-product algorithm[END_REF] [START_REF] Wiberg | Codes and decoding on general graphs[END_REF]. In addition, as it was proposed in [START_REF] Yeung | An iterative algorithm and low complexity hardware architecture for fast acquisition of long PN codes in UWB systems[END_REF] [START_REF] Principe | Rapid acquisition of gold codes and related sequences using iterative message passing on redundant graphical models[END_REF], the use of Redundant Graphical Model (RGM) improves significantly the decoder performance. It relies on the following property : if g(D) is the sequence's characteristic polynomial, it satisfies g(D 2 n ) = g(D) 2 n . This property is exploited to create additional parity check equations. Polynomial g n (D) = g(D 2 n ) also generates a parity check matrix H n similar to H. Theses matrices can be concatenated to create a larger parity check matrix H RGM [START_REF] Yeung | An iterative algorithm and low complexity hardware architecture for fast acquisition of long PN codes in UWB systems[END_REF]:

H RGM =      H 0 H 1 . . . H n RGM -1      (7)
where n RGM is the number of RGMs used for decoding. These RGMs increase the column weight of the parity check matrix, while keeping the row weight constant. If the decoding is successful (all parity check equations are satisfied), the soft decision output of the decoder is convert into a binary representation with a hard decision rule. Then, according to the Fibonacci representation (Fig. 4), the first r bits of the codeword are the content of the shift registers at initialization.

C. Determination of the initial state of the sequence s n

The decoder provides a vector of r bits representing the initial state of sequence si,j , denoted by A s. We want to find the initial state of sequence s n at the beginning of the frame:

A sn = (s n (0), • • • , s n (41))
T . This task is achieved in 2 steps.

The first step consists in finding the initial state of the sequence s decim , A s decim , which gives si,j by decimation by a factor 4: si,j (k) = s decim (4k) [START_REF] Arazi | Decimation of m-sequences leading to any desired phase shift[END_REF] shows that there is a fixed transposition matrix B between the state vectors of a m-sequence and its decimated by 2 version, and detailed the procedure to compute this matrix B. As a consequence:

A s decim = B 2 A s
The second step eventually provides vector A sn . It exploits the "shift and add" relation between the sequences s decim and s n :

s decim (k) = s n (k + i + j) ⊕ s n (k + i)
Using the state transition matrix G, defined by [START_REF] Peterson | Introduction to spreadspectrum communications[END_REF]:

G =          0 1 0 0 • • • 0 0 0 1 0 • • • 0 0 0 0 1 • • • 0 . . . . . . . . . 0 • • • . . . 0 0 0 0 • • • 1 1 g r-1 g r-2 g r-3 • • • g 1          (8)
we have:

A sn = (G i+j + G i ) -1 A s decim
These two steps are finally combined in a single transposition matrix T :

A sn = T A s T = (G i+j + G i ) -1 B 2
It is important to note that matrix T could be computed offline and stored in memory.

D. Flow chart of the algorithm

A flow chart of the algorithm is presented in Fig. 5. The decoding operation must be implemented at the chip rate, which requires a very high data rate decoding capability.

E. Robustness to a multipath channel

The algorithm is by construction robust to multipath whose delays are larger than a chip period. This is due to the chip level equalization of Eq. 4. A multipath component delayed by more than one chip will be scrambled by sequences c I and c Q and will be considered as an additional noise source by the decoder. In order to cope with multipath, the receiver implements a search window of length W chips. At each time q, the receiver runs the detection algorithm for time q to q + W -1. If a valid codeword is detected within this window, the detection is declared successful.

III. SIMULATION RESULTS

Let us first define the conventional synchronization hypothesis :

• H 0 : the receiver is not synchronized with the beginning of the frame.

• H 1 : the receiver is synchronized with the beginning of the frame. The performance of the algorithm are measured by the probabilities of correct detection P d , false alarm P f a and missed detection P m , defined as follows: Âxn is the estimated initial state of sequence x n , given by the decoder output. Performances are measured with the following simulation configurations:

• When measuring P m , the receiver is synchronized with the beginning of the frame (i.e. hypothesis H 1 is satisfied), while its is not synchronized when P F A is evaluated.

• The decoder implements a Self-Corrected Min-Sum (SCMS) message-passing algorithm [START_REF] Savin | Self-corrected min-sum decoding of ldpc codes[END_REF], which provides a very high data rates decoding capability, while performing close to the "optimal" Belief Propagation (Sum-Product) decoding. The decoder stops when either all the parity check equations are satisfied or the maximum number of iteration N iter is reached.

• The number of variable is M = 1500 at the decoder input. This value is selected so that it is smaller than the size of a Power Control Group (PCG) which contains 1536 chips [START_REF] Gpp | -A -Physical Layer Standard for cdma[END_REF]. A PCG is equivalent to a slot in the WCDMA system.

• The channel is AWGN. Simulations over 10 7 frames did not produce any false alarm. This means that P f a < 10 -6 with a good confidence level. Fig. 6 shows the probability of missed detection for a chip level MMSE, ZF or MRC equalizer. The number of RGMs is set to its largest possible value n RGM = 6 for M = 1500 variables and N iter = 100. P m is measured as a function of the Signal to Noise Ratio (SNR) at the input of the receiver. At P m = 0.1, the MMSE equalizer gives a gain of almost 1 dB and 3 dB with respect to ZF and MRC. Fig. 7 shows the probability of missed detection as a function of the number of RGMs for a receiver configuration with chip level MMSE equalizer. It is not worth increasing the number of RGMs above 4.

According to [START_REF] Padovani | The application of spread spectrum to PCS has become a reality: Reverse Link Performance of 1S-95 Based cellular systems[END_REF], the target Eb/No at a base station is around 7 dB for the 9.6 kbps service which corresponds to our simulation configuration (all the PCGs are transmitted). Since the combined coding-spreading gain equals 128, the target SNR at a base station is around -14 dB. According to Fig. 7, the algorithm implementing a chip level MMSE equalizer detects the scrambling code of an interfering user at a SNR of 2.5 dB. This corresponds to a very strong interferer since its received power is 16.5 dB larger than a desired user attached to the femtocell. However,the performance of the algorithm is not good enough to detect certain lower but still strong interferes. For instance, an interferer which power is 15 dB larger than the desired user will not be detected even it will have very damaging effects. The bad detection performance of the algorithm is due to the weight of the characteristic polynomial of the scrambling code generator. This effect has been well studied in the framework of fast correlation attacks [START_REF] Mihaljevic | A low-complexity and highperformance algorithm for the fast correlation attack[END_REF]. The weight equals 20 while it should be below 4 or 5 to give good decoding performance [START_REF] Canteaut | Improved fast correlation attacks using parity-check equations of weight 4 and 5[END_REF]. One way to improve performance is to generate a parity check matrix which row weight is below 4 or 5. This could be achieved by implementing the technique proposed in [START_REF] Chose | Fast correlation attacks: An algorithmic point of view[END_REF].

IV. CONCLUSION

An algorithm which performs a blind identification of the scrambling code of a reverse link CDMA2000 transmission has been presented. It is dedicated to the 64-ary orthogonal modulation used by the radio configurations 1 and 2 of the standard. The simulation results show that with a BP decoding algorithm it is possible to obtain a reliable estimation at a SNR of 2.5 dB in an AWGN channel. It is also noticeable that the false alarm probability is below 10 -6 for all the tested SNR values.

Fig. 3 .

 3 Fig. 3. Long scrambling code generation

Fig. 5 .

 5 Fig. 5. Flow chart of the algorithm

P

  d = P (Ic = 1 and Âxn = A xn |H 1 ) P f a = P (Ic = 1|H 0 ) P m = 1 -P d Ic is the indication function of the decoder: Ic = 1 if all parity check equations are satisfied 0 otherwise

Fig. 6 . 6 1

 66 Fig. 6. Influence of the equalizer -AWGN channel -N iter = 100, M = 1500 and n RGM = 6

Fig. 7 .

 7 Fig. 7.Sensitivity to the number of RGM -AWGN channel -MMSE equalizer , N iter = 100 and M = 1500