Sylvain Lazard
email: sylvain.lazard@loria.fr

Mariano Luis
email: luis@fceia.unr.edu.ar

Pe

Sylvain Petitjean
email: sylvain.petitjean@loria.fr

Intersecting Quadrics: An Efficient and Exact Implementation

Keywords: I.1.2 [Symbolic and Algebraic Manipulation]: Algorithms; I.3.5 [Computer Graphics]: Computational Geometry and Object Modeling Algorithms, Experimentation Intersection of quadrics, C++ implementation, experimental analysis

We present the first complete, exact and efficient C++ implementation of a method for parameterizing the intersection of two implicit quadrics with integer coefficients of arbitrary size. It is based on the near-optimal algorithm recently introduced by Dupont et al. [2].

Unlike existing implementations, it correctly identifies and parameterizes all the connected components of the intersection in all cases, returning parameterizations with rational functions whenever such parameterizations exist. In addition, the coefficient fields of the parameterizations are either minimal or involve one possibly unneeded square root.

We prove upper bounds on the size of the coefficients of the output parameterization and compare these bounds to observed values. We give other experimental results and present some examples.

INTRODUCTION

Computing an explicit representation of the intersection of two general quadrics (i.e., quadratic surfaces) is a fundamental problem in areas such as solid modeling, computational * Work done while this author was visiting LORIA (supported by the International Relations Delegation of INRIA).

Past work. Until recently, the only known general method for computing a parametric representation of the intersection between two arbitrary quadrics was that of J. Levin [7]. This method is based on an analysis of the pencil generated by the two quadrics, i.e. the set of linear combinations of the two quadrics.

Though useful for curve tracing, Levin's method has serious limitations. When the intersection is singular or reducible, a parameterization by rational functions is known to exist, but Levin's pencil method fails to find it and generates a parameterization that involves the square root of some polynomial. In addition, since it introduces algebraic numbers of very high degree (for instance in the computation of eigenvalues and eigenvectors), a correct implementation using exact arithmetic is essentially out of reach. In addition, when a floating point representation of numbers is used, the method may output results that are wrong (geometrically and topologically) and it may even fail to produce any parameterization at all and crash.

Over the years, Levin's seminal work has been extended and refined in several different directions. Wilf and Manor [START_REF] Wilf | Quadric-surface intersection curves: shape and structure[END_REF] use a classification of quadric intersections by the Segre characteristic (see [START_REF] Bromwich | Quadratic Forms and Their Classification by Means of Invariant Factors[END_REF]) to drive the parameterization of the intersection by the pencil method. Recently, Wang, Goldman and Tu [START_REF] Wang | Enhancing Levin's method for computing quadric-surface intersections[END_REF] further improved the method making it capable of computing structural information on the intersection and its various connected components and able to produce a parameterization by rational functions when such a parameterization exists. Whether the refined algorithm is numerically robust is open to question.

Another method of algebraic flavor was introduced by Farouki, Neff and O'Connor [START_REF] Farouki | Automatic parsing of degenerate quadric-surface intersections[END_REF] for parameterizing the intersection in degenerate situations. In such cases, using a combination of classical concepts (Segre characteristic) and algebraic tools (factorization of multivariate polynomials), the authors show that explicit information on the morphological type of the intersection curve can be reliably obtained. A notable feature of this method is that it can output an exact parameterization of the intersection in simple cases, when the input quadrics have rational coefficients. No im-plementation is reported however.

Rather than restricting the type of the intersection, others have sought to restrict the type of the input quadrics, taking advantage of the fact that geometric insights can then help compute the intersection curve [START_REF] Miller | Geometric algorithms for detecting and calculating all conic sections in the intersection of any two natural quadric surfaces[END_REF][START_REF] Shene | On the lower degree intersections of two natural quadrics[END_REF]. Specialized routines are devised to compute the intersection curve in each particular case. Even though such geometric approaches are numerically more stable than the algebraic ones, they are essentially limited to the class of so-called natural quadrics, i.e., the planes, right cones, circular cylinders and spheres.

Apart from [START_REF] Dupont | Near-optimal parameterization of the intersection of quadrics[END_REF], perhaps the most interesting of the known algorithms for computing an explicit representation of the intersection of two arbitrary quadrics is the method of Wang, Joe and Goldman [START_REF] Wang | Computing quadric surface intersections based on an analysis of plane cubic curves[END_REF]. This algebraic method is based on a birational mapping between the intersection curve and a plane cubic curve. The cubic curve is obtained by projection from a point lying on the intersection. Then the classification and parameterization of the intersection are obtained by invoking classical results on plane cubics. The authors claim that their algorithm is the first to produce a complete topological classification of the intersection (singularities, number and types of connected components, . . . ). Numerical robustness issues have however not been studied and the intersection may not be correctly classified. Also, the center of projection is currently computed using Levin's (enhanced) method: with floating point arithmetic, the center of projection will in general not exactly lie on the curve, which is another source of numerical instability.

Contributions.

In this paper, we present the first exact, robust, efficient and usable implementation of an algorithm for parameterizing the intersection of two arbitrary quadrics, given in implicit form, with integer coefficients. (Note that quadrics with rational or finite floating-point coefficients can be trivially converted to integer form.) This implementation is based on the parameterization method described in [START_REF] Dupont | Near-optimal parameterization of the intersection of quadrics[END_REF].

Precisely, our implementation has the following features: • it computes an exact parameterization of the intersection of two quadrics with integer coefficients of arbitrary size; • it correctly identifies, separates and parameterizes all the connected components of the intersection and gives all the relevant topological information; • it places no restriction of any kind on the type of the intersection or the type of the input quadrics; • the parameterization is rational when one exists; otherwise the intersection is a smooth quartic and the parameterization involves the square root of a polynomial; • the parameterization is either optimal in the degree of the extension of Z on which its coefficients are defined or, in a small number of well-identified cases, involves one extra possibly unnecessary square root; • the implementation is carefully designed so that the size of the coefficients is kept small; • it is fast and efficient and can routinely compute parameterizations of the intersection of quadrics with input coefficients having ten digits in less than 50 milliseconds on a mainstream PC.

Our C++ implementation can be queried via a web interface at http://www.loria.fr/isa/qi. The code will be released to educational and research audiences shortly.

The paper is organized as follows. In Section 2, we recall the main ideas of the parameterization algorithm we introduced in [START_REF] Dupont | Near-optimal parameterization of the intersection of quadrics[END_REF] and describe its implementation. In Section 3, we focus on two cases covering the two main parameterization philosophies used in our implementation, prove theoretical bounds on the size of the output coefficients and compare those bounds to observed values. We give more experimental results and performance evaluation in Section 4. Finally, we show the output produced by our implementation for some examples in Section 5, before concluding.

ALGORITHM AND IMPLEMENTATION

In this section, we give a brief presentation of the basic ideas underpinning our near-optimal parameterization method [START_REF] Dupont | Near-optimal parameterization of the intersection of quadrics[END_REF]. We then move on to a description of the main design choices we made to implement it.

Preliminaries

In what follows, all the matrices considered are 4 × 4 real matrices, unless otherwise specified. We call quadric associated with a symmetric matrix S the set

QS = {x ∈ P 3 | x T Sx = 0},
where P 3 = P 3 (R) denotes the real projective space of dimension 3 (x T Sx is quadratic and homogeneous in the coordinates of x). In the rest of this paper, points and parameterizations are assumed to live in projective space. Recall that a point of P 3 has four coordinates.

We define the inertia of S and QS as the pair

σS = (max (σ + , σ -), min (σ + , σ -)),
where σ + (resp. σ -) is the number of positive (resp. negative) eigenvalues of S. The rank of S is the sum σ + + σ -. Recall that Sylvester's Inertia Law asserts that the inertia of S (and thus the rank) is invariant by a real projective transformation [START_REF] Lam | The Algebraic Theory of Quadratic Forms[END_REF]. We call projective cones (or simply cones) the quadrics of rank 3. For the benefit of the reader, we recall that, in affine real space, quadrics of inertia (4, 0) are empty, quadrics of inertia (3, 1) are ellipsoids, hyperboloids of two sheets or elliptic paraboloids, and quadrics of inertia (2, 2) are hyperboloids of one sheet or hyperbolic paraboloids (see [START_REF] Dupont | Near-optimal parameterization of the intersection of quadrics[END_REF] for a complete characterization of affine quadrics). Also, quadrics of inertia (2, 1) are cones or cylinders. All the quadric surfaces except those of inertia (3, 1) are ruled surfaces.

Given two matrices S and T , let R(λ, µ) = λS + µT . The set {R(λ, µ) | (λ, µ) ∈ P 1 } is called the pencil of matrices generated by S and T . For the sake of simplicity, we sometimes write a member of the pencil R(λ) = λS -T, λ ∈ R ∪ {∞}. Associated to a pencil of matrices is a pencil of quadrics {Q R(λ,µ) | (λ, µ) ∈ P 1 }. Recall the classical result that the intersection of two distinct quadrics of a pencil is independent of the choice of the two quadrics.

The binary quartic form det R(λ, µ) is called the determinantal equation of the pencil. The quadrics of the pencil of rank less than or equal to 3 are exactly the quadrics Q R(λ,µ) such that det R(λ, µ) = 0.

Near-optimal parameterization algorithm

Main ideas

Let {Q R(λ,µ) | (λ, µ) ∈ P1 }, with R(λ, µ) = λS + µT , be a pencil of quadrics. The main idea of existing methods for parameterizing the intersection of two quadrics based on an analysis of their pencil (Levin's and derivatives) is as follows: find a quadric QR of some particularly simple form in the pencil generated by QS and QT (assume QR = QS), parameterize this quadric, plug the parameterization X in the equation of QS, solve the resulting equation X T SX = 0, and plug the result in X, finally giving the parameterization of the intersection.

The key to making this procedure work in practice is to find a quadric QR that is ruled and thus admits a parameterization that is linear in one of its parameters so that the equation X T SX = 0 has degree 2. Levin's main result was to prove that a pencil of quadrics always contains at least one "simple" ruled quadric [START_REF] Levin | A parametric algorithm for drawing pictures of solid objects composed of quadric surfaces[END_REF]. However, since such quadrics are found by first finding the zeros of the determinant of the upper left 3 × 3 submatrix of R(λ, µ), a cubic equation, and since cubic equations have generically no rational root (by Hilbert's Irreducibility Theorem), Levin's algorithm introduces non-rational numbers at an early stage and, in practice, floating-point arithmetic has to be used, resulting in numerical robustness problems.

The principal contribution of [START_REF] Dupont | Near-optimal parameterization of the intersection of quadrics[END_REF] was to show that, by a careful choice of the intermediate quadric QR, the appearance of algebraic numbers can be kept to a minimum. One major result is encapsulated in Theorem 3 of [START_REF] Dupont | Near-optimal parameterization of the intersection of quadrics[END_REF]: except when the intersection is reduced to two real points, the pencil contains at least one ruled quadric whose coefficients are rational and such a quadric can be easily computed. In addition, thanks to new worst-case optimal (in the number of square roots) parameterizations of ruled projective quadrics, we can always find such a rational ruled quadric QR with a parameterization involving only one square root. Some of the basic ingredients used in our algorithm or to infer information about the intersection are the Segre classification of pencils and its refinement over the reals (the Canonical Form Theorem for pairs of real symmetric matrices -see [START_REF] Uhlig | Simultaneous block diagonalization of two real symmetric matrices[END_REF]), a projective setting, ad hoc projective transformations to compute the canonical form of a projective quadric, and Sylvester's Inertia Law [START_REF] Lam | The Algebraic Theory of Quadratic Forms[END_REF].

The basic principles underlying the design of our implementation are as follows:

• compute structural information on the intersection and its various real components as early as possible; • use the structural information gathered to drive the parameterization process and make the right choices so that the output is optimal or near-optimal from the point of view of the degree of the extension of Z on which its coefficients are defined.

In our implementation we were interested not just in optimizing the number of square roots in the output but also in minimizing the size of the output coefficients. For this reason, the basic philosophy is to use as intermediate ruled quadric QR a quadric with rational coefficients of the smallest rank that we can easily find, the rationale being, for instance, that the parameterization of a cone involves coefficients of smaller asymptotic size than the coefficients of the parameterization of a quadric of inertia (2, 2). There are essentially two cases: (i) QR has rank 4; (ii) QR has rank 3 or less.

QR has rank 4

The main case where QR has rank 4 is when the intersection is a smooth quartic (we leave the other cases aside due to lack of space). In this situation, the quartic determinantal equation det R(λ) has no multiple root. It could well be that at least one of its simple roots is rational and that a QR with rank less than 4 could have been used, but checking this via the Rational Root Theorem can be very time consuming 1 . Since generically a degree-four equation has no rational root, we prefer instead to isolate the real zeros of the determinantal equation using an implementation of Uspensky's algorithm [START_REF] Rouillier | Efficient isolation of a polynomial real roots[END_REF]. We then take (at most two) rational test points λi outside the isolating intervals in the areas where det R(λ) > 0. If one of the quadrics R(λi) has inertia (4, 0), the intersection is empty (it is a complex smooth quartic), a consequence of Finsler's Theorem (see [START_REF] Dupont | Near-optimal parameterization of the intersection of quadrics[END_REF]). Otherwise, we proceed.

We now have a quadric R0 = R(λ0) of inertia (2, 2) and a range of values I = [a, b] such that λ0 ∈ I and det R(λ) > 0 for all λ ∈ I. In the worst case, the parameterization of QR involves two square roots [START_REF] Dupont | Near-optimal parameterization of the intersection of quadrics[END_REF]. We can improve this situation as follows. First, compute a point p0 on QR 0 . Approximate this point by a point p with integer coordinates (recall that p is a projective point). Find the quadric QR = Q R(λ 1 ) through p. If p is close enough to p0, then λ1 ∈ I and det R > 0. We thus have a quadric of inertia (2, 2) containing a point in P 3 (Z): such a quadric can be parameterized with at most one square root [START_REF] Dupont | Near-optimal parameterization of the intersection of quadrics[END_REF].

Plugging the parameterization X((u, v), (s, t)) of QR, with (u, v), (s, t) ∈ P 1 , in the equation of any other quadric of the pencil gives a bihomogeneous equation that has degree two in (u, v) and two in (s, t). Solving this equation for (s, t) in terms of (u, v) and replugging in the parameterization of QR gives a parameterization of the smooth quartic:

X(u, v) = X1(u, v) ± X2(u, v) p ∆(u, v),
where X1(u, v) (resp. X2(u, v)) is a vector of homogeneous polynomials of degree 3 (resp. 1) and ∆(u, v) is a homogeneous polynomial of degree 4. If δ = det R is a square, then all of these polynomials have rational coefficients and the parameterization is optimal in terms of the degree of the extension of Z on which it is defined. If δ is not a square, then we can only conclude that the parameterization is near-optimal : it might well be that there exists another quadric Q R of inertia [START_REF] Dupont | Near-optimal parameterization of the intersection of quadrics[END_REF][START_REF] Dupont | Near-optimal parameterization of the intersection of quadrics[END_REF] in the pencil such that det R is a square, implying that √ δ could have been avoided in the output (see Section 5.2 for an example). Finding such a quadric however implies, in general, finding a rational point on a hyperelliptic curve (see [START_REF] Dupont | Near-optimal parameterization of the intersection of quadrics[END_REF]), a problem known to be very hard.

QR has rank < 4

Though not generic, the situation where QR has rank < 4 happens quite often in practice since it covers in particular all the types of intersection corresponding, in the Segre characterization, to the determinantal equation having a single multiple root λ0. Indeed, in that case, the multiple root is both real (otherwise its complex conjugate would also be a multiple root of det R(λ)) and rational (otherwise its algebraic conjugate would also be a multiple root of det R(λ)). So the associated quadric QR = Q R(λ 0 ) has rational coefficients and has rank 3 or less.

The general philosophy for parameterizing the intersection is to parameterize QR, plug the parameterization in any other quadric of the pencil, and solve the resulting equation in the parameters. There are however many situations in which this procedure can be simplified by the fact that we can find a rational point on QR outside its singular locus and thus parameterize QR rationally, and that we know enough information on the intersection to greatly simplify the solving and factorization of the equation in the parameters.

Let us illustrate this on the example of an intersection consisting of a cubic and a line that are tangent. The determinantal equation in this case has a quadruple root corresponding to a cone QR of inertia (2, 1). By the above argument, QR has rational coefficients. So the vertex c of QR has rational coordinates. c is the point of tangency of the cubic and the line of the intersection. Assume QR = QS. The line of the intersection is necessarily rational (otherwise its conjugate would be in the intersection). This line can be found by intersecting the cone QR with the plane tangent to QS at c. Picking any point with rational coordinates on this line other than c gives a non-singular rational point on the cone. A projective cone having a rational point other than its singular locus can be rationally parameterized. Plugging this parameterization in QS gives an equation in the parameters of the cone which factors into two terms of total degree 1 and 3. Each factor can then be solved rationally for one parameter in terms of the other. The linear factor yields the line of the intersection and the cubic factor yields the cubic.

Implementation

Our implementation builds upon the LiDIA [START_REF]A C++ Library for Computational Number Theory[END_REF] and GMP [4] C/C++ libraries. LiDIA was originally developed for computational number theory purposes, but includes many types of simple parameterized and template classes that are useful for our application. Apart from simple linear algebra routines and algebraic operations on univariate polynomials, we use LiDIA's number theory package and its ability to manipulate vectors of polynomials, polynomials having other polynomials as coefficients, etc. On top of it, we have added our own data structures. We have compiled LiDIA so that it uses GMP multiprecision integer arithmetic. From now on, we refer to the multiprecision integers as bigints, following the terminology of LiDIA.

Our implementation consists of more than 17,000 lines of source code, which is essentially divided into the following chapters:

• data structures (1,500 lines): structures for intersections of quadrics, for components of the intersection, for homogeneous polynomials with bigint coefficients (coordinates of components), for homogeneous polynomials with bigint polynomials as coefficients, and basic operations on these structures. . . • elementary operations (2,000 lines): computing the inertia of a quadric of bigints, the coefficients of the determinantal equation, the gcd of the derivatives of the determinantal equation, the adjoint of a matrix, the singular space of a quadric, the intersection between two linear spaces, applying Descartes's Sign Rule, the Gauss decomposition of a quadratic form into a sum of squares, isolating the roots of a univariate polynomial using Uspensky's method, . . . • number theory and optimizations (1,500 lines): gcd optimizations of the bigint coefficients of a polynomial, a vector or a matrix, optimizations of the coefficients of pairs and triples of vectors, reparameterization of lines so that its representative points have small height, . . . • quadric parameterizations (2,000 lines): parameterization of a quadric of inertia (2, 2) with bigint coefficients going through a rational point, of a cone (resp. conic), of a cone (resp. conic) with a rational point, of a pair of planes, . . . • intersection parameterizations (9,000 lines): dedicated procedures for parameterizing the components of the intersection in all possible cases, i.e. when the determinantal equation has no multiple root (1,500 lines), one multiple root (3,000 lines), two multiple roots (1,500 lines) or when it vanishes identically (3,000 lines). • printing and debugging (1,000 lines): turning on debugging information with the DEBUG preprocessor directive, checking whether the computed parameterizations are correct, pretty printing the parameterizations, . . . .

HEIGHT OF OUTPUT COEFFICIENTS

In this section, we prove theoretical bounds on the height of the coefficients of the parameterizations computed by our intersection software. We do this for the two types of intersection already outlined in the previous section: smooth quartic and cubic and tangent line.

The smooth quartic case is important because it is the generic intersection situation (given two random quadrics, the intersection is a smooth quartic with probability 1) and because it is also the worst case from the point of view of the height of the coefficients involved. The cubic and tangent line case is used to validate a key design choice we made, which is to take the quadric with rational coefficients of lowest possible rank to parameterize the intersection. We compare what happens for this type of intersection when we use for QR a quadric of inertia (2, 2) or a quadric of inertia [START_REF] Dupont | Near-optimal parameterization of the intersection of quadrics[END_REF][START_REF] Bromwich | Quadratic Forms and Their Classification by Means of Invariant Factors[END_REF].

In what follows, the size of an integer e is log 10 |e| + 1 if e is not zero and 1 otherwise. The size of an algebraic number e1 + √ δ e2, where e1, e2, δ are integers, is the maximum of the sizes of e1, e2 and δ. The size of a matrix, a polynomial or a parameterization is the maximum of the sizes of its coefficients. The height of a quantity E (a point, polynomial, matrix or parameterization) in terms of the coefficients of E is the ratio of the size of E and the size of E . When E depends on the coefficients of the input matrices S and T , we sometimes drop the reference to S and T and simply talk about the height of E.

In the following subsections, we give asymptotic upper bounds on the heights of parameterizations. For the sake of simplicity, we refer to these bounds as the heights of the parameterizations.

Smooth quartic

We consider here the case where the intersection is a smooth quartic. Let QR be the quadric of inertia (2, 2) used to parameterize the intersection and p a point of P 3 (Z) on QR, as described in Section 2.2.2. We assume that the size of p is constant. This assumption is very realistic in practice; experimentally, we have observed that the coordinates of p are integers between -2 and 2 most of the time. We prove the following result.

Proposition 3.1. The parameterization of a smooth quartic

X(u, v) = X1(u, v) ± X2(u, v) p ∆(u, v)
is such that the height in terms of S and T is asymptotically at most 27 for the polynomials of X1, 8 for the polynomials of X2 and 38 for ∆(u, v).

Proof. We first show how the parameterization of QR is computed and bound the height of its coefficients.

Let M0 be the projective transformation sending the point p0 = (1, 0, 0, 0) T to the point p. Let R0 denote the quadric obtained from R through the projective transformation M0: R0 = M T 0 R M0. It follows from Sylvester's Inertia Law [START_REF] Lam | The Algebraic Theory of Quadratic Forms[END_REF] that R0 has the same inertia as R, i.e. [START_REF] Dupont | Near-optimal parameterization of the intersection of quadrics[END_REF][START_REF] Dupont | Near-optimal parameterization of the intersection of quadrics[END_REF]. Moreover, the point p0 belongs to QR 0 since M0 p0 = p.

Note that the size of M0 is the size of p, i.e. 1, and thus the height of R0 in terms of S and T is the same as the height of R. For clarity, we rename R0 and p0 by R and p, respectively.

Let x denote the vector (x, y, z, w) T . Let L be 1/2 times the differential of quadric QR at p (one can trivially show that the matrix of L is the first row of R) and let i be such that R[1, i] = 0 (such an i necessarily exists). We compute the polynomial division of QR = x T Rx by Lx with respect to the variable x[i]. The result of the division is

R[1, i] 2 (x T Rx) = (Lx) (L x) + A, (1)
where

L [ξ] = -R[i, i] R[1, ξ]+2 R[1, i] R[i, ξ] for ξ = 1, . .

. , 4 and

A = c[j] x[j] 2 + c[k] x[k] 2 + 2 c jk x[j] x[k]
where j and k are equal to the two values in {2, 3, 4} distinct from i, and

c[ξ] = R[ξ, ξ] R[i, 1] 2 + R[i, i] R[ξ, 1] 2 -2 R[ξ, 1] R[i, 1] R[i, ξ], ξ ∈ {j, k}, c jk = R[j, k] R[i, 1] 2 + R[j, 1] R[k, 1] R[i, i] -(R[j, 1] R[k, i] + R[k, 1] R[j, i]) 1].
We assume in the following that c[j] = 0 (if c[j] = 0 but c[k] = 0, we exchange the roles of j and k; otherwise the analysis is different but similar and we omit it here). For clarity we denote in the following c = c[j] and r = R [1, i].

We consider the projective transformation M such that, in the new projective frame, the quadric QR has equation (up to a factor)

x T M T RM x = 4x y + z 2 -cw 2 .
In accordance with Equation (1) we choose x = Lx, y = L x. We apply Gauss' decomposition of quadratic forms into sum of squares to A and set z

= c[j] x[j] + c jk x[k] and w = x[k]
. Precisely, we define M such that its adjoint has its first row equal to L, its second row equal to L , and the last two rows equal to zero except for the entry [3, j] equal to c[j], the entry [3, k] equal to c jk , and the entry [4, k] equal to 1. Straightforward computations show that the four columns of M can be simplified by the factors rc, r, 2r, and 2r 2 , respectively. We then get

x T M T RM x = r 2 c (4xy + z 2 -det(R) w 2 ).
(

) 2
If i, j, k are equal to 2, 3, 4 respectively, M is equal to

M = 0 B B B B B @ R[2,2] -c R[2,2] R[1,3]-r R[2,3] M [1,4] -2 r 0 -r R[1,3] M [2,4] 0 0 r 2 M [3,4] 0 0 0 r c 1 C C C C C A , M [1, 4] = r (R[1, 4] (R[2, 2] R[3, 3] -R[2, 3] 2 ) + R[3, 4](r R[2, 3] -R[2, 2] R[1, 3]) + R[2, 4](R[1, 3] R[2, 3] -r R[3, 3])), M [2, 4] = r (R[1, 4] (R[1, 3] R[2, 3] -r R[3, 3]) + R[1, 3] (r R[3, 4] -R[1, 3] R[2, 4])), M [3, 4] = r (-r 2 R[3, 4] -R[2, 2] R[1, 3]R[1, 4] + r (R[1, 3] R[2, 4] + R[1, 4] R[2, 3])).
We can easily parameterize the quadric of Equation ( 2) and the parameterization of the original QR is, with δ = det (R) and (u, v) and (s, t) in P 1 (R),

M0M `ut √ δ, sv √ δ, (us -tv) √ δ, us + tv ´T . (3)
To evaluate the height of this parameterization in terms of S and T , first note that the matrix R is the matrix λS + µT of the pencil such that (λ, µ) ∈ P 1 is solution of p T (λS + µT )p = 0. So (λ, µ) = (-p T T p, p T Sp) has height 1 and R = λS + µT has height 2 in terms of S and T , since the size of p is 1. Also, δ = det R has height 4 in terms of R and thus height 8 in terms of S and T . When δ is a square, √ δ has height 4 in terms of S and T .

The heights of the four columns of M are 1, 3, 2 and 4 respectively in terms of R, and so 2, 6, 4 and 8 respectively in terms of S and T . The worst case for the height of the coefficients of the parameterization of QR happens when √ δ is a square, which we assume for the rest of the proof. It then follows from (3) that the coordinates of the parameterization of QR are polynomials of the form

ρ1 ut + ρ2 sv + ρ3 us + ρ4 tv (4)
where the heights of ρ1, . . . , ρ4 are 6, 10, 8 and 8, respectively, in terms of S and T . Let hs = 10 and ht = 8 denote the maximum heights of the coefficients of s and t in (4), respectively.

When substituting the parameterization of QR into the equation of one of the initial quadrics (say S), we obtain an equation which can be written as averaged the results. The plots of Figure 1 show that the observed height of the coefficients converges to 38 (i.e., the theoretical height) when no gcd computation is performed for simplifying the output parameterization, and converges to 36 when gcd computation is performed. We ran experiments with inputs of size up to 10,000 and observed the same limit of 36 on the height of the coefficients when gcd computation is performed. We do not have any satisfactory explanation as of why the theoretical height is not reached in that case.

a s 2 + b st + c t 2 = 0, (5)

Cubic and tangent line

We now consider the case of an intersection consisting of a cubic and a tangent line. In this case, we can parameterize the intersection using an intermediate rational quadric QR of inertia either (2, 2) or (2, 1): the pencil contains an instance of both types of quadrics.

We prove the following theoretical bounds on the height of the coefficients of the parameterizations of the cubic and the line. Proposition 3.2. When a quadric QR of inertia (2, 2) is used to parameterize the intersection, the parameterizations of the cubic and the line have asymptotic height at most in terms of S and T .

Proof. In this situation, the determinant of the intermediate quadric QR is necessarily a square because QR contains a rational line (see [START_REF] Dupont | Near-optimal parameterization of the intersection of quadrics[END_REF]). So the bounds found in the proof of Proposition 3.1 apply. In particular the heights hs, ht, ha, h b , and hc of the coefficients of Equations ( 4) and ( 5) apply. Equation ( 5) factors into two terms, one of degree 0 and the other of degree 2 in, say, (u, v), and both linear in, say, (s, t); (5) can be written as Proof. We follow the algorithm outline given in Section 2.2.3 to determine the height of the output.

(αs + βt)(α s + β t) = a s 2 + b st + c t 2 = 0,
Let λ0 be the quadruple real root of the determinantal equation. The determinantal equation can be written as det R(λ) = γ(λ -λ0) 4 .

Since the coefficients of det R(λ) have height 4 in terms of the coefficients of S and T (they are 4 × 4 determinants), λ0 has height 1. So the coefficients of the cone QR = Q λ 0 S-T have height 2.

Finding the singular point of QR amounts to finding a point c ∈ P 3 (Z) in the kernel of R, i.e. such that Rc = 0. If we decompose R such that Ru is the upper left 3 × 3 matrix of R and r4 is the first three coordinates of the last column of R, c such that cu is the first three coordinates and c4 is the last and we assume det Ru = 0 (R has rank 3 so one of its 3 × 3 minors is non-zero), then c is found by solving Rucu = -c4r4.

A solution is thus c = (-R * u r4, det Ru), where R * u is the adjoint of Ru. The coefficients of R * u and r4 have height respectively 2 and 1 in terms of the coefficients of R. So the coordinates of c have height 3 in terms of the coefficients of R, i.e. 6 in terms of the coefficients of S and T .

Since the line of the intersection is the (double) intersection of QR and the tangent plane to QS at c, a point p on this line is any point satisfying

Rp = Sc. (6)
(Observe that if p is a solution, any a1p + a2c is also solution.) The right-hand side of ( 6) has height 6 + 1 = 7 in the coefficients of S and T . As above, one can assume that det Ru = 0 and there is a unique point p having zero as last coordinate. This point satisfies pu = R * u (Sc)u and the height of its coordinates is 4 + 7 = 11. Overall, the coefficients of the line (c, p) have height 11. Now, we need to parameterize a cone (QR) containing a rational point (p). First, we apply to R a projective transformation P sending the point (0, 0, 0, 1) T to c and the point (0, 0, 1, 0) T to p. We are left with the problem of parameterizing a cone going through the origin. Such a cone has equation a1x 2 + a2xy + a3y 2 + a4yz + a5xz = 0 [START_REF] Levin | A parametric algorithm for drawing pictures of solid objects composed of quadric surfaces[END_REF] and a parameterization is given by

0 B B @ a5 0 a4 0 0 a4 a5 0 -a1 -a3 -a2 0 0 0 0 1 1 C C A 0 B B @ u 2 v 2 uv s 1 C C A ,
where (u, v, s) lives in a real quasi-projective space (see [START_REF] Dupont | Near-optimal parameterization of the intersection of quadrics[END_REF]) and where a1, a2, a3 have height 2 and a4, a5 have height 11 + 2 = 13. Lifting the parameterization to the original space by multiplying by matrix P , we find a parameterization X(u, v, s) of QR such that the coefficients of u 2 , v 2 , uv have height 13 and the coefficients of s have height 6. Plugging X(u, v, s) in the equation of any other quadric of the pencil gives an equation in the parameters of the form: where b(u, v) is a polynomial of degree 2 and height 1 + 6 + 13 = 20 and c(u, v) a polynomial of degree 4 and height 2×13+1 = 27. Observe first that a = 0 since the singularity of the cone, which is a point of the intersection, is reached at (u, v) = (0, 0) and at this point s = 0 necessarily (because X(u, v, s) is a faithful parameterization of the cone). We also know that (8) has a linear factor corresponding to the line of the intersection and it is not too difficult to realize that this factor is a5 u + a4 v, where a4 and a5 are as in [START_REF] Levin | A parametric algorithm for drawing pictures of solid objects composed of quadric surfaces[END_REF]. So, after factoring out the linear term, (8) can be rewritten

as 2 + b(u, v)s + c(u, v) = 0, (8)
b (u, v)s + c (u, v) = 0,
where b (u, v) is a polynomial of degree 1 and height 20 -13 = 7 and c (u, v) is a polynomial of degree 3 and height 27 -13 = 14. We can solve this equation rationally for s and, multiplying to clear the denominators, we get a parameterization of the cubic of height 14 + 6 = 20.

The difference in the heights of the parameterizations underscored in the above two propositions is vindicated by some experiments we made. In the smooth quartic case, random examples could be generated by taking input quadrics with random coefficients. But we have to proceed differently here. We start with a canonical pair of quadrics intersecting in a cubic and a tangent line and apply to this pair a random transformation. More precisely, given a canonical pair S, T , four random bigints r1, r2, r3, r4 and a random real projective transformation P , we take as new input quadrics the quadrics of matrix: S = P T (r1S + r2T )P, T = P T (r3S + r4T )P.

If we take all the random coefficients in the range [-3 √ 10 s , 3 √ 10 s ], then the quadrics S and T have size s (the size of the canonical pair S, T can be neglected).

Figure 2 shows the height of the coefficients of the parameterization of the cubic when a quadric QR of inertia (2, 2) or (2, 1) is used (with mild optimizations turned on). The plots clearly show that the coefficients of the cubic are smaller when a cone is used to parameterize the intersection. The fact that the observed heights are, in the limit, so different from the theoretical bounds (8 instead of 22 when a cone is used) is possibly a consequence of the way S and T are generated: it certainly does not reflect a truly random distribution in the space of quadrics with integer coefficients of size s intersecting in a cubic and a tangent line.

Figure 3 further reinforces our choice of using a cone: the parameterizations have not just smaller coefficients, they are

EXPERIMENTAL RESULTS

We now report on some experimental results and findings from our implementation.

The experiments were made on a Dell Precision 360 with a 2.60 GHz Intel Pentium CPU. LiDIA, GMP and our own code were compiled with g++ 3.2.2.

Implementation versions

In what follows, we compare three versions of our implementation:

• unoptimized : nothing is done to simplify the coefficients either during the computations or in the parameterizations computed; • mildly optimized : some gcds are performed at an early stage (optimization of the coefficients and of the roots of the determinantal equation, optimization of the coordinates of singular and rational points, . . . ) to avoid hampering later calculations with unnecessarily big numbers; Let us finally mention that we tried a fourth version of our implementation where the extraction of the square factors is done by fully factoring the numbers (using the Elliptic Curve Method and the Quadratic Sieve implemented in LiDIA [START_REF]A C++ Library for Computational Number Theory[END_REF]). But this version is almost of no interest: for small input coefficients, the strongly optimized version already finds all the necessary factors, and for medium to large input coefficients, integer factoring becomes extremely time consuming.

•

Performance evaluation

Let us first discuss the impact of the MAXFACTOR variable on the output. Figure 4 shows that a value of 10 5 has a dramatic impact on computation time while all values less than 10 4 are acceptable. We have determined that the best compromise between efficiency and complexity of the output is obtained by setting MAXFACTOR to 10 3 , which we assume

now.

Figure 5 shows the evolution of the aggregate computation time in the smooth quartic case, which is the most computationally demanding case, with the three versions outlined above. What we infer from these plots is first that the computation times for the unoptimized and mildly optimized versions are very similar, while we observe (see Figure 1) a dramatic improvement in the height of the output coefficients with the mildly optimized version for reasonably small inputs. This explains our choice of putting the mild optimizations in the form of a preprocessor directive, not a binary argument: they might as well have been called mandatory optimizations.

A second lesson to be learned from Figures 1 and5 is that for a size of input ranging from roughly 5 to 60, the computation time is roughly 30% larger for the strongly optimized version than for the mildly optimized. At the same time, the height of the output is between 20% (input size of 5) and 5% (input size of 60) smaller. For large values of the input size (Figure 4), the difference in computation time between the mildly optimized and the strongly optimized versions drops to less than 10%, but not much is gained in terms of height of the output.

Another interesting piece of information inferred from Figure 1 is that the standard deviation of the height of the output coefficients is large for small input size in the strongly optimized version. This means that in the good cases the height of the output is dramatically smaller than the height in the mildly optimized case, and in the bad cases is similar to it.

Deciding to spend time on optimization essentially depends on the application. For most real-world applications, where the size of the input quadrics is small by construction, we believe optimizing is important: it should be kept in mind that the computed parameterizations are often the input to a later processing step (like in boundary evaluation) and limiting the growth of the coefficients at an early stage makes good sense.

A last comment that can be made looking at Figure 4 concerns the efficiency of our implementation. Indeed, those plots show that we can compute the parameterization of the intersection of two quadrics with coefficients having 400 digits in 1 second and 1,000 digits in 5 seconds (on average).

Efficiency can be measured in a different way. In Figure 6, we have plotted the total computation time, with the strongly optimized version, for a file containing 120 pairs of quadrics covering all intersection situations over the reals. The "random" quadrics were generated as in Section 3.2. For an input size s = 500, the total computation time is roughly 72 seconds, on average, for the 120 pairs of quadrics, i.e. 0.6 second per intersection. This should be compared to the 1.7 seconds on average needed to compute the intersection in the smooth quartic case for the same size of input (Figure 4). This difference is simply explained by the fact that very degenerate intersections (like when the determinantal equation vanishes identically, which represents 36 of the 120 quadrics in the file) are usually much faster to compute.

Our last word will be on memory consumption. Our implementation eats up very little memory. In the smooth quartic case, the total memory chunks allocated sum up to less than 64 kilobytes for input sizes up to 20. It takes input sizes of more than 700 digits to get to the 1 MB range of used memory.

EXAMPLES

We now give three examples of parameterizations computed by our algorithm. Other examples can be tested by querying our parameterization server.

Comparing our results with the parameterizations computed with other methods does not make much sense since our implementation is the first to output exact parameterizations in all cases. However, for the sake of illustration, our first two examples are taken from the paper describing the plane cubic curve method of Wang, Joe and Goldman [START_REF] Wang | Computing quadric surface intersections based on an analysis of plane cubic curves[END_REF]. the quadric QR of inertia (2, 2) used to parameterize the intersection is such that its determinant is not a square.

As explained in Section 2.2.2, the parameterization is thus only near-optimal in the sense that it is possible, though not necessary, that the square root can be avoided in the coefficients.

It turns out that in this particular example it can be avoided. Consider the cone QR corresponding to the rational root (-1, 21) of the determinantal equation:

QR : -QS + 21 QT = 2x 2 -y 2 -w 2 .
QR contains the obvious rational point (1, 1, 0, 1), which is not its singular point. This implies that it can be rationally parameterized. Plugging this parameterization in the equation of QS or QT gives a simple parameterization for C:

X(u, v) = 0 B B B B B @ u 2 + 2 v 2 2 uv u 2 -2 v 2 0 1 C C C C C A ± 0 B B B B B @ 0 0 0 1 1 C C C C C A √ 2 u 4 +4 u 2 v 2 +8 v 4 .

Example 3

Our last example illustrates the fact that our implementation is complete in the sense that it computes parameterizations in all possible cases. It concerns two quadrics intersecting in two tangent conics. The execution trace (with debugging information turned off) is given in Output 2.

CONCLUSION

We have presented a C++ implementation of an algorithm for parameterizing intersections of quadrics. The implementation is exact, efficient and covers all the possible cases of intersection. This implementation is based on the LiDIA library and uses the multiprecision integer arithmetic of GMP.

Future work will be devoted to understanding the gaps between predicted and observed values for the height of the coefficients of the parameterizations, to working out predicates and filters for making the code robust with floating point data (many classes and data structures have already been templated for a future use with floating point coefficients) and to porting our code to the CGAL geometry algorithms library.

  where a, b and c depend on (u, v). It follows from (4) that the heights of a, b, c are ha = 21 (= 2hs + 1), h b = 19 (= hs + ht + 1) and hc = 17 (= 2ht + 1), respectively. When substituting the solution (s = 2c, t = -b ± √ b 2 -4ac) into the parameterization (3) we obtain a parameterization of the smooth quartic in which each coordinate has the formα(u, v) ± β(u, v) p b 2 -4ac.The height of the coefficients of α is 27 (= hs +hc = ht +h b ), the height of the coefficients of β is 8 (= ht) and the height of the coefficients of ∆ is 38 (= 2h b = ha + hc).

Figure 1

 1 Figure1shows how the observed height of the coefficients of ∆(u, v) evolves as a function of the input size s for three different versions of our implementation (see Section 4 for the details). For each value of s in a set of samples between 0 and 60, we have generated random quadrics with coefficients in the range [-10 s , 10 s ], computed the height of the coefficients of the parameterization of the smooth quartic and

Figure 1 :

 1 Figure 1: Evolution of the height of ∆(u, v) (smooth quartic case) as a function of the size of the input, with the standard deviation displayed on the optimized plot.

Proposition 3 . 3 .

 33 where α, β are constants and α , β are polynomials in (u, v). Since αβ + βα = b, α and the coefficients of α have height at most h b . Similarly, ββ = c thus β and the coefficients of β have height at most hc. Substituting the solutions (s = β, t = -α) and (s = β , t = -α ) into the parameterization (3), we get parameterizations of the cubic and the line whose coefficients have height at most hs + hc = ht + h b = 27. When a quadric QR of inertia (2, 1) is used to parameterize the intersection, then asymptotically the parameterization of the line has height at most 11 and the parameterization of the cubic has height at most 20 in terms of S and T .

Figure 2 :

 2 Figure 2: Observed height of the parameterization of the cubic in the cubic and tangent line case.

Figure 3 :

 3 Figure 3: Computation time for the cubic and tangent line case.

  strongly optimized : mildly optimized, plus extraction of the square factors of some bigints (like in the smooth quartic case, where √ δ = √ det R can be replaced by b √ a if det R = a b 2 ) and gcd simplifications of the coefficients of the final parameterizations. For the extraction of the square factors of an integer n, the strongly optimized version finds all the prime factors of n up to min ( 3 √ n , MAXFACTOR), where MAXFACTOR is a predefined global variable.

Figure 4 :

 4 Figure 4: Evolution of execution time in the smooth quartic case as a function of the size of the input for very large input sizes.

Figure 5 :

 5 Figure 5: Evolution of execution time in the smooth quartic case as a function of the input size, with the standard deviation shown on the optimized plot.

Figure 6 :

 6 Figure 6: Computation time for 120 pairs of quadrics covering all intersection cases, with standard deviation.

>> quadric 1 :

 1 -4*x^2 -56*x*y -24*x*z -79*y^2 -116*y*z + 70*y*w -85*z^2 -20*z*w + 9*w^2 >> quadric 2: 6*x^2 + 84*x*y + 36*x*z + 45*y^2 + 160*y*z -210*y*w + 131*z^2 + 30*z*w -45*w^2 >> complex intersection: two tangent conics >> real intersection: two tangent conics >> parameterization of conic: [-39*u^2 + 443*u*v -7254*v^2, 3*u^2 -66*u*v + 1388*v^2, 6*u^2 -132*u*v + 701*v^2, -5*u^2 + 110*u*v -3005*v^2] parameterization of conic: [-39*u^2 + 443*u*v -4004*v^2, 3*u^2 -66*u*v + 1138*v^2, 6*u^2 -132*u*v + 201*v^2, -5*u^2 + 110*u*v -1205*v^2] >> time spent: 10 ms

If however one of the initial quadrics has rank 3, then it should be used to parameterize the intersection. Doing so results in a parameterization having the same algebraic complexity in the worst case, but of smaller coefficient size.

ACKNOWLEDGMENTS

The authors wish to acknowledge Laurent Dupont for a preliminary implementation of the parameterization algorithm in MuPAD, Guillaume Hanrot for his C implementation of Uspensky's algorithm, Daniel Lazard for his help in designing the parameterization algorithm, and Etienne Petitjean for his NetTask socket management tool which makes possible the querying of our parameterization software via a web interface.

Example 1

Our first example is Example 4 from [START_REF] Wang | Computing quadric surface intersections based on an analysis of plane cubic curves[END_REF]. The two quadrics are a quadric of inertia (2, 1) (elliptic cylinder) and a quadric of inertia (2, 2) (hyperboloid of one sheet). The curve of intersection C has implicit equation

A rendering of the intersection is given in Figure 7.a (made with Surf [START_REF]Surf : A Tool To Visualize Algebraic Curves and Surfaces[END_REF]).

In [START_REF] Wang | Computing quadric surface intersections based on an analysis of plane cubic curves[END_REF], the authors find the following parameterization for C: The authors report a computation error on this example (measured as the maximum distance from a sequence of sample points on the curve to the input quadrics) of order O(10 -7 ).

Our implementation outputs the following exact and simple result in less than 10 ms:

The polynomials involved in the parameterization are defined in Q[u, v], which means we are in the lucky case where the intermediate quadric of inertia (2, 2) found to parameterize the intersection has a square as determinant. So the parameterization obtained is optimal in the extension of Z on which its coefficients are defined.

Example 2

Our second example is Example 5 from [START_REF] Wang | Computing quadric surface intersections based on an analysis of plane cubic curves[END_REF]. It is the intersection of a sphere and an ellipsoid that are very similar (see Figure 7.b):

In [START_REF] Wang | Computing quadric surface intersections based on an analysis of plane cubic curves[END_REF], the authors compute the parameterization (9) with Our implementation gives the result displayed in Output 1. Since the polynomials of X(u, v) involve a square root √ 10,


Fatima Nasser

Zhong-Yang Li

Nadine Martin

Philippe Gueguen
email: pilippe.gueguen@obs.ujf-grenoble.fr

Non-biased estimation of modal parameters of multi-component ambient vibrations in high-rise buildings using a self-governed random decrement technique

The estimation of modal parameters of a structure by applying the random decrement technique has been introduced since late 1960's, and is widely used because of its simplicity and rapidity. Since this technique can be applied to mono-component signals only, a filtering process is a necessary preliminary step for getting the random decrement signature when applying on multi-component signals. But filtering is the most critical step with a high influence on the accuracy of the modal parameter estimation, especially when the signal is composed of very closely-spaced modes. In this paper, we propose to apply the random decrement technique directly on a multicomponent signal without the filtering step, in order to provide what we called a multimode signature. The idea is then to propose an adapted method to separate the different modes in this multi-mode signature. For that purpose the modal parameters are first roughly approximated in a non-optimal way and with a bias from a time-frequency distribution. Based on this initialization, a final non-biased estimation is achieved via maximum likelihood estimator. The results conducted in this paper, over simulated and real world ambient vibrations on high rise buildings, suggest positively that the proposed method is totally autonomous and is capable of yielding very good estimation results.

Introduction

In civil engineering, the study of the dynamic behaviour of the building plays a significant role in detecting its possible damages, and determining its safety and reliability. The dynamic properties of the building are directly related to the mass and the stiffness, and the modal parameters are the basic contributor in giving information about the inherent dynamic properties of the building. So that, knowing the modal parameters would provide sufficient information about these two physical properties. As such, it is of high importance to estimate these parameters with high precision.

However, the signals measured on the buildings under study are the response of these buildings to the ground excitation which is assumed to be white Gaussian noise. To estimate the modal parameters for such type of signals, the impulse response should be firstly estimated, and thus a deconvolution method is a must. Here comes the motivation of using the Random Decrement Technique (RDT). The RDT was developed by Henry Cole [START_REF] Cole | On-the-line Analysis of Random Vibrations[END_REF] in the late 1960's. Driven by the motive of online failure detection by monitoring the damping of a structure excited by white noise. Cole intuitively came up with a process of averaging the time segments with the same initial conditions to get a function proportional to the free decay. When a large enough number of selected samples from a very long time history conditions are averaged, the random component of response averages out to zero and what survives is only the free decay, or the impulse response, or what is known as the Random Decrement Signature (RDS) [START_REF] Cole | On-the-line Analysis of Random Vibrations[END_REF] . This technique was given a mathematical basis later on by Vandiver et al in 1982 [START_REF] Vandiver | A Mathematical Basis of the Random Decrement Vibration Signature Analysis Technique[END_REF] , who showed that for a practical case of a linear time-invariant system excited by a zeromean, stationary, Gaussian random process, the RDS of the output is shown to be proportional to the autocorrelation of the process [START_REF] Vandiver | A Mathematical Basis of the Random Decrement Vibration Signature Analysis Technique[END_REF] . The modal parameters are estimated using the RDS [START_REF] Asmuussen | Modal Analysis Based on the Random Decrement Technique: Application to Civil Engineering Structures[END_REF][START_REF] Ibrahim | Random Decrement Technique for Modal Identification of Structures[END_REF] . Thus in order to obtain credible modal parameter estimations, one would have to pay great attention on the estimation of the RDSs themselves. In practice, the RDT is applicable over mono-component signals only, thus a filtering process is a necessary preliminary step for getting the RDS when applying over multicomponent signals.

During the study of the RDT, it is observed that filtering is one of the most delicate and critical processes, and it is difficult to be performed firmly, especially when the multicomponent signal is composed of very closely spaced frequency modes. Firstly, we intended to tune the filtering bandwidth manually, as it is prevalent in most cases, e.g. [START_REF] Zubaydi | The Filtering Effect of Random Responses of Stiffened Plates on Their Random Decrement Signatures and Natural Frequencies[END_REF] , however, the manual tuning of the filtering process is found to be very time consuming with high influence on the modal parameter precision. Then in Surveillance 7 [START_REF] Nasser | Automatic Parameter Setting of Random Decrement Technique for the Estimation of Building Modal Parameters[END_REF] , we proposed to set the filter bandwidth automatically. This proposition has made the estimation of the modal parameters faster, easier, reliable, and independent of user interaction. But it doesn't however rule out the user know how completely, because there are still several factors which require user intervention, for e.g. the choise of the filter type and the filter order. From here, we thought that it would be great if this step could be totally skipped. Even though, the RDT has been widely investigated since its development [START_REF] Antoni | The Discrete-Time Random Decrement Technique: Closed-form Solutions for the Blind Identification of SIMO Systems[END_REF][START_REF] Yang | Application of the Random Decrement Technique in the Detection of an Induces Crack on an Offshore Platform Model, Computer Methods for Offshore Structures[END_REF] , however, it had been rarely carried out without filtering. This work responds to this shortage and proposes a novel approach to extend the applicability of the RDT for modal parameter estimation without the necessity of the filtering process. This makes the originality of this paper. The application of the RDT over multi-component signals without passing by the filtering process leads to the estimation of what we called multi-mode RDS. Working with the multi-mode RDS allowed the proposition of a new method for the modal parameter estimation, mainly the frequency and the damping ratio. The proposition states first to roughly approximate the modal parameters, using a nonoptimal way. Based on this initialization step, a final non-biased estimation is achieved from a maximum likelihood estimator.

Nevertheless, when dealing with the RDT, one should not turn a blind eye toward a very important aspect, which is the processing parameters of this technique, that play a significant role in the accuracy of the RDS extraction and thus the modal parameter estimation. Proceeding from this fact, an intense study of the influential parameters of the RDT, mainly on the damping ratio estimation, has been carried out through this work, and will be briefly presented throughout this paper. The rest of this paper is organized as follows: the problematic of this paper is formulated in section 2. Section 3 shows a brief introduction to the RDT and the estimation of the multi-mode RDS. The proposed method is then explored in section 4. Section 5 presents briefly the processing parameters of the RDT and their effect on the modal parameter estimation. Section 6 illustrates the obtained results over simulated and real world signals. Results are finally discussed in section 7 in terms of relevance and perspective.

Problem formulation

The signal under study in this paper is considered as

[ ] [ ] [ ] [ ] n e n P n + * = h n Y , (1)
The index n is interpreted as discrete time

, [ ] [ ] ∑ = = N k k n h n 1 h [ ] N k , 1 ∈ ∀
where N is the number of modes in the signal and being equal to 1 for single degree of freedom systems (SDOF), [ ] ) ( )

n e f n Dk n f k k k ω π π ξ sin 2 1 2 - - = k h , (2)
where

2 1 2 k Dk ξ π ω - = k f the pseudo pulsation. ) (t e
is an additive white Gaussian noise of zero mean and unknown variance. In order to estimate the modal parameters of such types of signals we have to deal with the impulse response separately, here comes the motivation behind using the RDT, since it is a method able to estimate a signature similar to the impulse response of the building known as the RDS. The way of estimating the RDS by using the RDT is explained in section 3.

Multi-mode Random Decrement Signature

The RDT is a method which transforms the stochastic process [ ] n Y , that is assumed to be stationary process, into RDS. The RDS is defined as the mean value of a stochastic process on a condition of the process itself [START_REF] Asmuussen | Modal Analysis Based on the Random Decrement Technique: Application to Civil Engineering Structures[END_REF][START_REF] Ibrahim | Random Decrement Technique for Modal Identification of Structures[END_REF] ,

[ ] [ ] [ ] [ ] n Y Y T n Y E D τ τ + = , (3)
An RDS is referred to

) (τ Y D
, where τ indicates the process time lag,

[ ] [ ] [ ] n Y T n Y E τ + is the mean value of [ ] τ + n Y
being the process from which the mean value is calculated, on [ ] n Y T being the process where the condition is fulfilled.

For a time series, the RDS can be estimated unbiased as an empirical mean

[ ] [ ] [ ] ) ( 1 ˆ1 ∑ + = = M i n y i Y i T n y M D τ τ , (4)
where

) ( ˆτ Y D
is the estimated RDS, [ ] n y is the measurement and M is the number of points in the process which fulfils the triggering condition

[ ] M i , 1 ∈ ∀
. The initial condition of the time segments in the averaging process are presented by the triggering condition at the time lag 0 = τ . Asmussen [START_REF] Asmuussen | Modal Analysis Based on the Random Decrement Technique: Application to Civil Engineering Structures[END_REF] has defined the different triggering conditions used as a specific formulation of the applied general triggering condition,

[ ] [ ] [ ] { } 2 1 2 1 , b n Y b a n Y a T n Y < ≤ < ≤ =  . ( 5
)
where 1 a ,

n Y  is the discrete time differentiation of [ ] n Y .
Asmussen [START_REF] Asmuussen | Modal Analysis Based on the Random Decrement Technique: Application to Civil Engineering Structures[END_REF] , as well, has presented the most frequently used triggering conditions as level crossing,

{ } a t Y T L t Y = = ) ( ) ( , positive point, { } 2 1 ) ( ) ( a t Y a T P t Y > ≤ = , local extrema, { } 0 ) ( , ) ( 2 1 ) ( = > ≤ = t Y a t Y a T E t Y 
, and zero-up crossing,

{ } 0 ) ( , ) ( 2 ) ( > = = t Y a t Y T Z t Y



. The number of triggering points controls the estimation time and the accuracy of the estimates. All the modes of the input signal are preserved in the estimated RDS, so that if the input is single-mode signal, then the estimated RDS will contain only a single mode, and vice versa. Since we proposed to apply the RDT directly over the multi-component ambient vibration signal under study, so the estimated RDS will be the multi-mode RDS that contains the same natural frequencies

[ ] [ ] ∑ = = N k k m s m s 1 with [ ] ) 1 2 cos( 2 2 0 k k m f k k m f e A m s k k ξ π ξ π - = - , ( 7
)
where m is interpreted as the time index of the RDS, with

N k ≤ ≤ 1
and k A 0 being the initial amplitude of each mode.

A new modal parameter estimation algorithm

Despite of their rapidity and simplicity, Fourier transform based approaches are hardly able to estimate and track the model parameters with no bias. Parameteric or modelbased estimators can give significant improvements in the time-frequency resolution at the expense of a higher computational complexity. Accordingly, a non-biased estimation is achieved by using the parametric Maximum Likelihood Estimator (MLE), this is a deterministic approach that is maximized by means of a stochastic technique based on the simulated annealing algorithm. However, good starting values of the parameters are helpful in reducing wanders among local minima. Therefore, we initialize the parameter vector of MLE by using the nonparametric spectral based approach.

Maximum Likelihood Estimation

This section describes a parametric spectral estimator method using MLE in order to find the optimum set of parameters that minimize the error function. The aim is to reduce the bias in the frequency resolution of the conventional methods by using the characteristics associated to the MLE. This will allow the computational of higher accurate damping ratios, and increasing the frequency resolution. The parameters of each mode of (7) form a vector,

[ ] [ ] T k k k k T k T k f A , , , ,..., 0 0 1 φ ξ = = θ θ θ , ( 8
)
where T indicates the transpose, and k 0 φ the initial phase of each mode. We consider the MLE of θ being equivalent to the minimization of a least square approach under the assumption of white Gaussian additive noise. Thus it results in the following non-linear equation

) ( min arg ) ( max arg ˆ4 4 θ θ θ θ θ LS MLE l l k k × × ℜ ∈ ℜ ∈ = = , (9)
[ ] [ ]

2 2 / 2 / 4 min arg ) ( ∑ - = - = ℜ ∈ × L L m LS m s m y l k θ θ , ( 10
)
where L is the sample number. Direct minimization of ( 10) is difficult due to the high non-linearity of the function and the parameter number. Classical optimization techniques such as gradient descent, Gauss-Newton and EM algorithm do not ensure convergence to the global minimum when local minima are numerous. This problem can be overcome with meta-heuristic approaches, and in particular, with the simulated annealing algorithm [START_REF] Martin | Frequency and amplitude tracking for short nonstationary and nonlinear signals[END_REF] . The simulated annealing technique is efficient when a desired global extremum is hidden in many local extrema. The simulated annealing has an analogy with thermodynamics where metal cools and anneals. In the same way, after an initialization of the parameter vector to estimate as will be explained in section 4.2, an iterative loop, controlled by a scalar referred to as a temperature, generates a new candidate of the vector that minimizes the cost function, namely the least square function

) ( θ LS l
. This candidate is accepted or not according to the Metropolis acceptance rule.

Spectral based estimation

In the aim of obtaining a rough estimation of the modal parameters, the classical, nonparameteric spectral based approach from time-frequency distribution will be used. This step will serve as an initialization step to provide a priliminary estimation for the MLE approach.

Accordingly, the usage of two windows centered at . The average of these two periodograms yields the Welch spectral estimator of the multi-mode RDS, from which for each detected frequency peak a related damping ratio can be calculated using the magnitude of the detected peak at the two different instants 1 t and 2 t . One of the drawbacks of the STFT approach is that the frequency resolution is limited due to the windowing of the RDS. The usage of the RDT with the multi-mode RDS has posed some problems which necessitates the intense study of the processing parameters of this technique, that play a significant role in the accuracy of the RDS extraction and thus the modal parameter estimation.

RDT Processing parameter study

For the efficient use of the RDT, the data processing parameters associated with the technique need to be understood. Some of these parameters play a significant role in accurate estimation of modal parameters; the parameter settings in this work are judged depending on their influence on the accuracy of the estimated damping ratio. Based on this judgment, the optimum parameter settings will be chosen for implementation of RDT over multi degree of freedom systems. On one hand, the length of the RDS, rds L (points) , and the length of the signal under study, sig L (points) , are the operational parameters of RDT that are manipulated according to the user requirements. Those parameters should be carefully chosen in order to study their effect on the damping ratio estimation. On the other hand, the sampling frequency of the signal under study, s F , is a characterizing parameter of the input, and it is taken into account in the behaviour of the algorithm, and the accuracy of the parameters to estimate. That are in turn, the natural frequency, k f , and the damping ratio, k ξ , of each mode k in the signal under study. Then the influence of 2 parameters are to be studied according to the estimation of 2 modal parameters, and taking into consideration the signal sampling frequency. When the natural frequency is normalized by the signal sampling frequency, both rds L and sig L can then be normalized in terms of their number of periods, such as,

) k rds rds R(f L N =
(periods) , and

) k sig sig R(f L N = (periods) , with k s k f F R(f = ) riod) (points/pe
. Subsequently, the parameters to be studied are rds N and sig N . Accordingly, those processing parameters and their impact on the accuracy of damping ratio estimate using RDT were evaluated firstly on SDOF systems. The conclusion of the processing parameter study is based on 20 independent noise realizations. sig N was chosen to range from 500 periods to 4000 periods . For should be chosen as the result of a compromise. On one hand, the signal length should be neither too long in order not to loose the time resolution of the signal, nor too short in order to assure that it contains a large number of averaged segments that yield the exponentially decaying form. On the other hand, the length of the signature should be large enough to assure that it is of exponentially decaying form, however, not too long in order to prevent the beating phenomena that appears with very long RDS. Noting that it should not be shorter than two periods, otherwise the decaying phenomena will be deteriorated.

From the results of this study that are summarized in fig. 2, it is always recommended to be in the green region which yields an acceptable estimation error of the modal parameters of % 30 ≤ error Estim . It is not recommended to be in the orange region that gives a high estimation error of

% 50 ≤ < error Estim 30%
, however, if it is unavoidable to be in that region it should be expected that the estimations might be less reliable, and for sure it is totally unrecommeded to be in the red region since it associates a very high, unacceptable estimation error of % 50 > error Estim .

Application on ambient vibration signals

Buildings are permanently excited by some natural solicitation sources, for example, the background seismic noise from the earth, the natural phenomena around the building, internal source like human steps, vihicles, and the different types of rotating machineries in the buildings, and the natural activity of the earth. These different excitations induce various types of vibrations in the building, which are referred to as ambient vibrations. This section is aimed at the discussion of the accuracy and efficiency of the proposed algorithm in this paper for the modal parameter estimation over simulated ambient vibration signals, and real-world ambient vibration signals recorded on real-world buildings. This discussion is guided by the experience gained from the intense study over the RDT processing parameters, explained in section 5, using the "optimum" parameter settings needed to yield reliable estimations.

Number of periods in the signal

Number of periods in RDS

Analysis of simulated ambient vibration signals

To study the relevancy of the proposed algorithm on modal parameter estimation, three simulated signals have been devised. These signals are of different configurations, as shown in Table 1, to assure the general functionality of the proposition. Based on the intense study of the processing parameters of the RDT briefed in section 5, sig N and rds N for all the three simulated signals were selected from the recommended green zone in fig. 1, as 5 and 3100 periods respectively. It is worth noting as well that the triggering condition used in this section was the positive point type, this was encouraged by the experimental study we've done over the RDT processing parameters, and by what was recommended by Assmusen [START_REF] Asmussen | Random Decrement: Identification of Structures Subjected to Ambient Excitation[END_REF] to analyse ambient data. The

. 2 , 1 . 2 (%) k ξ [ ] 7 , 5 , 3 , 1 [ ] 2
f ˆ(Hz) k ξ ˆ(%) k f ˆ(Hz) k ξ ˆ(%) k f ˆ(Hz) k ξ ˆ(%)
1 st mode As shown in the table, both natural frequencies and damping ratios were estimated with better precision for higher energy modes (modes in lower frequencies), except for at which the modes are very closely spaced in frequency. This result is expected when dealing with optimization methods like the simulated annealing, as it is very difficult to estimate the correct parameters of the weak energy modes with such types of methods.

k k k f f Energy ξ × ∝ 2 1 ) ( , (12)
As presented by equation (12), it is not surprising that for all the signals of this section, the estimations of the last mode is not as good as the other modes in lower frequencies, as its energy is much lower than the others which makes it indifferentiable from the noise. Nevertheless, all the measures of normalized errors indicate that on average the proposed model exhibit very good estimation performance.

Analysis of real-world ambient vibration signals

The signal of this part is a real-world ambient vibration signal, recorded at the top of the 13 th floor city hall building of Grenoble (France), fig. 3, using multiple sensors placed on several stories to measure simultaneously the vibrations in three directions, longitudinal, transverse and vertical. This signal is sampled at 200 Hz.

For this signal we study the analysis uniquely over the data in the longitudinal direction.

Based on the information provided to us by the geophysicists, we will be focusing on the first three modes represented respectively as, the first longitudinal at 1.15 Hz, the first transverse at 1.2 Hz, and the first torsion mode at 1.44 Hz. The results reported in Table 3 confirm that the proposed algorithm does to a large extent translate into very acceptable, yet very good modal parameter estimation performance even over signals with very closely spaced frequency modes which has not been done before.

Conclusions

Unlike the classical random decrement technique, for which the filtering process is an essential preliminary step for estimating the random decrement signature when applying on multi-component signals, this paper demonstrates how to use the Random Decrement Technique without filtering for modal parameter estimations of buildings subjected to ambient excitations. This is done by proposing a new method for modal parameter estimations, which allows dealing with multi-mode Random Decrement Signatures instead. The resulting multi-mode RDS are then used to identify the modal parameters of the signal under study by using a non-parametric method from a time-frequency distribution with bias on the frequency estimations, as an initialization step. This step is then followed by a final non-biased parametric method using the maximum likelihood estimator for an accurate, and robust modal parameter estimations. The proposed algorithm helped in making the RDT totally autonomous and more friendly user, and the estimation of the modal parameters faster, easier, reliable, and totally independent of user interaction.

In the future we aim to propose a modal parameter modelization method which can estimate directly the number of modes without any initialization step.

  ambient vibration of the building, * indicates the convolution between the impulse response [ ] n hof that building and the seismic noise ) (t P that is assumed to be white Gaussian noise of zero mean and unknown variance. The impulse response is characterized by its natural frequency )

2 b

 2 are the triggering levels, and [ ]

  response of the signal under study presented in equation(2). We assume that the model of the multi-mode RDS is defined as [ time process consisting of a deterministic multi-component process [ ] m s , embedded in the residual of the RDT algorithm which can be considered as an additive white Gaussian noise [ ] m e RDS with zero mean and unknown variance,

  rds L being the length of the estimated multi-mode RDS as shown in fig. 1. Thus two periodograms at two time instants

Figure 1 .

 1 Figure 1. The estimated multi-mode RDS (left), the periodograms calculated using the STFT with a window length equal to

  range from 3 to 12 cycles; along with a damping ratio ranges between 1% to 5% where values of real world signals can be found. The SDOF systems analyzed experimentally also have estimated different signal to noise ratio (SNR) for the additive noise;

Figure 2 .

 2 Figure 2. Good regions (green), intermediate regions (orange) and bad regions (red) when setting the number of periods of the RDS according to the signal

Figure 3 .Table 3 .

 33 Figure 3. Grenoble city hall, France.

Table 1 . Three simulated signals 1 Sig , 2 Sig and 3 Sig with different configurations 1 Sig 2 Sig 3

 1

	Sig

Table 2 . The estimated modal parameters of 1 Sig , 2 Sig and 3 Sig along with their normalized error reported in parenthesis 1 Sig 2 Sig 3 Sig k

 2

Table 2

 2 summarizes the results of the modal parameter estimation for the three simulated signals. Parameter estimations shown in the table are the natural frequencies and the damping ratios with normalized errors reported in parantheses.

Sig

Acknowledgements

This work has been supported by French Research National Agency (ANR) through RISKNAT program (project URBASIS ANR-09-RISK-009).

Number of periods in the signal Number of periods in RDS