Yurii Fedorovich Smirnov

Maurice R Kibler

In memoriam

Keywords: numbers: 21.10.-k, 21.60.-n, 24.10.-i, 02.20.-a

describe some particular facets of the contributions of the late Professor Smirnov in theoretical physics and mathematical physics. They also relate some personal

I. INTRODUCTION

Yurii Fedorovich Smirnov passed away in 2008. He was a famous physicist with interests in nuclear physics, atomic and molecular physics, condensed matter physics and mathematical physics. More generally, he was at the origin of deep achievements around symmetry methods in physics.

I have not the shoulders to carry the weight of all the fields in which Yurii was recognized as a superb researcher. It is enough to say that he contributed to many domains of mathematical physics (e.g., finite groups embedded in compact or locally compact groups, Lie groups and Lie algebras, quantum groups, special functions) and theoretical physics (e.g., nuclear, atomic and molecular physics, crystal-and ligand-field theory). Let me mention, among other fields, that he obtained alone and with collaborators very important results in the theory of clustering in nuclear systems, in projection operator techniques, in the theory of heavy ion collisions, and in the so-called J-matrix formalism for quantum scattering theory. As another major contribution, he proposed a method, the so-called (e,2e) method, an analog of the (p,2p) method of nuclear physics.

It is not the purpose of these notes to give an extensive list and analyse the numerous papers by Yurii. I shall rather focus on some particular facets of his works corresponding to the above mentioned fields. I shall also devote a part of these notes to some more personal reminiscences related to some of his other fields of interest (crystal-field theory, nonbijective transformations, and the use of quantum groups in nuclear physics). Yurii Smirnov signed and/or co-signed eleven books and more than 250 scientific articles (only a few part of them shall be mentionned here). He translated into Russian several scientific books. As an example, he translated a book on the harmonic oscillator written by Marcos Moshinsky in 1969, precisely the book which was a starting point for their common book, on the same subject, published in 1996 [START_REF] Moshinsky | The Harmonic Oscillator in Modern Physics[END_REF]. He was a member of the editorial board of several journals and a councillor of the scientific councils of the Skobeltsyn Institute of Nuclear Physics and of the Chemistry Department of Moscow State University as well as of the Institute for Theoretical and Experimental Physics (ITEP) in Moscow.

III. A FIRST SEMINAL WORK

Following the pioneer work of Talmi [START_REF] Talmi | [END_REF], Marcos Moshinsky and Yurii F. Smirnov were independently interested in the description of pairs of nucleons in a harmonic-oscillator po-tential. In 1959, Moshinsky developed a formalism to connect the wave functions in two different coordinate systems for two particles (with identical masses) in a harmonic-oscillator potential [3]. In this formalism, any two-particle wave function |n 1 ℓ 1 , n 2 ℓ 2 , λµ , expressed in coordinates with respect to the origin of the harmonic-oscillator potential, is a linear combination of wave functions |nℓ, NL, λµ , expressed in relative and centre-of-mass coordinates of the two particles. The so-called transformation brackets nℓ, NL, λ|n 1 ℓ 1 , n 2 ℓ 2 , λ make it possible to pass from one coordinate system to the other. Moshinsky gave an explicit expression of these coefficients in the case n 1 = n 2 = 0 and derived recurrence relations that can be used to obtain the coefficients for n 1 = 0 and n 2 = 0 from those for n 1 = n 2 = 0 [3]. Along this vein, Brody and Moshinsky published extensive tables of transformation brackets [START_REF] Brody | Tables of Transformation Brackets for Nuclear Shell-Model Calculations[END_REF].

At the end of the fiftees, Smirnov worked out a parent problem, viz. the calculation of the Talmi coefficients for unequal mass nucleons, and he gave a solution for the case n 1 = 0 and n 2 = 0 [START_REF] Yu | [END_REF]. (Indeed, the transformation brackets and the Talmi coefficients are connected via a double Clebsch-Gordan transformation.) The coefficients nℓ, NL, λ|n 1 ℓ 1 , n 2 ℓ 2 , λ , called transformation brackets by Moshinsky and total Talmi coefficients by Smirnov, are now referred to as Moshinsky-Smirnov coefficients. In this respect, both the Moshinsky-Smirnov coefficients and the Talmi coefficients were revisited at the end of the seventies in terms of generating functions in the framework of the approaches by J.S. Schwinger and by V.

Bargmann of the harmonic-oscillator bases [6]. (The work by M. Hage Hassan [6], who prepared his State Doctorate thesis at the Institut de Physique Nucléaire de Lyon, constitutes a very deep and original approach to the Talmi coefficients and Moshinsky-Smirnov coefficients.) It should be noted that the transformation brackets or Moshinsky-Smirnov coefficients are also of importance for atoms and molecules as shown by Marcos and Yurii in their book [START_REF] Moshinsky | The Harmonic Oscillator in Modern Physics[END_REF] written during the time Yurii was an invited professor at the Instituto de Física of the Universidad Nacional Autónoma de México.

The harmonic oscillator is a central ingredient in numerous works by Smirnov and Moshinsky. Many applications of the nonrelativistic and relativistic harmonic oscillators to modern physics (from molecules, atoms, and nuclei to quarks) were pedagogically exposed in the book by Marcos and Yurii [START_REF] Moshinsky | The Harmonic Oscillator in Modern Physics[END_REF] with a special attention paid to the many-body problem (in the Hartree-Fock approximation and from the point of view of unitary symmetry), the nuclear collective motion, and the interacting boson model. Their common interests were not limited to applications based on the harmonic oscillator. Let us simply mention that both of them were interested in group theoretical methods and symmetry methods in physics and also contributed to several fields of mathematical physics as, for instance, state labelling problem, special functions, generating functions, and nonbijective canonical transformations.

We shall continue in Section IV with some contributions of Yurii mainly in nuclear physics with some excursions in some other domains.

IV. CONTRIBUTIONS IN NUCLEAR PHYSICS AND RELATED AREAS

A. Nuclear structure models Yu.F. Smirnov and some of his collaborators contributed significantly to the theory of clustering in nuclei (see the monographs [START_REF] Neudatchin | Nucleon Clusters in Light Nuclei[END_REF][START_REF] Nemets | sky, Nucleon Clusters in Atomic Nuclei and Many-Nucleon Transfer Reactions[END_REF]). His researches in this domain started, as seen above, with the derivation of transformation brackets [START_REF] Yu | [END_REF]. Then, it was necessary to elaborate a formalism for deriving many-nucleon fractional parentage coefficients (FPCs) including one, two, and sometimes three shells. This problem was solved by Yu.F. Smirnov and co-workers via a factorization of each FPC into a space part and a spin-isospin one [START_REF] Kukulin | [END_REF].

The coincidence of the spin-isospin part of the FPCs with the isoscalar factor part of the Clebsh-Gordan coefficients of the group SU(4) was demonstrated. Calculation technique of cluster characteristics of oscillator shell-model systems (spectroscopic amplitudes, factors, etc.) as a whole was first developed in [10,11]. For 1p-shell nuclei, a complete determination of these characteristics was presented in [12].

For lightest nuclei Yu.F. Smirnov with co-workers built the translationally-invariant shell model which is free of nonphysical oscillations of the center of mass of a nucleus [13].

Cluster model was proposed by Yu.F. Smirnov with co-authors for the description of electromagnetic form-factors of lightest nuclei [14].

The formalism of FPCs for many-boson systems, as nuclei described by the interacting boson model in the SU(3) limit, was built by Yu.F. Smirnov and co-authors in [15].

In connection with shell model calculations, Yu.F. Smirnov and K.V. Shitikova [START_REF] Yu | Physics of Elementary Particles and Atomic Nuclei[END_REF] gave a significant contribution to the hyperspherical function method. In the same direction, he also developed some group theoretical aspects of the generalized hyperspherical function method [START_REF] Filippov | Microskopic Theory of Collective Excitations of Atomic Nuclei[END_REF]. As a further important contribution in the field of mathematical physics, Smirnov developed, in collaboration with R.M. Asherova and V.N. Tolstoy, the method of extremal projection operators for simple Lie groups [START_REF] Asherova | [END_REF] and for semisimple complex Lie algebras [19] (this powerful and universal method was generalized by V.N. Tolstoy to all finite-dimensional contragredient superalgebras, to affine Kac-Moody algebras and superalgebras, and finally to the q-deformations of these algebras and superalgebras). Yu.F. Smirnov and collaborators applied the extremal projection operator method to several groups [20]- [27]; in particular, he elaborated a formalism for the Clebsch-Gordan coefficients and the 3nj-symbols of SU (3) and of the q-deformations of SU(2) and SU(3).

B. Theory of nuclear reactions

The theory of cluster quasi-elastic knock-out nuclear reactions was developed by Yu.F. Smirnov and co-authors in [28]. The role of cluster de-excitation in this process was demonstrated in [29].

Yu.F. Smirnov and Yu.M. Tchuvil'sky pointed out by making use of group-theory methods a general property of heavy projectile -heavy target nucleus channel, namely, systems exhibiting great distinctions between the wavefunction of the channel and the wavefunction of the respective compound. As a result, the concept of the co-called structural forbidness of heavy fragmentation (and fusion) was introduced into the theory of heavy ion collisions [30].

Calculation methods concerning the theory of alpha decay of giant isoscalar quadrupole resonances of nuclei were presented in the paper [31].

C. Quark degrees of freedom of nuclear systems Independently of V.A. Matveev and P. Sorba [32] -the authors of the concept of hidden color of the six-quark system -Yu.F. Smirnov and Yu.M. Tchuvil'sky demonstrated that if the deuteron would consist of two three-quark nucleons, then its wavefunction should contain components with hidden nucleon isobars (∆, etc.) [33]. The method of calculation of the statistical weights of these components in the deuteron wavefunction was worked out in [33].

Yu.F. Smirnov was an active participant of the team who developed methods of calculation relative to the dynamics of the six-quark system [34]. (The paper [34] contains a proof of the connection between the isoscalar factors of the Clebsch-Gordan coefficients of the unitary group SU(n) and the Clebsch-Gordan coefficients of the companion permutation group S n .) This team also emphasized the role of color-magnetic attraction and as a consequence of quark configuration mixing in this system [35].

D. Composite particle interaction

Yu.F. Smirnov and co-authors proposed an approach for treating the interaction between composite particles through the use of deep attractive potentials with redundant states [36].

He was the principal researcher of the group which made the first calculations of the ground state of the three-body (3α) system with an interaction of such a type [37].

E. The J-matrix formalism

The J-matrix formalism, a formalism utilizing L 2 bases in quantum scattering theory, was suggested in the mid-seventies by E.J. Heller, H.A. Yamani, W.P. Reinhardt, and L.

Fishman, people coming from the atomic physics community. The same formalism utilizing the oscillator basis was independently rediscovered in 1980 by G.F. Filippov, a good friend of Yurii. Then, Yu.F. Smirnov become very interested in this new quantum scattering formalism and performed a lot of works by developing and using it in various applications. In particular, (with Yu.I. Nechaev) he developed exact mathematical grounds of the J-matrix approach with the help of an oscillator basis [38], studied (together with J.M. Bang, A.I. Mazur, A.M. Shirokov, and S.A. Zaytsev) its relations with the well-known R-matrix and P -matrix formalisms in scattering theory [39], extended (together with A.M. Shirokov and S.A. Zaitsev) the J-matrix formalism to the case of the so-called true few-body scattering [40] and to the case of relativistic systems described by the Dirac equation. The first applications of the J-matrix formalism (in the oscillator basis) to the nuclear shell model was presented in a series of papers, published by Yu.F. Smirnov with A.I. Mazur and V.A. Knyr, devoted to the hypernucleus production reactions. Another application of this formalism, an application to systems decaying via three-body channels, was published by Yu.F. Smirnov with T.Ya. Mikhelashvili and A.M. Shirokov [41]. Later Yu.F. Smirnov with Yu.A. Lurie and A.M. Shirokov studied in detail, by means of this approach, three-body decays of loosely-bound 11 Li and 6 He nuclei within cluster models. As further interesting extensions of this approach, the three-body J-matrix formalism was used to calculate the three-body S-matrix and a direct numerical calculation of the S-matrix poles made it possible to improve essentially variational calculation results for the binding energies, ground state rms radii, electromagnetic transitions, etc. It is worth noting that the first J-matrix calculation of three-body S-matrix poles was performed by Yurii Smirnov together with A.M. Shirokov and L.Ya. Stotland in the studies of atomic He and H -ion. A review of some of the results obtained using the J-matrix approach by Yurii Smirnov and collaborators can be found in [START_REF] Yu | Harmonic oscillator representation in the theory of scattering and nuclear reactions[END_REF].

The J-matrix formalism is based on the solutions of three-term recurrence relations. This brought naturally Yurii Smirnov to study general properties of eigenenergies and eigenvectors for three-term recurrence relations (equivalent to second-order finite difference equations).

Together with P.A. Braun, A.M. Shirokov, and S.K. Suslov he developed an approximate method for analyzing the spectrum structure and eigenvector properties of three-term recurrence relations by replacing them by second-order ordinary differential equations of the Schrödinger type, a replacement which makes it possible to develop quantum intuition. This approach was used by Yu.F. Smirnov and collaborators in various problems, as for instance in the analysis of general properties of Clebsch-Gordan and Racah coefficients for SU [START_REF] Talmi | [END_REF] and SU(1, 1) groups [START_REF] Yu | [END_REF] as well as in the study of level clustering in the high-J spectra of non-rigid spherical top molecules [44] and of partly-filled shell ions in crystalline fields [45].

Interesting results were also obtained by Yurii Smirnov and collaborators in the study of exact solutions of three-term recurrence relations. In particular, together with A.M. Shirokov and N.A. Smirnova, he discovered a parameter symmetry of the interaction boson model [46] which was studied by him later (with O. Costaños, A. Frank, A.M. Shirokov, and N.A. Smirnova) in more complicated versions of the interacting boson model and in some other algebraic models.

F. The (e,2e) method

In collaboration with V.G. Neudatchin, in 1967-1969, Yu.F. Smirnov proposed a new method for the experimental investigation of the electronic structure of atoms, molecules, and solids. This method, based on the measurement of the process of quasi-elastic knock-out of electron by high-energy (in atomic scale) electron with coincidences registration, is called the (e,2e) method. It is the natural analog of the (p,2p) method popular in nuclear physics.

The proposed method was successfully applied in many laboratories around the world and provided data on the structure of many-electron systems (see for example [47] for a review about this major contribution). [START_REF] Sviridov | Spektroskopiya Kristallov (Akad. Nauk SSSR[END_REF]. At that time, the references for mathematical aspects of crystaland ligand-field theory were based on works by Y. Tanabe, S. Sugano, and H. Kamimura from Japan [START_REF] Tanabe | [END_REF], J.S. Griffith from England [START_REF] Griffith | The Theory of Transition-Metal Ions[END_REF], and Tang Au-chin and his collaborators from China [START_REF] Au-Chin | Theoretical Method of the Ligand Field Theory[END_REF] (see also some contributions by the present author [START_REF] Kibler | [END_REF]). The two abovementioned books by Smirnov and his colleagues shed some new light on the mathematical analysis of spectroscopic and magnetic properties of partly filled shell ions in molecular and crystal surroundings. In particular, special emphasis was put on the derivation of the Wigner-Racah algebra of a finite group of molecular and crystallographic interest from that of the group SO(3) ∼ SU(2)/Z 2 .

My second (indirect) contact with Yurii is related to an invitation to participate in the fifth workshop on Symmetry Methods in Physics organized by Yurii F. Smirnov and Raya M.

Asherova in Obninsk in July 1991. Unfortunately, I did not get my visa on time so that my participation was reduced to a paper in the proceedings of the workshop edited by Smirnov and Asherova [START_REF] Kibler | Symmetry adaptation and two-photon spectroscopy of ions in molecular or solid-state finite symmetry[END_REF].

In the beginning of the 1990's, I had a chance to get in touch with another facet of Yurii's work. In 1989, a Russian speaking student from Switzerland, C. Campigotto, spent one year in the group of Prof. Smirnov. He started working on the so-called Kustaanheimo-Stiefel transformation, an R 4 → R 3 transformation associated with the Hopf fibration S 3 → S 2 with compact fiber S 1 . Such a nonbijective transformation makes it possible to connect the Kepler-Coulomb system in R 3 to the isotropic harmonic oscillator in R 4 . (More precisely, the Kustaanheimo-Stiefel transformation allows to pass from a four-dimensional harmonic oscillator subjected to a constrain to the three-dimensional hydrogen atom, see for instance [START_REF] Kibler | [END_REF]). Then, Campigotto (well-prepared by Smirnov and his team, especially A.M. Shirokov and V.N. Tolstoy) came to Lyon to prepare a French doctorate thesis [56]. He defended his thesis in 1993 with G.S. Pogosyan (representing Yu.F. Smirnov) as a member of the jury.

A fourth opportunity to get involved with Yurii came from our mutual interest in quantum groups and in nuclear and atomic spectroscopy. I meet him for the first time in Dubna in 1992. We then started a collaboration (partly with R.M. Asherova) on qand qp-boson calculus in the framework of Hopf algebras associated with the Lie algebras su(2) and su(1, 1) [START_REF] Yu | Some aspects of q-boson calculus[END_REF]. In addition, we pursued a group theoretical study of the Coulomb energy averaged over the nℓ N -atomic states with a definite spin [START_REF] Kibler | [END_REF]. We also had fruitful exchanges in nuclear physics. Indeed, Prof. Smirnov and his colleagues D. Bonatsos (from Greece), S.B. Drenska, P.P. Raychev and R.P. Roussev (all from Bulgaria) developed a model based on a oneparameter deformation of SU(2) for dealing with rotational bands of deformed nuclei and rotational spectra of molecules [59] (see also [60]). Along the same line, a student of mine, R. Barbier, developed in his thesis a two-parameter deformation of SU(2) with application to superdeformed nuclei in mass region A ∼ 130 -150 and A ∼ 190 [61]. It was a real pleasure to receive Yurii in Lyon on the occasion of the defence of the Barbier thesis in 1995. Indeed, from 1992 to 1995, Yurii made four stays in Lyon (one with his wife Rita and one with his daughter Tatyana) and we jointly participated in several meetings, one in Clausthal in Germany (organised by H.-D. Doebner, V.K. Dobrev, and A.G. Ushveridze) and two in Bregenz in Austria (organised by B. Gruber and M. Ramek).

VI. CLOSING

Yurii had many students (some of them are now famous physicists), many collaborators in his country and abroad, and had an influence on many scientists. He was also an exceptional teacher. It was very pleasant, profitable and inspiring to be taught by Prof. Smirnov and/or to discuss with him. I personally greatly benefited from discussions with Yurii Smirnov.

Yurii Fedorovich Smirnov will remain an example for many of us. We will remember the exceptional qualities of the man as a scientist, as a teacher and as a generous person. Yurii, we shall not forget you.

  II. YURII FEDOROVICH SMIRNOV Yurii Fedorovich Smirnov was a Russian physicist. He was born in the city of Il'inskoe (Yaroslavl' region, Russia) in 1935. He graduated from Moscow State University. Then, he completed his Doctorate thesis at the same university under the guidance of Yurii M. Shirokov and benefited from fruitful contacts with other distinguished physicists as for instance Yakov A. Smorodinsky. He achieved his career in the (Skobeltsyn) Institute of Nuclear Physics and in the Physics Department of (Lomonosov) Moscow State University with many stays abroad. The last fifteen years of his life were shared between Moscow and Mexico city where he was a visiting Professor and got a Professor position successively at the Instituto de Física and, then, at the Instituto de Ciencias Nucleares of the U.N.A.M. (he spent almost 11 years in Mexico). He received prestigious awards: the K.D. Sinel'nikov Prize of the Ukrainian Academy of Sciences in 1982 and the M.V. Lomonosov Prize in 2002. He was a member of the Academy of Sciences of Mexico.

V

  . SOME PERSONAL REMINISCENCES My first contact with the work of Yurii Smirnov goes back to 1978 when my colleague J. Patera showed me, on the occasion of a NATO Advanced Study Institute organised in Canada by J.C. Donini, a beautiful book written by D.T. Sviridov and Yu.F. Smirnov [48]. This book dealt with the spectroscopy of d N ions in inhomogeneous electric fields (a part of a disciplinary domain known as crystal-and ligand-field theory in condensed matter physics and explored via the theory of level splitting from a theoretical point of view). In 1979, B.I. Zhilinskiȋ, while visiting Dijon and Lyon in France in the framework of an exchange programme between USSR and France, provided me with another interesting book, dealing with f N ions in crystalline fields, written by D.T. Sviridov, Yu.F. Smirnov and V.N. Tolstoy

Acknowledgments Part of these notes were presented at the 13th International Conference on Symmetry Methods in Physics (SYMPHYS-XIII) organized in memory of Prof. Yurii Fedorovich Smirnov by the Bogoliubov Laboratory of Theoretical Physics of the Joint Institute for Nuclear Research and the International Center for Advanced Studies at Yerevan State University (the conference was held in Dubna, Russia, 6-9 July 2009). I acknowledge the Organizing Committee of SYMPHYS-XIII, especially George S. Pogosyan, for their kind invitation to participate to this interesting conference. The part not presented at SYMPHYS-XIII were written on the basis of notes communicated by Vladimir G. Neudatchin, Andrey M. Shirokov, Yurii M. Tchuvil'sky, and Valeriy N. Tolstoy. I am very much indebted to them for their invaluable help in the preparation of this paper. Thanks are also due to Corrado Campigotto, Mehdi Hage Hassan, and Rutwig Campoamor-Stursberg for useful correspondence.


Pierre-Antoine Adragna
email: pierre-antoine.adragna@univ-savoie.fr

Maurice Pillet
email: maurice.pillet@univ-savoie.fr

Fabien Formosa
email: fabien.formosa@univ-savoie.frserge.samper@univ-savoie.fr

Serge Samper

Inertial tolerancing and capability indices in an assembly production

Keywords: assemblage, tolérancement, inertiel, indice, capabilité assembly, inertial, tolerancing, capability, index

Le tolérancement traditionnel considère la conformité d'un lot lorsque le lot de production satisfait les spécifications. La caractéristique est considérée pour elle-même mais pas en fonction de son incidence dans l'assemblage dont elle fait partie. Le tolérancement inertiel propose une autre alternative de tolérancement afin de garantir l'assemblage final. L'inertie I= δ 2 + σ 2 n'est pas tolérancée par un intervalle de tolérance mais par un scalaire représentant l'inertie maximale que la caractéristique ne doit pas dépasser. Nous détaillerons comment calculer ces tolérances inertielles selon deux cas de figures, si on vise à respecter une inertie ou un indice de capabilité Cpk sur la caractéristique résultante dans le cas particulier de tolérances uniformes et le cas général de tolérances non uniformes. Un exemple servira à comparer les différentes méthodes de tolérancement.

ABSTRACT. Traditional tolerancing considers the conformity of a batch when the batch satisfies the specifications. The characteristic is considered for itself and not according to its incidence in the assembly. Inertial tolerancing proposes another alternative of tolerancing in order to guarantee the final assembly characteristic. The inertia I = σ 2 + δ 2 is not toleranced by a tolerance interval but by a scalar representing the maximum inertia that the characteristic should not exceed. We detail how to calculate the inertial tolerances according to two cases, one aims to guarantee an inertia of the assembly characteristic the other a tolerance interval on the assembly characteristic by a Cpk capability index, in the particular but common case of uniform tolerances or more general with non uniform tolerances. An example will be detailed to show the results of the different tolerancing methods.

Introduction

The traditional tolerancing considers the conformity of a batch when the batch satisfies the specifications. The characteristic is considered for itself and not regarding its incidence on the final assembly resultant. It has been showed that inertial tolerancing (I = σ 2 + δ 2 which is no more based on a [Min Max] interval but on the Taguchi loss function) proposes another tolerancing method to guarantee the final assembly while allowing larger variability in the case of centered production.

This paper proposes a method to calculate the inertial tolerances of the components of a 1D mechanical assembly chain in two cases:

-we want to guarantee an inertial tolerance on the final assembly, -we want to guarantee a minimum of the Cpk index on the tolerance interval [Min; Max] of the assembly characteristic.

Different cases are considered, even the general case where the components tolerances are not uniformly distributed. A comparison with the traditional tolerancing will show the difference on the allowed variability of components. The two considered cases have different hypothesis of application. We will discuss on the choice of using the first or the second approach.

An industrial case of application will be treated as an example.

Inertial tolerancing, capability indices and conformity

The aim of tolerancing is to determine an acceptation criterion on the components characteristics x i to guarantee the quality of the assembly resultant Y. In the case of a good design, when the x characteristic is produced on the target, the quality is optimal. As x gets an offset from the target, the function of the assembly will be more sensitive to the conditions of use and the environment, and can lead to a non-satisfaction of the customer. Taguchi demonstrated that the financial loss associated to an offset from the target is proportional to the square of this offcentering L = K(T-X)² . (Pillet et al, 2001) shows that in the case of a batch, the financial loss associated is L = K(σ² + δ²). Then he defines

I x = ∑ i=1 n ( ) x -Target 2 n = σ x 2 + δ x 2 [1]
This function is called Inertia by analogy to the mechanical inertia. Here I x represents the inertia of the x characteristic, σ x is the standard deviation of the batch distribution and δ x is the offset of the mean to the target.

To qualify the capability of a process with the inertial tolerancing, two capability indices have been defined:

Cp = I Max σ Batch [2]
which indicates the capability of a centered process.

Cpi = I Max I Batch [3]
which indicates the capability considering the process off-centering.

Compared to the traditional tolerancing, the proposed approach of the inertial tolerancing is quite different. The aim is no more to guarantee a rate of parts out of tolerance, but to guarantee the centering of components around the target in order to guarantee the quality of the assembly. The reflection is no more based on the proportions out of tolerances but on the inertias of the components, the normality of the batch distribution is no more a necessary criterion.

Tolerancing of an assembly

For the tolerancing of assembly systems, the problem consists of finding the elementary characteristics x i of the components in order to obtain a final characteristic Y satisfying the functional condition of the assembled product for the customers needs. As a general rule, it is possible to approximate the system behavior around the target by a linearization at the first order. The final characteristic behavior can be expressed by the following relation:

Y = α 0 + ∑ i=1 n α i .x i [4]
Where α i is the influence coefficient of the component i on the assembly resultant Y, α 0 is the target value of Y and n is the number of components in the assembly. For the computation of the components tolerances, the general case will be considered where tolerances are not uniformly distributed with the use a difficulty coefficient β i ≥ 1, also called feasibility coefficient. The simplified but common case will also be presented, where the tolerances are uniformly distributed β i = 1, and the incidence coefficients are all equal α i = 1 .

Review of the traditional tolerancing methods

Before the presentation of the inertial tolerancing, here is a brief reminds of the traditional tolerancing methods of assembly systems. Three commonly used traditional methods are presented.

Worst of cases tolerancing

In this case, one considers that the final characteristic of the assembly will be respected in any cases of assembly. In the general case where tolerances are nonuniformly distributed, the β i coefficients will be used for the components. The tolerance expression of a components is R xi = β i . R x , the assembly resultant is then:

R Y = ∑ i=1 n ( ) | | α i . β i . R x [5] R xi = β i . R Y ∑ i=1 n ( ) | | α i . β i [6]
Where R Y represents the tolerance interval of the functional condition of the assembly, and R xi is the tolerance interval of the components. Tolerances can be distributed regarding different methods by changing the β i coefficients [START_REF] Graves | Tolerance Analysis Taylored to your organization[END_REF]:

-Uniform distribution of the tolerance, -Considering the tolerances of standard components, or conception rules, -Proportional to the square of the nominal length, -Considering the process capabilities, In the case of uniform distribution of the tolerance β i = 1 , and same incidences of the components α i = 1 , the relation [5] becomes then:

R Y = ∑ i=1 n R xi = n . R xi [7] R xi = R Y n [8]
The well known inconvenient of this method is the high price of production due to the tightened tolerances of the components.

Statistical tolerancing

The statistical tolerancing has been developed to consider the low probability of having several characteristics in limit of their tolerances simultaneously (Shewhart, 1931) [START_REF] Evans | Statistical tolerancing: the state of the art[END_REF]. Under the hypothesis that the x i variables are independent with a standard deviation R xi , the equation [4] gives the following relation:

σ Y = ∑ i=1 n (
)

α i 2 . σ xi 2 [9]
In the general case of non uniform distribution, one consider R xi = β i . σ xi , then:

R Y = ∑ i=1 n ( ) α i 2 . β i 2 . σ xi 2 [10] R xi = β i . R Y ∑ i=1 n ( ) α i 2 . β i 2 [11]
With tolerances proportional to the standard deviation [START_REF] Chase | A survey of research in the application of tolerance analysis to the design of mechanical assemblies[END_REF], in the case of a uniform repartition, and same incidence of the components α i = 1 , the relation [10] becomes then:

R Y = ∑ i=1 n R xi 2 = n . R xi 2 [12] R xi = R Y n [13]
In this tolerancing method, the basic hypothesis is to consider the centering of all characteristics on their target values. The inconvenient of this method is that it does not guarantee the conformity of the assembly characteristic in all configurations of the components. It can be possible that the components respect their tolerance intervals, but their off-centering from the target lead to non-conformity on the final condition in its tolerance interval.

Inflated statistical tolerancing:

Several methods are proposed in order to reduce the negative aspect of the statistical tolerancing. A proposed method is the inflated statistical tolerancing [START_REF] Graves | How to reduce costs using a tolerance analysis formula tailored to your organization[END_REF] [START_REF] Graves | Tolerance Analysis Taylored to your organization[END_REF]. This method consists of using the inflated coefficient in the tolerancing of the components based on the statistical method.

R Y = f . ∑ i=1 n ( ) α i 2 . β i 2 . R xi 2 [14] R xi = β i . R Y f . ∑ i=1 n ( ) α i 2 . β i 2 [15]
Where f represents the inflated coefficient generally chosen around f = 1,5 to 1,6. In the case of a uniform distribution of the tolerances and same incidences of the components, one has the following relation:

R Y = f . ∑ i=1 n R xi 2 [16] R xi = R Y f . n [17]
In the case where f = 1, one finds the statistical tolerancing method. In the case where f = n , one finds the worst of cases tolerancing method. A discussion (Graves et al, 2000) on different situations leads to choose different f values. Graves proposes an interesting approach to optimize the f coefficient following different capability indices. Although this is an improved method, it is possible to find a situation where the final assembly characteristic will not well be respected.

Inertial tolerancing of an assembly guarantying an inertia on the resultant

The following results come from [START_REF] Pillet | Inertial Tolerancing in the case of assembled products[END_REF]. These results are reminded in order to be used further. The calculations are done under the hypothesis of a uniform distribution of the components tolerances.

σ Y = ∑ i=1 n α i 2 . σ xi 2 [18] δ Y = ∑ i=1 n α i . δ xi [19]
The inertia of the resultant characteristic Y is defined by the relation [1]. Replacing relations [18] and [19], one have the inertia of the resultant in function of the components inertias and off-centering:

I Y = ∑ i=1 n α i 2 . I xi 2 + 2 . ∑ i=1 n α i . α j . δ xi . δ xj [20]
The first part of the equation corresponds to the addition of the squared inertias. The double product corresponds to the case where all off-centering of the components are on the same side. A discussion is necessary to treat different hypothesis.

Hypothesis 1: Worst of cases

This hypothesis of components in their worst of cases considers that the component inertia is only due to its off-centering from the target,

I xi = δ xi 2 = | | δ xi ,
the relation [20] becomes then:

I Y = ∑ i=1 n α i 2 . δ xi 2 + 2 . ∑ i=1 n α i . α j . δ xi . δ xj [21]
In the case where all the α i = 1, and the tolerances are uniformly distributed I xi = δ xi , one have the following relations:

I Y = n 2 . I x 2 = n . I x [22] I x = I Y n [23]
In the general case where the α i coefficients are non-equal and the tolerances are not uniformly distributed I xi = β i . δ x , one has the following relation:

I xi = β i . I Y ∑ i=1 n | | α i . β i [24]

Hypothesis 2: random distribution of the averages

This hypothesis is close to the consideration made for the traditional statistical tolerancing. In this case, the double product ∑ i=1 n α i . α j . δ xi . δ xj is null, the relation

[20] becomes then:

I Y = ∑ i=1 n α i 2 . I xi 2 [25]
In the case where α i = 1, one thus have:

I Y = n . I x 2
[26]

I x = I Y n [27]
In the general case of non-uniform distribution of tolerances, one considers I xi = β i .σ xi . From relation [25] one obtains:

I Y = ∑ i=1 n (
)

α i 2 . β i 2 . σ xi 2 [28]
thus the components inertia can be

I xi = β i . I Y ∑ i=1 n ( ) α i 2 . β i 2 [29] 3.2.3. Hypothesis 3: off-centering of δ = k . σ of all components.
This hypothesis considers that all components have systematic off-centering equals to δ = k . σ . In this case the component inertia is:

I x = σ x 2 + (k . σ x ) 2 = σ x . 1 + k 2 [30] then σ x = I x 1 + k 2 [31]
and

δ x = I x . 1 - 1 1 + k 2 [32]
Equation [20] becomes then:

I Y = ∑ i=1 n α i 2 . I xi 2 + 2 . ∑ i=1 n α i . α j . ⎝ ⎛ ⎠ ⎞ 1 - 1 1 + k 2 . I xi . I xj [33]
In the case where all α i = 1, one obtains then:

I Y = n.I x 2 + n.( n -1 ). ⎝ ⎛ ⎠ ⎞ 1 - 1 1 + k 2 .I x 2 [34] thus I x = I Y n. ⎝ ⎛ ⎠ ⎞ n.k 2 + 1 1 + k 2 [35]

Hypothesis 4: off-centering of m components out of n

In this hypothesis, the designer determines the number of characteristics that can have a systematic off-centering. Under these conditions, the double product is reduced. In the case where all α i = 1 and tolerances are uniformly distributed, one have:

I Y = n.I x 2 + m.(m -1). ⎝ ⎛ ⎠ ⎞ 1 - 1 1 + k 2 .I x 2 [36] thus I x = I Y .( ) 1 + k 2 n.( ) 1 + k 2 + m.k 2 .(m -1)
[37]

Inertial tolerancing guarantying a Cpk index on the assembly resultant: the corrected inertial tolerancing

In the case where the aim is to guarantee a Cpk index on the assembly resultant, the conformity of this final characteristic is considered regarding a tolerance interval. In the most common case where the target is centered in the tolerance interval Target = LT + UT 2 where LT is the lower tolerance and UT is the upper tolerance and R Y = UT -LT 2 is the length of the tolerance interval. The capability index Cpk is defined by:

Cpk = Min ⎝ ⎜ ⎛ ⎠ ⎟ ⎞ R Y 2 -( ) Target -x 3.σ Y , R Y 2 -( ) x -Target 3.σ Y [38]
Which is equivalent to

Cpk = R Y 2 -δ Y 3.σ Y [39]
Replacing with the equations [18] and [19], one has the following relation:

Cpk = R Y 2 -∑ i=1 n α i . δ xi 3 . ∑ i=1 n α i 2 . σ xi 2 [40]
In order to define the components inertia in function of the tolerance interval to guarantee on the assembly resultant, one considers that the components are centered. In this case, the components inertia is only due to their standard deviation:

I xi = σ xi 2 + 0 = σ xi .
The resultant is also centered, from where I Y = σ Y , one have then:

Cp = R Y 6 . σ Y [41] σ Y = I Y = R Y 6 . Cp [42]
In the general case of non-uniform distribution, one considers I xi = β i . σ xi . In order to differentiate the Cp index to a prediction index, it will be renamed as the ICC coefficient (Inertial Corrected Coefficient). As the hypothesis of centered production is considered, one has I Y = σ Y . From the relation [29] one obtains the following tolerancing relation:

I xi = β i . R Y 6 . ICC . ∑ i=1 n ( ) α i 2 . β i 2 [43]
In the case of a uniform repartition of the tolerances and all α i = 1, one obtains then:

I xi = R Y 6 . ICC . n [44]

1. Variations of the Cpk index in function of the components off-centering

Now it is possible to define the components inertia in function of the tolerance interval of the assembly resultant, let us interest in the variations of the Cpk index while components are in limit of their inertial tolerances. This hypothesis allows expressing the standard deviation in function of the off-centering, one has the expression of the batch inertia given by [1], I xi = σ xi 2 + δ xi 2 thus the components deviation can be expressed as follows:

σ xi = I xi 2 -δ xi 2 [45]
Which is replaced in [39]. Let us also consider that all δ xi = δ , one obtains then:

Cpk = R Y 2 -n.| | δ 3 . n.( ) I xi 2 -δ 2 [46]
The component inertias are defined by [44], then:

Cpk = R Y 2 -n.| | δ 3 . R Y 2 36 . ICC 2 -n.δ 2 [47]
The Cpk function is pair, the domain of study can then be reduced to δ > 0 , this allows to suppress the absolute value function on δ. The derivate function of Cpk at δ is:

∂Cpk ∂δ = ⎝ ⎛ ⎠ ⎞ R Y 2 -n. δ .n.δ 3 . ⎝ ⎛ ⎠ ⎞ R Y 2 36 . ICC 2 -n.δ 2 3/2 - n 3 . R Y 2 36 . ICC 2 -n.δ 2 [48]
The variations of the Cpk function are studied thanks to the sign of its derivate function continuously defined on [0; I xi [ The function has a minimum for:

δ Min = R Y 18.ICC 2 [49]
And the value of Cpk is:

Cpk Min = ICC 2 - n 9 [50]
These results are interesting because they allow to find the minimum of the ICC coefficient value to guarantee a minimum of the Cpk index in function of the number of components in the assembly. Thus it is also possible to know how many components can compose the assembly in order to guarantee a minimum of the resultant Cpk index with a given value of the ICC coefficient:

ICC = Cpk 2 + n 9 [51] n = 9.( ) ICC 2 -Cpk 2 [52]
NOTE. -The results found considering the hypothesis that all component offcentering are equal δ xi = δ correspond to those found with the gradient method. It consists of calculating the Δf gradient of the Cpk function depending on n variable δ xi and its Hf Hessian matrix, and then to find the point for which the gradient is null and the eigen values of the Hessian matrix are all strictly positive in order to have a local minimum. In our case, the local minimum of the Cpk is found for all δ xi = R Y 18.ICC 2 . A demonstration by recurrence on n and i proves it but will not be exposed in this paper.

2. Variations of the Cpk function regarding the number of components and their off-centering

We will see the influence of the off-centering and the number of components on the Cpk index value. Figure 1 confirms the previous results on the Cpk variations. An assembly system which components tolerancing has been done with ICC = 1 or 1,5 has a Cpk which is minimum when all δ xi = 0,056 or δ xi = 0,025 respectively.

Figure1. Variations of Cpk indices with different configurations of the number of components and ICC coefficient values in function of the components off-centering in limit of their inertial tolerances.

The Cpk minimum can be calculated for n < 9.ICC 2 , from relation [50]. This remark can also be observed on the following table where ICC = 1 and 1,5, and n = 1 to 12. For an ICC index chosen to ICC = 1, the minimum value of the Cpk index can be calculated up to n = 9, for n = 12 it is not possible to evaluate the Cpk minimum. With a limit study, it can be showed that this minimum tends towardinfinity.

Cpk, R Y = 1 n = 3 n = 6 n = 9 n = 12 ICC = 1 0,816 0,577 0,000 -infinity ICC = 1,5 1,384 1,258 1,118 0,957

Table 2. Different values of the minimum of the Cpk index for different number of components n and different ICC coefficient values.

A link can be observed between the evolution of the minimum of the Cpk index and the ICC coefficient value. With an ICC coefficient ICC = 1, it is impossible to guarantee a Cpk index, Cpk > 1 due to this worst configuration of all δ xi , even for n = 1 when δ ≠ 0. For an ICC coefficient ICC = 1,5, it is possible to guarantee a Cpk index Cpk > 1,1 up to n = 9 components in the assembly. It is nearly possible to guarantee Cpk = 1 for 12 components.

Let us see the variations of the ICC coefficient value for the calculation of the inertial tolerances in function of the Cpk index value to guarantee, and the number of components in the assembly.

3. ICC values to guarantee a minimum of the Cpk index in function of the number of components.

From relation [51], one can choose the ICC coefficient value in order to guarantee a minimum of the Cpk index of the assembly in its tolerance interval.

Pp index garantying a given

Validation of the abacus in the case of the generalized tolerancing

Previous results are obtained in the case of a particular inertial tolerancing with uniform distribution of the tolerances and all α i = 1. What are the results in the general case of a non-uniformly distributed tolerancing and all incidences are not equal? In the general case of non-uniform repartition of the tolerances, the inertia is defined by the relation [43] Replacing [43] and [45] in the equation of the Cpk defined by [40], we obtain the following expression of the Cpk:

Cpk = R Y 2 -∑ i=1 n α i . | | δ xi 3 . R Y 2 36 . ICC 2 -∑ i=1 n α i 2 . δ xi 2 [53]
One can feel that with a variable change δ = α i . δ xi , the minimum of the Cpk can be found for the general case for δ xi = R Y 18 . ICC 2 . α i , in this case Cpk = ICC 2 -n 9 as in the particular case. One can conclude that the abacus presented in figure 2 is also applicable in the case of generalized inertial tolerancing.

An application case

The example will be used to compare the results of the different tolerancing methods.

The target value of the gap is 1mm with a tolerance interval of 1mm, thus +/-0,5mm. The tolerance distribution is uniform. The following relation defines the assembly relation:

Y = X 1 -X 2 -X 3 -X 4 -X 5
[54]

Application of the different tolerancing methods

The inertial tolerancing of the two-presented cases (guarantying an inertia or an interval) will be compared to the traditional tolerancing methods. In the case of application, the tolerance interval R Y = 1, the number of components is n = 5. In the case of an inertial tolerancing guarantying a capability index Cpk = 1 on the assembly resultant, the ICC coefficient to take into account for the inertial tolerances calculation is:

ICC = Cpk 2 + n 9 = 1,25 [55] X 1 X 2 X 3 X 4 X 5
Mini on the condition (LT): 0.5 Maxi on the condition (UT) : 1.5

Traditional tolerancing

Inertial tolerancing

Worst of case: The maximum batch dispersion in the case of the traditional tolerancing is useful to compare the allowed dispersion on the components with the two methods. For inertial tolerancing, when the batch is centered, the maximum batch dispersion equals to the inertial tolerance I xi . It is then possible to compare the maximum dispersion allowed for centered components with the different tolerancing methods.

R xi = R Y n R xi = 0,

Discussion on the case of application

One can see that the maximum authorized dispersions for the traditional tolerancing by worst of cases or statistical are similar to those of the inertial tolerancing guarantying an inertia by worst of cases σ xiMax = I xi = 0,033, or statistical σ xiMax = I xi = 0,075 respectively. But the maximum allowed dispersion by the inflated statistical tolerancing is σ xiMax = 0,05, which is lower than the maximum authorized dispersion by the corrected inertial tolerancing guarantying a minimum of the Cpk index Cpk = 1 which is I xi = 0,06.

Before concluding that this last method guarantying a Cpk index gives larger dispersion than the inflated statistical tolerancing, let us compare the influence of the number of components on this dispersion difference. First, the inflated coefficient f = 1,5 can be compared to the ICC coefficient ICC = 1,25 that is lower. Choosing a ICC coefficient ICC = 1,5 is the same as guarantying a Cpk index Cpk = 1,30 for an assembly of 5 components, from relation [52], but this also guarantee a Cpk index Cpk = 1 for an assembly composed of up to 11 components from relation [50], 11,25 = n =9.( ) ICC 2 -Cpk 2 . It is possible to conclude that for an assembly of less than 12 components, the inertial tolerancing guarantying a Cpk index Cpk = 1 is larger than the inflated statistical tolerancing. Over 12 components, the inflated statistical tolerancing gives larger dispersions than the inertial tolerancing guarantying a Cpk index, but the inertial tolerancing allows to guarantee a Cpk index on the final assembly characteristic, that cannot do the inflated statistical tolerancing.

The comparison between traditional statistical tolerancing and inertial tolerancing thanks to the authorized dispersion is dangerous; this comparison can be done only for centered batches.

Conclusion

In this paper has been proposed a corrected inertial tolerancing method, which aims to guarantee the conformity of the assembly resultant in a tolerance interval thanks to a Cpk index. The ICC coefficient value to take into account for the tolerances calculation is found regarding the number of components in the assembly, and the minimum value of the Cpk index to guarantee on the assembly resultant.

The application of this tolerancing method allows to apply the inertial tolerancing method but with a link to the traditional tolerancing. That is to say that this method allows to guarantee a Cpk index on the assembly resultant while applying the inertial tolerancing based on a statistical approach.

Figure 2 .

 2 Figure 2. Value of the ICC coefficient to take into account in order to guarantee different values of the Cpk index in function of the number of tolerance.

Figure 2

 2 Figure2is an abacus of the inertial tolerancing guarantying a Cpk index on the assembly resultant in the case of a uniform distribution of the tolerances and same incidences of the components α i = 1 . Two relations have to be considered while dealing with the particular but common case of a uniform distribution of the tolerances for an assembly where the components have the same incidences, the relations [44], I xi = R Y 6 . ICC . n , and [51], ICC = Cpk 2 + n 9 .

.

  One obtains the following table of variations.

	δ		0	δ =	R Y 18.ICC 2	I xi =	R Y 6.ICC. n
	Cpk		ICC		ICC 2 -	9 n	+ Infinity
	∂Cpk ∂δ	-	R Y 2.n.ICC 2	_	0	+	+ Infinity

Table 1 .

 1 Variations of the assembly Cpk function depending on the component offcentering δ for a given number of components toleranced by the corrected inertial tolerancing.

Table 2 .

 2 Comparison of different traditional and inertial tolerancing methods. For the traditional tolerancing, the tolerance interval of the components is R

xi , and in the case of a centered batch, the maximum batch dispersion is σ xiMax = R xi 6 . For the inertial tolerancing, the components inertias are expressed by σ xiMax = I xi .


Dima Batenkov
email: dima.batenkov@weizmann.ac.il

Niv Sarig
email: niv.sarig@weizmann.ac.il

Yosef Yomdin
email: yosef.yomdin@weizmann.ac.il

An "algebraic" reconstruction of piecewise-smooth functions from integral measurements

Introduction

This paper presents some results on a well-known problem in Algebraic Signal Sampling and in other areas of applied mathematics: reconstruction of piecewise-smooth functions from their integral measurements (like moments, Fourier coefficients, Radon transform, etc.). Our results concern reconstruction (from the moments) of signals in two specific classes: linear combinations of shifts of a given function, and "piecewise D-finite functions" which satisfy on each continuity interval a linear differential equation with polynomial coefficients. Let us start with some general remarks and a conjecture. It is well known that the error in the best approximation of a C k -function f by an N -th degree Fourier polynomial is of order C N k . The same holds for algebraic polynomial approximation and for other basic approximation tools. However, for f with singularities, in particular, with discontinuities, the error is much larger: its order is only C √ N . Considering the so-called Kolmogorow N -width of families of signals with moving discontinuities one can show that any linear approximation method provides the same order of error, if we do not fix a priori the discontinuities' position (see [START_REF] Ettinger | Linear versus non-linear acqusition of step-functions[END_REF], Theorem 2.10). Another manifestation of the same problem is the "Gibbs effect" -a relatively strong oscillation of the approximating function near the discontinuities. Practically important signals usually do have discontinuities, so the above feature of linear representation methods presents a serious problem in signal reconstruction. In particular, it visibly appears near the edges of images compressed by JPEG, as well as in the noise and low resolution of the CT and MRI images. Recent non-linear reconstruction methods, in particular, Compressed Sensing ( [START_REF] Candes | Compressive sampling[END_REF][START_REF] Donoho | Compressed sensing[END_REF]) and Algebraic Sampling ( [START_REF] Dragotti | Sampling Moments and Reconstructing Signals of Finite Rate of Innovation: Shannon Meets Strang-Fix[END_REF][START_REF] Maravic | Exact Sampling Results for Some Classes of Parametric Non-Bandlimited 2-D Signals[END_REF][START_REF] Prandoni | Approximation and compression of piecewise smooth functions[END_REF][START_REF] Elad | Shape from moments-an estimation theory perspective[END_REF][START_REF] Gustafsson | Reconstructing planar domains from their moments[END_REF]), address this problem in many cases. Both approaches appeal to an a priori information on the character of the signals to be reconstructed, assuming their "simplicity" in one or another sense. Compressed sensing assumes only a sparse representation in a certain (wavelets) basis, and thus it presents a rather general and "universal" approach. Algebraic Sampling usually requires more specific a priori assumptions on the structure of the signals, but it promises a better reconstruction accuracy. In fact, we believe that ultimately the Algebraic Sampling approach has a potential to reconstruct "simple signals with singularities" as good as smooth ones. In par-ticular, the results of [START_REF] Eckhoff | Accurate reconstructions of functions of finite regularity from truncated Fourier series expansions[END_REF][START_REF] Kvernadze | Locating discontinuities of a bounded function by the partial sums of its Fourier series[END_REF][START_REF] Gelb | Detection of edges in spectral data II. Nonlinear enhancement[END_REF][START_REF] Tadmor | High resolution methods for time dependent problems with piecewise smooth solutions[END_REF][START_REF] Prandoni | Approximation and compression of piecewise smooth functions[END_REF] strongly support (also apparently do not accurately formulate and prove) the following conjecture: There is a non-linear algebraic procedure reconstructing any signal in a class of piecewise C k -functions (of one or several variables) from its first N Fourier coefficients, with the overall accuracy of order C N k . This includes the discontinuities' positions, as well as the smooth pieces over the continuity domains. At present there are many approaches available to a robust detection of discontinuities from Fourier data (see [START_REF] Gelb | Detection of edges in spectral data II. Nonlinear enhancement[END_REF][START_REF] Eckhoff | Accurate reconstructions of functions of finite regularity from truncated Fourier series expansions[END_REF][START_REF] Kvernadze | Locating discontinuities of a bounded function by the partial sums of its Fourier series[END_REF] and references therein). The remaining problem seems to be an accurate estimate of the accuracy of the solution of the nonlinear systems arising. Our results below can be considered, in particular, as a step in this direction. On the other hand, they have been motivated by the results in [START_REF] Dragotti | Sampling Moments and Reconstructing Signals of Finite Rate of Innovation: Shannon Meets Strang-Fix[END_REF][START_REF] Maravic | Exact Sampling Results for Some Classes of Parametric Non-Bandlimited 2-D Signals[END_REF][START_REF] Prandoni | Approximation and compression of piecewise smooth functions[END_REF], and in [START_REF] Gustafsson | Reconstructing planar domains from their moments[END_REF][START_REF] Elad | Shape from moments-an estimation theory perspective[END_REF].

Linear combinations of shifts of a given function

Reconstruction of this class of signals from sampling has been described in [START_REF] Dragotti | Sampling Moments and Reconstructing Signals of Finite Rate of Innovation: Shannon Meets Strang-Fix[END_REF][START_REF] Maravic | Exact Sampling Results for Some Classes of Parametric Non-Bandlimited 2-D Signals[END_REF]. We study a rather similar problem of reconstruction from the moments. Our method is based on the following approach: we construct convolution kernels dual to the monomials. Applying these kernels, we get a Prony-type system of equations on the shifts and amplitudes.

Let us restate a general reconstruction problem, as it appears in our specific setting. We want to reconstruct signals of the form

F (x) = N i=1 j,l a i,j,l f (l) i (x + x j ) (1)
where the f i 's are known functions of x = (x 1 , . . . , x d ), and the form (1) of the signal is known a priori. The parameters a i,j,l , x j = (x j 1 , . . . , x j d ) are to be found from a finite number of "measurements", i.e. of linear (usually integral) functionals like polynomial moments, Fourier moments, shifted kernels, evaluation over some grid and more. In this paper we consider only linear combinations of shifts of one known function f (although the method of "convolution dual" can be extended to several shifted functions and their derivatives -see [START_REF] Sarig | Reconstruction of "Simple" Signals from Integral Measurements[END_REF]). First we consider general integral "measurements" and then restrict ourselves to the moments and Fourier coefficients. In what follows

x = (x 1 , . . . , x d ), t = (t 1 , . . . , t d ), j is a scalar index, while k = (k 1 , . . . , k d ), i = (i 1 , . . . , i d ) and n = (n 1 , . . . , n d ) are multi-indices. Partial ordering of multi- indices is given by k ≤ k ⇔ k p ≤ k p , p = 1, . . . , d. So we have F (x) = s j=1 a j f (x + x j ). (2)
Let the measurements μ k (F ) be given by μ k (F ) = F (t)ϕ k (t)dt, for a certain (multi)-sequence of functions ϕ k (t), k ≥ 0 = (0, . . . , 0). Given f and ϕ = {ϕ k (t)}, k ≥ 0 we now try to find certain "triangular" linear combinations

ψ k (t) = 0≤i≤k C i,k ϕ i (t) (3)
forming, in a sense, some "f -convolution dual" functions (similar to a bi-orthogonal set of function) with respect to the system ϕ k (t). More accurately, we require that

f (t + x)ψ k (t) = ϕ k (x). (4)
We shall call a sequence ψ = {ψ k (t)} satisfying (3), ( 4) f -convolution dual to ϕ. Below we find convolution dual systems to the usual and exponential monomials. We consider a general problem of finding convolution dual sequences to a given sequence of measurements as an important step in the reconstruction problem. Notice that it can be generalized by dropping the requirement of a specific representation (3):

ψ k (t) = k i=0 C i,k ϕ i (t).
Instead we can require only that f (t)ψ k (t) be expressible in terms of the measurements sequence μ k . Also ϕ k in (4) can be replaced by another a priori chosen sequence η k . This problem leads, in particular, to certain functional equations, satisfied by polynomials and exponents (as well as exponential polynomials and some kinds of elliptic functions). Now we have the following result:

Theorem 1. Let a sequence ψ = ψ k (t) be f -convolution dual to ϕ. Define M k by M k = 0≤i≤k C i,k μ i .
Then the parameters a j and x j in (2) satisfy the following system of equations ("generalized Prony system"):

s j=1 a j ϕ k (x j ) = M k , k = 0, . . . . (5)
Proof We have

M k = 0≤i≤k C i,k μ i = F (t) 0≤i≤k C i,k ϕ i (t)dt = F (t)ψ k (t) = s j=1 a j f (t + x j )ψ k (t)dt = s j=1 a j ϕ k (x j ).
In specific examples we can find the minimal number of equations in (5) necessary to uniquely reconstruct the parameters a j and x j in (2).

Reconstruction from moments

We are given a finite number of moments of a signal F as in [START_REF] Candes | Compressive sampling[END_REF] in the form

m n = F (t)t n dt. (6)
So here

ϕ n (x) = x n 1 1 • • • x n d d for each multi-index n = (n 1 , . . . , n d ).
We look for the dual functions ψ n satisfying the convolution equation

f (t + x)ψ n (t)dt = x n (7)
for each multi-index n. To solve this equation we apply Fourier transform to both sides of [START_REF] Ettinger | Linear versus non-linear acqusition of step-functions[END_REF]. Assuming that f (ω) ∈ C ∞ (R d ), f (0) = 0 we find (see [START_REF] Sarig | Reconstruction of "Simple" Signals from Integral Measurements[END_REF]) that there is a unique solution to [START_REF] Ettinger | Linear versus non-linear acqusition of step-functions[END_REF] provided by

ϕ n (x) = k≤n C n,k x k , (8)
where

C n,k = 1 ( √ 2π) d n k (-i) n+k ∂ n-k ∂ω n-k ω=0 1 f (ω) .
This calculation is symbolic and works for more general cases. The actual calculation in our polynomial case is done using straightforward matrix calculations. We set the generalized polynomial moments as

M n = k≤n C n,k m k (9)
and obtain, as in Theorem 1, the following system of equations:

s j=1 a j (x j ) n = M n , n ≥ 0. (10)
This system can be solved explicitly in a standard way (see, for example, [START_REF] Nikishin | Rational Approximations and Orthogonality[END_REF][START_REF] Dragotti | Sampling Moments and Reconstructing Signals of Finite Rate of Innovation: Shannon Meets Strang-Fix[END_REF][START_REF] Sarig | Signal Acquisition from Measurements via Non-Linear Models[END_REF]). In one-dimensional case it goes as follows (see [START_REF] Nikishin | Rational Approximations and Orthogonality[END_REF]): from [START_REF] Kisunko | Cauchy type integrals and a D-moment problem[END_REF] we get that for z = (z 1 , . . . , z d ) the generalized moments generating function (d = 1 yet, notice that the formulas are still multidimensional)

I(z) = n∈N d M n z n = s j=1 a j d l=1 1 1 -x j l z l (11)
is a rational function. Hence its Taylor coefficients satisfy linear recurrence relation, which can be reconstructed through a linear system with the Hankel-type matrix formed by an appropriate number of the moments M n 's. This is, essentially, a procedure of the diagonal Padé approximation for I(z) (see [START_REF] Nikishin | Rational Approximations and Orthogonality[END_REF]). The parameters a j , x j are finally reconstructed as the poles and the residues of I(z).

For several variables, although the formulas are the same as above, the generalization of the solution of the Prony system is more involved and should be addressed separately.

In one dimensional case with the derivatives f (l) included we have

F (x) = s j=1 r l=0 a j,l f (l) (x + x j ). (12)
The corresponding moment-generating function in this case takes the form

I(z) = s j=1 r l=0 l q=0 l q (-1) q+l a j,l /(x j ) l (1 -x j z) q+1 . ( 13
)
which is still a rational function (d-dimensional case with derivatives is similar). We would like to stress that in this case the dual polynomials ψ k are not changed and they are given as in [START_REF] Gelb | Detection of edges in spectral data II. Nonlinear enhancement[END_REF]. Therefore also the formula for the generalized moments M n is the same as in (9).

Fourier case

In the same manner as in section 2.1 we now choose ϕ k (x) = e ikx . We get immediately ψ k (x) = 1 f (k) e -ikx . Indeed,

f (t + x)ψ k (t)dt = f (t + x) 1 f (k) e ikt dt = f (k) f (k) e -ikx = ϕ -k (x). ( 14
)
Here the triangular system of equations ( 3) is actually not triangular any more but still since

ψ k (x) = 1 f (k) ϕ -k (x)
we can express the generalized moments through the original ones via

M k = 1 f (k) μ -k [F ]
. Now exactly as before we can find a generalized Prony system in the form

1 f (k) μ -k [F ] = M k = j a j e -ikx j = j a j ρ k j ( 15
)
where ρ j = e -ixj . In this case we get a rational exponential generating function and we can find its poles and residues on the unit complex circle as we did in the polynomial case.

Further extensions

The approach above can be extended in the following directions: 1. Reconstruction of signals built from several functions or with the addition of dilations also can be investigated (a perturbation approach where the dilations are approximately 1 is studied in [START_REF] Sarig | Signal Acquisition from Measurements via Non-Linear Models[END_REF]). 2. Further study of "convolution duality" can significantly extend the class of signals and measurements allowing for a closed -form signal reconstruction.

Reconstruction of piecewise D-finite functions from moments

Let g : [a, b] → R consist of K+1 "pieces" g 0 , . . . g K with K ≥ 0 jump points a = ξ 0 < ξ 1 . . . < ξ K < ξ K+1 = b
Furthermore, let g satisfy on each continuity interval some linear homogeneous differential equation with polynomial coefficients: D g n ≡ 0, n = 0, . . . , K where

D = N j=0 kj i=0 a i,j x i d j dx j (a ij ∈ R) (16)
Each g n may be therefore written as a linear combination of functions {u i } N i=1 which are a basis for the space N D = {f : D f ≡ 0}:

g n (x) = N i=1 α i,n u i (x), n = 0, 1, . . . , K (17)
We term such functions g "piecewise D-finite". Many real-world signals may be represented as piecewise Dfinite functions, in particular: polynomials, trigonometric functions, algebraic functions. The sequence {m k = m k (g)} is given by the usual moments

m k (g) = b a x k g(x)dx
We subsequently formulate the following Piecewise D-finite Reconstruction Problem. Given N, {k i }, K, a, b and the moment sequence {m k } of a piecewise D-finite function g, reconstruct all the parameters {a i,j }, {ξ i }, {α i,n }.

Below we state some results (see [START_REF] Batenkov | Moment inversion problem for piecewise D-finite functions[END_REF] for detailed proofs) which provide explicit algebraic connections between the above parameters and the measurements {m k }.

The first two theorems assume a single continuity interval (compare [START_REF] Kisunko | Cauchy type integrals and a D-moment problem[END_REF]).

Theorem 2. Let K = 0 and D g ≡ 0 with D given by [START_REF] Sarig | Reconstruction of "Simple" Signals from Integral Measurements[END_REF].

Then the moment sequence {m k (g)} satisfies a linear recurrence relation

(E -a I) N (E -b I) N • N j=0 k j i=0 a i,j Π (i,j) (k, E) m k = 0 ( 18
)
where E is the discrete forward shift operator and Π (i,j) (k, E) are monomials in E whose coefficients are polynomials in k: Π (i,j) (k, E) = (-1) j (i+k)! (i+k-j)! E i-j . Theorem 3. Denote

E(E) def = (E -a I) N (E -b I) N , v (i,j) k def = E(E) • Π (i,j) (k, E) m k , h j (z) def = ∞ k=0 v (0,j) k z k , G j (x) def = E(x) d j dx j g(x)
Assume the conditions of Theorem 2. Then [START_REF] Batenkov | Moment inversion problem for piecewise D-finite functions[END_REF] The vector of the coefficients a = (a i,j ) satisfies a linear homogeneous system

Ha =       v (0,0) 0 v (1,0) 0 . . . v (kN ,N ) 0 v (0,0) 1 v (1,0) 1 . . . v (k N ,N ) 1 . . . . . . . . . . . . v (0,0) M v (1,0) M . . . v (k N ,N ) M            a 0,0 a 1,0 . . . a k N ,N      = 0 (19) for all M ∈ N. (2) v (i,j) k = m i+k (G j (x))
. Consequently, h j (z) is the moment generating function of G j (x).

(3) Denote p j (x) = kj i=0 a i,j x i . Then the functions Φ = {1, h 0 (z), . . . h N (z)} are polynomially dependent:

N j=0 h j (z) z max kj p j (z -1 ) = Q(z) where Q(z) is a polynomial with deg Q < max k j .
The system of polynomials {z max kj p j (z -1 )} is called the Padé-Hermite form for Φ.

To handle the piecewise case, we represent the jump discontinuities by the step function H(x) def = 0 x < 0 1 x ≥ 0 and write g as a distribution

g(x) = g 0 + K n=1 g n (x)H(x -ξ n ) (20
α N,K         =      m 0 m 1 . . . m M      (22)
The above results can be combined as follows to provide a solution of the Reconstruction Problem: The constants M and M determine the minimal required size of the corresponding linear systems (19) and ( 22) in order for all the solutions of these systems to be also solutions of the original problem. It can be shown that:

1. There exists no uniform bound on M without any additional information on the nature of the solutions. Explicit bounds may be obtained for simple function classes such as piecewise polynomials of bounded degrees or real algebraic functions.

) Theorem 4 .Proposition 5 .

 45 Let K > 0 and let g be as in (20) with operator D annihilating every piece g n . Then the operatorD def = K n=1 (x -ξ i ) N I • D (21)annihilates the entire g as a distribution. Consequently, conclusions of Theorems 2 and 3 hold with D replaced by D as in (21). Let K ≥ 0 and {u i } N i=1 be a basis for the space N D , where D annihilates every piece of g. Assume[START_REF] Tadmor | High resolution methods for time dependent problems with piecewise smooth solutions[END_REF] and denote c n i,k = ξ n+1 ξnx k u i (x) for n = 0, . . . , K. A straightforward computation gives ∀ M ∈ N:

  (a) Let N, {k i }, K, a, b and {m k (g)} be given. If K > 0, replace D (still unknown) with D according to (21). (b) Build the matrix H as in (19). Solve Ha = 0 and obtain the operator D * = D a which annihilates g. (c) If K > 0, factor out all the common roots of the polynomial coefficients of D * with multiplicity N . These are the locations of the jump points {ξ n }. The remaining part is the operator D † which annihilates every g n . (d) By now D † and {ξ n } are known. So compute the basis for N D † and solve (22).

For every specific D, an explicit bound M = M (D) may be computed for the system (22).The above algorithm has been tested on exact reconstruction of piecewise polynomials, piecewise sinusoids and rational functions.


Ayush Bhandari
email: ayushbhandari@ntu.edu.sg

Pina Marziliano

Sampling of Sparse Signals in Fractional Fourier Domain

In this paper, we formulate the problem of sampling sparse signals in fractional Fourier domain. The fractional Fourier transform (FrFT) can be seen as a generalization of the classical Fourier transform. Extension of Shannon's sampling theorem to the class of signals which are fractional bandlimited shows its association to a Nyquist-like bound. Thus proving that signals that have a non-bandlimited representation in FrFT domain cannot be sampled. We prove that under suitable conditions, it is possible to sample sparse (in time) signals by using the Finite Rate of Innovation (FRI) signal model. In particular, we propose a uniform sampling and reconstruction procedure for a periodic stream of Diracs, which have a nonbandlimited representation in FrFT domain. This generalizes the FRI sampling and reconstruction scheme in the Fourier domain to the FrFT domain.

Introduction

Shannon's sampling theorem [START_REF] Shannon | Communications in the presence of noise[END_REF] provides access to the digital world. Our understanding of this sampling theorem together with the reconstruction formula is solely based on the frequency content of the signal of interest. This is where the indispensable Fourier transform comes into the picture. Almeida [START_REF] Almeida | The fractional Fourier transform and time-frequency representations[END_REF] introduced the fractional Fourier transform or the FrFT-a generalization of the Fourier transform-to the signal processing community in 1994. The generalization of the Fourier transform by FrFT has several interesting consequences from the signal processing perspective. For instance, non-bandlimited signals in the Fourier domain can still have a compactly supported representation in FrFT domain [START_REF] Xia | On bandlimited signals with fractional Fourier transform[END_REF], when dealing with non stationary distortions, the FrFT based filters can perform better than Fourier domain based filters (in sense of mean square error) [START_REF] Kutay | Optimal filtering in fractional Fourier domains[END_REF] etc. To give the reader an idea about the growing popularity of FrFT, it would be worth mentioning that on at least eight occasions including, [START_REF] Xia | On bandlimited signals with fractional Fourier transform[END_REF][START_REF] Zayed | On the relationship between the Fourier and fractional Fourier transforms[END_REF][START_REF] Erseghe | Unified fractional Fourier transform and sampling theorem[END_REF][START_REF] Zayed | New sampling formulae for the fractional Fourier transform[END_REF][START_REF] García | Orthogonal sampling formulas: A unified approach[END_REF][START_REF] Candan | Sampling and series expansion theorems for fractional Fourier and other transforms[END_REF][START_REF] Torres | Sampling theorem for fractional bandlimited signals: A self-contained proof. application to digital holography[END_REF][START_REF] Tao | Sampling and sampling rate conversion of band limited signals in the fractional Fourier transform domain[END_REF]], Shannon's sampling theorem [START_REF] Shannon | Communications in the presence of noise[END_REF][START_REF] Unser | Sampling-50 years after Shannon[END_REF] was independently extended to the class of fractional bandlimited signals. In [START_REF] Ozaktas | Introduction to the fractional Fourier transform and its applications[END_REF], the FrFT of a signal or a function, say x(t), is defined by

x θ (ω) = FrFT{x(t)} = x(t)K θ (t, ω)dt (1)
where

K θ (t, ω) def =      1-j cot θ 2π e j t 2 +ω 2 2 cot θ-jωt csc θ , θ = pπ δ(t -ω), θ = 2pπ δ(t + ω), θ + π = 2pπ (2)
is the transformation kernel, parametrized by the fractional order θ ∈ R and p is some integer. The FrFT of a time-frequency representation e.g. Gabor Transform results in rotation of the plane by the fractional order of the FrFT [START_REF] Almeida | The fractional Fourier transform and time-frequency representations[END_REF]. Thus, we denote fractional order by θ and from now on, we will use fractional order and angle interchangeably. The inverse-FrFT with respect to angle θ is the FrFT at angle -θ, given by,

x(t) = ∞ -∞ x θ (ω) K -θ (t, ω) dω. (3)
Whenever θ = π/2, (1) collapses to the classical Fourier transform definition. A direct consequence of the generalization of the Fourier transform by the FrFT results in a modification in the idea of bandlimitedness. Its impact is visible in the change that manifests in Shannon's sampling theorem for fractional bandlimited signals [START_REF] Tao | Sampling and sampling rate conversion of band limited signals in the fractional Fourier transform domain[END_REF], which is stated in Theorem 1.

Theorem 1 (Shannon-FrFT). Let x(t) be a continuoustime signal. If the spectrum of x(t), i.e. x θ (ω) is fractional bandlimited to ω m which means, x θ (ω) = 0, when |ω| > ω m , then x(t) is completely determined by giving its ordinates at a series of equidistant points spaced T = π ωm sin θ seconds apart. This theorem has an equivalence to the Shannon's sampling theorem for θ = π/2. The reconstruction formula for fractional bandlimited signals is given in [START_REF] Tao | Sampling and sampling rate conversion of band limited signals in the fractional Fourier transform domain[END_REF],

x(t) = λ * θ (t) n∈Z λ θ (nT ) x(nT )sinc ((t -nT ) ω m csc θ) (4) where λ θ (•) def = e j(•) 2 cot θ 2
is a domain independent chirp modulation function and the '*' in the superscript denotes complex conjugation. If x (t) is the approximation of x(t), then x (t)x(t) 2 = 0 when ω m ωs 2 sin θ-the Nyquist rate for FrFT-where ω s = 2π/T is the sampling frequency. Note that all the aforementioned results are equivalent to Shannon's sampling theorem with respect to Fourier domain for θ = π/2. Theorem 1 (for FrFT) has a striking similarity with the Shannon's sampling theorem (for FT), in that, sampling non-bandlimited signals is impossible. Consider Dirac's delta function or δ(t). Using (2), we have,

δ θ (ω) = FrFT {δ (t)} = 1-j cot θ 2π λ θ (ω) (5)
which is a non-bandlimited function (and least sparse when compared to the time-domain counterpart) and thus, Theorem 1 fails to answer the following question: If x(t) is a fractional non-bandlimited signal, then, how can we sample and reconstruct such a signal? To make this statement clear, we introduce the fractional convolution operator, which is denoted by ' * θ * θ * θ '. Accordingly, filtering x(t) by a filter, h(t), in 'fractional sense' 1 is equivalent to [START_REF] Kraniauskas | Method for defining a class of fractional operations[END_REF],

x(t) * θ h(t) = 1-j cot θ 2π λ * θ (t)•([x(t)λ θ (t)] * [h(t)λ θ (t)]) (6)
where ' * ' denotes the usual convolution operator. In light of this definition, we wish to address the problem of recovering parsimonious x(t) from the samples of its filtered version, i.e., y(nT

) = x(t) * θ h(t)| t=nT , n ∈ Z.
This problem has a natural/strong link with that of sparse sampling [START_REF] Marziliano | Sampling innovations[END_REF]16,[START_REF] Blu | Sparse sampling of signal innovations[END_REF]. The Heisenberg-Gabor uncertainty principle for the FrFT [START_REF] Shinde | An uncertainty principle for real signals in the fractional Fourier transform domain[END_REF] (a generalization of the Fourier duality) asserts that the product of spreads of x θ (ω) and x(t) has a lower bound which is proportional to sin 2 θ 4 (assuming that x = 1). This implies that sparsity in one domain will lead to loss of compact support in canonically conjugate domain. Our contribution in this article is to propose a sampling and reconstruction scheme for signals which have a sparse representation in time domain and whose fractional spectrum is non-bandlimited. We model our sparse signal as a continuous periodic stream of Diracs which is being observed by an acquisition device which deploys a sincbased filter. The paper is organized as follows: We assume that the reader is familiar with basic ideas outlined in [START_REF] Unser | Sampling-50 years after Shannon[END_REF]16,[START_REF] Blu | Sparse sampling of signal innovations[END_REF]. In Section II, we introduce our sparse signal model and the definition of the fractional Fourier series (FrFS). Using these as preliminaries, in Section III, we derive an equivalent representation of our signal in FrFT domain. In Section IV, we discuss the sampling theorem and its completeness and Section V is the conclusion.

Preliminaries

Sparse Signal Model

We model our sparse signal as a periodic stream of K Diracs, i.e.

x(t) = K-1 k=0 c k n∈Z δ(t -t k -nτ ) (7)
1 We adhere to this modified definition of convolution operator as it inherits the fractional Fourier duality property, in that,

FrFT {x(t) * θ h(t)} = λ * θ (ω) • x θ (ω) h θ (ω), which does not hold for the FrFT of x(t) * h(t) unless θ = π 2 .
with period τ , weights

{c k } K-1 k=0 arbitrary shifts, {t k } K-1 k=0 ⊂ [0, τ ).
In sense of [16], the signal has 2K degrees of freedom per period and the rate of innovation being ρ = 2K τ . From now on, the signal x(t) will denote the stream of Diracs.

Fractional Fourier Series (FrFS)

Periodic signals can be expanded in FrFT domain as a fractional Fourier series or FrFS [START_REF] Pei | Fractional Fourier series expansion for finite signals and dual extension to discrete-time fractional Fourier transform[END_REF]. The FrFS of a periodic signal, say x(t), can be written as,

x(t) = m∈Z x θ [m]Φ θ (m, t) (8)
where,

Φ * θ (m, t) = sin θ -j cos θ τ e j t 2 +(2πm sin θ/τ ) 2 2 cot θ-j2πmt/τ
constitutes the basis for FrFS expansion for a τ -periodic x(t). The FrFS coefficients are given by,

x θ [m] = τ x(t)Φ * θ (m, t)dt = x, Φ θ (m, •) (9)
where τ denotes the integral width and a, b = a(t)b * (t)dt denotes the inner product. The well-known Fourier series (FS) is just a special case of FrFS for θ = π 2 .

Stream of Diracs in Fractional Fourier Domain

In Fourier analysis, the Poisson summation formula (PSF) plays an important role. It is a well-known fact that a stream of Diracs (Dirac comb) in time-domain is another stream of Diracs in Fourier domain. In this subsection, we will derive the equivalent representation of Dirac comb in FrFT domain. This can be seen as a generalization of the Poisson summation formula for Dirac comb in FrFT domain.

Theorem 2. Let n∈Z δ(t -nτ ) be a Dirac comb, then n∈Z δ (t -nτ ) FrFT ←→ 1 τ 2π 1-j cot θ k∈Z δ θ [kω 0 sin θ] e -j t 2 + (kω 0 sin θ) 2 2 cot θ+jkω0t
where

ω 0 = 2π τ . Proof. Let s(t) def = n∈Z δ(t -nτ ).
The proof is done by expanding s(t) in FrFS basis or,

s(t) = k∈Z s, Φ θ s θ [k] Φ θ (k, t). (10)
The coefficients of this expansion are given by,

s θ [k] (9) = s, Φ θ (k, t) = κ (θ) √ τ t0+τ t0 s (t) Φ * θ (k, t) dt, ∀t 0 ∈ R = κ (θ) √ τ τ /2 -τ /2 δ (t) e j(t 2 +(kω0 sin θ) 2 /2) cot θ-jkω0t dt (since s(t + τ ) = s(t) and s(t) = δ (t) , t ∈ -τ 2 , τ 2 ) = κ (θ) √ τ e j((kω0 sin θ) 2 /2) cot θ (5) = κ (θ) √ τ 2π 1 -j cot θ δ θ [kω 0 sin θ] (11)
where κ (θ) = √ sin θj cos θ. Back substitution of ( 11) in ( 10) results in,

s(t) = 1 τ 2π 1-j cot θ × k∈Z δ θ [kω 0 sin θ] e -j t 2 + (kω 0 sin θ) 2 2 cot θ+jkω0t
.

This concludes the proof.

For sake of convenience, we will assume that the constant

1-j cot θ 2π
has been absorbed in τ . Note that at θ = π 2 , s(t) = 1 τ k∈Z e jkω0t which is the result of applying the PSF on s(t) in Fourier domain. Our immediate goal now is to derive the FrFS equivalent of x(t) in [START_REF] Zayed | New sampling formulae for the fractional Fourier transform[END_REF]. Since x(t) is a linear combination of some s(t) delayed by some time shift t k , it will be useful to recall shift property of FrFT [START_REF] Almeida | The fractional Fourier transform and time-frequency representations[END_REF] which states that,

FrFT {s (t -t k )} = s θ (ω -t k cos θ) e j 1 2 t 2 k sin θ cos θ-jωt k sin θ . (12)
Therefore, call

x(t) = K-1 k=0 c k • s k (t)
where s k (t) is the time-shifted version of s(t) with shift parameter t k . Using Theorem 2 and the shift-property of FrFT, we have,

s k (t) = n∈Z δ(t -t k -nτ ) (8) = m∈Z FrFT{ δ(t -t k )} | ω=mω0 sin θ Φ θ (m, t) (12) = 1 τ m∈Z e j cot θ 2 (t 2 k -t 2 )+jmω0(t-tk)
PSF for Dirac Comb in FrFT .

Having obtained the FrFT-version of s k (t), we can write,

x(t) = K-1 k=0 c k • n∈Z δ(t -t k -nτ ) = K-1 k=0 c k m∈Z e j cot θ 2 (t 2 k -t 2 )+jmω0(t-tk) = e -jt 2 cot θ 2 m∈Z 1 τ K-1 k=0 c k e j cot θ 2 (t 2 k )-jmω0tk p[m]
e j 2πm τ t .

Note that x(t) is non-bandlimited, however, it can be completely described by the knowledge of p[m] which in turn can be expanded as a linear combination of K complex exponentials.

Sampling and Reconstruction of Sparse Signals in Fractional Fourier Domain

We assume that a sinc-based kernel is used to prefilter x(t). In particular, we let the sampling kernel to be ϕ n (t) = e -j cot θ 2 t 2 sinc(t -nT ). Integer translates of ϕ n (t) form an orthonormal basis and the FrFT of ϕ(t)(= ϕ 0 (t)) is given by ϕ θ (ω) = with finite rate of innovation ρ = 2K τ . Let the sampling kernel/prefilter ϕ(t) be an ideal low-pass filter which has fractional bandwidth [-Bπ, Bπ], where B is chosen such that B ≥ ρ. If the filtered version of x(t), i.e. y(t) = x(t) * θ ϕ(-t) is sampled uniformly at locations t = nT , n = 0, . . . , N -1 then the samples,

1-j cot θ 2π e -j cot θ 2 ω 2 rect(ω/2π
y (nT ) = x(t) * θ ϕ (-t)| t=nT , n = 0, . . . , N -1,
y (nT ) = x(t) * θ ϕ (-t)| t=nT , n = 0, . . . , N -1 = λ * θ (nT ) τ m∈Z p[m] × e j 2πm τ t , (B csc θ) sinc ((B csc θ) (t -nT )) .
The inner product in the above step is further simplified using the Fourier integral,

e j 2πm τ t , (B csc θ) sinc ((B csc θ) (t -nT )) = rect( m Bτ csc θ )e j 2πm τ (nT ) .
We can therefore conclude that,

y (nT ) = λ * θ (nT ) τ m∈Z p[m] rect( m Bτ csc θ )e j 2πm τ (nT ) = λ * θ (nT ) τ M θ m=-M θ p[m]e j 2πm τ (
nT ) , n = 0, . . . , N -1 where

M θ = Bτ csc θ 2 .
Signal reconstruction from its samples: is equivalent to finding a corresponding vector A that forms a null space of a suitable submatrix of p[m] i.e. P (2M θ -K+1)×(K+1) -which is essentially the set Null(P) = A ∈ R K+1 : P • A = 0 . For details of this computation, the reader is referred to (cf. Pg. 1427, [16]). Figure 1 shows the layout of this algorithm.

Call p[m] = K-1 k=0 a k u m k -a linear combination of K-complex ex- ponentials, u k = λ * π/2 √ ω 0 t k with weights a k = c k • λ θ (t k ). The problem of calculating {a k } K-1 k=0 and {u k } K-1 k=0 is based on finding a suitable polyno- mial A(z) = K-1 k=0 1 -u k z -1

Conclusion

We presented a scheme for sampling and reconstruction of sparse signals in fractional Fourier domain. A direct consequence of modeling our signal of interest as a Finite Rate of Innovation signal, is that, the outcome bears an acute resemblance with the results previously derived, for the Fourier domain case. This simplifies the problem to the extent that reconstruction strategy remains unchanged and as we have shown, one can obtain the precise locations and amplitudes of the stream of Diracs using the annihilating filter method. Since time and frequency domains are special cases of the FrFT domain, it turns out that the number of values (M θ ) required for exact reconstruction of time domain signal depends on the chirp rate of transformation, i.e. θ.

2 .

 2 are a sufficient characterization of x(t), provided that N ≥ 2M θ + 1 and M θ = Bτ csc θ Proof. Using the following FrFT pair,1-j cot θ 2π λ * θ (ω) • rect( ω 2πB )FrFT ←→ (B csc θ) λ * θ (t) sinc (Bt csc θ) we define our sampling kernel as,ϕ B (t -nT ) = λ * θ (t) ϕ (B csc θ (t -nT ))which is compactly supported over [-Bπ, Bπ]. Prefiltering and sampling x(t) results in,

Figure 1 :

 1 Figure 1: Sampling and reconstruction of periodic stream of Diracs in FrFT domain.

  ). In light of the definition in[START_REF] Erseghe | Unified fractional Fourier transform and sampling theorem[END_REF], prefiltering the input signal x(t) with the kernel/low-pass filter ϕ(-t) and sampling can be written as, y (nT ) = x(t) * θ ϕ(-t)| t=nT . The main result is in the form of the following theorem.

	Theorem 3. Let x(t) be a τ -periodic stream of Diracs
	weighted by coefficients {c k }	K-1 k=0 and locations {t k }	K-1 k=0

  whose inverse ztransform yields the annihilating filter coefficients, A[m] which annihilate p[m]. In matrix notation, finding A[m]


Thierry Blu
email: thierry.blu@m4x.org

The Generalized Annihilation Property A Tool For Solving Finite Rate of Innovation Problems

We describe a property satisfied by a class of nonlinear systems of equations that are of the form F(Ω)X = Y.

Here F(Ω) is a matrix that depends on an unknown Kdimensional vector Ω, X is an unknown K-dimensional vector and Y is a vector of N ≥ K) given measurements. Such equations are encountered in superresolution or sparse signal recovery problems known as "Finite Rate of Innovation" signal reconstruction. We show how this property allows to solve explicitly for the unknowns Ω and X by a direct, non-iterative, algorithm that involves the resolution of two linear systems of equations and the extraction of the roots of a polynomial and give examples of problems where this type of solutions has been found useful.

Introduction

We consider the signal resulting from the convolution between a window ϕ(t) and the sum of K Diracs with amplitude x k located at time t k . Given the N uniform samples y n (T = sampling step )

y n = K k=1 x k ϕ(nT -t k )
where n = 1, 2, . . . , N, (

then FRI problems (see [START_REF] Vetterli | Sampling signals with finite rate of innovation[END_REF][START_REF] Blu | Sparse sampling of signal innovations[END_REF]) consist in retrieving the parameters t k and x k . Solving such problems is conceptually interesting because it shows how to break the standard Nyquist-Shannon bandlimitation rule for the exact reconstruction of signals from their uniform samples [START_REF] Shannon | A mathematical theory of communication[END_REF].

The system of consistent equations (1) can be expressed under the generic form of a nonlinear problem as shown in Fig. 1 (see next page), where the parameters Ω = [ω 1 , ω 2 , . . . ω K ] are related unambiguously to the unknowns t k 's. Because of the variety of settings adapted to this general approach, it happens to be necessary to distinguish between the parameters ω k -which we shall call "abstract parameters"-and the locations t k : typically, the ω k 's will be the zeros of some polynomial and from these ω k 's, we will be able to retrieve the t k 's using a functional relation of the form ω k = λ(t k ) for some invertible function λ(t).

At first sight, solving such a nonlinear system of equations is a daunting task. Fortunately, if the matrix F(Ω) satisfies a property that we shall call "Generalized Annihilation Property" (GAP), this reduces to solving two linear systems of equations sandwiching a nonlinear step that amounts to polynomial root extraction in practical cases.

The filters ϕ(t) that satisfy the GAP are thus especially interesting, since the related FRI problems enjoy a straight non-iterative solution.

The Generalized Annihilation Property (GAP)

We carry on with the previously identified general nonlinear problem, namely

F(Ω) X = Y, (3)
where the unknowns are Ω = [ω 1 , ω 2 , . . . ω K ] and X = [x 1 , x 2 , . . . x K ], and where the measurements are Y = [y 1 , y 2 , . . . y N ]. This system is said to satisfy the Generalized Annihilation Property whenever there exist K + 1 constant matrices, A k , and K + 1 scalar functions of Ω, h k (Ω), such that we have the identity

K k=0 h k (Ω) A k F(Ω) = 0. (4)
for any vector of parameters Ω. By right multiplying with X, the above equation implies that any solution Ω of ( 3) is also a solution of the (generalized) annihilation equation

K k=0 h k (Ω)A k Y = 0. (5)
This equation can be expressed in a matrix form AH = 0 where the unknown is

H = [ h 0 (Ω), h 1 (Ω), . . . , h K (Ω) ] T and the matrix A = A 0 Y, A 1 Y, . . . , A K Y .
Thus, in order to solve [START_REF] Shannon | A mathematical theory of communication[END_REF] for Ω and X, the idea consists in finding the scalar coefficients h k (Ω) that satisfy [START_REF] Stoica | Introduction to Spectral Analysis[END_REF], then retrieving ω 1 , ω 2 , . . . , ω K from the knowledge of h k (Ω), and finally finding X such that F(Ω) X = Y. Without elaborating on the conditions that make this solution unique, a

     ϕ(T -t 1 ) ϕ(T -t 2 ) • • • ϕ(T -t K ) ϕ(2T -t 1 ) ϕ(2T -t 2 ) • • • ϕ(2T -t K ) . . . . . . . . . ϕ(N T -t 1 ) ϕ(N T -t 2 ) • • • ϕ(N T -t K )      F(Ω)      x 1 x 2 . . . x K      X =      y 1 y 2 . . . y N      Y (2)
Figure 1: Algebraic equivalent of the consistency equations [START_REF] Vetterli | Sampling signals with finite rate of innovation[END_REF].

minimal requirement is that the matrices A k have at least K rows.

In the simple case where the h k (Ω)'s are related to the ω k 's through a polynomial relation

K k=0 h k (Ω)z -k = K k=1 (1 -ω k z -1 ), (6)
solving (3) boils down to a three-step algorithm that can be summarized as follows:

1. Compute a solution H = [ 1, h 1 , . . . , h K-1 , h K ] T of A 0 Y, A 1 Y, . . . , A K Y H = 0; 2. Compute the roots ω k of the z-transform H(z) = K k=0 h k z -k ; 3. Compute a solution X of F(Ω) X = Y.
Example-Spectral estimation problems boil down to a nonlinear problem of the form (3) involving the Vandermonde matrix:

F(Ω) =      ω 1 ω 2 • • • ω K ω 2 1 ω 2 2 • • • ω 2 K . . . . . . ω N 1 ω N 2 • • • ω N K     
where the frequencies to retrieve, f k , are related to ω k through ω k = e j2πf k . This problem satisfies the GAP for band-diagonal matrices A k which are more precisely given by:

A k = 0 N -K,k I N -K 0 N -K,K-k ,
where 0 m,n is the m × n zero matrix and I n is the n × n identity matrix. A minimal-yet not sufficient-condition for the unicity of the solution is N ≥ 2K. Since the A k can be seen as shifting operators by k samples, the annihilation equation is analogous to a filtering equation-with an annihilating filter. The annihilation algorithm is then equivalent to Prony's method [START_REF] Prony | Essai expérimental et analytique[END_REF]. Of course, spectral estimation in the presence of noise has been addressed by numerous researchers since the 1970's [5, 6, 7, 8, 9, 10, 11].

Some GAP Kernels

The GAP is actually shared by many interesting filters that can be used in sampling schemes, resulting in easily solvable FRI problems. Among them, the first ones to be identified were the periodized sinc, the infinite (i-e., not periodized) sinc and the Gaussian kernels [START_REF] Vetterli | Sampling signals with finite rate of innovation[END_REF]. Even more interestingly, recent research indicates that this property may somewhat be related to the Strang-Fix conditions which makes a very intriguing connection with approximation theory [START_REF] Dragotti | Sampling moments and reconstructing signals of finite rate of innovation: Shannon meets Strang-Fix[END_REF], and considerably broadens the class of FRI-admissible kernels. In all cases investigated so far, the scalar coefficients h k (Ω) satisfy (6).

Periodized sinc (Dirichlet) filter

Solving the FRI problem in the case of a periodic stream of Diracs is equivalent to considering [START_REF] Vetterli | Sampling signals with finite rate of innovation[END_REF] where ϕ is a periodized sinc kernel, e.g., a Dirichlet kernel

ϕ(t) = k ′ ∈Z sinc(B(t -k ′ τ )) = sin(πBt) Bτ sin(πt/τ )
where τ is the period of the Dirac stream and B some bandwith (chosen so that Bτ is an odd integer) [START_REF] Blu | Sparse sampling of signal innovations[END_REF]. This problem can be reformulated using the annihilation equation (4) by defining the following annihilation matrices

A k = 0 Bτ -K,k I Bτ -K 0 Bτ -K,Bτ -k W
where W = [e -j2πmn/N ] for |m| ≤ ⌊Bτ /2⌋ and 1 ≤ n ≤ N , is the N -DFT submatrix of size Bτ × N . Then, the abstract parameters ω k are related to the locations t k through ω k = e -j2πt k /τ . This kernel has been found useful for the estimation of UWB channels [START_REF] Maravić | Lowsampling rate UWB channel characterization and synchronization[END_REF] and for image superresolution [START_REF] Baboulaz | Exact feature extraction using finite rate of innovation principles with an application to image super-resolution[END_REF].

Infinite sinc filter

The filter ϕ(t) is given by ϕ(t) = sinc Bt with B = 1/T . When ϕ * x (t) is sampled uniformly at frequency B, the nonlinear system of equations satisfies the GAP. The abstract parameters ω k are related to the locations t k through ω k = t k and the annihilation matrices are given by

A k =        K K K K-1 • • • K 0 0 • • • 0 0 . . . . . . . . . . . . . . . . . . . . . . . . . . . 0 0 • • • • • • K K K K-1 • • • K 0        ×           1 0 • • • • • • 0 0 2 k . . . . . . . . . . . . 3 k . . . . . . . . . . . . . . . 0 0 • • • • • • 0 N k          

Gaussian filter

The filter ϕ(t) is given by ϕ(t) = exp(-t 2 /σ 2 ). When ϕ * x (t) is sampled uniformly at frequency T -1 , the nonlinear system of equations satisfies the GAP. The abstract parameters ω k are related to the locations t k through ω k = exp(2t k T /σ 2 ) and the annihilation matrices are given by

A k = 0 N -K,k I N -K 0 N -K,K-k ×         e T 2 σ 2 0 • • • • • • 0 0 e (2T ) 2 σ 2 . . . . . . . . . . . . . . . 0 0 • • • • • • 0 e (N T ) 2 σ 2        
A version of this solution (actually, for a Gabor kernel) was used in Optical Coherence Tomography, showing the possibility to resolve slices of a microscopic sample below the coherence length of the illuminating reference light [START_REF] Blu | A new highresolution processing method for the deconvolution of optical coherence tomography signals[END_REF].

Finite Support Strang-Fix filters

Through linear combinations of its shifts, the finite support filter ϕ(t) is assumed to reconstruct polynomials up to some degree L -1 (standard Strang-Fix condition [START_REF] Strang | A Fourier analysis of the finite element variational method[END_REF]) or exponentials e a l t where a l -a 0 is linear with l = 0, 1, . . . , L -1. More precisely, in the standard Strang-Fix case, we denote by c l,n the coefficients such that n∈Z c l,n ϕ(nT -t) = t l where l = 0, 1, . . . , L -1, by T the sampling step, and by [0, S] the support of ϕ(t).

Then, the abstract parameters ω k are related to the locations t k through ω k = t k and the annihilation matrices are given by

A k =      c k,1 c k,2 • • • c k,N c k-1,1 c k-1,2 • • • c k-1,N . . . . . . . . . c k-L+1,1 c k-L+1,2 • • • c k-L+1,N     
.

Additionally, there is a constraint on the minimal number of samples N for the GAP to hold, which is that N be larger than ⌈(S + max k {t k })/T ⌉.

Conclusion

We have shown how to unify the different techniques used in FRI signal reconstruction through an algebraic property that we call the Generalized Annihilation property. In essence, this property allows to solve nonlinear system of equations within two noniterative steps. We hope that this property can be used to solve other FRI problems (i.e, with new kernels) in particular in dimensions higher than 1 (for instance, like in [START_REF] Kandaswamy | Analytic sensing: reconstructing pointwise sources from boundary Laplace measurements[END_REF]), and maybe to solve other types of problems not directly related to sampling.