Member, IEEE Philippe Benabes

Accurate time-domain simulation of continuous-time sigma-delta modulators

Keywords: sigma-delta, simulation, analog-to-digital conversion, continuous-time

In this paper we present a methodology for the simulation of continuous-time sigma-delta converters. This method, based on a fixed-step algorithm, permits not only a time-domain simulation of the modulator output but also the simulation of intermediary signals. The method is based on the discretization of the continuous-time models and the use of a discrete simulator such as Simulink, which is more efficient than an analog simulator. By using filters with a sampling frequency higher than the modulator output frequency, the model can simulate input signals with a bandwidth which is higher than half the modulator sampling frequency. The transformation is exact in terms of Noise Transfer Function and asymptotically exact in terms of Signal Transfer Function (the Transfer Function from the modulator input to each stage filter output rapidly tends to the continuous-time model transfer function when the number of steps increases).

I. INTRODUCTION

Sigma-Delta (Σ∆) circuits are very attractive analog-todigital converters [START_REF] Schreier | Understanding Delta-Sigma Data Converters[END_REF] as they achieve high accuracy with few critical analog components. The speed and integrator bandwidth limitations of switched-capacitor Discrete-Time (DT) modulators can be mitigated by the use of Continuous-Time (CT) modulators [START_REF] Benabes | A methodology for designing continuous-time sigma-delta modulators[END_REF] [START_REF] Ortmanns | Continuous-Time Sigma-Delta A/D Conversion[END_REF]. The latter are not easy to integrate, but possess one key advantage over their discretetime competitors: no sampling is performed within the filter itself, allowing central frequencies in the Gigahertz range [START_REF] Naderi | On the Design of Undersampling Continuous-Time Band-Pass Delta-Sigma Modulators for Gigahertz Frequency A/D Conversion[END_REF]. An implicit anti-aliasing is performed by the analog filter [START_REF] Keller | On the Implicit Anti-Aliasing Feature of Continuous-Time Cascaded SigmaDelta Modulators[END_REF], thus reducing the constraints at the input of the converter. On the other hand, CT circuits are more difficult to design and to simulate than DT circuits. Classical design methodologies usually start with high-level design and simulations before transistor-level simulations [START_REF] Yantchev | Propagation characteristics of the fundamental symmetric Lamb wave in thin aluminum nitride membranes with infinite gratings[END_REF] [START_REF] Francken | A High-Level Simulation and Synthesis Environment for Σ∆ Modulators[END_REF]. The behavior of electronic components can be modeled at high level for example to analyze the impacts of environmental and process variations such as in [START_REF] Suarez | Behavioral Modeling Methods for Switched-Capacitor Σ∆ Modulators[END_REF]. When the high level-simulations are performed by an analog simulator, they require a very long computational time. An analytical integration of the differential equations can be performed such as in [START_REF] Guo | Efficient Look-Up-Table-Based Modeling for Robust Design of Σ∆ ADCs[END_REF], but this method is exact only if the input signal can be expressed in an analytical form. An alternative but equivalent methodology consists in using an equivalent DT model of CT modulators [START_REF] Benabes | A methodology for designing continuous-time sigma-delta modulators[END_REF]. Such a method supposes that the input signal of the CT filters has a particular shape (for example, a sampled and held signal) and

The author is with the Department of Signal Processing and Electronics, Supélec, 91190 Gif/Yvette, FRANCE (e-mail: philippe.benabes@supelec.fr).

thus the CT differential equations are solved for this particular signal. The output of the filter is calculated only at the A.D.C sampling times. However, the input bandwidth is limited to half the sampling frequency. As can be seen in [START_REF] Vogels | Efficient Time-Domain Simulation of Continuous-Time Sigma Delta A/D Converters Using Analytical Integration[END_REF], we can use a state space representation of the filters, or alternatively we can use direct s transfer functions as seen in [START_REF] Keller | A method for the Discrete-Time simulation of Continuous-time Sigma-Delta Modulators[END_REF]. In this paper, an extension of the simulation method of CT modulators based on Oversampled Discrete-Time (ODT) models [START_REF] Sobot | Tunable Continuous-Time Bandpass Σ∆ Modulators With Fractional Delays[END_REF] is presented. With this method, each sampling period is divided into a fixed number of steps. The outputs of the filters are then calculated for each sub-sampling time. The advantage of DT simulations over CT algorithms is that the outputs of the filters are calculated only at each sub-sampling time. CT algorithms (especially those with a variable step) calculate more samples than required, thus increasing the simulation time.

This transformation is exact in terms of Noise Transfer Function (NTF) and asymptotically exact in terms of Signal Transfer Function (STF). The STF of the model rapidly tends to the STF of the CT model when the number of steps increases. Simulations of the response of signals with a bandwidth higher than half the sampling frequency is then possible. Furthermore, an estimation of the ideal step-size derived from the bandwidth of the input signal will be performed.

This paper is structured as follows: Section II describes both, the synthesis and the analysis method for continuous filter modulators. An application of the simulation of a high order modulator with extra loop delay is illustrated in Section III. The extension to the simulation of intermediary signals is explained in section IV. Finally, concluding remarks are given in section V.

II. SYNTHESIS AND SIMULATION METHODS OF CONTINUOUS-TIME MODULATORS

A. Equivalency between Continuous-Time and Discrete-Time Filters

The relationship between the CT g(s) and the DT filter transfer function f (z) of a Σ∆ modulator can be expressed using the well known formula [START_REF] Keller | A method for the Discrete-Time simulation of Continuous-time Sigma-Delta Modulators[END_REF], [START_REF] Benabes | Fixed-Step Simulation of Continuous-Time Sigma Delta Modulators[END_REF],

f (z) = (1 -z -1 )Z t=nTs L -1 g(s)e -dTss s ( 1
)
where d is the loop delay between the A.D.C input and the D.A.C output (normalized to the sampling period), Z stands for the Z transform, L -1 is the inverse Laplace Transform, and T s is the sampling period. This formula ensures the equivalency of the Noise Transfer Function (NTF) between the CT modulator shown in Fig. 1 and the DT equivalent topology shown in Fig. 2. This formula is valid for square feedback signals but can be extended to any signal shape, see [START_REF] Gao | Excess loop delay effects in continuous-time delta-sigma modulators and the compensation solution[END_REF] or [START_REF] Oliaei | Design of continuous-time sigma-delta modulators with arbitrary feedback waveform[END_REF]. Formula (1) can be inverted and used for the synthesis of a CT modulator, i.e. in order to obtain the CT filter transfer function that ensures that the NTF of both modulators will be the same. When the delay d is non zero, (1) has no solution as a degree of freedom is missing, making it impossible to map the transfer function. An optimization method such as in [START_REF] Aboushady | Systematic Approach for Discrete-Time to Continuous-Time Transformation of sigma-delta Modulators[END_REF] can be used to optimize the coefficients of the CT filters in order to maintain the stability and the performance of the loop. As explained in [START_REF] Loeda | On the Design of High-Performance Wide-Band Continuous-Time Sigma-Delta Converters Using Numerical Optimization[END_REF], in order to ensure the equivalency, a feedback term can also be added between D.A.C output and A.D.C input. All feedback techniques are summarized in [START_REF] Yahia | A new technique to compensate the influence of the feedback DAC delay in continuous-time bandpass converters[END_REF].

With such a feedback term, (1) can now be written as:

f (z) = (1 -z -1 )Z t=nTs L -1 g(s)e -dTss s + T (z) (2)
Using this methodology, the excess loop-delay [START_REF] Benabes | Fixed-Step Simulation of Continuous-Time Sigma Delta Modulators[END_REF] is no longer a hindrance, but becomes simply a parameter that can be adjusted [START_REF] Keller | A Comparative Study on Excess Loop Delay Compensation Techniques for Continuous-Time Sigma-Delta Modulators[END_REF], creating a compromise between the gain margin of the loop and the feasibility of the feedback D.A.C. The proposed method in [START_REF] Yahia | Bandpass sigma-delta modulators synthesis with high loop delay[END_REF] consists in decomposing f (z) into rational fractions, and converting each term using precalculated solutions of (2). To obtain these solutions, we used Maple c software, as shown in annex A. The degree of T (z) is the integer part of d plus one [START_REF] Benabes | Synthesis and analysis of sigma-delta modulators employing continuous-time filters[END_REF].

Table I gives solutions of (2) for several transfer functions and delays between 0 and 2. For the third-order, solutions may be found in [START_REF] Sobot | Tunable Continuous-Time Bandpass Σ∆ Modulators With Fractional Delays[END_REF] and [START_REF] Benabes | Influence of the feedback D.A.C delay on a continuous-time bandpass sigma-delta converter[END_REF]. They may also be achieved when using the Maple code given in annex A.

B. Signal transfer function

Assuming that the input signal is a band-limited signal (limited to the half of the sampling frequency), and that the quantizer can be modeled by an additive white noise, the signal transfer function of the DT topology can be expressed as:

ST F DT (ϕ) = f (e 2jπ ϕ F S ) 1 + f (e 2jπ ϕ F S ) (3)
where F s is the sampling frequency, and ϕ the signal frequency. The signal transfer function of the CT topology is [START_REF] Kuo | Digital control systems. 2nd Edition[END_REF]:

ST F CT (ϕ) = g(2jπϕ) 1 + f (e 2jπ ϕ F S ) (4)
f being related to g by equation [START_REF] Schreier | Understanding Delta-Sigma Data Converters[END_REF]. It is clear that even if the DT topology can model the NTF of the modulator correctly, the STF obtained by (3) is different from (4). For example, with a first-order lowpass modulator, where

f (z) = z -1
1-z -1 , and g(s) = FS s , the DT signal transfer function is a simple delay:

ST F DT (z) = z -1 ⇒ ST F DT (ϕ) = e -2jπ ϕ F S . (5)
This transfer function clearly differs from the CT transfer function which is:

ST F CT (ϕ) = 1 -e 2jπ ϕ F S 2jπ ϕ FS . (6)
As a consequence, when we deal with non constant signals, the DT model cannot be used to perform exact time-domain simulations of the behavior of a CT modulator. Furthermore, the DT model is unable to model an input signal with a frequency higher than half the sampling frequency. In order to enhance the signal-transfer function and remove the frequency limitation, we propose using an oversampled model of the DT modulator (ODT)

C. Oversampled model of a sigma delta modulator

Let us now consider the oversampled model of a Σ∆ modulator shown in Fig. 3. The sampling frequency of the A.D.C is still F s , but the digital filter (F ) runs now at kF s . The feedback signal is held during k samples. This is a special case of multi-rate systems as seen in [START_REF] Benabes | Fast simulation of bandpass Continuous-Time Sigma delta Modulators[END_REF]. In order to simplify the notations, the Z variable denotes functions running at frequency kF s , while the z variable denotes a function running at frequency F s . A delay r is introduced at the output of the zero-order hold, in order to model the delay d between A.D.C input and D.A.C output. r will be chosen as the integer part of k.d (r = ⌊k.d⌋, where ⌊⌋ denotes the integer part).

Next we consider the transfer function between the A.D.C output and its input. We calculate the response of the filter Y * [n] to a discrete impulse at the A.D.C output in the case of both: CT modulator, and ODT modulator.

In the case of the CT modulator:

g(s) d f (z) T (z) 1 Tss 0 ≤ d ≤ 1 z -1 1-z -1 dz -1 1 Tss 1 ≤ d ≤ 2 z -1 1-z -1 z -1 + (d -1)z -2 1 Tss-a 0 ≤ d ≤ 1 e -ad e a -1 a z -1 1-e a z -1 1-e -ad a z -1 1 Tss-a 1 ≤ d ≤ 2 e -ad e a -1 a z -1 1-e a z -1 e -ad (e a -1) a z -1 + 1-e -a(d-1) a z -2 1 (Tss) 2 0 ≤ d ≤ 1 1-2d 2 z -1 1-z -1 + z -1 1-z -1 2 -d 2 2 z -1 1 (Tss) 2 1 ≤ d ≤ 2 1-2d 2 z -1 1-z -1 + z -1 1-z -1 2 1-2d 2 z -1 -(d-1) 2 2 z -2 TABLE I CT TO DT TRANSFORMATION TABLE Fig. 3. Oversampled Σ∆ modulator Y * [n] = L -1 t=nTs [b(s)g(s)] with b(s) = 1-e -Tss s e -dTss (7)
In the case of the ODT modulator:

Y * [n] = Z -1 N =kn [B(Z)F (Z)] with B(Z) = Z -r (1-Z -k ) k(1-Z -1 ) (8)
These two responses should be equal for each A.D.C sampling time at frequency F s . A sufficient condition is that the response to a discrete impulse of F (Z) is the same, for each sampling time at kF s , as the response of g(s) to an impulse with a width equal to T s /k delayed by d-r/k. Indeed, an impulse with a width of T s is the sum of k impulses with a width of T s /k delayed by multiples of T s /k.

This condition can be expressed as:

F (Z) = (1 -Z -1 )Z t=nTs/k L -1 g(s)e -(d-r/k)Tss s (9)
We denote D = kd -r. It can be clearly seen from its definition that 0 ≤ D < 1. (9) can be rewritten as:

F (Z) = (1 -Z -1 )Z t=nTs/k L -1 g(s)e -D Ts k s s (10)
If we replace D by d and T s /k by T s , this equation is the same as [START_REF] Schreier | Understanding Delta-Sigma Data Converters[END_REF]. Table II gives solutions of (10), i.e. shows how to obtain F (z) from g(s) for first-order and second-order filters.

The transformation will be denoted as ODTT, Oversampled Discrete-Time Transform, and we will write it as:

g(s) F (Z) 1 Ts k s Z -1 1-Z -1 -DZ -1 1 Ts k s-a e -ad e a -1 a Z -1 1-e a Z -1 -1-e -aD a Z -1 1 ( Ts k s) 2 1-2D 2 Z -1 1-Z -1 + Z -1 1-Z -1 2 + D 2 2 Z -1 1 ( Ts k s-a) 2 e -aD (1+aD-e a (a(1-D)-1))Z -1 a 2 (1-e a Z -1 ) + (e a -1)e a(1-D) Z -2 a 2 (1-e a Z -1 ) 2 + 1-(1+aD)e -aD a 2 Z -1 TABLE II ODT TRANSFORMATION TABLE F (Z) = ODT T (g(s), k, D)
where k is the oversampling ratio and D the difference between the real delay and the modeled delay normalized to T s /k. We can see that in case of k = 1, this transformation is the same as the transformation of (1), i.e. f (z) = ODT T (g(s), 1, d).

III. APPLICATION TO THE SIMULATION OF HIGH-ORDER MODULATORS WITH LOOP DELAY

The methodology proposed in section II will be used to simulate the time-domain response of a CT Σ∆ modulator. In this section we try to obtain the output signal of the modulator. We will evaluate the error committed by this model in terms of Signal Transfer Function (STF).

A. Simulation of the output of a modulator

A CT Σ∆ modulator can always be expressed as Fig. 4, where g x (s) is the transfer function between the modulator input and the A.D.C input, and g q (s) is the transfer function between the D.A.C output and the A.D.C input.

With the definition r = ⌊kd⌋, and D = kd -r, we use the ODTT transformation as:

F x (Z) = ODT T (g x (s), k, 0) (11)
F q (Z) = ODT T (g q (s), k, D) (12)
The ODT topology is shown in Fig. 5 If the input signal was sampled and held at frequency kF s , then both modulators would be completely equivalent. Let us take as an example a classical second-order modulator. The DT prototype is given in Fig. 6.

In this modulator,

f (z) = z -1 1 -z -1 + 1 2 z -1 1 -z -1 2 . ( 13
)
A CT equivalent topology can be found using table I. With a loop delay d between 0 and 1, we obtain:

g(s) = 1 2(T s s) 2 + 2d + 3 4T s s , (14)
T (z) = d(3 + d) 4 z -1 . ( 15
)
The prototype of the CT modulator is given in Fig. 7. This modulator is equivalent to the Fig. 6 modulator in terms of NTF, but the STF will have to be evaluated separately. It should be noted that the proposed methodology does not give the signal coefficients x 1 and x 2 . They will be calculated using STF considerations.

In order to simulate the CT modulator, we transform the CT modulator into an ODT topology. This operation takes into account the feedback signal and the input signal of the modulator. For this example, we choose d = 0.5. With this value,

g q (s) = 1 2(T s s) 2 + 1 T s s , (16)
and

g x (s) = x 1 1 2(T s s) 2 + x 2 1 T s s . ( 17
)
In order to simplify the example, the oversampling ratio k will be chosen as an even number. As a result, r = k/2, and D = 0. Using [START_REF] Keller | A method for the Discrete-Time simulation of Continuous-time Sigma-Delta Modulators[END_REF], and ( 12), then

F q (Z) = 1 + 4k 4k 1 k Z -1 1 -Z -1 + 1 2 1 k Z -1 1 -Z -1 2 , (18)
F x (Z) = x 1 + 4kx 2 4k 1 k Z -1 1 -Z -1 + x 1 2 1 k Z -1 1 -Z -1 2 . ( 19
)
The resulting ODT topology is shown in Fig. 8. In order to verify the theory, we have simulated the behavior of the second order DT modulator shown in Fig. 6, the equivalent CT modulator shown in Fig. 7, and the ODT model shown in Fig. 8, using Simulink software. The chosen OSR is equal to k = 8. Only the ten first output samples were plotted so as to make the figures more readable. In a first simulation, the input signal is a constant signal equal to 0 and the A.D.C inputs are compared in Fig. 9. DacM ax is the D.A.C maximum output voltage. It can be seen that these signals are equal for each sample time nT s . Furthermore, the inputs of the comparators of the CT and the ODT modulator are equal at each sub-sampling time nT s /k.

We have then applied a low frequency sine input signal sampled and held at frequency kF s and verified in Fig. 10 that the comparator inputs are still equal for each sub-sampling time for both the CT model and the ODT model. It can be verified that the DT model does not give the right response.

At last, we have applied a low frequency CT sine input signal and verified in Fig. 11 that the comparator inputs differ increasingly for the CT and ODT case. However, the fact that the STF of the CT modulator is correctly approached by the one of the ODT modulator at low frequencies will be verified in the next section.

B. Signal transfer function evaluation

Assuming that the quantizer can be modeled by an additive white noise, the signal transfer function of the CT topology can be expressed as:

ST F CT (ϕ) = g x (2jπϕ) 1 + f (e 2jπ ϕ Fs ) (20)
The STF can be used even for frequencies higher than F s /2, knowing that a signal at frequency ϕ is aliased into a signal at frequency ϕ -mF s at the modulator output (where -Fs 2 < ϕ -mF s < Fs 2 ). The signal transfer function of the ODT topology can be expressed as:

ST F ODT (ϕ) = F x (e 2jπ ϕ kFs ) 1 + f (e 2jπ ϕ Fs ) , (21)
with -kFs 2 < ϕ < kFs 2 . The error introduced by using the ODT model instead of the CT model can be evaluated using the ratio:

err = ST F ODT (ϕ) ST F CT (ϕ) = F x (e 2jπ ϕ kFs ) g x (2jπϕ) (22)
we denote φ = ϕ kFs . When ϕ varies from -kFs 2 to kFs 2 , then φ varies from -0.5 to 0.5, and

err(φ) = F x (e 2jπφ ) g x (2jπφkF s ) . ( 23
)
This error depends on the order of the filter, k, and its poles. In the case of integrators, the error does not depend on k.

For a first-order integrator,

err(φ) = 2jπφ e 2jπφ -1 , ( 24
)
for a second-order integrator,

err(φ) = (2jπφ) 2 e 2jπφ + 1 2 (e 2jπφ -1) 2 , ( 25
)
and for a third-order integrator,

err(φ) = (2jπφ) 3 e 2jπφ 2 + 4e 2jπφ + 1
6 (e 2jπφ -1)

3 . ( 26
)
These results were obtained from table II.

A pre-distortion filter may be introduced in order to compensate this STF distortion. This filter will not flatten the error from -kFs 2 to kFs 2 , but it will work efficiently from -4 to kFs 4 . The filter must also be a symmetrical FIR filter in order to maintain a linear phase correction (simple delay). Fig. 12 compares the STF modulus before and after correction for a ST F distortion with and without correction for a second-order integrator first-order integrator and Fig. 13 for a second-order integrator (ϕ is denoted as phi on the figures).

The most simple usable FIR (2 nd order) can be expressed as, for a n th order integrator by

C n (z) = λ n + (2 -λ n )Z -1 + λ n Z -2 (27)
with λ 1 ≈ 0.045, λ 2 ≈ -0.045, λ 3 ≈ -0.045.

It can be verified that the coefficient k has not appeared in the previous expressions.

C. Example of a second-order modulator

As an application, we can consider the modulator shown in Fig. 7. An equivalent bandpass example can be found in [START_REF] Benabes | Influence of the feedback D.A.C delay on a continuous-time bandpass sigma-delta converter[END_REF]. The model shown in Fig. 8 must be changed in order to introduce the STF correction filters, resulting in Fig. 14 topology. The STF of this modulator has been evaluated and compared with the one obtained for the ODT model in Fig. 15. The curves are given for the case x 1 = 1 and x 2 = 0. In order to evaluate the efficiency of our methodology, this STF could be compared with the one that would be obtained by making a bilinear transform of the filters of the CT modulator. It was shown in [START_REF] Sobot | Tunable Continuous-Time Bandpass Σ∆ Modulators With Fractional Delays[END_REF] that this topology is not strictly equivalent in terms of NTF. In the case of bandpass modulators, there would be a shift in the central frequency of the modulator [START_REF] Benabes | Influence of the feedback D.A.C delay on a continuous-time bandpass sigma-delta converter[END_REF]. In order to obtain a good NTF approximation, an oversampling ratio of at least 64 or 128 would be required. Using the ODT model, the NTF is correct for any value of k.

Table III compares the simulation times of a second-order lowpass modulator for 100000 output samples in the simulink environment. The first model is the one of Fig. 6 with two CT integrators. The algorithm used to solve the equations is ODE45. The DT model is the model of Fig. 7. The third and fourth models (bilinear with k=8 and k=128) are obtained by replacing the CT integrators of Fig. 7 by integrators obtained by a bilinear transform. The last model (ODT) is the one of Fig. 8. The fastest method is the one that uses the DT model (k = 1) but it does not provide accurate results in terms of STF. The bilinear method provides good results for k = 128 but implies long simulation times. Our methodology provides good results in terms of NTF and STF with a reasonable simulation time compared to the simulation time required by a continuous algorithm. This kind of results can be verified in [START_REF] Vogels | Efficient Time-Domain Simulation of Continuous-Time Sigma Delta A/D Converters Using Analytical Integration[END_REF]. The advantage of our methodology is to be able to simulate the response to signals with large bandwidth by adjusting k.

IV. INTERMEDIARY SIGNALS SIMULATION

A. A new ODT model taking into account intermediary signals

The goal of the previous section was to simulate the output of the modulator as accurately as possible. This methodology can be extended to the simulation of intermediary signals, especially the outputs of each filter. Indeed, the previous model provides a correct output value of the modulator in the case of a constant or sampled input signal, but the output of intermediary filters does not fit the real output of the CT filters. The methodology proposed in section III is now extended to intermediary signals. If the input signal was a sample and held signal at frequency kF s , then the methodology would provide the exact intermediary values at every sub-sampling time. For a CT input signal, the model gives an approximation of the outputs of all filters. Note that in this paper, we consider only single-loop topologies, but the methodology could be extended to any kind of topology.

Fig. 16 shows the topology of a n th order CT modulator: We denote the transfer functions between D.A.C output and the m th filter output as g m q . From Fig. 16, it can be seen that

g 1 q (s) = k 1 G 1 (s) g 2 q (s) = G 2 (s)(k 2 + k 1 G 1 (s)) ... g n q (s) = G n (s)(k n + G n-1 (s)(k n-1 + G n-2 (s)(...))) (28)
We also denote the transfer function between the modulator input and the m th filter output as g m

x . It can also be seen that

g 1 x (s) = x 1 G 1 (s) g 2 x (s) = G 2 (s)(x 2 + x 1 G 1 (s)) ... g n x (s) = G n (s)(x n + G n-1 (s)(x n-1 + G n-2 (s)(...))) (29
) Each transfer function can be transformed into its ODT equivalent by

F m x = ODT T (g m x , k, 0), m = 1..n F m q = ODT T (g m q , k, D), m = 1..n (30)
The output of F m q would be equal to the output of filter G m if it was connected only to the modulator feedback through k m . The output of F m

x would be an approximation of the output of filter G m if it was connected only to the modulator input through x m .

The resulting topology is given in Fig. 17. This topology provides the output of each filter separately. It is, of course, functionally correct. However, the quantity of performed calculations is not optimal as each CT filter output requires a separate DT filter. . As a consequence F m x and F m q have the same denominator, and the latter divides the denominator of F m+1 x and F m+1 q . F m

x and F m q can be expressed as

F m x (Z) = T m x (Z) m p=1 Dp(Z) F m q (Z) = T m q (Z) m p=1 Dp(Z) (31)
In the case of lowpass modulators, the terms D p (Z) are first-order filters; in the case of bandpass modulators, they are second-order filters.

The proposed topology used for the simulation is shown in Fig. 18. The integrators are replaced by their denominators and the coefficients are replaced by the numerators of the filters.

The equivalency is ensured if:

F m x (Z) = N m x (Z)+N m F (Z) ( N m-1 x (Z)+N m-1 F (Z)(...) ) D m-1 (Z) Dm(Z) F m q (Z) = N m q (Z)+N m F (Z) ( N m-1 q (Z)+N m-1 F (Z)(...) ) D m-1 (Z) Dm(Z) (32)
equations will be written recursively:

T m x (Z) m p=1 Dp(Z) = 1 Dm(Z)    N m F (Z) T m-1 x (Z) m-1 p=1 Dp(Z) + N m x (Z)     T m q (Z) m p=1 Dp(Z) = 1 Dm(Z)     N m F (Z) T m-1 q (Z) m-1 p=1 Dp(Z) + N m q (Z)     (33) resulting in: T m x (Z) = N m F (Z)T m-1 x (Z) + N m x (Z) m-1 p=1 D p (Z) T m q (Z) = N m F (Z)T m-1 q (Z) + N m q (Z) m-1 p=1 D p (Z) (34)
In case m = 1 (first stage), the equations (34) will result in:

T 1 x (Z) = N 1 x (Z) T 1 q (Z) = N 1 q (Z) (35)
Equations ( 34) are linear equations where the unknown parameters are the coefficients of the polynomials N m F , N m x , N m q . If the delay r is chosen correctly so that 0 < D = kd -r < 1, then the order of T m x (Z) and T m q (Z) is the same as the order of m-1 p=1 D p (Z), and the constant term is zero (this can be seen from table II). It is also clear from (34) that the constant term of N m F , N m x and N m q is zero. The degree of the resulting solution depends on the degree of the numerators and denominators. It should be noticed that this system always provides a solution. Indeed, increasing the degree of N m F , N m x and N m q by one adds three unknown parameters and only two equations. Consequently, there is a minimal degree for which the number of unknown parameters is at least equal to the number of equations. Our goal is to obtain the minimum order solution of these equations.

If we take into account that the system (34) usually has more equations than unknown terms, then the proposed algorithm can be expressed as follows:

• choose the minimum order for N m F , N m x and N m q , • solve the system (34) in a mean square sense,

• verify if the solutions provided by the above are valid,

• increase the order of N m F , N m x and N m q by one and start again if the equation system is not solved. The previous loop will end when the number of unknown parameters is equal to the number of equations. In the case of lowpass modulators, the degree of N m F , N m x and N m q will be m-1 and the degree of N 0

x and N 0 q is 1. It can be noticed that for orders higher than 2, the terms N m F , N m x and N m q are no longer simple coefficients as it is the case in usual models but when m > 2 they are polynomials. This is why a third-order modulator will be considered in the next example. In the case of bandpass modulators the degree of N m F , N m x and N m q is usually 2 for the first three stages. It will be 3 for the last stage.

B. Intermediary signal transfer functions evaluation

As in section III-B, it is possible to define a transfer function between the modulator input and the outputs of the filters. We consider a sine input signal at frequency ϕ. This signal is aliased into a signal at frequency ϕ -mF s at the modulator output. Due to the sampling function, the feedback signal contains components at frequencies ϕ -mF s + kF s . We consider in this evaluation the component for which k = m, i.e. a component at frequency ϕ. The output of filter k is due to 2 terms: one coming from from the input signal x and another coming from the output signal q

For the CT model (Fig. 16), the global transfer function from input to filter k output, defined previously, can be expressed as:

T F CT (m, ϕ) = g m x (2jπϕ) -g m q (2jπϕ)b(2jπϕ)ST F CT (ϕ) ( 36
)
where g m x , g m q , and ST F CT (ϕ) have been defined respectively in (29), (28), and [START_REF] Yahia | Bandpass sigma-delta modulators synthesis with high loop delay[END_REF], and b(s) is the transfer function of the D.A.C defined in [START_REF] Francken | A High-Level Simulation and Synthesis Environment for Σ∆ Modulators[END_REF].

On the other hand, this transfer function for the ODT modulator can be expressed as:

T F ODT (m, ϕ) = F m x (e 2jπϕ ) -F m q (e 2jπϕ )B(e 2jπϕ )ST F ODT (ϕ) (37)
where F m x , F m q , and ST F ODT (ϕ) have been defined respectively in (30) and ( 21), and B(Z) is the transfer function of the D.A.C defined in [START_REF] Suarez | Behavioral Modeling Methods for Switched-Capacitor Σ∆ Modulators[END_REF].

C. Application to the simulation of a low pass modulator

We consider a classical third-order lowpass modulator. Its topology is given in Fig. 19 with g 1 = 0.182, g 2 = 0.46, g3 = 1.495, t 1 = 0.664. These values can be easily found using, for example, the tools of Schreier [START_REF] Schreier | Understanding Delta-Sigma Data Converters[END_REF]. The chosen D.A.C delay is equal to half the sampling period (dT s = 0.5), and k is equal to 8. With these values, r = 4 and D = 0. By applying the section IV-A methodology, the equivalent ODT model is given by Fig. 20 with

N 1 x (Z) = Z -1 N 1 q (Z) = Z -1 N 2 F (Z) = Z -1 N 2 x (Z) = 1.13e -2 Z -1 N 2 q (Z) = 1.011Z -1 N 3 F (Z) = 1.5Z -1 -0.5Z -2 N 3 x (Z) = 2.18e -4 Z -1 + 3.27e -4 Z -2 N 3 q (Z) = 1.029Z -1 -0.028Z -2 (38)
If we had used a classical bilinear transformation of the filters, coefficients N 2

x and N 3 x would have been zero and coefficients N 2 q , N 3 q and N 3 F would have been Z -1 . However if we increase the value of k, then N 2

x and N 3 x tend towards zero and the other coefficients tend towards Z -1 . Fig. 21 compares the shape of the outputs of each filter obtained by a simulation with Simulink of the CT model and the ODT model. The input signal was a sampled and held signal at kF s . It can be seen that the ODT model provides the real value of each filter output for each sub-sampling time.

In addition to section III-B, it is possible to evaluate the transfer function between the modulator input and each filter output. The definition was given in section IV-B: a sine input signal at frequency ϕ is applied at the modulator input and we consider the component at this same frequency at each filter output. As in [START_REF] Sobot | Tunable Continuous-Time Bandpass Σ∆ Modulators With Fractional Delays[END_REF], it can be seen that taking k equal to 4 times the ratio between the input signal bandwidth and F s /2 provides accurate results for the simulation of intermediary signals.

The simulations times are almost the same as the ones obtained in table III. Simulating 100000 samples of a secondorder modulator requires 7 seconds (it was 6 seconds when intermediary signals were not expected).

D. Extension to non-idealities simulations

The ODT model can be extended to the simulation of non ideal behaviors. Some of them can be directly deduced from the ODT model while others need some mathematical developments. To start with, discrete-time noise sources can be easily added at the input or the output of all filters. This can Fig. 20. Equivalent ODT modulator help evaluate, for example, the effect of the high frequency noise components that are aliased in the signal band. The spectrum shape of the noise sources has to be derived from the CT noise models. The effect of the clock jitter [START_REF] Shoaei | Continuous-Time Delta-Sigma A/D Converters for High Speed Applications[END_REF] could also be included in this model. However, this would require a mathematical extension of the methodologies proposed, for example, in [START_REF] Reddy | Fundamental Limitations of Continuous-Time DeltaSigma Modulators Due to Clock Jitter[END_REF]. It is also possible to model the real transfer function of the amplifiers used in the modulator, as long as they can be expressed by their Laplace transform. This can be easily done by replacing the ideal filters G 1 to G n in Fig. 16 by their real transfer functions. Some non linearities can also be introduced, by including, for example, a non-linear transfer function at the output of the amplifiers. Nevertheless, such a model cannot simulate correctly, for example, non linear capacitances in an integrator. In any case, high level models can never replace a transistor-level simulation where all nonlinearities and parasitical elements are considered.

V. CONCLUSION

A methodology for time-domain simulations of continuoustime modulators was proposed. This methodology is based on a fixed step discretization of each output sample. Compared with variable step methods, the fixed step method ensures an exact equivalency in terms of noise transfer function, even for large simulation steps. The optimal step value is chosen from the bandwidth of the input signal. This allows for the simulation of the response to input signals with a bandwidth larger than half the modulator sampling frequency. This methodology can also help simulate the behavior of all the modulator outputs of the filters on top of the modulator output signal. Simulations are very fast as they use discretized equations. STF considerations have shown that the ODT method describes the behavior of a CT modulator better than a classical transformation method such as the bilinear transform.

Fig. 1 .Fig. 2 .

 12 Fig. 1. CT Σ∆ modulator

Fig. 4 .

 4 Fig. 4. CT Σ∆ modulator

Fig. 7 .

 7 Fig. 7. 2 nd order CT Σ∆ modulator

Fig. 8 .

 8 Fig. 8. Equivalent 2 nd order ODT Σ∆ modulator

Fig. 9 .

 9 Fig. 9. Simulink simulation of A.D.C input with zero input signal

Fig. 11 .

 11 Fig. 11. Simulink simulation of A.D.C input with CT sine input signal

Fig. 12 .

 12 Fig. 12.ST F distortion with and without correction for a first-order integrator

Fig. 14 .

 14 Fig. 14. ODT model of a second-order modulator including signal predistortion filter

Fig. 15 .

 15 Fig. 15. Comparison of STF obtained with ODT and CT model

Fig. 16 .

 16 Fig. 16. General CT topology

Fig. 17 .

 17 Fig. 17. ODT topology providing intermediary outputs

  has been transformed in order to be expressed as a classical multi-feedback topology. As the poles are conserved by the ODT transformation, F m x and F m q have the same poles. Furthermore, these poles are included in the poles of F m+1x and F m+1 q

Fig. 18 .

 18 Fig. 18. ODT topology for simulation

Fig. 19 .

 19 Fig. 19. Classical third-order lowpass modulator

Fig. 21 .

 21 Fig. 21. Intermediary signals shape

Fig. 22 .

 22 Fig. 22. Transfer function between input and first integrator output

  y=[zeros(1,or) x0(lngsolx+1:lngsolx+lngsoly)]; z=[zeros(1,or) x0(lngsolx+lngsoly+1:end)];

APPENDIX A MAPLE CODE TO SOLVE (1)

The following code was used to solve [START_REF] Schreier | Understanding Delta-Sigma Data Converters[END_REF] and to fill in table I. In this code, which deals with first-order terms ( 1 s-a ), the delay d is between 0 and 1. Larger delays can be considered by changing the sixth line of the code.

APPENDIX B MATLAB CODE TO SOLVE (34)

The following code can be used to solve (34). This code was tested only in the case when the size of s1 was the same as the size of s2. ------------------------------------------


Zhiguo Song

member, IEEE Caroline Lelandais-Perrault

member, IEEE Philippe Bénabes

Synthesis of Complex Subband Hybrid Filter Banks A/D Converters using Adaptive Filters

Hybrid Filter Banks (HFB) A/D converters (ADC) may be an attractive solution for future cognitive radio application. In this paper, using complex adaptive finite impulse response (FIR) filter whose coefficients are initialized with zero, their optimal values are achieved by iterative adjustments. This synthesis method is validated for a complex subband HFB by simulation in the time-domain.

I. INTRODUCTION

As one of the key technologies for the cognitive radio, A/D converters (ADCs) are expected to convert alternatively a wide band with a low resolution (in sensing mode) or a narrow bandwidth with a higher resolution (in communication mode). Hybrid Filter Banks (HFB) ADCs may be an appropriate solution for these requirements [START_REF] Lelandais-Perrault | Wide-band, band pass and versatile Hybrid Filter Bank A/D, conversion for software radio[END_REF] [START_REF] Lelandais-Perrault | Hybrid Filter Bank A/D conversion systems applied to future telecommunication scenarios[END_REF]. Being parallel architecture, HFB ADCs are composed of an analog frequency decomposition (analysis bank), a conversion-stage (ADC bank) and a digital reconstruction (synthesis bank). With a given ADC sampling rate, classical HFBs enlarge the conversion bandwidth. The resolution of conversion is related to the digital reconstruction complexity and the resolution of each ADC [START_REF] Petraglia | High speed A/D conversion incorporating a QMF bank[END_REF], [START_REF] Velazquez | Design of hybrid filter bank for analog/digital conversion[END_REF], [START_REF] Löwenborg | A class of twochannel hybrid analog/digital filter banks[END_REF]. Furthermore, when a narrow band conversion is needed, the digital part may be reconfigured to reconstruct this subband only, which defines subband HFB ADCs proposed in [START_REF] Asemani | Subband architecture for Hybrid Filter Bank A/D converters[END_REF].

The synthesis methods of the HFBs proposed in [START_REF] Velazquez | Design of hybrid filter bank for analog/digital conversion[END_REF] and [START_REF] Löwenborg | A class of twochannel hybrid analog/digital filter banks[END_REF] synthesize jointly the analog part and digital part. These methods require high-precision analog components to obtain a given transfer function. Therefore, their hardware implementation becomes more difficult and more expensive with the scaling of integrated circuit technologies [START_REF] Lelandais-Perrault | Wide-band, band pass and versatile Hybrid Filter Bank A/D, conversion for software radio[END_REF]. Relaxing the constraints upon the realization of the analog components, and assuming that the knowledges of the analysis frequency response are obtained, the synthesis filters can be achieved by calculating the inverse FFT of the ideal synthesis filter frequency response [START_REF] Velazquez | Hybrid filter banks for analog/digital conversion[END_REF]. With the same hypothesis, [START_REF] Lelandais-Perrault | Wide-band, band pass and versatile Hybrid Filter Bank A/D, conversion for software radio[END_REF] proposed a least mean square global approximation (LMSGA) synthesis method that minimizes a criterion describing the near prefect reconstruction of the classical HFB. This method is applied also for synthesizing the subband HFB [START_REF] Asemani | Subband architecture for Hybrid Filter Bank A/D converters[END_REF]. We proposed in [START_REF] Song | Adaptive Equalization for Calibration of Subband Hybrid Filter Banks A/D Converters[END_REF] an adaptive method to compensate analog errors digitally for a subband HFB. Initializing with the precalculated coefficients by the LMSGA method, although the analog errors degrade dramatically the performance of the subband HFB, the iterative adjustments of the synthesis filters can recover a good resolution finally.

Real part

Imaginary part ADC MT Hi(s) ui(n) MT ADC H0(s) u0(n) ADC MT HM-1(s) uM-1(n) Analysis bank x(t)
Synthesis bank of the subband m

xQm(n) xIm(n) + + FI m,0 (z) FI m,i (z) FI m,M -1 (z) FQ m,0 (z) FQ m,i (z) FQ m,M -1 (z)
H(s) F (z) dm(n) ŝm(n) x(t) xm(n) s(t) em(n) u(n)
Fig. 2. Principle of synthesis for complex subband HFB ADC: adaptive channel equalization (phase 1). Once the synthesis process is achieved, the subband HFB is ready to convert the extern wideband signal x(t) (phase 2). The think line designates the signal in complex form.

In this paper, we propose the adaptive equalization [START_REF] Widrow | The complex LMS algorithm[END_REF] [10] as a synthesis method for the subband HFB. The digital filters are all initialized to zero. Their coefficients are adjusted iteratively to reach an optimal reconstruction. Furthermore, as complex signals are widely used in digital processing for radiocommunication, we propose here an architecture (Figure 1) that reconstruct complex digital signals which are more convenient for the downstream treatments.

Figure 2 presents the principle of the proposed synthesis method in the M -chennal complex subband HFB context. Two linked signals are provided by the generator: the test signal s(t) and the desired signal d m (n). s(t) and d m (n) are respectively analog wideband signal and complex digital signal of the subband m ∈ [0, M -1]. The error e m (n) is used to adjust the synthesis filter coefficients in the subband m. The synthesis process will be achieved until the mean-square error (MSE) of e m (n) is small enough to meet the reconstruction accuracy requirements. Thus, the subband HFB is ready to convert the wideband signals x(t). The quantization errors are neglected here.

In section II, the synthesis method using adaptive channel equalization is theoretically justified for a complex subband HFB ADC. In section III, initializing the synthesis filters coefficients to zero, an eight-channel complex adaptive subband HFB ADC is simulated in the time-domain. A wideband signal different from the test signal s(t) is converted after the synthesis processing. Then, the performances of the complex subband HFB are discussed.

II. ADAPTIVE EQUALIZATION FOR SYNTHESIZING A

COMPLEX SUBBAND HFB ADC For each subband, the synthesis process is similar, so we can just focus on the reconstruction of one subband m.

A. Test signal s(t)

The test signal is a critical point to realize the proposed synthesis method. As a wideband signal, s(t) should cover spectrally all of the band to convert, i.e. its spectrum should be "rich" enough. Furthermore, it is necessary to avoid the additional interferences generated by the frequencies apart the passband. Filtered by bandpass filters with high selectivity, zero-mean white noise signal may satisfy the needs of the test signal for synthesizing the complex subband HFB.

B. Desired signal d m (n)

The desired signals provide the references for reconstructing digitally the subband signals of the test signal. Denoting S(jΩ ) the Fourier transform of s(t), Ω is the pulsation and its normalization is given by w = Ω T , where T is the Nyquist period of sampling. Using complex modulated passband filter G m (e jω ), S(jΩ) is split into M subbands which bandwidth is 2π M (Figure 3). Let define G m (e jω ) in two parts:

G m (e jω ) and G m (e jω ), where m ∈ [0, M 2 -1] and m ∈ [ M 2 , M -1]. So, G m (e jω ) can be presented as follows ⎧ ⎨ ⎩ G m (e jω ) = 1; R = (n + 2m ) π M ≤ ω ≤ (n + 2m + 2) π M G m (e jω ) = 0; ω / ∈ R (1)
⎧ ⎨ ⎩ Gm (e jω ) = 1; R = -(n + 2(m -M 2 )) π M ≥ ω ≥ -(n + 2(m -M 2 ) + 2) π M Gm (e jω ) = 0; ω / ∈ R (2)
where n ∈ Z. After being sampled at the rate 1 MT , the selected subband signals can be founded in the bass frequency. According to (1) [START_REF] Lelandais-Perrault | Hybrid Filter Bank A/D conversion systems applied to future telecommunication scenarios[END_REF] and the extraction of the subband signal illustrated in Figure 3, the desired signal of the subband m can be given by (3) for n is odd

D m (e jω ) = ⎧ ⎨ ⎩ 1 M S(e j( ω M +2m π M +(n+1) π M ) ) 1 M S(e j( ω M +(2m +1) π M +(n+1) π M ) ) (3)
¡ ¡ ¡ ¡ ¡ ¢ ¢ ¡¢ ¢ ¢ ¡¢ ¢ ¢ ¡¢ ¢ ¢ ¡¢ ¢ ¢ ¡¢ ¢ ¢ ¡¢ ¢ ¢ ¢ £ £ ¤£ ¤£ ¤£ £ £ ¤£ ¤£ ¤£ £ £ ¤£ ¤£ ¤£ £ £ ¤£ ¤£ ¤£ £ £ ¤£ ¤£ ¤£ £ £ ¤£ ¤£ ¤£ £ £ £ £ £ £ £ £ £ £ ¥ ¥ ¤¥ ¤¥ ¤¥ ¥ ¥ ¤¥ ¤¥ ¤¥ ¥ ¥ ¤¥ ¤¥ ¤¥ ¥ ¥ ¤¥ ¤¥ ¤¥ ¥ ¥ ¤¥ ¤¥ ¤¥ ¥ ¥ ¤¥ ¤¥ ¤¥ ¥ ¥ ¤¥ ¤¥ ¤¥ ¥ ¥ ¥ ¥ ¥ ¦ ¦ ¤¦ ¤¦ ¤¦ ¤¦ ¦ ¦ ¤¦ ¤¦ ¤¦ ¤¦ ¦ ¦ ¤¦ ¤¦ ¤¦ ¤¦ ¦ ¦ ¤¦ ¤¦ ¤¦ ¤¦ ¦ ¦ ¤¦ ¤¦ ¤¦ ¤¦ ¦ ¦ ¤¦ ¤¦ ¤¦ ¤¦ ¦ ¦ ¦ ¦ ¦ ¦ § § ¤ § ¤ § ¤ § ¤ § § § ¤ § ¤ § ¤ § ¤ § § § ¤ § ¤ § ¤ § ¤ § § § ¤ § ¤ § ¤ § ¤ § § § ¤ § ¤ § ¤ § ¤ § § § ¤ § ¤ § ¤ § ¤ § § § § § § § ¡¡¡ ¡¡ ¡¡ ¡¡ ¡¨© ¡© ¡© © ¡© ¡© © ¡© ¡© © ¡© ¡© © ¡© ¡© © © © ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ 0 GM 2 -1 (e jω ) GM 2 (e jω ) GM-1(e jω ) Gm (e jω )
G0(e jω ) rad/s When n is even, D m (e jω ) is given by

S(jΩ)Gm (e jω ) (n + M ) π M Dm (e jω ) Dm (e jω ) n is even 1 n is odd rad/s rad/s 0 -(n + M ) π M (n + 2m + 2) π M (n + 2m ) π M n π M -n π M 0 -π π -π π
D m (e jω ) = ⎧ ⎨ ⎩ 1 M S(e j( ω M +(2m +1) π M +(n+2) π M ) ) 1 M S(e j( ω M +2m π M +n π M ) ) (4)
D m (e jω ) can be achieved by using the same method for extracting D m (e jω ). Hence, the digital complex desired signals are ready for synthesizing the subband HFB adaptively.

C. Synthesis of the complex subband HFB

The wideband signal of test s(t) is decomposed by the analog filters into M channels, then undersampled at the rate of 1 MT to be u(n) denoting the inputs of the synthesis filters, u(n) is given by

u(n) = [u 0 (n), u 1 (n), • • • , u M -1 (n)] T , (5)
where v T means the transpose of the vector v, and u i (n) is denoted by

u i (n) = [u i (n), u i (n -1), • • • , u i (n -L + 1)] ( 6
)
where L is the length of the synthesis filter. Let the vectors fm,i = [ fm,i (0), fm,i (1), • • • , fm,i (L -1)] denote the complex synthesis filters coefficients in subband m channel i, which will be adjusted adaptively to minimize the reconstruction error. Let fm denote all the synthesis filter coefficients in the subband m

fm = [ fm,0 , fm,1 , • • • , fm,M-1 ] (7)
Let e m (n) denote the difference between the desired signal and the estimated subband signal ŝm (n). ŝm and e m (n) are defined respectively as follows

ŝm (n) = f * m u(n) (8) e m (n) = d m (n) -ŝm (n) (9)
Let define J m the chosen cost function MSE, its dependence on the complex synthesis coefficients is shown as follows

J m = E(|e m (n)| 2 ) = E(|d m (n)| 2 ) -f * m E(u(n)d * m (n)) - fm E(u(n)d m (n)) + f * m E u(n) T u(n) f T m (10)
The expectation E(u(n)d m (n)) is equal to the crosscorrelation between u(n) and d m (n), noted p m . The expectation E(u(n) T u(n)) is equal to the auto-correlation of the vector u(n), noted R uu . So, [START_REF] Haykin | Adaptive Filter Theory[END_REF] can be written by

J m = σ 2 dm -f * m p * m -f * m p m + f * m R uu f T m (11)
Since R uu is a matrix nonsingular, its inverse matrix R -1 uu exists. Minimizing J m in terms of fm , the optimal synthesis coefficients can be achieved by solving the Wiener-Hopf equations [START_REF] Haykin | Adaptive Filter Theory[END_REF] and given by

f opt m = R -1 uu p m ( 12
)
The calculation of f opt m with ( 12) is obviously too costly for inversing the large matrices. Thus, the complex LMS algorithm is chosen here due to its simplicity [START_REF] Widrow | The complex LMS algorithm[END_REF] [10], which updates iteratively the synthesis filter coefficients fm for approximating f opt m the (Figure 4) with the formulas as follows

fIm (n + 1) = fIm (n) + μe Im (n)u(n) (13) fQm (n + 1) = fQm (n) + μe Qm (n)u(n) ( 14
)
where μ is named step-size. The condition of stability of long complex LMS FIR filters is shown below [START_REF] Haykin | Adaptive Filter Theory[END_REF]

0 < μ < 1 MLE[u 2 i (n)] , (15)
where E u 2 i (n) denotes the expectation of the power spectral density of the synthesis bank inputs.

+ + + + FI m,0 (z) FI m,i (z) FI m,M -1 (z) FQ m,0 (z) FQ m,i (z) FQ m,M -1 (z) eQ m (n) - + eI m (n) + - dI m (n) dQ m (n) ŝQm(n) u0(n) ui(n) uM-1(n) ŝIm(n)
Fig. [START_REF] Velazquez | Design of hybrid filter bank for analog/digital conversion[END_REF]. Synthesis of the complex subband HFB ADC by using complex LMS algorithm.

D. Performance of the complex subband HFB ADC

Figure 1 shows the architecture of the complex subband HFB, which reconstruction function [START_REF] Song | Adaptive Equalization for Calibration of Subband Hybrid Filter Banks A/D Converters[END_REF] can be presented in frequency-domain by

T m,p (e jω ) = M -1 i=0 H i (j(Ω - 2πp MT ))F m,i (e jω ) (16)
where the index p ∈ [-(M -1), M -1]. Theoretically, the conditions for perfect reconstruction of the subband are shown as follows

⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ T m,0 (e jω ) = e -jωρ , p = 0, ρ ∈ R + * M -1 -(M -1)
T m,p (e jω ) = 0, p = 0

where T m,0 (e jω ) is defined as the distortion function, which is expected to a pure delay of ρ; the sum of T m,p =0 (e jω ) is defined as the aliasing function, which is expected to be 0. Hence, the performance characteristic of the subband HFB for reconstruction of the subband m can be defined here by the Distortion-to-Aliasing Ratio (DAR)

DAR dB (ω) = 20log Tm,m (e jω ) ce -jωρ i =m T m,i (e jω ) 2 (18)

III. SIMULATION RESULTS

A. Conditions of simulation

Considering the wideband signal to convert is in the frequency band [0.85 1.25]GHz. The analog signal is split into an analog bandpass filter bank which consists in eight pure Gm-LC resonators [START_REF] Lelandais-Perrault | Wide-band, band pass and versatile Hybrid Filter Bank A/D, conversion for software radio[END_REF], whose central frequencies and quality factors are respectively 2πΩ i and Q i . The transfer function of the analysis filter of each channel is given by

H i (s) = Ωi Qi s s 2 + Ωi Qi s + Ω 2 i (19)
According [START_REF] Lelandais-Perrault | Hybrid Filter Bank A/D conversion systems applied to future telecommunication scenarios[END_REF], the dynamic range of the radiocommunication application like GSM is defined to be 82 dB i.e a resolution of 14 bit. This can be considered as the goal of reconstruction accuracy. Each synthesis subband filters of the subband m are complex FIR filters whose coefficients are initialized by zero. The process of synthesis of the complex subband HFB is simulated in the time-domain with Matlab using C MEX-functions.

A pseudo random sequence is filtered by the complex subband anti-aliasing filters, their sum forms the test signal s(t). The signals of each subband m is sampled at the frequency 1 MT . Thus, complex desired signal d m (n) is obtained. Here, the purpose of the simulations is to synthesize the subband 0. In order to reduce the aliasing terms, it's efficiency to allocate a small part of frequency borders (low and high frequencies) at each subband spectrum as guard band (GB) [START_REF] Asemani | Subband architecture for Hybrid Filter Bank A/D converters[END_REF], which equals 7% of the subband bandwidth in this paper. Figure 5 shows the spectrums of s(t) and d 0 (n).

B. Simulation results of synthesis

The synthesis process is stopped as soon as the MSE of e 0 (n) is small enough, i.e. the performance of the subband HFB satisfies the application requirements. Here, the stop condition is defined by dynamic range of 82 dB, which corresponds to a resolution of 14 bit. In the end of the synthesis process, the complex synthesis coefficients may be considered as the optimal coefficients noted by f opt 0 . In Figure 6, the above part shows the convergence of the coefficients VS the iteration number. The performance of the subband HFB is calculated in the time-domain by desired signal to error ratio, which noted SNR conventionally. The convergence of the SNR is shown in the below part, which progresses for 14.6 bit. With f opt 0 , the performances are calculated in the frequencydomain. The frequency responses of the synthesis filters in the subband 0 are presented in the above part of Figure 7. According to (18), the performance DAR is calculated to be 90 dB, which is presented in the below part.

C. Simulation results of validation

After being synthesized with good performance, the complex subband HFB is ready to convert the others signals differed from s(t) (Figure 1). In order to validate this aspect, fixing f opt 0 , the comb sinus signals are applied as x(t). The simulation result shows that the resolution of conversion of the subband 0 reaches 14.5 bit.

IV. CONCLUSION

In this paper, negating the quantization errors, a synthesis method using adaptive filters for a complex subband HFB ADC is proposed. The optimal synthesis coefficients can be approached iteratively by using the complex LMS algorithm. This method is easy for implementing in a digital circuit. The resolution of conversion for one subband can reach 14 bit, this meets severals needs of cognitive radio.

Fig. 1 .

 1 Fig. 1. Continuous-time M -channel complex subband HFB ADC, T denotes the Nyquist sampling period of input wideband signal. The subscripts denote the subband number and the channel number respectively, m, i ∈ [0, M -1].

Fig. 3 .

 3 Fig. 3. Above: the passband filters for selecting the subband signals, where wheren ∈ Z, m ∈ [0, M 2 -1], m ∈ [ M 2 , M -1]. Below: extraction of each subband desired signal D m (e jω ) from S(jΩ).

Fig. 5 .

 5 Fig. 5. Specturum of the test signal s(t) and the desired signal d 0 (n).

Fig. 6 .

 6 Fig. 6. Above: the convergence of the coefficient indexed by 64 VS the iteration number. Below: the evolution of the SNR presented in bit for the reconstruction of the subband 0.

Fig. 7 .

 7 Fig.7. Above: the frequency responses of the complex filters in the subband 0. Below: frequency response of complex subband HFB ADC.


J Juillard

G Arndt

E Colinet

Enhancing the frequency stability of a NEMS oscillator with electrostatic and mechanical nonlinearities

We show how to achieve large-amplitude oscillation and good frequency stability in a nonlinear NEMS oscillator by using a softening nonlinearity to locally balance the hardening behaviour of the Duffing nonlinearity for a given motion amplitude. The case when the softening behaviour is obtained via electrostatic biasing is theoretically investigated, under parallel-plate assumptions, for a fully-capacitive pulse-actuated resonant structure. This is validated by simulated data, in the case of a resonant clamped-clamped beam.

Motivations and results

), where 0  is the natural pulsation of the resonant structure, A is the oscillation amplitude, normalized with respect to the gap and  and  are positive design-dependent coefficients respectively representing the strength of the Duffing hardening and that of the electrostatic softening. From (1), we establish that:

5 / 2 2 1 0                  opt A A dA d . (2)
Thus, provided the structure is designed so that   2  , there exists an optimal oscillation amplitude around which small perturbations on A have little impact on the value of  . In order to ensure a large SNR, opt A should be made as large as possible, but safely maintained below pi A , for which electrostatic pull-in [START_REF] Elata | On the static and dynamic response of electrostatic actuators[END_REF] occurs (Fig. 34). Similar analyses can be made in the case of flexible structures; although they are not tractable analytically, they lead to results that are qualitatively similar to (1-2) (Fig. 5). These theoretical results are compared with simulations, in the case of the resonant structure developed in the ANR-funded M&NEMS project (Fig. 6). This design approach remains qualitatively valid for several other actuation schemes (PLL, etc.).

Word count: 496

Corresponding author: Jérôme Juillard, Dept. SSE, Supélec, 3 rue Joliot-Curie, 91192 Gif-sur-Yvette Cedex, France, Tel: +33-1-69-85-14-08, Fax: +33-1-69-85-14-29, E-mail: jerome.juillard@supelec.fr. Parallel-plate structure

V b V gnd V b /2 - + Pulse generation + + + b(t) + V c Q(t) 3 a     2 0 2 1         Q j j H tot + +   a a V c  , Duffing NL Microbeam Feedback NL   c elec V a f ,

Control electronics

  2 / 3 2 2 0 1 1 4 3 1 A A         (1)   0 1 1 1 2 2 2     pi pi A A   (3) Clamped-clamped beam         2 2 2 2 2 2 0 1 1 4 4 3 1 A A A K A A E A            (4)   0 1 1 1 2 / 3 2 2     pi pi A A   (5)
Fig. 3 -For a pulse-actuated parallel-plate structure, it is possible to find analytical expressions for the optimal oscillation amplitude (2) and for the pull-in amplitude, which is governed by ( 3). ( 3) reduces to a third -degree polynomial in Api 2 . The pull-in amplitude corresponds to a maximum of the potential energy of the system. For a clamped-clamped beam, the expression of the closed-loop oscillation pulsation contains elliptic integrals (4): Aopt can only be determined numerically. Each point of the curve corresponds to a nonlinear transient simulation of the closed-loop system, with a pulse amplitude varying between 10mV and 290mV. The relative error is defined as the ratio of the estimated standard deviation of the (amplitude-dependent) pulsation to its estimated mean. It is minimal close to the expected oscillation amplitude (Fig. 5). Note that, for small pulse amplitudes, the detected signal is drowned in the measurement noise, resulting in poor frequency stability.

Fig. 1 -

 1 Fig. 1 -Principle of closed-loop pulse-actuation. The moving structure (in green) is the midpoint of a biased capacitive half-bridge. The amplifier's output voltage is proportional to the charge of the moving electrode. Short voltage pulses are fed back to the central electrode when the amplifier's output voltage crosses zero [2-3], so that oscillation builds up. Measurement noise b(t) is responsible for frequency instability.

Fig. 2 -

 2 Fig.2-Equivalent block diagram of an electrostatically-actuated, stress-stiffened NEMS (shaded area) with its nonlinear feedback electronics and actuation scheme. The system may be looked at as a linear system with one complex feedback nonlinearity. Resonant NEMS usually have large quality factors, making it possible to study them with describing function analysis.

Fig. 4 -

 4 Fig. 4 -Map of Aopt/Api vs.  and for a parallel-plate resonant structure. There exists no optimal amplitude below the lower boundary and Aopt>Api above the upper boundary.

Fig. 5 -

 5 Fig. 5 -Plot of Aopt and Api vs.  for the microbeam of the M&NEMS project (0=2.2 ×10 7rad.s -1 , =0.72, Vpi=30V). The pull-in voltage Vpi is defined as the bias voltage for which the center position becomes unstable. =0.16 corresponds to Vb=12V.

Fig. 6 -

 6 Fig.6-Simulated frequency stability of the M&NEMS structure under a 12V bias. The measurement noise is assumed white and has a uniform distribution. Each point of the curve corresponds to a nonlinear transient simulation of the closed-loop system, with a pulse amplitude varying between 10mV and 290mV. The relative error is defined as the ratio of the estimated standard deviation of the (amplitude-dependent) pulsation to its estimated mean. It is minimal close to the expected oscillation amplitude (Fig.5). Note that, for small pulse amplitudes, the detected signal is drowned in the measurement noise, resulting in poor frequency stability.


J Juillard
email: jerome.juillard@supelec.fr

K Jafaridinani
email: kian.jafaridinani@supelec.fr

E Colinet
email: eric.colinet@cea.fr

In this paper, we present a new weighted least-squares (WLS) approach for parameter estimation based on binary data. Two WLS criteria are studied. We show that these two criteria do not have the same asymptotical behavior although they are closely related. Particularly, in the presence of noise, one of the criteria used for determining the system parameters provides an appropriate estimation, whereas the other one leads to an underestimation of the system parameters. These asymptotical results are illustrated by simulations in Gaussian and non-Gaussian contexts.

INTRODUCTION

In this paper, we give some theoretical results concerning the asymptotic bias of least-square approaches to parameter estimation methods based on binary measurements. This work is based on [START_REF] Colinet | A weighted least-squares approach to parameter estimation problems based on binary measurements[END_REF]. It is originally motivated by the need to add integrated low-cost self-test features to micro fabricated devices, such as MEMS and NEMS. Even though there exists a wide range of applications where identification methods based on binary observations are necessary or desirable [START_REF] Wang | System identification using binary sensors[END_REF], the focus is brought here on the test of microelectronic devices.

Context and state of the art

It is well-known that, as characteristic dimensions become smaller, the dispersions afflicting electronic devices tend to become larger. Typical sources of dispersion are variations in the fabrication process or changes in the operating conditions, such as temperature and ambient pressure. Imperfect knowledge of semiconductor physics or of nano-mechanics is also a large source of uncertainty. As a consequence, it is usually impossible to guarantee a priori that a given device will function properly. The tests that are run after fabrication, under different operating conditions, ensure that only suitable devices are commercialized. However, these tests are very costly: a typical figure is that one third of the cost of a micro fabricated device is due to testing [START_REF] Charlot | Electrically induced stimuli for MEMS self-test[END_REF]. An alternate solution is to integrate low-cost self-test (and selftuning) features, such as parameter estimation routines, in each device, so that it can adapt to changing conditions. However, most parameter estimation methods [START_REF] Ljung | System identification -theory for the user[END_REF][START_REF] Walter | Identification of parametric models from experimental data[END_REF] do not easily scale down to the micro-or nano-world. This is because these methods rely on high-resolution digital measurements of the system output. Their integration requires the implementation of highresolution analogue-to-digital converters (ADCs) and, thus, results in longer design times, larger silicon areas and increased costs. Our objective is then to develop a parameter estimation method that relies on very low-resolution (ideally binary) measurements, in order to keep the added cost as small as possible.

In the field of micro-electronics, this issue has been addressed by [START_REF] Negreiros | Ultimate low cost analog BIST[END_REF]: the proposed method consists in using a white Gaussian input to excite the unknown (linear) system and to estimate the power spectral density (PSD) of the binary output, which is assumed to be stationary. From this estimated PSD, it is possible to analytically derive the modulus of the transfer function of the unknown system. However, one cannot use this approach to obtain any information concerning the phase of the transfer function. Moreover, the generation of a white Gaussian input requires a high-resolution digital-to-analogue converter (DAC), whose implementation can be costly. These issues are partially solved in [START_REF] Juillard | Initialization of the BIMBO self-test method using binary inputs and outputs[END_REF], where a white Bernoulli input is used as an exciting signal (thus requiring only a 1-bit DAC). Provided the mixing properties of the linear system are good (i.e. its impulse response does not vanish too quickly), it is possible to establish an analytical relationship between the cross-covariance function of the binary inputs and outputs and the impulse response of the unknown system and, thus, to identify the system.

Outside the context of micro-electronics, the most significant contributions come from [START_REF] Wigren | Adaptive filtering using quantized output measurements[END_REF] and Wang and his co-workers [START_REF] Wang | Joint identification of plant rational models and noise distribution functions using binaryvalued observations[END_REF][START_REF] Wang | System identification using binary sensors[END_REF][START_REF] Zhao | Identification of Wiener systems with binary-valued output observations[END_REF]. In [START_REF] Wang | System identification using binary sensors[END_REF], Wang. introduced a new method for estimating parameters from binary (or quantized) data. The unknown system is excited by a periodic signal and, as in [START_REF] Rafaljowicz | Linear systems identification from random threshold binary data[END_REF], the threshold of the quantizer is randomly specified by a partially known dithering signal. However, the hypotheses on the quantizer's input and on the random threshold are much less strict than in [START_REF] Rafaljowicz | Linear systems identification from random threshold binary data[END_REF]: the main constraint is that the cumulative distribution function (cdf) of the threshold must be invertible. Although it is usually possible to generate continuous dithering signals, this constraint prevents their approach from being used with

Consistency of weighted least-square estimators for parameter estimation problems based on binary measurements

discrete-valued dithering signals, which are very easy to generate in the context of microelectronic applications. Under these hypotheses, the parameter estimation problem reduces to solving a linear system. This approach is generalized in [START_REF] Wang | Joint identification of plant rational models and noise distribution functions using binaryvalued observations[END_REF], where it is shown that the cdf does not have to be known a priori: it can be estimated along with the parameters of the system. This work has been extended from finite impulse response (FIR) systems to infinite impulse response (IIR) systems and to nonlinear Wiener systems [START_REF] Wang | Joint identification of plant rational models and noise distribution functions using binaryvalued observations[END_REF][START_REF] Zhao | Identification of Wiener systems with binary-valued output observations[END_REF]. Finally, Wigren has developed a least mean squares (LMS) approach to the problem of online parameter estimation from binary observations [START_REF] Wigren | Adaptive filtering using quantized output measurements[END_REF]. This method is based on an approximation of the quantizer, which makes it possible to define an approximate gradient of the least-squares criterion (this is developed in section II). Under some hypotheses [START_REF] Ljung | Analysis of recursive stochastic algorithms[END_REF], it is possible to guarantee the asymptotical convergence of this method to the nominal parameters (or at least to the boundary of the model set).

Each of these methods has its merits and weaknesses. It is our opinion that Wigren's approach is the one that fits better in the context of test of electronics system: as opposed to Wang's or Rafajlowicz's methods, it does not require a varying threshold (dithering signal), which can only be implemented with a high-resolution DAC or a surfaceconsuming analog component. Furthermore, it can be used online, without specific input design, as well as offline.

However, it does rely on an approximation of the quantizerand, thus, it is not an exact method, stricto sensu, as opposed to [START_REF] Colinet | A weighted least-squares approach to parameter estimation problems based on binary measurements[END_REF]. Moreover, convergence is not guaranteed in the presence of measurement noise at the quantizer's input.

The method introduced in [START_REF] Colinet | A weighted least-squares approach to parameter estimation problems based on binary measurements[END_REF] bears some likeness to Wigren's, in the sense that it is based on the minimization of a least-squares criterion. Its main distinguishing feature is that no approximation of the quantizer is made (hence, it is not based on a "pseudogradient"): rather, it relies on a weighted least-squares (WLS) criterion where the parameter-dependent weights are chosen in order to smooth out the discontinuities. The asymptotical convergence of this approach to the nominal system parameters can be guaranteed, even in the presence of measurement noise, provided the signal at the quantizer's input is Gaussian and centred. Note that these hypotheses used in [START_REF] Colinet | A weighted least-squares approach to parameter estimation problems based on binary measurements[END_REF] only provide sufficient conditions for the asymptotical convergence and should not be considered as restrictive.

Outline of the article

In this paper, we present some new asymptotical results concerning WLS criteria in the presence of measurement noise at the quantizer's input. The structure of the article is the following. In section II, the notations are introduced and the WLS criteria presented in [START_REF] Colinet | A weighted least-squares approach to parameter estimation problems based on binary measurements[END_REF]) is re-introduced, along with their asymptotical properties. In section III, we study an alternative WLS approach: we prove that this approach is as efficient as the one proposed in [START_REF] Colinet | A weighted least-squares approach to parameter estimation problems based on binary measurements[END_REF] in the absence of measurement noise. On the other hand, we show that this other approach leads to a systematic error in the presence of measurement noise. Some examples are given in section IV. Section V contains some concluding remarks.

PRELIMINARIES

Framework and notations

Let us consider a discrete-time invariant linear system H .

We assume H has a finite impulse response of length L , i.e. the impulse response can be represented by a column vector

( ) L k k 1 = = θ θ
. Let l u be the known scalar value of the system input at time l and l y be the (scalar) value of the system output, so that:

θ φ T l l y = , (1)
where

( ) l L l k k l u 1 + - = = φ
is the (column) vector of observations at time l . Let l d be an additive dithering signal and l b an additive noise at the quantizer's input. Let x m and 2 x σ denote the first-and second-order moments of any signal x . The system output is measured via a 1-bit ADC so that only the sign ( )

l l z S s =
of the system output is known, where

( ) ( )    - = ≥ = otherwise , 1 0 if , 1 x S x x S , (2)
and

l l l l b d y z + + = .
We are interested in finding an estimate θ ˆ of θ , based on N observations of l u , l s and, if need be, l d . It is assumed that l u , l d and l b are stationary and uncorrelated and that l z is Gaussian and centred. The estimated quantities are denoted by a hat (e.g.

l l l d y z + =
). In order to simplify some expressions, l d and l b are assumed to be white and centred and θ ˆ has length L .

These notations are summed up in Fig. 1.

WLS approach for parameter estimation based on binary data

In [START_REF] Colinet | A weighted least-squares approach to parameter estimation problems based on binary measurements[END_REF], we show that the problem of parameter estimation from binary measurement can be treated by minimizing WLS criteria of the form:

( )( ) ( ) ∑ ∑ = = - = N l l N l l l l N f z f s s z f J 1 1 2 ˆ4 1 , (3)
+ ℜ ∈ ∀ > x x f (in practice, ( ) 2 x x f = is often used). It is simple to show that ( ) θ N f J is
continuously differentiable everywhere, except where:

( ) 0 1 = ∑ = N l l z f . (4)
Since ( ) . f is positive, (4) reduces to:

0 θ φ φ =           +           N T N T d d M M 1 1 ˆ, (5)
which is a system of N equations with L ( N L < ) unknowns. Provided l u has a density and 0

/ ≠ ∃ l d l
, it can be proved that (5) has no solution with probability 1. When

0 = l d , N l ≤ ≤ 1 , (4)
= ⇒ = N f N J J , (6)
where N J 0 is the (discontinuous) non-weighted least-squares criterion. In [START_REF] Colinet | A weighted least-squares approach to parameter estimation problems based on binary measurements[END_REF], the case when

( ) 2 x x f
= is studied in depth and the properties of the corresponding WLS criterion, N J 1 , are established in a probabilistic framework, when N goes to infinity. In particular, it is shown that N J 1 converges surely to:

( )

      - - = ∞ 2 1 1 acos 1 r r r J π , ( 7
)
where r is the correlation coefficient of z and z ˆ. Under the assumptions made in section II.1, the expression of r is:

( ) θ R θ θ R θ θ R θ θ 2 2 2 2 u T d u T d b u T d r + + + + = σ σ σ σ , ( 8
)
where u R is the L L × covariance matrix of u . By straightforward differentiation, one can then prove the following properties.

Property 1 -Provided u R has full rank, the minimum of

∞ 1
J is reached for:

- θ θ = ˆ if 0 2 ≠ d σ . - 0 , ˆ> ∀ = λ λθ θ if 0 2 = d σ . Property 2 - ( ) θ 1∞ J
is convex in a neighbourhood of its minimum.

Note that these properties hold regardless of the presence or absence of noise at the quantizer's input. Note also that it is not possible to identify the system's gain if 0 2 = d σ . The dithering signal acts as reference amplitude, as in [START_REF] Wang | System identification using binary sensors[END_REF].

AN ALTERNATIVE WLS APPROACH

An alternative approach to the problem of parameter estimation based on binary measurements is to define a criterion of the form:

( )( ) ∑ = - = N l l l l N f s s z f N G 1 2 41 , (9)
where ( ) . f has the same properties as in section II.2, in order to ensure continuous differentiability of N f G . In the case when ( ) 2 x x f = , the corresponding criterion, N G 1 , lends itself to analysis in much the same way as N J 1 . In particular, if l z and l z ˆ are stationary:

( ) ( ) 2 2 1 1 Ê 4 1 lim l l l N N s s z G G - = = ∞ → ∆ ∞ , (10)
where ( ) . E denotes the expected value. Under the assumption that l z and l z ˆ are Gaussian and centred, we have:

( ) ( ) ( ) ( ) ( ) ( )         + = - = ∫ ∫ ∫ ∫ ∫ ∫ ∞ ∞ - ∞ - ∞ ∞ ∞ - ∞ ∞ - ∞ 0 0 2 0 0 2 2 2 1 , , 4 1 , 4 1 dz z d z z p z dz z d z z p z dz z d z z p z S z S z G , ( 11
)
where ( )

0 0.5 1 1.5 2 2.5 3 3.5 0 0.2 0.4 0.6 0.8 1 σ b λ Minimization of J 1 N , σ d = 1 Minimization of G 1 N , σ d = 1 Minimization of J 1 N , σ d = 2 Minimization of G 1 N , σ d = 2 θ θ θ θ = [1 1] m b = m u = m d = 0 σ u = σ d = 1 0 0.5 1 1.5 2 2.5 3 3.5 -0.2 0 0.2 0.4 0.6 0.8 1 σ b λ Influence of dither distribution function Minimization of J 1 N with d n Bernoulli Minimization of G 1 N with d n Bernoulli Minimization of J 1 N with d n uniform Minimization of G 1 N with d n uniform θ θ θ θ = [1 1] m b = m u = m d = 0 σ u = σ d = 1
( ) 2 ˆ2 2 2 2 2 1 2 2 1 2 1 exp , ˆr z z r z z r z z p z z z z z z -                 - + - - = σ πσ σ σ σ σ
, and thus:

( ) ∞ ∞ =       - - = 1 2 2 2 1 1 acos J r r r G z z σ π σ . ( 12
)
The optimal θ ˆ is found by differentiating (12). This yields:

0 θ θ R 0 θ = ∂ ∂ + ⇔ = ∂ ∂ ∞ ∞ ∞ 2 ˆ1 2 1 1 r dr dJ J G z u σ , ( 13
)
where

2 1 1 2 r dr dJ - - = ∞ π , ( 14
) and θ R θ R θ 2 ˆu u z z z r r - = ∂ ∂ σ σ σ . ( 15
)
Thus, ( 13) is equivalent to:

( )

0 θ R θ R = - - u z z u r r σ σ 2 1 âcos . ( 16
)
If u R has full rank, this reduces to:

( )

θ θ z z r r σ σ 2 1 âcos - = . ( 17
)
Thus the optimal θ ˆis (at least) collinear to θ . Two cases must now be distinguished, depending on whether a dithering signal is present at the quantizer's input or not.

3.1 Case 1:

0 2 ≠ d σ Suppose now that 1 = r : in this case, ( )
0 1 acos 2 = - = r r
and ( 16) is verified. In other words, any vector θ ˆ such that 1 = r minimizes ∞ 1 G . Now, it is clear that r can be equal to 1 if and only if:

     = = θ θ ˆ0 2 b σ . ( 18
)
Thus the following property holds.

Property 3 Provided R has full rank, the minimum of

∞ 1 G is reached for θ θ = ˆ if 0 2 ≠ d σ and 0 2 = b σ . Consider now that 0 2 ≠ b σ
and, as a consequence, that 1 ≠ r . ( 17) can be rewritten as:

( )

θ θ z z r r σ σ 2 acos 1 ˆ-= . ( 19
) Let ( ) ( ) z z z z r g r r σ σ σ σ λ 2 acos 1 = - = . ( 20
)
Taking the square of (20) and multiplying by 2 z σ yields: We know that ( )

( ) ( ) ( ) ( ) ( ) ( ) ( ) 2 2 2 2 2 2 2 2 2 2 2 2 1 d u T d b u T d u T d b r g r g r g σ σ σ λ λ σ σ σ λ = - + + ⇔ T d b d r g r g 2 2 2 1- + + = σ σ σ λ (22) 0 0.5 1 1.5 2 2.5 3 3.5 0 0.2 0.4 0.6 0.8 1 σ b λ Influence of m u Minimization of J 1 N Minimization of G 1 N θ θ θ θ = [1 1] m u = 0.6 m b = m d = 0 σ u = σ d = 1 0 0.5 1 1.5 2 2.5 3 3.5 0 0.2 0.4 0.6 0.8 1 σ b λ Influence of m b Minimization of J 1 N Minimization of G 1 N θ θ θ θ = [1 1] m b = 1.5 m u = m d = 0 σ u = 1,σ d = 2
1 , 1 0 < ∀ < ≤ r

SIMULATION RESULTS

In this section, we compare the results obtained with WLS criteria N G 1 and N J 1 defined in the previous sub-sections. We , depending on which criterion is tested, and P is the number of experiments. P is chosen so that the standard deviation of λ ~ is small. Note that the angle between p θ ˆ and θ must be small for (24) to give a meaningful estimation of λ .

Since L N >> , this should indeed be the case: as a measure of precaution, one can verify that the correlation coefficient between p θ ˆ and θ is very close to 1, for all experiments.

The results are shown in Fig. 2. The value of λ estimated with N J 1 is always very close to 1, regardless of the value of d σ or b σ : this is as expected from the asymptotic analysis (Property 1) and shows that it is possible to properly estimate the system parameters regardless of the level of noise at the quantizer input, when N J 1 is selected as the criterion.

Conversely, when N G 1 is used, λ is equal to 1 only when 0 = Several other scenarios can be tested in which the constraints on the distributions of l u , l d and l b are made less strict. In Fig. 3,l u and l b are Gaussian and centred, whereas l d is uniform or Bernoulli. This test shows that, when b σ increases, λ becomes larger than 1 when N J 1 is used. In Fig. 4,l d and l b are Gaussian and centred, whereas l u is

Gaussian with nonzero mean: simulations show that λ becomes smaller than 1 when N J 1 is used. Finally, in Fig. 5, l u and l d are Gaussian and centred, whereas l b is Gaussian with nonzero mean. However, this does not seem to have any influence on λ , which stays very close to 1 when N J 1 is used. In these last three scenarios, the qualitative behaviour of λ when N G 1 is used does not differ much from the case when all signals are Gaussian and centred. One can also notice that the bias is always larger when N G 1 is used.

These simulations tend to show that the hypotheses made for obtaining Properties 1-6 may somehow be relaxed when 0 =

CONCLUSION

In this paper, we present a new WLS approach for parameter estimation based on binary data: this approach is an alternative to the one presented in [START_REF] Colinet | A weighted least-squares approach to parameter estimation problems based on binary measurements[END_REF]. However, we have shown that these two criteria, although they are closely related, do not have the same asymptotical behavior. In particular, we have established that, in the presence of noise, using N G 1 for determining the system parameters leads to an underestimation of θ , even though the direction of the vector is properly determined. These asymptotical results are confirmed by simulations performed in a Gaussian context. Other scenarios, where not all signals are Gaussian or centred, have also been tested: these show that the hypotheses made for obtaining the theoretical results can somehow be relaxed.

  Fig. 1. Framework and notations

Fig. 2 .

 2 Fig. 2. λ ˆ vs. b σ when l u , l b and l d are Gaussian and centred.

Fig. 4 .

 4 Fig. 4. λ ˆ vs. b σ when l d and l b are Gaussian and centred and l u is Gaussian with nonzero mean.

  the presence of measurement noise at the quantizer's input, one cannot hope to achieve a consistent estimation of the nominal parameters θ by minimizing ∞

  distributions with zero mean: u σ is set equal to 1, whereas d σ and b σ may vary. In order to satisfy the asymptotic conditions, one must impose L N >> (Colinet and Juillard, 2008). For each experiment, θ ˆ is estimated with N G 1 and N J 1 , thanks to a BFGS quasi-Newton method, implemented by the fminunc function of Matlab. Each experiment is repeated a large number of times, in order to estimate λ ~, the expected value of λ , for different values of d σ and b σ . More precisely, for given d σ and b σ ,

  the estimated value of λ decreases, which illustrates Property 4. This shows that this criterion is not suitable for parameter estimation in a noisy context. Moreover, λ is not only a function of b σ but also of d σ : as can be expected, the larger d σ is with respect to b σ , the closer λ is to 1.

Preprints of the 15th IFAC Symposium on System IdentificationSaint-Malo, France, July 6-8, 2009


Mohammed Nabil El

Guillaume Bouleux
email: guillaume.bouleux@univ-st-etienne.fr

Rémy Boyer
email: remy.boyer@lss.supelec.fr

Sylvie Marcos
email: marcos@lss.supelec.fr

SEQUENTIAL ESTIMATION OF THE RANGE AND THE BEARING USING THE ZERO-FORCING MUSIC APPROACH

In this paper, we consider the range and bearing estimation of near-field narrow-band sources from noisy data observed across a passive sensor array. For some difficult scenarios as for correlated and largely spaced sources at low SNRs, or correlated and closely spaced sources, the Near FieLd (NFL) version of the MUltiple SIgnal Classification (MUSIC) algorithm is no longer reliable. In this paper, we adapt and extend the sequential Zero-Forcing MUSIC (ZF-MUSIC) algorithm, which avoids the delicate search of multiple maxima, to the NFL context. In order to compare the NFL ZF-MUSIC with the Second-Order Statistics Weighted Prediction (SOS-WP) algorithm, we assumed an uniform sampling in the spatial domain. However, the proposed algorithm can be exploited for general array geometries. Finally numerical simulations show that the variance of the proposed algorithm achieves the Cramér-Rao Bound (CRB) in difficult scenarios and for sufficient Signal to Noise Ratio (SNR).

INTRODUCTION

The localization of radiating sources by an array of passive sensors arises in several applications, such as radar, sonar, seismology, geophysics, radio astronomy and oceanography. A variety of algorithms [START_REF] Schmidt | Multiple emitter location and signal parameter estimation[END_REF] has been proposed to solve the estimation problem of the bearing of far-field narrowband sources. In this scenario, the propagating waves can be considered to have a planar wavefront when they reach the sensor array. However, when the sources are located in the near-field region, the waves impinging on the sensors can no longer be assumed to have a planar wavefront. Thus the techniques based on the far-field assumption show unsatisfactory performance. A good approximation of the nonlinear propagation delay can be done by a second order Taylor expansion, called the Fresnel approximation. Many localization methods for far-field sources can be extended to the NFL context. Among all the existing methods [START_REF] Van Der Veen | Subspace based signal analysis using Singular Value Decomposition[END_REF], we focus on sequential estimation schemes. There exists essentially two families of sequential estimation schemes: [START_REF] Schmidt | Multiple emitter location and signal parameter estimation[END_REF] based on the Fourier transform, such as the RELAX approach [START_REF] Li | Angle and waveform estimation via RELAX[END_REF] -this algorithm uses a relaxation based method to minimise the maximum likelihood function-and (2) based on a high-resolution scheme such as the Sequential MUSIC This project is funded by both region Île de France and Digeteo Research Park.

algorithms [START_REF] Oh | A sequential estimation approach for performance improvement of eigenstructurebased methods in array processing[END_REF]. Those algorithms adopt a step by step deflation procedure to reduce the signal subspace dimension. At each step, the observed signal is then projected onto the null space of the previously estimated bearings.The zeroforcing approach, called ZF-MUSIC and introduced recently in the far-field case [START_REF] Bouleux | Zero-forcing based sequential MUSIC algorithm[END_REF], is based on a different principle. It is a sequential high-resolution MUSIC approach which leaves unchanged the signal subspace but scales appropriately the MUSIC pseudo-spetrum. In this paper, we present the extended method of the sequential ZF-MUSIC algorithm in the NFL context. Results and simulations show that the proposed algorithm has a lower variance than the standard NFL MU-SIC [START_REF] Huang | Near-field multiple source localization by passive sensor array[END_REF] and the Second-Order Statistics Weighted Prediction (SOS-WP) [START_REF] Grosicki | A weighted linear prediction method for near-field source localization[END_REF] algorithms. Furthermore, the variance of the NFL ZF-MUSIC achieves the CRB in difficult scenarios and for sufficient SNR.

PROBLEM FORMULATION

Consider M radiating near-field and narrowband sources observed by an Uniform Linear Array1 (ULA) of N sensors with interelement spacing d. The signal at the output of the (n + 1) th sensor can be described by

x n (t) = M ∑ l=1 s l (t)e jτ nl + v n (t), t = 1, . . . , T, n = 0, . . . , N -1,
where T (T > N) is the number of snapshots, s l (t) is the l th radiating signal, v n (t) is the additive complex white Gaussian noise, which is assumed to be temporally and spatially uncorrelated, of density probability function N (0, σ 2 ) where σ is a positive real number and τ nl denotes the time delay associated with the signal propagation time from the l th source to the (n + 1) th sensor. Consequently τ nl can be expressed as [START_REF] Huang | Near-field multiple source localization by passive sensor array[END_REF]

τ nl = 2πr l λ 1 + n 2 d 2 r 2 l - 2nd sin θ l r l -1 , (1)
where r l and θ l are the range and the bearing of the l th source, and λ is the signal wavelength. A good approximation of the time delay τ nl , which corresponds to the so-called Fresnel region, i.e., 0.62(d 3 (N -1) 3 /λ ) 1/2 < r l < 2d 2 (N -1) 2 /λ , consists of the second-order Taylor expansion, which is given by [START_REF] Grosicki | A weighted linear prediction method for near-field source localization[END_REF]

τ nl ≈ ω l n + φ l n 2 + O d 2 r 2 l , (2)
where O d 2 r 2 l represents terms of order larger or equal to

d 2 r 2 l
, and ω l and φ l are the so-called electric angles which are given by ω l = -2π d λ sin(θ l ), and φ l = π d 2 λ r l cos 2 (θ l ). Then the signal can be expressed as

x n (t) = M ∑ l=1 s l (t)e j(ω l n+φ l n 2 ) + v n (t). (3)
Rather than estimating the range and bearing directly, we will estimate the electric angles. Once the estimation of the electric angles is done, we can deduce easily the range and the bearing according to [START_REF] Grosicki | A weighted linear prediction method for near-field source localization[END_REF]:

θ l = -arcsin ω l λ 2πd , (4)
r l = πd 2 λ φ l cos 2 (θ l ). (5)
Consequently the received signal vector x(t) can be expressed as x(t) = As(t) + n(t), where

x(t) = [x 1 (t) . . . x N (t)] T , s(t) = [s 1 (t) . . . s M (t)] T and v(t) = [v 1 (t) . . . v N (t)]
T denote the received signal vector, the source signal vector and the noise vector, respectively. The array manifold matrix is given by

A = [a(ω 1 , φ 1 ) . . . a(ω M , φ M )], (6)
where the elements of the steering vector are given by

[a(ω l , φ l )] n+1 = e j(ω l n+φ l n 2 ) , n = 0, . . . , N -1. (7)
The problem here addressed is then the estimation of the couple {ω l , φ l }, l = 1, . . . , M given the noisy data X = [x(1) x(2) . . . x(T )].

NFL ZERO-FORCING SEQUENTIAL MUSIC ALGORITHM

The sequential ZF-MUSIC algorithm [START_REF] Bouleux | Zero-forcing based sequential MUSIC algorithm[END_REF] can be modified to estimate the electric angles. Let RX be the sample spatial covariance of the noisy observation. According to the data model (3), we have

RX = 1 T XX H , (8)
and we assume that lim

T →∞ RX = E[XX H ] = AR S A H + σ 2 I,
where R S is the source covariance matrix. We sort the eigenvalues λ i of the matrix RX in descending order

λ 1 ≥ • • • ≥ λ M ≥ λ M+1 ≥ • • • ≥ λ N , then the expression (8) can be writ- ten as RX = N ∑ n=1 λ n u n u H n , (9)
where u n is the eigenvector associated with the eigenvalue λ n . The noise-subspace is then spanned by U = [u M+1 . . . u N ] and the noise projector is defined by

Π ⊥ = UU H = ∑ N n=M+1 u n u H n .
Consequently, the estimate of the electric angles of the l th source are defined as the maximizers of the spectral form of the sequential NFL ZF-MUSIC given by

( ωl , φl ) = arg max (ω,φ ) f (N) l (ω, φ ) a(ω, φ ) H Π ⊥ a(ω, φ ) , (10)
where l ∈ [1 . . . M] and the zero-forcing function is a quadratic function defined by

f (N) l (ω, φ ) = a N (ω, φ ) H P ⊥ l a N (ω, φ ), (11)
where a N (ω, φ ) = 1 √ N a(ω, φ ) and P ⊥ l = I -P l is the orthogonal projector onto the space spanned by the l -1 previously estimated electric angles ( ω1 , φ1 ), . . . , ( ωl-1 , φl-1 ), where

P l = I A l (A H l A l ) -1 A H l for l = 2, . . . , M and
A l = [a N ( ω1 , φ1 ) . . . a N ( ωl-1 , φl-1 )] for l = 2, . . . , M
The search of the maximizers of (10) can be done by defining a set of discrete points spanning the electric angles intervals of interest and by performing a two dimensional search. However, according to (4) and ( 5), the Fresnel region, must be taken into account. Then the sequential NFL ZF-MUSIC algorithm can be described as follow:

Init. Apply the spectral form of the sequential NFL ZF-MUSIC, with P ⊥ 1 = I, i.e., f

(N) 1 (ω, φ ) = 1.
Then, compute the projector P ⊥ 2 = Ia N ( ω1 , φ1 )a N ( ω1 , φ1 ) H . Loop. For l ∈ [2, . . . , M], compute the zero-forcing function according to (11) and solve the criterion (10).

For a fast computation, instead of performing a costly 2-D search, we can discretize our 2-D search with a large gap according to the Fresnel region. The previously estimated electric angles can be used as initial points. We then refine our estimation by using an unconstrained nonlinear optimization, like the fminsearch matlab function which is a simplex direct search method [8]. Finally, Fig. 4 shows that increasing the number of sensors leads to decreasing the main lobe width. Furthermore, we notice that the peak is more selective for the second pulsation.

SIMULATIONS

In this section, the NFL ZF-MUSIC algorithm is compared with the standard NFL MUSIC algorithm and with the SOS- WP (Second-Order Statistics Weighted Prediction) introduced by Grosicki et al. [START_REF] Grosicki | A weighted linear prediction method for near-field source localization[END_REF]. The context of these simulations is an Uniform and Linear Array (ULA) of 12 sensors spaced by a quarter wavelength. Thus, according to (4) and ( 5), the Fresnel region is delimited by: |ω| < 1.5708 rad and 0 < φ < 0.0694 rad. Two sources are assumed to be nearfield narrowband complex circular Gaussian sequences with zero mean and variance, σ 2 s , equal to one. Consequently, the SNR and the source covariance matrix are defined, respectively, as SNR[dB] = 10 log 10 ( 1 σ 2 ), and

R S = 1 γ γ * 1 ( 12
)
where |γ| determines the degree of correlation between these two sources. where |γ| determines the degree of correlation between these two sources. The simulations were run with T = 15 snapshots and N trial = 100 Monte Carlo trials. For each scheme, the Mean Squared Error (MSE) for the first electric angle is evaluated by

MSE(ω) = 1 2N trial [∑ N trial i=1 ( ω1 (i) -ω 1 ) 2 + ∑ N trial i=1 ( ω2 (i) -ω 2 ) 2 ]
, where ω1 (i) and ω2 (i) represent the estimates of ω 1 and ω 2 repectively. The MSE for the second pulsation is evaluated in the same way. As a comparison, we have also plotted the corresponding Cramér-Rao Bounds (CRB) [START_REF] Grosicki | A weighted linear prediction method for near-field source localization[END_REF].

In the first case, we consider two largely spaced sources which belong to the Fresnel zone, i.e., (ω 1 , φ 1 ) = (0.5, 0.015) rad and (ω 2 , φ 2 ) = (1.2, 0.06) rad. According to Fig. 5 and Fig. 6, where we have considered different degrees of correlation between the two largely spaced sources, the SOS-WP algorithm cannot estimate the electric angles with a good accuracy even with high SNRs. The NFL ZF-MUSIC and the standard NFL MUSIC algorithms are equivalent for SNRs higher than 20 dB. However, we note that the standard NFL MUSIC algorithm is in good agreement with the CRB from about 15 dB, unlike the sequential NFL ZF-MUSIC algorithm which is more accurate at low SNRs and continues to be in very good agreement with the CRB from 0 dB and 5 dB when the sources are moderately and highly correlated, respectively.

In the second case, we consider two closely spaced sources which belong to the Fresnel zone, ie., (ω 1 , φ 1 ) = (0.78, 0.019) rad and (ω 2 , φ 2 ) = (0.85, 0.063) rad. According to Fig. 7 and Fig. 8, where we have considered different values of γ for the two closely spaced sources, the standard NFL MUSIC and the SOS-WP algorithms are no longer reliable. However, the sequential NFL ZF-MUSIC algorithm is in good agreement with the CRB of the first electric angle from 5 dB and 10 dB for moderately and highly correlated sources, respectively. It is the same case for the CRB of the second pulsation from 5 dB and 15 dB for moderately and highly correlated sources, respectively.

CONCLUSION

A sequential algorithm for the estimation of the range and bearing of near-field sources has been proposed. The algorithm is an extension of the far-field zero-forcing MUSIC algorithm presented recently. Without any deflation of the signal subspace, the MUSIC criterion is scaled by an appropriate function which removes the previously estimated ranges and bearings. The proposed algorithm has been tested in some critical scenarios, e.g., for highly correlated sources, and/or for closely spaced sources, the results have shown that the variance of the NFL ZF-MUSIC algorithm achieves the CRB for sufficient SNR and outperforms the standard NFL MUSIC and the SOS-WP algorithms.

Fig. 1 .

 1 (a) and Fig.1.(b) represent the plot of the zeroforcing function for one source, (ω, φ ) = (1.2, 0.05) rad, and then for two sources, (ω 1 , φ 1 ) = (0.5, 0.015) rad and (ω 2 , φ 2 ) = (1.4, 0.06) rad, respectively. From Fig.1.(a), we note that the zero-forcing function is equal to zero for the considered source without affecting too much the other area, and the zero-forcing function is almost equal to one everywhere else. Fig.1.(b) shows that the latter function has two reversed peaks equal to zero for (ω 1 , φ 1 ) and (ω 2 , φ 2 ), this represents the cancellation of the zero-forcing function for the two considered electric angles. As the previous scheme, this latter function is almost equal to one everywhere else. Consequently, Fig.1suggests that the zero-forcing function f (N) l (ω, φ ) is equal to zero for all the previously estimated electric angles without affecting too much the other area.

Figure 1 :

 1 Figure 1: Zero-forcing function for N = 12 and ; (a) one source, (b) two sources.

Figure 2 :

 2 Figure 2: The longitudinal and cross section of the normalised NFL MUSIC spectrum and the normalised NFL ZF MUSIC spectrum at : (a)ω 1 = 1.4 rad, (b)ω 2 = 0.5 rad, (c)φ 1 = 0.015 rad, (d)φ 2 = 0.06 rad.

Figure 3 :

 3 Figure 3: The longitudinal and cross section of the normalised NFL MUSIC spectrum and the normalised NFL ZF MUSIC spectrum at : (a)ω 1 = 1.4 rad, (b)ω 2 = 1.1 rad, (c)φ 1 = 0.015 rad, (d)φ 2 = 0.05 rad.

Figure 4 :

 4 Figure 4: The width of the zero-forcing function corresponding to an attenuation of 3 and 6 dB for (ω, φ ) = (1.3, 0.05) rad and different values of N ; (a) the width of the main lobe of the longitudinal section at φ = 0.05 rad, (b) the width of the main lobe of the cross section at ω = 1.3 rad.

Figure 5 :

 5 Figure 5: MSE(ω rad) vs. SNR[dB], Largely spaced sources (ω 1 = 0.5, φ 1 = 0.015) rad and (ω 2 = 1.2, φ 2 = 0.06) rad, with 12 sensors and 15 snapshots ; (a) γ = 0, (b) γ = 0.3, (c) γ = 0.5, (d) γ = 0.8.

Figure 6 :

 6 Figure 6: MSE(φ rad) vs. SNR[dB], Largely spaced sources (ω 1 = 0.5, φ 1 = 0.015) rad and (ω 2 = 1.2, φ 2 = 0.06) rad, with 12 sensors and 15 snapshots ; (a) γ = 0, (b) γ = 0.3, (c) γ = 0.5, (d) γ = 0.8.

Figure 7 :

 7 Figure 7: MSE(ω rad) vs. SNR[dB], Closely spaced sources (ω 1 , φ 1 ) = (0.78, 0.019) rad and (ω 2 , φ 2 ) = (0.85, 0.063) rad, with 12 sensors and 15 snapshots ; (a) γ = 0, (b) γ = 0.3, (c) γ = 0.5, (d) γ = 0.8.

In order to compare the NFL ZF-MUSIC with the SOS-WP algorithm[START_REF] Grosicki | A weighted linear prediction method for near-field source localization[END_REF], we assumed an uniform sampling in the spatial domain. However, this assumption can be easily relaxed to more general geometries.