Jinjun Xiong

Kingho Tam
email: ktam@ee.ucla.edu.

Lei He

Buffer Insertion Considering Process Variation *

A comprehensive probabilistic methodology is proposed to solve the buffer insertion problem with the consideration of process variations. In contrast to a recent work, we point out, for the first time, that the correlation between the required arrival time and the downstream loading capacitance must be considered in order to solve the problem "correctly". We develop an efficient bottom-up recursive algorithm to calculate the joint probability density function that accurately captures the above correlation, and propose effective pruning rules to exclude probabilistically inferior solutions. We verify our buffer insertion using timing analysis with both device and interconnect variations, and show that compared to the conventional buffer insertion algorithm using nominal device and interconnect parameters, our new buffer insertion methodology can reduce the probability of timing violation by up to 30%.

INTRODUCTION

Ultra deep submicron process exhibits significant variations [START_REF] Borkar | Parameter variations and impacts on circuits and microarchitecture[END_REF] [START_REF] Bowman | Impact of die-to-die and within-die parameter fluctuations on the maximum clock frequency distribution[END_REF]. One of the causes is due to critical dimension being scaled quicker than the development of its controlling process technology [START_REF] Visweswariah | Death, taxes and failing chips[END_REF]. The variability of physical dimensions such as the effective channel length of a transistor is therefore proportionately increasing. As a result of these effects, the performance of manufactured circuits differs significantly from what a circuit simulator predicates under nominal circuit parameters.

Studies on process variations have been mainly focused on variability modeling and statistical timing analysis (STA). For example, [START_REF] Jess | Statistical timing for parametric yield prediction of digital integrated circuits[END_REF] proposed three approaches to compute the probability density functions (p.d.f.) that describe the timing distribution of the circuits. [START_REF] Agarwal | Computation and refinement of statistical bounds on circuit delay[END_REF] advocates to estimate circuit timing variations by computing bounds. However, there are very few works on design optimization that consider process variation effects. For example, essential to high performance circuit design, buffer insertion has been studied extensively in literature e.g., [START_REF] Van Ginneken | Buffer placement in distributed RC-tree networks for minimal Elmore delay[END_REF] and [START_REF] Alpert | Wire segmenting for improved buffer insertion[END_REF]. Yet, none of them has considered the effect of process variations. To the best of our knowledge, the only work that entails such a flavor is by [START_REF] Khandelwal | A probabilistic approach to buffer insertion[END_REF], where a buffer insertion problem that considers the effect of wire length variation has been formulated. However, as we shall point out in Section 3, an unrealistic assumption is made in [START_REF] Khandelwal | A probabilistic approach to buffer insertion[END_REF] in an attempt to solve such a problem.

The major contribution of this paper is two-fold: (1) we formulate a buffer insertion problem with the consideration of process variations for both devices and interconnects; and

(2) we propose a comprehensive probabilistic methodology to solve the above problem efficiently. Preliminary experiment results have shown that by using our methodology, the timing constraints can be satisfied with a much higher probability than the conventional design under nominal values [START_REF] Van Ginneken | Buffer placement in distributed RC-tree networks for minimal Elmore delay[END_REF].

The rest of the paper is organized as follow. Section 2 gives the preliminary and modeling aspects of this work. We formulate the problem for buffer insertion considering process variation (BIPV) in Section 3. Section 4 discusses the details of the algorithm and the methodology for solving the BIPV problem. We present experiment results in Section 5 and conclude the paper with discussion of our future work in Section 6.

PRELIMINARY AND MODELING

Buffer Insertion Preliminary

For simplicity of presentation, we follow the same argument as [START_REF] Van Ginneken | Buffer placement in distributed RC-tree networks for minimal Elmore delay[END_REF] by assuming that the routing tree is given as a binary routing tree and the legal buffer positions (nodes) are directly after the branching points of the tree1 . For a given buffered routing tree, we associate every legal buffer position t in the tree with two numbers: the input loading capacitance (or downstream loading capacitance) Lt and the required arrival time Tt. Denoting c and r as interconnect's unit length capacitance and sheet resistance, respectively, we model each interconnect segment in the routing tree with length l i as a π model, where the resistance is given by r × li, and the capacitance is given by c × li. For a given technology, we associate each buffer from the library with three numbers: the input capacitance C b , the output resistance R b and the intrinsic delay T b . Under the Elmore delay model, the L t and T t can be computed as follows.

If node t is obtained by adding a wire of length li at its direct downstream node n, then

Lt = Ln + c • li (1)
T t = T n -r • l i • L n - 1 2 • r • c • l 2 i . (2)
If node t is obtained by adding a buffer at its direct down-stream node n, then

Lt = C b (3) Tt = Tn -T b -R b • Ln. (4)
If node t is obtained by merging two nodes m and n, then

Lt = Ln + Lm (5) Tt = min(Tn, Tm). (6)
It has been proved in [START_REF] Van Ginneken | Buffer placement in distributed RC-tree networks for minimal Elmore delay[END_REF] that the buffer insertion problem, without considering process variation, can be solved optimally via dynamic programming. Moreover, by properly defining the dominance relationship (or pruning rule) between two solutions, i.e., solution (L1, T1) dominates solution (L 2 , T 2 ) if condition L 1 < L 2 and T 1 > T 2 are satisfied, [START_REF] Van Ginneken | Buffer placement in distributed RC-tree networks for minimal Elmore delay[END_REF] proved that by keeping only dominating solutions at every node, the dynamic programming approach can solve the problem in polynomial time without loosing optimality.

Device and Interconnect Variations

Device variation and interconnect variation are considered as two typical effects of process variations in this work. Because of the process variations in gate thickness, doping density and channel length, buffers' C b , R b , and T b deviate from their respective nominal values. We model such an effect by describing them as random variables that are characterized by a joint probability density function (j.p.d.f.) g(C b , R b , T b ), whose domain is given by ΩC b b ,T b . Due to variations in masking, the manufactured interconnect width and space are also different from their nominal values. We model such an effect by assuming that interconnect's c and r are random variables and that they can be described by the j.p.d.f. as h(c, r), whose domain is given by Ωc,r. As the mechanisms of inducing variations on devices and interconnect are different, it is reasonable to assume that the device and interconnect variations are mutually independent.

BIPV PROBLEM FORMULATION

Formulation 1. Buffer Insertion considering Process Variation (BIPV) Problem: Given the topology of a routing tree with parasitic capacitance and resistance, required arrival times and loading capacitances specified at all sinks, determine the placement of buffers in the routing tree such that the probability of the required arrival time at the root meeting the design specification is maximized with the consideration of process variations for both interconnect and devices, and as a secondary objective, the number of buffers used are minimized.

Even though the conventional buffer insertion problem can be solved optimally via dynamic programming [START_REF] Van Ginneken | Buffer placement in distributed RC-tree networks for minimal Elmore delay[END_REF], solving BIPV problem "optimally" becomes much complicated. The difficulty lies in that:

• Lt and Tt are no longer fixed values, but two random variables. How do we compare two random solutions to tell that one is better than the other?

• Lt and Tt are interdependent, as both are complicated functions of all downstream random variations. Therefore, j.p.d.f., instead of their respective p.d.f., should be used to correctly characterize Lt and Tt.

Obviously, the conventional deterministic dominance relationship between two solutions is no longer applicable. A straight forward extension is to use the mean value as a metric to compare two random solutions. Despite the simplicity of this metric, however, by using the following small example, we show that this simple parametric scheme do not suffice to provide a good guideline for designers to distinguish different design alternatives.

Suppose we are given four buffer insertion solutions, and for simplicity, we assume that all solutions have the same loading capacitance L with different required arrival time Tt, which is described by a discrete probability mass function. The four solutions are: (A) (L, T 1 ) with P(T 1 = 100) = 1; (B) (L, T2) with P(T2 = 200) = 0.5 and P(T2 = 0) = 0.5; (C) (L, T3) with P(T3 = 1000) = 0.1 and P(T3 = 0) = 0.9; (D) (L, T4) with P(T4 = 200) = 0.9 and P(T4 = -800) = 0.1.

It is obvious that the above four solutions have very different characteristics and designers definitely want to distinguish them for a given design. However, after some simple calculations, we find that all four solutions result in the same mean value of 100 for the required arrival time Tt. Therefore, if the metric is purely based on the mean values, we cannot distinguish the above four designs at all! It may seem tempting for people to argue that they can extend the mean value metric (first order moment) by including variance (second order moment) or even higher order moments as additional metrics to compare solutions. But it is easy to construct cases such that none of the above simple metrics would work in order to distinguish different solutions with variations 2 .

Another common practice that people tend to adopt is to ignore the dependence between Lt and Tt. For example, in [START_REF] Khandelwal | A probabilistic approach to buffer insertion[END_REF] where the wire length is the only source of variation3 , it is assumed that there is no correlation between L t and T t , so that the pruning rules can be carried out based on simple multiplication of their respective p.d.f.s. However, this assumption does not hold in general. For example, if a downstream node n is connected to an upstream node t through an interconnect of length l, then from (1) and ( 2), we can clearly see that any increase in l by δ increases L t by cδ while decreases T t by rδ + 1 2 rc(2δl + δ 2 ). Treating T t and Ct as two independent parameters ignores the fact that increasing l in fact collectively decreases Tt and increases Ct, thus causing the solution to degrade much more quickly than before. Therefore, the j.p.d.f. must be used to correctly describe L t and T t 's interdependence.

BIPV ALGORITHMS

Algorithm Overview

We follow the same dynamic programming paradigm as [START_REF] Van Ginneken | Buffer placement in distributed RC-tree networks for minimal Elmore delay[END_REF] to solve the BIPV problem. The overall algorithm (BIPV-ALG) is shown in the top block of Fig. 1. We first traverse the routing tree bottom-up once and build a set of dominant BIPV solutions for all sub-trees. After we reach the root, we pick an optimal solution from the set of kept dominant solutions, optimizing for the required arrival time and the number of buffers inserted. We then back-track the chosen optimal solution to determine the solution for each sub-tree recursively. The key part to this algorithm is the bottom-up traversal of the routing tree. Therefore, we also present the detailed bottom-up algorithm in the bottom block of Fig. 1.

BIPV-ALG(t)

Input: root of the routing tree t.

Output: Solution to BIPV problem. Z = BIPV-BOTTOM-UP-ALG(t);

Z * = PICK-BEST-SOL(Z); BACK-TRACK-SOL(Z * ); BIPV-BOTTOM-UP-ALG(t) /*

Computation of J.P.D.F.

According to (1) and (2), or (3) and ( 4), we know that the upstream variations do not change the distribution of Lt and Tt, as Lt and Tt are only functions of node t's downstream random variations. Therefore, L t and T t 's distributions are independent of their upstream random variations. Moreover, if (Lt1, Tt1) and (Lt2, Tt2) share some (or do not share any) common down-stream paths, then Lt1, Tt1, Lt2, and Tt2 are mutually correlated (or independent).

Even thought (Lt, Tt) depend on their downstream random variations in a complicated way, the way we compute (L t , T t ) is via the recursive equations according to either (1) and (2), or (3) and (4), or ( 5) and [START_REF] Van Ginneken | Buffer placement in distributed RC-tree networks for minimal Elmore delay[END_REF]. Therefore, we develop an efficient algorithm to compute the j.p.d.f. of (Lt, Tt) in a recursive fashion while we traverse the routing tree. This is particularly useful in the context of dynamic programming, because only incremental computation is necessary at each node. To develop the recursive computation formula, we employ the multivariate transformation technique [START_REF] Hogg | Introduction to Mathematical Statistics[END_REF]. We denote the j.p.d.f. of node t's direct downstream nodes as fL n ,T n with its domain given by ΩL n ,T n .

J.P.D.F. after Adding a Wire

To compute the j.p.d.f. of (L t , T t ) obtained from ( 1) and (2) via the multivariate transformation technique, we introduce two new random variables X and Y as follows:

X = c; ( 7 ) Y = r; ( 8 )
Lt = Ln + c • li; ( 9 )
Tt = Tn -r • li • Ln - 1 2 • r • c • l 2 i . (10)
By transforming variables, we obtain c = X;

( 1 1 ) r = Y ;

( 1 2 )

L n = L t -X • l i ; ( 1 3
)

Tn = Tt + Y • li • Lt - 1 2 • Y • X • l 2 i . (14)
It is easy to verify that the mapping between (c, r, Ln, Tn) and (X, Y, L t , T t ) is a one-to-one mapping. Therefore, the Jacobian is given by

J = ∂c ∂X ∂c ∂Y ∂c ∂L t ∂c ∂T t ∂r ∂X ∂r ∂Y ∂r ∂L t ∂r ∂T t ∂L n ∂X ∂L n ∂Y ∂L n ∂L t ∂L n ∂T t ∂T n ∂X ∂T n ∂Y ∂T n ∂L t ∂T n ∂T t = 1 0 0 0 0 1 0 0 -li 0 1 0 1 2 • Y • l 2 i li • Lt -1 2 • X • l 2 i Y • li 1 = 1. (15)
As c and r are independent of Ln and Tn, the j.p.d.f. of (c, r, Ln, Tn) is given by: fc,r,L n ,Tn = h(c, r) • fL n ,Tn . Therefore, we obtain the j.p.d.f. of (X, Y, Lt, Tt) as follows:

fX,Y,L t,Tt = |J| • fc,r,L n,Tn (c, r, Ln, Tn)

( 1 6 ) = h(c, r) • fL n,Tn (Ln, Tn) = h(x, y) • fL n ,T n (Lt -X • li, Tt + Y • li • Lt - X • Y • l 2 i 2 )
Then the j.p.d.f. of (L t , T t ) is obtained by

fL t ,Tt = ¡ Ω X,Y (L t ,T t ) fX,Y,L t ,Tt dXdY, (17)
where Ω X,Y (L t , T t ) is the domain for X, Y in terms of L t and Tt. Knowing c and r's domain Ωc,r and (Ln, Tn)'s domain ΩL n,T n, we can deduce ΩX,Y (Lt, Tt) according to ( 7), (8) ( 9) and (10).

J.P.D.F. after Adding a Buffer

To compute the j.p.d.f. of (Lt, Tt) obtained from (3) and (4) via the multivariate transformation technique, we introduce three new random variables X, Y and Z as follows:

X = L n ; ( 1 8 ) Y = R b ; ( 1 9 ) Z = T b ;
( 2 0 )

Lt = Cb; ( 2 1 ) Tt = Tn -Tb -Rb • Ln. ( 22
)
By transforming variables, we obtain

Cb = Lt; ( 2 3 ) R b = Y ; ( 2 4 ) T b = Z;
( 2 5 ) Ln = X;

( 2 6 )

Tn = Tt + Z + Y • X. (27)
It is easy to show that the mapping between (C b , R b , T b , Ln, Tn) and (X, Y, Z, L t , T t ) is a one-to-one mapping. Therefore, the Jacobian is given by

J = ∂C b ∂X ∂C b ∂Y ∂C b ∂Z ∂C b ∂L t ∂C b ∂T t ∂R b ∂X ∂R b ∂Y ∂R b ∂Z ∂R b ∂Lt ∂R b ∂Tt ∂T b ∂X ∂T b ∂Y ∂T b ∂Z ∂T b ∂L t ∂T b ∂T t ∂L n ∂X ∂L n ∂Y ∂L n ∂Z ∂L n ∂L t ∂L n ∂T t ∂T n ∂X ∂T n ∂Y ∂T n ∂Z ∂T n ∂L t ∂T n ∂T t = 0 0 0 1 0 0 1 0 0 0 0 0 1 0 0 1 0 0 0 0 Y X 1 0 1 = -1. ( 28
)
As C b , R b , and T b are independent of Ln and Tn, the j.p.d.f. of (C b , R b , T b , L n , T n ) is given by: f C b ,R b ,T b ,Ln,Tn =g(C b , R b , T b ) • f Ln,Tn . Therefore, we obtain the j.p.d.f. of (X, Y, Z, Lt, Tt) as follows:

fX,Y,Z,L t ,T t = |J| • fC b ,R b ,T b ,Ln,Tn = g(C b , R b , T b ) • fL n ,T n (29) = g(Lt, Y, Z) • fL n ,T n (X, Tt + Z + Y • X).
Then j.p.d.f. of (L t , T t ) is obtained by

fL t ,T t = ¡ Ω X,Y,Z (Lt,Tt)
fX,Y,Z,L t ,T t dXdY dZ, (30)

where Ω X,Y,Z (X, Y, Z) is the new domain for X, Y and Z in terms of L t and T t . Knowing C b , R b and T b 's domain ΩC b ,R b ,T b and (Ln, Tn)'s domain ΩL n ,T n, we can deduce ΩX,Y,Z(X, Y, Z) according to (18), ( 19), (20), ( 21) and (22).

J.P.D.F. after Merging Two Solutions

To compute the j.p.d.f. for (Lt, Tt) after merging two solutions according to ( 5) and ( 6), we introduce two new random variables as follows:

X = L m (31) Y = T n + T m (32) Lt = Ln + Lm (33) Tt = min(Tn, Tm) ( 3 4 )
Because of the min-function, there is no one-to-one mapping relation between (Lm, Tm, Ln, Tn) and (X, Y, Lt, Tt), thus we cannot use the multivariate transformation technique to compute the j.p.d.f. of (L t , T t ) directly. But we note that the original domain ΩL m,Tm,Ln ,Tn for (Lm, Tm, Ln, Tn) can be divided into two disjoint sub-domains Ω1 and Ω2, where Ω1 is the sub-domain with Tm ≤ Tn and Ω2 is the other subdomain with Tm > Tn. For each disjointed sub-domain, we can show that there exists a one-to-one mapping between (L m , T m , L n , T n ) and (X, Y, L t , T t ), thus we can compute the j.p.d.f. for each sub-domain by using the multivariate transformation technique. Then by combining the two subdomains' j.p.d.f., we can compute the j.p.d.f. for the whole domain [START_REF] Hogg | Introduction to Mathematical Statistics[END_REF]. We first find the j.p.d.f. for sub-domain Ω1 with Tm ≤ Tn. We have:

Tt = min(Tm, Tn) = Tm (35)
After transformation of variables, we have

Lm = X; ( 3 6 ) L n = L t -X; ( 3 7 ) T m = T t ; ( 3 8 ) T n = Y -T t . ( 39
)
The Jacobian is given by

J = ∂Lm ∂X
= 1 0 0 0 -1 0 1 0 0 0 0 1 0 1 0 -1 = 1. ( 40
)
As (Lm, Tm) is independent of (Ln, Tn), the j.p.d.f. of (Lm, Tm, Ln, Tn) is given by: fL m ,T m ,L n ,T n =fL m ,T m • fL n ,T n , Therefore, we can obtain the j.p.d.f. of (X, Y, Lt, Tt) in sub-domain Ω1 as follows:

fΩ 1 = |J | • fL m ,T m ,L n ,T n = fL m ,T m • fL n ,T n = fL m ,T m (X, Tt) • fL n ,T n (Lt -X, Y -Tt). (41)
Similarly, for the sub-domain Ω2 with Tm > Tn, We have:

fΩ 2 = |J| • fL m,Tm,Ln ,Tn = fL m,Tm (X, Y -Tt) • fL n ,Tn (Lt -X, Tt). ( 42
)
Combining ( 41) and ( 42) together, we obtain the The j.p.d.f. for (X, Y, L t , T t ) in the whole domain Ω X,Y,Lt,Tt as follows [START_REF] Hogg | Introduction to Mathematical Statistics[END_REF]:

f X,Y,Lt,Tt = f Ω 1 (X, Y, L t , T t ) + f Ω 2 (X, Y, L t , T t ) (43) = f Lm ,Tm (X, T t ) • f Ln,Tn (L t -X, Y -T t ) + fL m ,Tm (X, Y -Tt) • fL n ,Tn (Lt -X, Tt)
Then the j.p.d.f. of (L t , T t ) is obtained by

f Lt,Tt (L t , T t ) = ¡ Ω X,Y (Lt,Tt) f X,Y,Lt,Tt dXdY, ( 44
)
where ΩX,Y (Lt, Tt) is the new domain for (X, Y ) in terms of L t and T t . Knowing (L m , T m )'s domain Ω Lm,Tm and (L n , T n )'s domain Ω Ln,T n , we can deduce Ω X,Y (L t , T t ) from (31),(32), ( 33) and (34).

J.P.D.F. Initialization

To carry on the above recursive computation of j.p.d.f., we have to set up the initial conditions starting from sinks. As sink's L s and T s are constants and buffers are always inserted at sink's upstream node with a wire between them, (1) and ( 2) should be used for j.p.d.f. initialization. After transformation of variables, we have:

c = (Lt -Ls)/li; ( 4 5 ) r = 2(Ts -Tt) (Lt + Ls)li . (46)
Therefore, the Jacobian is given by

J = ∂c ∂Lt ∂c ∂Tt ∂r ∂Lt ∂r ∂Tt = 1 0 -2(Ts-Tt) (Lt+Ls) 2 -2 (L t +L s )l i = -2 (Lt + Ls)li . (47)
As the j.p.d.f. for (c, r) is already known as h(c, r), we can obtain the j.p.d.f. of (L t , T t ) as follows:

f Lt,Tt = |J| • h(c, r) = 2 (L t + L s )l i • h( L t -L s l i , 2(T s -T t ) (L t + L s )l i ) (48)

Pruning Rules

A straight-forward way to extend the conventional dominance relationship between two solutions in the presence of process variation is as follows: solution (L1, T1) is said to dominate solution (L2, T2) if condition P (L1 ≤ L2) = 1 and P (T 1 ≥ T 2 ) = 1 are satisfied. In other words, solution (L 1 , T1) always results in a larger required arrival time but with a less loading capacitance when compared to solution (L2, T2). The physical interpretation of this criterion is well understood. However, there are two problems when it comes to practical implementation. First, for a continuous j.p.d.f. for two random variables (L t and T t ), the domain is usually defined over the whole feasible region: i.e., 0 ≤ Lt ≤ ∞ and 0 ≤ Tt ≤ ∞. Therefore, it is almost impossible to satisfy the conditions of P (T1 ≥ T2) = 1 and P (L1 ≤ L2) = 1 for any two given solutions. Second, assuming the first problem can be solved, there is no guarantee that the number of solutions after pruning will increase polynomially. Chances are that it will, most likely, grow exponentially, and this has been experimentally confirmed by [START_REF] Khandelwal | A probabilistic approach to buffer insertion[END_REF]. In the following, we propose a new set of pruning rules that not only solves the above problems, but also supports designers' intuition.

We recognize that for designers, there is always a design goal in their minds when they compare different design alternatives. For the same example as shown in Section 3, if the arrival time for a design is required to be 100, then obviously solution (A) is the best choice, as it satisfies the required arrival time constraints without any uncertainty. On the other hand, if the arrival time is required to be 200, then solution (D) would be the best choice, as this gives the highest probability to satisfy the design requirement.

Based upon the above observations, we give the following definition of dominance relationship between two solutions, which is closely related to designers' willingness to accept uncertainty for a given design. Recall that the (100α) th percentile of a p.d.f. f (x) is a number πα such that the area under f (x) to the left of of πα is α [START_REF] Hogg | Probability and Statistical Inference[END_REF]. That is,

α = ¡ π α -∞ f (x)dx. ( 49
)
In other words, πα gives a measure of designers preference for certainty in choosing the design parameter x in the presence of variations, such that the final design would have x less than πα with (100α)% certainty.

π α l β π u α π u β π l L1 L2 T2 T1 (a) (b) (1) (1)
( Suppose designers choose πα l and πα u as Lt's two percentiles with 0 ≤ α l < αu ≤ 1, and π β l and π βu as Tt's two percentiles with 0 ≤ β l < βu ≤ 1, which reflect designers preference for certainty in choosing different solutions (L t , T t ) to the BIPV problem. Then solution (L 1 , T 1 ) is said to dominate solution (L 2 , T 2 ) if the following conditions are satisfied:

π (1) αu < π (2) α l (50) π

(1)

β l > π (2) βu (51)
Another way to look at this dominance relationship is that solution (L1, T1) has a high probability of producing solutions with a larger required arrival time and a smaller loading capacitance. A graphical interpretation of this dominance relationship is shown in Fig. 2. Knowing α l , α u , β l , and β u , we can compute π α l , π αu , π β l , and π β u according to (49), which requires us to know the p.d.f. of Lt and Tt, respectively. As we already know the j.p.d.f. of (Lt, Tt) from section 4.2, the p.d.f. of Lt and Tt can be computed as two marginal p.d.f.'s of the j.p.d.f., respectively [START_REF] Hogg | Probability and Statistical Inference[END_REF].

EXPERIMENT

A preliminary version of the above BIPV algorithm is implemented in C++ on a Linux machine. Because there is no existing work that can take into account the accurate computation of j.p.d.f. for (L t , T t ) in the context of buffer insertion, direct comparison of our work with existing works is not possible. The only work that considered process variation for buffer insertion is [START_REF] Khandelwal | A probabilistic approach to buffer insertion[END_REF], which does not consider either the interconnect parasitic variation or the device variation. Moreover, it is assumed that all random variables are independent, which is not correct in general. Therefore, in the following, we compare our BIPV design with the conventional design under nominal design parameters (DUN) similar to [START_REF] Van Ginneken | Buffer placement in distributed RC-tree networks for minimal Elmore delay[END_REF]. For a given technology, foundry needs to characterize the j.p.d.f. for both devices and interconnects with consideration of process variations. Without loss of generality, we assume that the j.p.d.f. of both interconnect and devices are available to us as user inputs in the following. Normal distribution has been assumed for all test cases, and all parasitic values are adapted from [START_REF] Alpert | Wire segmenting for improved buffer insertion[END_REF] as shown in Table 1. Different combinations of α l , αu, β l , and βu values have been tested for pruning. Because the choice of these values only reflects designers' preference to certainty and does not significantly affect our conclusions, we only report experiment results with αl = βl = 0.2 and αu = βu = 0.8 in the following.

Global interconnects with multiple terminals are used in our experiment. For simplicity, every test case is assumed to have a uniform required arrival time specified at different sink terminals. To make the test case challenging, we set the required arrival time at sinks to be half of the critical path delay without buffers. Results from five test cases are reported in Table 2. Column 2 gives the total number of terminal pins for test cases that are typical for global interconnects. The range of interconnect length between two adjacent Steiner points in the routing tree is reported in column 3. Column 4 gives the required arrival time specified at all sink terminals. Column 5 and 6 are the results measured at the root of the buffered routing trees from the conventional DUN design, while column 7 and 8 are the results from our BIPV design. Column 5 and 7 are the probabilities of the designed buffered routing tree that would violate the required timing constraints. The number in column 6 and 8 are the total number of buffers inserted.

According to Table 2, it clearly shows that the conventional DUN design exhibits very large timing violation when process variations are considered, while our process variation aware BIPV design consistently results in smaller timing violation than the DUN design. Moreover, when the test case becomes larger, the improvement of BIPV design over the DUN design is more significant. For example, for test case 4, the conventional DUN design gives 43% timing violation, while our BIPV design results in only 13% violation, and the relative reduction is 30%. In terms of the total number of buffers inserted in the final buffered routing tree, the BIPV design uses sightly more buffers than the conventional DUN design, but there is never more than one extra buffer in all experiments. This further confirms that the BIPV algorithm inserts buffers at locations such that the probability of achieving timing constraints is optimized.

CONCLUSION AND FUTURE WORK

We have proposed a probabilistic methodology to solve the problem of buffer insertion considering process variations for both devices and interconnects. A very recent work solved the buffer insertion with wire length uncertainty by assuming that there is no correlation between the variations of required arrival time and downstream loading capacitance. We have shown that this assumption does not hold in general and that the j.p.d.f. should be employed to accurately capture the above correlation. We have developed an efficient algorithm to calculate the j.p.d.f. recursively. Our preliminary experiment results have shown that compared to the conventional buffer insertion under the nominal parameters, the proposed variation aware buffer insertion can reduce the probability of timing violation by up to 30%.

Note that the proposed methodology is very general and can be extended to consider other process variation effects. For example, we have assumed that there is no correlation between device and interconnect variations, and that wires not sharing downstream paths are independent. However, in a modern process, the variation between devices and interconnects may be interdependent, and wires nearby may also exhibit correlation due to chemical-mechanical planarization (CMP). One of our future work is to extend the proposed methodology to handle more general process variations, like across-chip-line-width variation (ACLV). We will also study the process variation aware routing topology generation and device sizing problems in the future.

Figure 1 :

 1 Figure 1: BIPV Algorithm.

Figure 2 :

 2 Figure 2: Graphic interpretation of dominance relationship between (L1, T1) and (L2, T2), where (a) refers to (50) and (b) refers to (51).

Table 1 : Device and interconnect variations.

 1

	Parameter	Mean	Standard Deviation
	Buffer input capacitance	0.04 pF	0.01 pF
	Buffer output resistance	180 Ω	30 Ω
	Buffer intrinsic delay	20 ps	5 ps
	Unit wire capacitance	0.08 fF/µm	0.02 fF/µm
	Unit wire resistance	0.08 Ω/µm	0.02 Ω/µm

Table 2 : Experimental result.

 2

Note that the methodology to be presented in this work does not depend on these assumptions. Proceedings of the Design, Automation and Test in Europe Conference and Exhibition (DATE'05) 1530-1591/05 $

20.00 IEEE

Readers are encouraged to verified that the last two solutions in the above small example in fact have the same variance.

[START_REF] Visweswariah | Death, taxes and failing chips[END_REF] We do not model the interconnect length variation in this work as that is not a typical process variation. Proceedings of the Design, Automation and Test in Europe Conference and Exhibition (DATE'05) 1530-1591/05 $ 20.00 IEEE

Proceedings of the Design, Automation and Test in Europe Conference and Exhibition (DATE'05) 1530-1591/05 $ 20.00 IEEE

Proceedings of the Design, Automation and Test in Europe Conference and Exhibition (DATE'05)

* This work is partially supported by NSF CAREER award CCR-0093273, a UC MICRO grant sponsored by Fujitsu Laboratories of America, Intel and LSI Logic, and a Faculty Partner Award by IBM.


Baohua Wang
email: baohuaw@eecs.umich.edu

Pinaki Mazumder

EM Wave Coupling Noise Modeling Based on Chebyshev Approximation and Exact Moment Formulation

This paper presents a new mathematical approach to modeling EM wave coupling noise so that it can be easily integrated into chip-level noise analysis tools. The new method employs Chebyshev approximation technique to model the distributed sources arising in the Telegrapher's equations due to EM wave coupling. A uniform plane wave illumination metric is provided to determine the order of approximation. Closed-form formulas for the noise transfer functions' moments are derived. By utilizing the formulated moments, reduced order models can be efficiently obtained to generate the induced noise caused by EM wave illumination. The accuracy of the proposed method is verified by Hspice simulation.

Introduction

Signal integrity is of paramount importance for reliable operation of very deep submicron CMOS VLSI chips. Noise analysis is performed by designers to verify the stability and timing specifications of the chip [START_REF] Kumar | Interconnect and noise immunity design for the pentium 4 processor[END_REF], [START_REF] Shepard | Harmony: static noise analysis of deep submicron digital integrated circuits[END_REF]. Previously, IR drop, Ldi/dt noise, coupling noise induced by signal switching were mainly considered [START_REF] Chen | Power supply noise analysis methodology for deep-submicron vlsi chip design[END_REF], [START_REF] Zhou | Timing analysis with crosstalk is a fixpoint on a complete lattice[END_REF]. However, interconnect noises induced by external electromagnetic wave illumination must also be addressed. For example, PC board can carry high-frequency signals through the long metal wires. They often behave as antenna by illuminating the pins of the mounted VLSI chip by systemgenerated EM waves. The induced noise incident on the chip pins can propagate inside the chip to severely affect the signal integrity of the internal circuit. EM effect, therefore, can manifest itself, e.g. as timing failures, due to impairment of signal integrity.

The analysis of the EM wave coupling to long wires is traditionally based on full-wave field solvers or SPICElike simulators [START_REF] Erdin | A spice model for incident field coupling to lossy multiconductor transmission lines[END_REF], [START_REF] Tkatchenko | High-frequency electromagnetic field coupling to long terminated lines[END_REF]. However, such type of approaches are generally unsuitable for chip-level noise analysis due to their prohibitive computational costs. This paper presents a new approach to modeling the external EM wave coupling in signal traces between two circuit components, e.g., buffer/repeaters inside a chip, or I/O ports on PCB. The proposed approach efficiently generates the induced noise waveforms that should be used as noise sources into chiplevel noise analysis flows to verify the signal integrity of the chip illuminated by external EM waves.

One difficulty in the modeling of EM wave coupling noise is the handling of distributed sources in the Telegrapher's equations describing EM wave coupling, since they correspond to an infinite number of lumped sources in the resultant equivalent TL (transmission line) circuit. To handle these sources, the proposed approach employs Chebyshev approximation technique because of its fast convergence and minmax property. After the order of approximation is determined by using the proposed uniform plane wave illumination metric, actual field values will be used for further computations. Therefore, the usual assumption of uniform plane wave illumination in analyzing EM wave coupling is not required any more.

Reduced order models (ROM's) are proposed in the paper to generate noise waveforms. ROM's have been derived for RC, RLC, TL circuits [START_REF] Pillage | Asymptotic waveform evaluation for timing analysis[END_REF], [START_REF] Odabasioglu | prima: passive reduced-order interconnect macromodeling algorithm[END_REF], [START_REF] Wang | The chebyshev expansion based passive model for distributed interconnect networks[END_REF]. They either require efficient moment generation, which is, however, lacked for the TL circuit considering EM wave coupling, or do not provide closed-form representations. This paper presents closed-form representations of the transfer functions' moments for the used TL circuit. Therefore, ROM's can be quickly obtained to generate noise waveforms.

The rest of the paper is organized as follows. Section 2 introduces the modeling EM wave coupling with long interconnects, driver circuit model, and its noise transfer functions. Section 3 introduces the modeling of distributed sources using Chebyshev approximation, and the proposed uniform plane wave illumination metric to determine the order of approximation. Section 4 presents the derived closedform moments' representations, and ROM's of noise transfer functions. Section 5 presents the experimental results.

EM Wave Coupling Noise Modeling

Telegrapher's Equations with EM Coupling

The behavior of the long interconnect (TL) being illuminated by EM waves, is governed by the Telegrapher's equations, with distributed source terms having been incorporated [START_REF] Paul | Analysis of Multiconductor Transmission lines[END_REF]. For a TL of length d and height h, locating in the xz plane given by y = 0, with x axis as its longitudinal direction, the governing equations in s-domain, when the line is illuminated by EM wave, is written as:

d dx V (s, x) + [R(x) + sL(x)] I(s, x) = V f (s, x) (1a) d dx I(s, x) + [G(x) + sC(x)]V (s, x) = I f (s, x) (1b)
V f (s, x) and I f (s, x) denote source terms describing the incident EM wave coupling. They are given as:

V f (s, x) = E i x,0,h (s) -E i x,0,0 (s) • e x - ∂V i f (s, x) ∂x (2a) I f (s, x) = -[G(x) + sC(x)]V i f (s, x). ( 2b
)
V i f (s, x) is the defined incident voltage, which is given, using the incident electric field E i of the EM wave, as:

V i f (s, x) = h 0 E i x,0,z (s) • e z dz.
In the above, e x and e z are unit vectors in the x direction and z direction, respectively.

Equivalent Circuit Description

Let V (s, x) in (1) be represented as the difference of two terms, i.e. V (s, x) = V s f (s, x) -V i f (s, x), where V s f (s, x) is the scattered voltage. (1) is reformulated as

d dx V s f (s, x) + [R(x) + sL(x)] I(s, x) = V s (s, x) (3a) d dx I(s, x) + [G(x) + sC(x)]V s f (s, x) = 0, (3b)
where V s (s, x) is the distributed voltage source, given by

V s (s, x) = E i x,0,h (s) -E i x,0,0 (s) • e x . ( 4
)
Fig. 1 shows the modeling of driver circuit under EM wave illumination based on [START_REF] Kumar | Interconnect and noise immunity design for the pentium 4 processor[END_REF]. The corresponding circuit model used in noise analysis is shown in Fig. 1(b), where effective driving resistance and parasitic capacitance are used.

) 0 , (t v i f ) , ( d t v i f ) , ( x t v s L C (a) Driver circuit using (3). ) 0 , (t v i f ) , ( d t v i f ) , ( x t v s L C D C s R i V (b) Driver circuit Model.

EM Wave Coupling Noise Transfer Functions

Using the principle of linear superposition, the output voltage at loading capacitor C L is represented as the summation of four voltages in s-domain, i.e.

V C L (s) = V i (s)H 1 (s) +V i f (s, 0)H 2 (s)- V i f (s, d)H 3 (s) +V dist (s), (5)
where V dist (s) denotes the voltage at the output produced by the distributed voltage sources, given by

V dist (s) = d 0 V s (s, x)H dist (s, x)dx. ( 6
)
Noise transfer functions H 1 (s), H 2 (s), H 3 (s) and H dist (s, x) are given as follows:

H 1 (s) = Z 0 Z L H A 1 cosh γ d + H B 1 sinh γ d (7a) H dist (s, x) = Z L (Z 0 cosh γ x + Z D sinh γ x ) H A dist cosh γ d + H B dist sinh γ d (7b) H 2 (s) = H dist (s, 0), H 3 (s) = H dist (s, d) (7c)
In the above,

r d = γ • d, γ x = γ • x,
and

H A 1 = Z 0 (R s + Z L + R s Y CD Z L ) H B 1 = Z 2 0 (1 + R s Y CD ) + R s Z L H A dist = Z 0 (Z D + Z L ), H B dist = Z 2 0 + Z D Z L Z L = 1 sC L ,Y CD = sC D , Z D = 1 1/R s + sC D Note that γ = (R 0 + sL 0 )(G 0 + sC 0 )
, which is the propagation constant of the TL, and Z 0 = (R 0 +sL 0 ) (G 0 +sC 0 ) , which is the characteristic impedance of the TL. The PUL (per unit length) parameters of the TL are R 0 , L 0 , G 0 and C 0 .

Handling of Distributed Sources Using

Chebyshev Approximation

Using Chebyshev Approximation

To model voltage V dist (s) generated by those distributed sources, Chebyshev polynomials are used to approximate [START_REF] Pillage | Asymptotic waveform evaluation for timing analysis[END_REF]. V s (s, x)H dist (s, x) in ( 6) is represented using Chebyshev polynomials up to order N as:

V s (s, x)H dist (s, x) ≈ N ∑ k=0 t k (s)T k (2x/d -1). ( 8
)
T k (•) is the k-th order Chebyshev polynomial. The coefficients t k (s) is given by

t k (s) =            N ∑ j=0 V s (s, x j ) N + 1 H dist (s, x j ) k = 0 N ∑ j=0 2V s (s, x j ) N + 1 H dist (s, x j )T k (y j ) k > 0 and y j = cos 2 j + 1 2N + 2 π , x j = d 2 (1 + y j ). (9)
Based on ( 6) and ( 8), Chebyshev approximation of V dist (s) is derived as

V dist (s)= d N+1 N ∑ j=0 V s (s, x j )H dist (s, x j ) f ( j, N), ( 10
)
where

f ( j, N) = 1 + ∑ N k=1 cos k(2 j+1) 2N+2 π 1+cos kπ 1-k 2 .

Uniform Plane Wave Illumination Metric

To determine the order of Chebyshev polynomials required in general EM wave coupling situation, the uniform plane wave illumination case is analyzed at first. In this situation, the electrical field of the incident EM wave is in the following form: v is the phase velocity of the wave in the medium, e f is the unit vector representing the direction of electric field, and e k is the unit vector representing the EM wave propagation direction. By denoting the projection of E i in xz plane as E i xz , β is then defined as the angle between E i and E i xz , θ is defined as the angle between x axis and E i xz , and α is defined as the angle between propagation direction and vector

E i x,y,z (s) = E 0 e f e -
E i × E i xz .
Using the defined uniform plane wave in [START_REF] Zhou | Timing analysis with crosstalk is a fixpoint on a complete lattice[END_REF], and using (4), V s (s, x) is given as

V s (s, x)=E 0 cos β cos θe -se kx x/v e -se kz h/v -1 . ( 12
)
e kx and e kz are the x and z components of unit vector e k respectively. From ( 6) and ( 12), the analytical representation

of V dist (s) is derived as V dist (s) = Z L E 0 cos β cos θ(C + -C -) e -se kz h/v -1 H A dist cosh γ d + H B dist sinh γ d ( 13
)
where ) with driving and loading conditions varied, the Chebyshev approximation method converges very quickly as N becomes larger. When N = 2, the Chebyshev approximations of V dist (s) match the analytical results well up to 6 GHz, the frequency of incident EM wave. The use of N = 3 matches the frequency up to 14 GHz, and the use of N = 5 matches a frequency range of 20 GHz.

C ± = (Z 0 ± Z D ) e ±(r∓se kx /v)d -1 2(γ ∓ se kx /v)
In actual EM wave illumination situation, the incident EM field distribution is arbitrary. Uniform plane wave is employed at first to estimate the required order of Chebyshev approximation. By inserting the PUL parameters of the used interconnection techniques into (13) and [START_REF] Wang | The chebyshev expansion based passive model for distributed interconnect networks[END_REF], N is determined by judging the error between the two calculated results around the maximum EM wave frequency in the actual illumination environment. Next, the actual field values are used to calculate the time-domain noise waveform V dist (t) according to [START_REF] Wang | The chebyshev expansion based passive model for distributed interconnect networks[END_REF].

Noise Transfer Function Modeling

Chebyshev approximation is proposed for modeling distributed sources in the last section, and the noise transfer functions at Chebyshev points are required. In this section, a set of concise moment formulas are proposed for the modeling of these noise transfer functions, so that reduced order models can be used to generate noise waveforms efficiently.

Moments and ROM of H 1 (s)

H 1 (s) is reformulated as H 1 (s) = 1 H C 1 cosh γ d + H D 1 sinh γ d γ d , ( 14
)
where

H C 1 = H A 1 /(Z 0 Z L ) = ∑ ∞ i=0 a C i s i = 1 + s(τ L s + τ D s )
and

H D 1 = H B 1 γ d /(Z 0 Z L ) = ∞ ∑ i=0 a D i s i = R s G d +(τ L d +τ d s )s+(L d C L +τ L d τ D s )s 2 +L d C L τ D s s 3
In the above, W d is the total W for the TL, i.e. W d = W 0 • d, for W ∈ {R, L, G,C}, and τ b a is a time constant defined as

τ b a = R a C b . By using the property cosh x = ∞ ∑ n=0 x 2n (2n)! , sinh x x = ∞ ∑ n=0 x 2n (2n + 1)! (15) H 1 (s) is rewritten as H 1 (s)= 1 ∑ ∞ n=0 (R d +sL d ) n (G d +sC d ) n H C 1 (2n)! + H D 1 (2n+1)! = 1 ∑ ∞ k=0 q k s k (16)
For the most general case that G 0 = 0, i.e. G d = 0, to determine the coefficient q k in (16), the following facts can be utilized:

1. the term s k can be formed by using a C i s i and a D i s i in H C 1 and H D 1 to multiply a term s k-i from (R d + sL d ) n (sC d ) n in (16).

the term s

k-i in (R d + sL d ) n (sC d ) n can be formed by letting n goes from k-i 2
to ki, and using term

(sC d ) n multiply the term s k-i-n in (R d + sL d ) n .
Hence, q k is derived as By (17), the approximation of H 1 (s) can be easily made. In one way, the series in the denominator of ( 16) is truncated to a finite order N, and a rational polynomial in the form of

q k = min(3,k) ∑ i=0 k-i ∑ n= k-i 2 n k -i -n •C n d L k-i-n d R 2n+i-k d (2n + 1)a C i + a D i (2n + 1)! (17)
1 ∑ N k=0 q k s k is used to model H 1 (s).
The other way is to use general Padé approximation procedure with the order of nominator polynomial chosen in a suitable manner. By (17), the moments of H 1 (s) can be generated efficiently for using in the general Padé approximation procedure.

For the completeness of formulation, q k , when G d = 0, is derived from (16) as well. It is given in the following:

q k = min(3,k) ∑ i=0    k-i ∑ n= k-i 2 n ∑ j=k-i-n n j n k -i -j •R n-j d L j d G n-k+i+ j d C k-i-j d (2n + 1)a C i + a D i (2n + 1)! + ∞ ∑ n=k-i k-i ∑ j=0 n j n k -i -j R n-j d L j d •G n-k+i+ j d C k-i-j d (2n + 1)a C i + a D i (2n + 1)! ( 18
)
Fig. 3 compares |H 1 ( f )| using analytical formula (7a) and using the truncated order representation of H 1 (s), the first approach. Here the TL has the same parameters as in Fig. 2. C D = 0 and the driver's resistance is varied from 400 Ω to 10 Ω. In both of driving conditions, the use of N = 18 matches the analytical results in the whole 10 GHz frequency range. The use of N = 9 matches up to 5 GHz. When N = 9 is chosen, it is suitable to use the general Padé approximation procedure to produce, for example, a 5-pole approximation of H 1 (s), relying on the known q k s.

Using ROM's by matching moments for TL circuits is a convenient approach to address TL effects, for example, [START_REF] Tanji | Closed-form expressions of distributed rlc interconnects for analysis of on-chip inductance effects[END_REF] presents a set of formulas for such a purpose. However, the limitation in [START_REF] Tanji | Closed-form expressions of distributed rlc interconnects for analysis of on-chip inductance effects[END_REF] is that low order moments need be approximated sometimes, which may reduce the accuracy of lower order models. The second limitation is that no formulations for higher order moments are available, therefore, no generalization can be made. This paper, however, gives closed-form moments' representations for the TL driver circuit by (17). The closed-form formulas for arbitrary order transfer function moments are very convenient in deriving ROM's or other analytical metrics for the TL driver circuit.

Moments and ROM of H dist (s)

H dist (s, x) is reformulated as H dist (s, x)= H E dist cosh γ x +H F dist sinh γ x γ x H C dist cosh γ d +H D dist sinh γ d γ d = ∑ ∞ m=0 p m s m ∑ ∞ k=0 q k s k (19)
where

H C dist = (1 + sτ D s )H A dist /(Z 0 Z L ) = H C 1 H D dist = (1 + sτ D s )H B dist /(Z 0 Z L ) = H D 1 H E dist = ∞ ∑ i=0 b E i s i = 1 + sτ D s H F dist = ∞ ∑ i=0 b F i s i = R s G x + sτ x s
Similarly, W x is used to denote the lumped parameter for the TL of length x, i.e. W x = W 0 x, for W ∈ {R, L, G,C}.

Comparing ( 16) and ( 19), the two q k s are the same. Following the similar way as deriving q k , the closed-form representation of p m , for G 0 = 0, is derived as

p m = min(1,m) ∑ i=0 m-i ∑ n= m-i 2 n m -i -n •C n x L m-i-n x R 2n+i-m x (2n + 1)b E i + b F i (2n + 1)! (20)
Similarly, using W x to replace W d in (18), for W ∈ {R, L, G,C}, and using b E i and b F i to replace a C i and a D i in (18), the representation for p m , when G 0 is not 0, can be obtained as well.

Fig. 4 shows |H dist ( f )| using the analytical formula (7b), and using the truncated series representation in (19). The same TL as in Fig. 2 is used. C D = 0, and the excitation point x is set at the half length point of the TL, i.e. x = d/2. For the purpose of demonstration, the numerator polynomial is truncated to N -1, and the denominator polynomial is truncated to N. The results show the correctness of moment function (20). Since p m and q k can be easily obtained from the derived formulas, and it is straightforward to obtain the moments for H dist (s, x) from p m and q k , the standard Padé approximation procedure can be employed to generate ROM's representing the function, H dist (s, x).

With the closed-form moment functions (17) and (20), either through direct truncation or using the general Padé procedure, rational functions can be formed, and the noise transfer function can be represented by several poles as below:

H(s) = W ∑ i=1 d i s -p i , ( 21
)
where p i is the i-th pole obtained from the rational form of noise transfer function, and d i is the residue. The step response, response to pulse source input, or other type of responses can be obtained conveniently by using (21).

Experimental Results

To demonstrate the proposed approach, here the experimental results are presented for a driver circuit, which drives a TL of 20mm length with PUL parameters: R 0 = 60 Ω, L 0 = 0.6 µH, and C 0 = 61 pF. At the other end of the TL, the load is modeled as a 300 fF capacitor. The modeled effective resistance R s is 200 Ω. The line has a separation of 0.35 mm with the ground. The EM wave illumination situation has a maximum transient electrical field strength of 2 kilovolts per meter. The bandwidth of the incident EM wave is around 10 GHz.

Fig. 5 illustrates the induced EM wave coupling noise in the circuit. The lower graph in Fig. 5(a) shows the transient incident EM field values around the half length point of the TL, which looks like a Gaussian shape waveform and has an amplitude of 2 KV/m and a spread of around 0.4 ns. However, the resulting distributed source (E i x,0,h (t) -E i x,0,0 ) • e x has a far less amplitude, which is about 100 times smaller than the amplitude of that location's incident field strength. (E i x,0,h (t) -E i x,0,0 (t)) • dx are used to connect adjacent two sections. In the figure, the noise waveform is only produced by the distributed sources, excluding V i f (t, 0) and V i f (t, d), since modeling the distributed sources is the major task which raises the difficulty. The noise waveform induced by distributed sources has an amplitude of around 0.2 V in the tested case. It causes around 20 ps delay variation. Since those pulses are as wide as 400 ps, logic faults can also occur if the incident EM wave becomes stronger. For example 8 KV/m field will produce a pulse around 0.4 ns duration and as strong as 0.9 V, which can cause the circuit operating in instable region and finally the logic state may be changed. Fig. 6 shows the results of the proposed modeling approach compared with Hspice simulation. Fig. 6(a) demonstrates the accuracy of Chebyshev approximation in the modeling of distributed sources. N = 3 is chosen as the approximation order, using the uniform plane wave metric and considering the bandwidth of the incident wave is around 10 GHz. In the diagram, the transient noise waveforms using Hspice simulation of twenty-section TL model and Chebyshev approximation with an order of 3 are very close to each other. It shows the effectiveness and accuracy of Chebyshev approximation approach. Fig. 6(b) compares the obtained noise waveforms in the first cycle, by combining Chebyshev expansion and Padé approximation based on the closed-form moment functions, and by using Hspice 20 TLs simulation. Note that in Fig. 6(a), exact transfer functions are used together with Chebyshev approximation. It has demonstrated the correctness of those derived formulas, with which ROM's of various degree of accuracies can be generated quickly.

Conclusion

This paper presents a novel modeling technique for considering EM wave coupling noise into conventional VLSI noise analysis flows. The proposed Chebyshev approximation approach for processing the distributed sources is shown as efficient and accurate. The paper gives the closedform representations of arbitrary order moments for the noise transfer functions. The derived closed-form representations can compute the moments exactly and quickly, thereby the proposed method will be useful in many design steps such as repeater insertion, timing models, and so on, to produce ROM's or other analytical metrics.

Figure 1 .

 1 Figure 1. Equivalent driver circuit under EM wave illumination.

  s v e k •(x,y,z) T = (cos β cos θ, sin β, cos β sin θ) e k = (cos α sin θsin α sin β cos θ, sin α cos β, cos α cos θsin α sin β sin θ)

  Rs=100Ω, C L = 5pF

Figure 2 .

 2 Figure 2. Comparisons between analytical results and chebyshev approximations.

Fig. 2

 2 Fig. 2 compares | V dist ( f ) E 0 | using Chebyshev approximations based on (10) and using formula (13). For the tested TL circuit (R 0 = 120 Ω/m, L 0 = 350 nH/m, C 0 = 120 pF/m, G 0 = 0, d = 15 mm, C D = 0 and α = β = θ = π 4) with driving and loading conditions varied, the Chebyshev approximation method converges very quickly as N becomes larger. When N = 2, the Chebyshev approximations of V dist (s) match the analytical results well up to 6 GHz, the frequency of incident EM wave. The use of N = 3 matches the frequency up to 14 GHz, and the use of N = 5 matches a frequency range of 20 GHz.In actual EM wave illumination situation, the incident EM field distribution is arbitrary. Uniform plane wave is employed at first to estimate the required order of Chebyshev approximation. By inserting the PUL parameters of the used interconnection techniques into (13) and[START_REF] Wang | The chebyshev expansion based passive model for distributed interconnect networks[END_REF], N is determined by judging the error between the two calculated results around the maximum EM wave frequency in the actual illumination environment. Next, the actual field values are used to calculate the time-domain noise waveform V dist (t) according to[START_REF] Wang | The chebyshev expansion based passive model for distributed interconnect networks[END_REF].

  Rs=10Ω, C L = 500fF

Figure 3 .

 3 Figure 3. Comparisons between analytical results and Padé approximations of H 1 (s).

  Rs=10Ω, C L = 500fF

Figure 4 .

 4 Figure 4. Comparisons between analytical results and Padé approximations of H dist (s).

Fig. 5 (

 5 b) shows the results from Hspice simulations only. To model the driver circuit into the Spice netlist, the TL is broken into 20 sections with equal lengths. Then lumped voltage sources Ê d/20 0

Figure 5 .

 5 Figure 5. Illustration of Field Stimulus and EM coupling noise.

Figure 6 .

 6 Figure 6. Verification of Modeling Accuracy.

Proceedings of the Design, Automation and Test in Europe Conference and Exhibition (DATE'05) 1530-1591/05 $ 20.00 IEEE


Laurent Dumas
email: dumas@math.uvsq.fr

Mohammed El Rhabi
email: elrhabi@gmail.com

Gilles Rochefort
email: grochefort@realeyes3d.com

An Evolutionary Approach For Blind Deconvolution Of Barcode Images With Nonuniform Illumination

This paper deals with a joint nonuniform illumination estimation and blind deconvolution for barcode signals by using evolutionary algorithms. Indeed, such optimization problems are highly non convex and a robust method is needed in case of noisy and/or blurred signals and nonuniform illumination. Here, we present the construction of a genetic algorithm combining discrete and continuous optimization which is successfully applied to decode real images with very strong noise and blur.

I. INTRODUCTION

This paper is devoted to the problem of blind deconvolution of linear barcode signals. Linear barcodes are one of the oldest technology related to automatic identification and data capture. Typical barcodes symbologies are UPC/EAN/JAN and are one variation of over 250. Linear barcodes are graphical representations of a finite sequences of digits. Information is encoded in the widths (lines) and the spacings of parallel lines. Up to a choosen symbology, each sequence of digits admits a unique graphical representation. So, retrieving digits from an acquired barcode require accurate estimation of the true lenght bars and spaces pattern. In order to help decoding, some symbology may include special start/stop group of bars and a validation checksum digit.

Ideally, the barcode is given as a binary one-dimensional signal when scanned by a laser or a binary image when scanned by a camera (1 represents white space and 0 indicates a black bar).

Linear barcodes are optimized to be read by a laser scanner. Typically UPC, EAN, JAN symbologies have been widely adopted in industrial applications (supermarket checkout system, labelling railroad cars, ...). Nowadays, every article has its own printed barcode. One may imagine using barcodes for its own purpose, for example price comparison, allergen database (gluten, peanut ...). Unfortunately, it is not usual to have its own personal barcode readers.

So, computer vision based barcode recognition could be an interesting alternative. Here, the problem could be separated in two steps: localization and decoding. Computer approach increases the performance of both steps, barcode localization and decoding (see [START_REF] Matas | AObject recognition using a tag[END_REF], [START_REF] Muniz | A robust software barcode reader using the Hough transform[END_REF]). In the computer vision approaches, low cost cameras (webcam, cameraphone, ...), without autofocusing or macro mode, could introduce some distortions and noise artifacts. As a consequence, decoding performances depend on the distance from the camera to the surface where the barcode appears ( [START_REF] Esedoglu | Blind deconvolution of barcode signals[END_REF]).

Anyway, acquired images are often blurred and distorted by various factors, including speckle noise, ambient light. Then, a way, to decode barcodes included in an image, is first to restore it. We would like to retrieve a barcode image with blur, noise and also distorted by nonuniform illumination, that has been acquired by a low cost camera (see an example on Figure 1). To reduce this problem, in [START_REF] Sanner | Ambient illumination bar code reader[END_REF], the author proposes some additional hardware device to overcome this drawback, but not convenient in the context of a mobile phone with a digital camera [START_REF] Shams | Bar code recognition in highly distorted and low resolution images[END_REF]. To consider the effect of nonuniform illumination, several methods have been proposed, among them the illumination cone method, a spherical harmonic based representation [START_REF] Zhang | Face recognition under variable lighting using harmonic image exemplar[END_REF], quotient image based approaches [START_REF] Wang | Generalized Quotient Image[END_REF] and correlation based method [START_REF] Savvides | Corefaces -Robust Shift Invariant PCA based Correlation Filter for illumination Tolerant Face Recognition[END_REF]. The performance of these methods are poor and many of them require either knowledge of light source or a large number of training data which is not tractable in real world applications.

Generally speaking, an image u could be seen as the product of a reflectance R and the illumination effect I (see [START_REF] Horn | Robot Vision[END_REF]). This model has been used in [START_REF] Kim | Joint nonuniform illumination estimation and deblurring for bar code signals Optic Express[END_REF] to deblur barcode signals under nonuniform illumination. The blurring operator is a Gaussian function with unknown standard deviation. The solution is estimated by a penalized non linear least square objective function, based on a proper parametrization of a linear barcode and nonuniform illumination. The solution is then estimated by the well-known gradient based method.

In this paper, we propose a novel method for blind deconvolution of barcode signals in the out-of-focus case, based on a genetic algorithm taking into account the highly non convexity of the problem and the fact that the target solution is a binary image.

The paper is organized as follows : section II gives a mathematical definition of the barcode and presents the problem. Section III proposes a solution to the problem by means of a genetic algorithm, describes the method and shows how to implement it. Finally, some numerical results illustrate this work in Section IV.

II. THE MATHEMATICAL MODEL

A. A mathematical definition for a barcode

Here, we give a general definition of a barcode in the bidimensional case. We suppose that f (x, y) is a 2D continuous function representing the intensity value of the barcode at location (x, y) :

f (x, y) = 1, if (x, y) ∈ B k 0, else (1)
where B k represents the region bounded by the k th bar. Assume that f d (i, j) is a discrete function acquired from f (x, y) such that :

f d (i, j) = j+1 j i+1 i f (x, y) dxdy (2)
where (i, j) = (0, 0) , (0, 1) , ..., (1, 0) , ..., (N, M ). f d (i, j) is the area of the region in pixel (i, j) bounded by the bar passing through that pixel, each pixel is considered as a square region.

B. The blurring model

Let u ∈ IR N M be an original N × M gray-scale image, and K ∈ IR N M ×N M represents a blurring (or convolution) operator, I ∈ IR N M the non uniform illumination, n ∈ IR N M the additive noise, and u 0 ∈ IR N M an observation which satisfies the relationship:

u 0 = IKu + n,
In general, the blurring operator K is supposed shift-invariant, and takes the following form:

K = k u, ( 3
)
where k is the point spreading function (PSF). In many real cases, the blurring operator K is unknown. Usually, we assume that K could be represented by few parameters p , these parameters are a priori unknown. For example, in the outof-focus case, the PSF could be approximated as :

k(x, y) = 1 πr 2 1 1 Br (x, y), (4)
where B r is the ball with the center 0 and radius r > 0 and 1 1 Br (x, y) defines its indicator function.

The non uniform illumination, I, is also unknown. Generally, it is assumed smooth and could be modelled by using B-spline functions that are uniformly located in the spatial domain ( [START_REF] Kim | Joint nonuniform illumination estimation and deblurring for bar code signals Optic Express[END_REF]):

I(x, y) = Lx i=1 Ly j=1 I ij β n (x -l i ) β n (y -l j ) , (5)
where β n (t) is a n-order B-spline function centered at zero, L x (resp. L y ) is the number of the B-spline functions, I ij the B-spline coefficients and l i , i ∈ {1, ..., L x } (resp. l j , j ∈ {1, ..., L y }) denote uniformly distributed center locations of the B-spline functions.

C. The inverse problem

If we suppose that n is a gaussian noise, the joint estimation of the blur kernel K, the illumination I and the restored image u is an ill-posed problem. An approximation of K (or r equivalently), I and u could be obtained by solving the following non convex optimization problem:

min r,I,u E(r, I, u) = Ω |IK(r)u -u 0 | 2 dxdy (6)
In order to control the noise, a regularization term of the type λ Ω ϕ(|∇u|)dxdy must be added to the cost function E. Yet, the choice of the function ϕ must be relevant and is often a difficult task. Moreover, the use of a deterministic descent method, such as steepest descent, is only able to seek for a local minima of the cost function, which is unfortunately highly non linear.

III. SOLVING THE INVERSE PROBLEM

As seen above, the solution of the inverse problem depends on the regularization term. The method of genetic algorithms is used here to overcome this drawback and to solve the inverse problem on a robust way with respect to noise and blur.

A. The evolutionary algorithm

Genetic algorithms are global optimization methods directly inspired from the Darwinian theory of evolution of species ( [START_REF] Goldberg | Genetic Algorithms in Search, Optimization, and Machine Learning[END_REF]). They consist in following the evolution of a certain number N p of possible solutions, also called population. To each element (or individual) x i ∈ O of the population is affected a fitness value inversely proportional to J(x i ), in case of a minimization problem for the cost function J. The population is regenerated N g times by using three stochastic principles called selection, crossover and mutation, that mimics the biological law of the 'survival of the fittest'.

They have show their efficiency in many applicative fields, in engineering science or in medicine, cite for instance, car shape optimization [START_REF] Dumas | Hybrid method for aerodynamic shape optimization in automotive industry[END_REF], turbine shape optimization [START_REF] Druez | A fully Adaptive Hybrid Optimization of aircraft engine blades[END_REF], pacemaker optimization [START_REF] Dumas | How Genetic Algorithms can improve a pacemaker efficiency[END_REF], to mention some results obtained by one of the author.

The genetic algorithm that is used here acts in the following way: at each generation, Np 2 couples are selected by using a roulette wheel process with respective parts based on the fitness rank of each individual in the population. To each selected couple, the Darwinian principles, namely crossover and mutation, are then successively applied with a respective probability p c and p m . A one-elitism principle is added in order to be sure to keep in the population the best element of the previous generation.

B. The cost function

The cost function that has been used is different from the one written in (6) because the regularization term has been removed and the L 2 norm is replaced by the L 4 norm:

min r,I,u J(r, I, u) = Ω |IK(r)u -u 0 | 4 dxdy. (7)
Here, The choice of the L 4 norm is not crucial but it has shown better convergence properties.

C. The search space

One key feature of genetic algorithms is their ability to deal, either with discrete or continuous search spaces, or even with both of them. It is precisely the case here as the search space for u is of finite type whereas the search space for I and r is of continuous type:

(i) The search space for u In this work, a particular symbology of barcode has been studied, namely the EAN13 type, which is currently the most used one. It encodes 13 characters by using the width of 30 vertical bars. It is worth noticing that the chosen approach can easily adapt to any other type of barcode. As shown in Figure 2, the signal u to be found is discretized in 61 integer values, gathered in a vector U :

U = (U 1 , ..., U 61 ) ∈ IN 61
which corresponds to the interval lengths (in terms of pixels) of each successive sequence of white and black colors (where u is respectively chosen equal to 0 and 1). Note in particular that U 1 and U 61 play a particular role bethey correspond to the left, respectively right, unknown margin at each side of the barcode in the image. Note also that the constraint

61 i=1 U i = M ( 8
)
where M is the width of the image has to be fulfilled.

(ii) The search space for r and I It is a given interval of IR + for r whereas it is an hypercube included in [0, 1] Lx , for I. The latter corresponds to the ordinates of L x points used for the reconstruction of I by cubic splines interpolation (see Figure 3 for an example of a nonuniform illumination obtained from an interpolation with 5 points).

Fig. 3. The continous search space for I: example of a nonuniform illumination

D. The darwinian principles

The darwinian principles, namely crossover and mutation, play a crucial role in the convergence properties of genetic algorithms and have to be chosen with care. (i) The darwinian principles for the discrete variable U The crossover principle between two vectors U and V in ∈ IN 61 is inspired from the original crossover in the case of chromosomes: it consists in creating two offsprings by exchanging a given number, m, of randomly chosen consecutive sequences of each vector (m-point crossover). As the obtained offsprings do not generally satisfy the constraint (8), a renormalization principle is added which consists in increasing (or reducing) the value of randomly chosen components from them. Actually, this last process acts as the mutation operator.

Assume for instance that m = 1 and denoting by U and V the two parents:

U = (U 1 , ..., U 61 ) V = (V 1 , ..., V 61 )
then the two offsprings are respectively:

U = (U 1 , ..., U k , V k+1 ..., V 61 ) V = (V 1 , ..., V k , U k+1 ..., U 61 )
for a random k ∈ {1, ..., 61}. Then two offsprings are mutated with the iterated process

U i = U i ± 1 for a random i ∈ {1, ..., 61} until 61 i=1 U i = M (ii)
The darwinian principles for the continuous variables For the continuous variables I and r, it consists in a barycentric combination for crossover and a non uniform mutation. For instance for the blur radius, considering two parents r 1 and r 2 , the crossover and the mutation respectively write as:

(r 1 , r 2 ) = (αr 1 + (1 -α)r 2 , αr 2 + (1 -α)r 1 )
and

r i = r i + β(r max -r i ) b (or r i = r i -β(r i -r min ) b with probability 0.5)
where α and β are random numbers between 0 and 1 and b ∈ [START_REF] Matas | AObject recognition using a tag[END_REF][START_REF] Shams | Bar code recognition in highly distorted and low resolution images[END_REF].

IV. RESULTS

Four examples of results are presented here. The first two examples are validation test cases: they consist in determining in the absence of noise, the continuous (respectively discrete variables) knowing the other ones. The third example is a complete deconvolution case of an image without knowing any information (illumination, blur, position of the barcode in the image, etc...) and with a strong noise. Finally the fourth example is an example of a barcode decoding from a real image, namely the one depicted in Figure 1.

Example 1: determination of blur radius and nonuniform illumination of an image with a known barcode, no noise

In this first example, the only unknowns are the blur radius and the nonuniform illumination function of an image of a known barcode. An evolutionary algorithm with 200 individuals and 150 generations is always able to approximate very accurately these two quantities. Actually, a quasi perfect agreement is achieved for the illumination function, as well as for the blur radius, in any configuration showing the robustness of the approach. For example, a comparison is made on Figure 4 between the original illumination (full line), the initial best illumination (dotted line) and the final best one after the optimization process (mixed dotted line, almost superimposed to the original illumination). Example 2: determination of a barcode characteristics with a known blur, no illumination, no noise

The objective is here to show that the evolutionary algorithm acting only on the discrete variables, that is the vector U ∈ IN 61 , is able to find an exact barcode corresponding to a given vector U 0 ∈ IN 61 . In this case, the blur radius r = 3 is fixed and the unknowns are the 61 coordinates of U . With a population number of 200 individuals and a generation number around 300, the algorithm always achieves a perfect convergence, even starting with a very bad initial guess (see Figure 5 for an example). As it can be observed, it is an example of a signal, with a strong noise and blur and with a nonuniform illumination. It is representative of the robustness and accuracy of our algorithm for deblurring and denoising. Obviously, the same kind of results can be reproduced with very similar signal.

The evolutionary algorithm presented in the previous section is applied with the following parameters: The convergence history of the optimization process is also depicted in Figure 8 . It shows that the evolutionary algorithm almost performs to achieve a 'perfect' solution, without any a priori information, such as the position of the barcode in the image and starting with a very far initial guess.

A statistical study on a large number of barcodes has also been done, showing that this approach can give positive results, that is a decodable barcode, on a large range of similar signals, even with a very strong noise, a strong blur and with nonuniform illumination.

Example 4: decoding of a real barcode image with blur, noise and nonuniform illumination:

The last example presents the decoding of a real barcode image of poor quality, namely the one depicted in Figure 1 (also visible in Figure 10) and taken with a Nokia N70. The code on the image is undistinguishable because of strong blur. Moreover, the image is also distorted by a strong nonuniform illumination. By using the same parameters for the genetic algorithm than for the previous example, the obtained best signal is depicted in Figure 9 (red, dotted line) and compared with the observed signal (blue, full line). The matching is almost excellent between both signals. Moreover, the obtained signal, also depicted in Figure 10, has been successfully decoded and is associated with the EAN13 code 4025515825135.

V. CONCLUSION

A robust method for blind deconvolution of barcode signals in the presence of blur, noise and with nonuniform illumination is presented here. Based on the resolution of the associated inverse problem with an evolutionary algorithm on a mixed search space, it allows to decode, without any information and in a robust and reproducible way, a very noisy and blurred barcode image. In particular, this approach has been successfully applied to decode barcodes images taken from a camera of very poor quality.

Fig. 1 .

 1 Fig. 1. A real barcode image with blur, noise and non uniform illumination

Fig. 2 .

 2 Fig. 2. The discrete search space for u: example of a barcode representation

Fig. 4 .

 4 Fig. 4. The exact illumination, the final best illumination (almost superimposed) and the initial best illumination (dotted line)

Fig. 5 .

 5 Fig. 5. The exact signal, the final best signal (almost superimposed) and the initial best signal (dotted line)

Fig. 6 .

 6 Fig. 6. The observed signal u 0

(Fig. 7 .

 7 Fig. 7. Comparison of the observed barcode, the initial and final best barcode.

Fig. 8 .

 8 Fig. 8. The convergence history

Fig. 9 .

 9 Fig. 9. The observed signal (full line) and the final best signal (dotted line)

  Fig. 10. Comparison of the observed barcode image and the final best barcode