Essay Review

Histories of algorithms: Past, present and future

• Maarten Bullynck • Université Paris 8, France
The book under review, first published some twenty years ago, is the product of a group of French historians of mathematics: Jean-Luc Chabert, Evelyn Barbin, Michel Guillemot, Anne Michel-Pajus, Jacques Borowczyk, Ahmed Djebbar et Jean-Claude Martzloff. An English translation by Chris Weeks was published by Springer in 1999 (Chabert et al., 1999) and a second French edition is now available and will be reviewed here. It was at the time of its first appearance-and is stillthe sole book-length study of the history of mathematical algorithms. Therefore this second corrected and extended edition of the French original is most welcome.

The original edition of Histoire d'algorithmes was published at a time of developing interest in the history of algorithms within the mathematics community and it was, in part, conceived for use in the mathematics classroom. The historical range of the book is impressive; from Ancient Mesopotamia and Egypt to modern computer developments, this volume endeavours to forge a history of mathematics through highlighting the role of algorithms. The book contains 15 chapters and closes with short biographies of the authors discussed and a general index. Each chapter opens with a general overview that introduces the reader to the chapter's topic. The topic is illustrated and developed through a series of excerpts taken (or translated into French) from the original, historical documents. The text is always followed by a rewriting in terms of elementary modern mathematics and by some analysis and posthistoire of the algorithm. While the main text of the book has been largely kept as it was in the first edition, this more recent posthistoire has been updated for this second edition.

The book, one finds, falls naturally into two parts. The first (chapters 1-5), comparative in nature, focuses on the structure of an algorithm at different times in different cultures and traditions, principally, though not exclusively, pre-modern and non-Western. The second part remains largely within the realm of modern Western mathematics (17th-20th centuries) and presents a more or less linear story of how certain kinds of algorithms used in numerical analysis were developed and generalized.

The comparative first five chapters cover arithmetical operations (Ch. 1), magic squares (Ch. 2), false position (Ch. 3), Euclid's algorithm (Ch. 4) and the calculation of π (Ch. 5). Topics in the second numerical analysis part of the book include Newton's method (Ch. 6), approximate numerical solutions to equations (Ch. 7), systems of linear equations (Ch. 9), interpolation (Ch. 10), quadratures (Ch. 11), differential equations (Ch. 12), approximation of functions (Ch. 13), and acceleration of convergence (Ch. 14). Interspersed with these topics from numerical analysis are a chapter on algorithms in number theory (Ch. 8) and a chapter on the development of a mathematical definition of algorithm in the nineteen-thirties by Gödel,Post,Church,Turing and Kleene (Ch. 15).

Finally, this second edition features a new final chapter, "Epilogue. Ecriture, temps, hasard". This "epilogue" discusses some recent topics not covered in the first edition of the book, such as programming languages, time complexity, probabilistic algorithms, pseudo-random generators and quantum algorithms. The choice of topics in this last chapter reflects the interests of the authors in recent developments in computer science. The treatment is quite summary and the chapter as a whole is rather an afterthought, lacking the coherence of the rest of the book. It does not feature excerpts of original papers, except for Edouard Lucas' description of the Tower of Hanoi problem. In particular, it would have profited from some of the literature on the history of computing, such as Knuth and Pardo (1977) and Wexelblatt (1981) on programming languages or Mahoney (1997) on theories of computation, not to mention more recent work.

In this second French edition, many errors found during the translation of the English edition or pointed out by reviewers have been corrected. Nevertheless, the book still contains quite a number of inaccuracies that more careful editing might have prevented. These are mostly due the use of secondary sources that have not been double-checked. For example, the Liber Abaci was first published in 1202 not in 1228 (p. 290); the ENIAC was not an analogue (!) machine, but rather one of the first digital computers (p. 436); Lambert's article dates from 1768, not 1761 (p. 188); BNF is the abbreviation of Backus Normal Form, only later did the acronym come to mean Backus Naur Form (p. 539), etc.

The bibliography remains rather patchy and is inconsistent in style throughout the chapters (details of reprints are not given, page numbers are often missing, first names of author are not always given, etc.). In this second edition a restricted number of more recent bibliographic references have been added, though no effort has been made to systematically update the bibliography. Besides the references to primary technical sources the book only rarely references the secondary and historical literature. Also, in some cases, the authors use older editions and translations of ancient texts (such as Colebrooke's Bhaskara translation or Peyrard's Euclid edition), while better, modern editions are available.

Nonetheless, despite these problems, this is an excellent source-book, rich in texts and themes, eminently suitable for use in the classroom. The reproduction or translation of the original texts greatly helps an appreciation of the historical texture and dimension of algorithms and may in many cases be an excellent point of departure for adding a historical dimension to mathematical lectures.

And yet there remains a sense of frustration with this book today, in that it no longer reflects the state of the art in the historiography of algorithms. Adapting some of the historical concepts, tools and other innovations of the last twenty years would have added considerably more depth and improved the treatment of some topics in this second edition of Histoire d'algorithmes. Though it is understandable that the book was not entirely reworked for this new edition, it is unfortunate that the occasion to engage with more recent work in a foreword or afterword has been missed. We would like to make up for this omission by giving an overview of what has been happening in the history of algorithms in the last 20 years and point to some work still in progress.

The history of algorithms in computer science. At its initial appearance in 1994, Histoire d'algorithmes was part of a rekindled interest in algorithms within the historiography of mathematics, although the origin of the study of algorithms as a mathematical and historical object reaches back to the decades between 1960 and 1980 when computer science was trying to establish itself as an independent academic discipline. 1 It is in this process of disciplinarisation that we find the first serious historical reflections on algorithms and other basic concepts of computer science. In this respect it is informative to take a look at two important contributions to the history of algorithms that were written by computer scientists in the early 1970s. Both contributions still serve, albeit in different ways, as points of reference. Donald E. Knuth was one of the people who had been promoting the concept of algorithm as a basic notion in computer science with his classic three-volume The Art of Computer Programming (1968Programming ( -1973)). During the writing of this he seized the occasion of a special twentieth-anniversary issue of the Communications of the Association for Computing Machinery to publish his pioneering paper "Ancient Babylonian Algorithms" ( Knuth, 1972). The opening phrase of the article explicitly states his motivation: "One of the ways to help make computer science respectable is to show that it is deeply rooted in history, not just a short-lived phenomenon" (p. 671). Drawing almost exclusively on Otto Neugebauer's work, Knuth analyses some Babylonian mathematical texts with the help of his modern understanding of computer arithmetic and programming. Knuth especially foregrounds the fact that a certain "algebraic" (we would rather say, procedural) format existed for expressing algorithms, that they used floating point arithmetic and that memory operations and even a kind of assignment are clearly distinguished. Knuth exclaims: "The above procedure reads surprisingly like a program for a stack machine like the Burroughs B5500" (p. 674), where one needs to know that the B5500 machine was a quite advanced computer for its time that had been designed with a highlevel programming language in mind. As Knuth (1970) had shown two years before in an analysis of John von Neumann's first computer program (dated 1946), the development of a formal device for expressing algorithms and memory management during a computation is highly non-trivial. 2 Sharing in the disciplinarisation vogue for computer science, that same year Herman H. Goldstine published his book The Computer from Pascal to von Neumann ( Goldstine, 1972). Goldstine had trained as a mathematician, then had specialised in ballistics and become an officer in the Ballistics Research Laboratory at the Aberdeen Proving Ground during World War II. In that function he was the liaison officer with the team at the Moore School of Electrical Engineering, University of Pennsylvania, that was building the ENIAC, one of the world's first computers. Goldstine was responsible for recruiting John von Neumann for the group and, after the war, continued his collaboration with von Neumann, co-authoring various important founding documents in computing. Reflecting his own experience, Goldstine's book prominently featured applied mathematics (Part II, chapter 4: "Numerical Analysis"). This contains the bold-and historically inaccurate-statement that "the state of numerical mathematics stayed pretty much the same as Gauss left it until World War II" (p. 287).

Five years later, Goldstine published a book-length elaboration of that statement, A History of Numerical Analysis ( Goldstine, 1977). In the introduction, Goldstine too thanks Otto Neugebauer for "help and inspiration in the preparation of this book". Describing theories and algorithms now belonging to the toolkit of a numerical analyst, Goldstine's story opens with Napier's logarithms and goes on in a linear progression until Gauss (and a final chapter quickly runs over Jacobi, Cauchy and Hermite). Bringing in the great names from the history of mathematics, Goldstine constructs an illustrious pedigree for the new discipline of numerical analysis, without even mentioning the more practical and down-to-earth problems of manual calculation or the effects of the birth of the digital electronic computer. For many (including the authors of the book under review) A History of Numerical Analysis remains the classic source-book on the history of algorithms used in numerical analysis.

Both Knuth's Babylonian paper and Goldstine's books attempted to connect the history of the young discipline with an older tradition in the history of mathematics, the 'great tradition' of Göttingen for which Neugebauer stood. Where Goldstine's book, adapting itself to a mathematical mold, can be regarded as the (pre)history of a mathematical discipline, viz. numerical analysis, Knuth's paper, in contrast, brought the tools and concepts of computer science into a classic part of the history of mathematics.

In the 1980s and 1990s, around the time that computer science departments were being established at many universities and the notion of "algorithm" was slowly becoming a topic within the secondary-school curriculum, some professional historians of mathematics started addressing algorithms as well. To some extent independently of the developments in computer science, these historians had also discovered the algorithm as an historical mathematical object in its own right.

The history of algorithms in ancient cultures. The concept of "algorithm" became an analytical tool of historiographic research and innovation because it allowed to connect with many other tendencies in the renewal of the historiography of mathematics since the 1990ies. Leaving a perspective centered on the Western World and on "modern" mathematics that privileged demonstration, abstraction and innovation, historians of mathematics have started looking at the practices, contexts, traditions and cultures of mathematics and have tried to consider mathematical practices as one among many structured rational practices.

Focussing on algorithms has proven very fruitful in the historiography of Ancient non-Western mathematics. Among the first to draw attention to ways of algorithmic thinking are Jim Ritter in his work on Egyptian and Babylonian mathematics and Karine Chemla in her work on ancient Chinese mathematics.

In order to bring out both the differences and similarities between structures found in Egyptian and Babylonian procedural texts, Ritter, 1989 and Ritter, 2004 uses a simple symbolic rewriting system (and sometimes a schematic flowchart) to construct a context that may "have an historical basis while others serve a heuristic purpose for the modern student" (Ritter, 2004, p. 178). This formal analysis helps to follow the gradual complexification of examples and the re-use of already known "subroutines", it also maps the operations of writing down and recalling intermediate results. It is an analytical tool that brings out algorithmic structures in ancient texts and may be used to compare them. Another aspect should be insisted upon, such an analysis is not constricted to mathematical texts, but can be applied to other formally tailored texts. Ritter's work brings out the formal similarity of rational structures occurring in a variety of cuneiform texts, not only mathematical, but also medical or juridical texts. This algorithmic or procedural style is thus transversal to (modernday) disciplinary boundaries.

The focus of Chemla's work has been the structure of procedural mathematical texts as they codify and build ever more general or abstract arguments (e.g. Chemla, 1987a, Chemla, 1987b, Chemla, 1991, Chemla, 1992and Chemla and Shushun, 2004). In many mathematical traditions the dynamic interaction between problems and procedures is a central part of mathematical activity. On the one hand, a procedure is slowly explained and generalised through its application in a number of problems. On the other hand, a sequence of problems gradually mobilises combinations of subprocedures explained before, and a given problem may even be solved by different combinations of subprocedures. Procedures thus acquire a certain generality. For instance in Liu Hiu's comments on the classic Nine Chapters on mathematical procedures this generality is stressed throughout by the introduction of an "abstract" vocabulary that brings out the ties that exist between the several chapters. Thus, the comments introduce a level of generality and ways of verification that, although they do not equate with the abstraction of classic demonstration, are fundamental devices of mathematics too.

Ritter's and Chemla's work have led to other studies in the fields of Egyptian, Chinese, Indian and Babylonian mathematics. See, e.g., the work of A. Imhausen, 2002and Imhausen, 2003, A. Bréard, 2000and Bréard, 2008, A. Keller (2006), C. Proust (2007). 3 There has been a warning that reading ancient mathematics with algorithmic lenses risks of running into the same anachronistic misreadings as reading it with algebraic glasses (Høyrup, 2008). Most researchers, however, have avoided that pitfall because they used a procedural reading as an analytical, nearly linguistic tool that helps to foreground particular aspects of old texts, not as a goal per se. This work has also fruitfully challenged and innovated paradigms governing the history of ancient Greek mathematics. In (Ritter and Vitrac, 1998) two styles in Greek mathematics were discerned: a demonstrative (e.g. Euclid) and an algorithmic one (e.g. in Heron and Diophantus). While the first style has been classically given higher status than the second one by historians, the algorithmic tradition should also be regarded as a constituent part of Greek mathematics. Looking at Greek mathematics from this perspective brings it in closer contact with "oriental" traditions in mathematics such as the Egyptian or Babylonian styles of mathematics that display a more algorithmic character. 4 Building on the idea of procedural texts and their development through examples, Bernard and Christianidis (2011) have shown how the composition of Diophantus' Arithmetic progresses example by example from simple to complex with an accumulative introduction of new devices of invention, mirroring the rhetorical construction of a text as used in classic Greek education.

As the quoted studies abundantly show, an important characteristic of the algorithmic style is the use of sequences of (numerically computed) problems and examples to explain and gradually develop an algorithm. Example or problem series are in fact a recurrent feature in nearly all mathematical traditions and they properly belong to the main body of mathematical knowledge. 5 In particular, these series have a long half-life, they are often copied, recycled and recombinated 6 and their transmission and "recycling" is an important topic of modern historiography. Only in late 18th century Europe was a (pedagogic) separation enforced between the textbook and theory on the one hand and the series of problems and examples on the other hand, making the latter seem like a minor matter.

These sequences or series are of course part of an educational context in which a procedure is learned through repeated exercise, but where a procedure also becomes more flexible in the hands of a pupil who proceeds through more difficult, intricate or compounded problems. It is also an expression of the master's knowledge and skill at developing series of problems that go through all useful variations and complexities, demonstrating the possibilities of combination and modification, while hinting at links with other topics. Such series of problems can show the various forms of a procedure, the possible extensions of a procedure, its combinations with other procedures, etc. It can help to broaden or generalise procedures and the concepts involved, or it can specialise to intricate particular cases where a number of procedures have to be carefully fitted and tied together.

The history of algorithms in 20th century mathematics. Contrary to Goldstine's often reproduced account in A History of Numerical Analysis, the first formations of numerical analysis as a discipline come long before the end of World War II, John von Neumann's contributions or the birth of the digital computer. Goldstine's narrow focus on the 16th to 18th century and on the greater names in mathematics has been criticised by Folta (1988), and, more recently, by Tournès (2014). Both pointed out that computational practices as now codified in numerical analysis may be found throughout the centuries, but very often not in "pure" mathematics, rather in fields of practical, mixed or applied mathematics. A history of numerical analysis has to take this aspect into account. As a discipline, Tournès argues, numerical analysis was instead to emerge in the late 19th century and to have a further momentum of disciplinarisation after World War I. 7 Indeed, the first World War and its computational needs, especially in ballistics, seems to have stimulated a first systematic review of the mathematical literature in order to find efficient methods of computation ( Gluchoff, 2011 andAubin andGoldstein, 2014). This led to the publication of the first manuals on practical calculation during the Interwar period.

Newer historiography has tried to integrate the immediate practical context in which many of the computational techniques were first developed into the mathematical narrative. In the modern period, there are often forms of 'mathematical workplaces' involved, to borrow an expression from (Epple, 2000). In these workplaces a group of persons, who may or may not be trained mathematicians, take up, apply and use specific parts of mathematical knowledge in a technological environment. The "technological" environment may vary from crafts or schools that practise account keeping or surveying to highly organised industrial or military complexes of research and development. There exist a number of interesting studies on the mathematical practitioners and human computors in such institutionally and socially defined contexts as the observatory (Wepster, 2010 andAubin et al., 2010), the computation bureau (Croarken, 1990 andGrier, 2005), K. Pearson's Biostatistics Lab (J. Barrow-Green in Aubin and Goldstein, 2014) or on German computational work in aerodynamics (Epple et al., 2005).

Many traditions of computation are intrinsically intertwined with the use and production of tables and/or instruments. Tables or instruments are part of the materiality of calculation, as such, they shape the way mathematical procedures are thought up and transformed. The material aspect is often hard to reconstruct because it frequently relies on an implicit, tacit know-how that is not codified into textbook or manual but rather passed along orally or by experience. Although the book under review refers to this materiality in its subtitle, "from the pebblestone to the computer chip" ("du caillou à la puce"), it is generally absent in the book's treatment of algorithms. The pivotal role of tables in performing Babylonian and Egyptian arithmetics is only mentioned in passing (pp. 13-25), the discussion of interpolation methods leaves little space for the production and use of large (trigonometric, logarithmic etc.) tables is (p. 355f.), and the digital computer hardly appears at all. Recent studies, however, have focussed on the materiality of and the tools used in computing and mathematics.

In 2003 an edited volume (Campbell-Kelly et al., 2003) gave a first introduction to the histories of mathematical tables. Among the topics treated were the industrial division of labour in de Prony's or Babbage's table projects, the Committees of mathematical tables around 1900 and modern day spreadsheets. The book rather introduced the question and hardly exhausted it, besides its focus was mainly on the English-speaking world in the modern period. A more ambitious approach was formulated in a project funded by the French ANR and directed by D. Tournès who had already been working on the development of tractional instruments and nomographs in 19th and 20th century mathematics (Tournès, 2000 andTournès, 2009). The project History of Numerical Tables (HTN) wants to investigate the role of numerical tables in various traditions of mathematics. Since 2008 over twenty researchers from various countries have regularly met to prepare a book on numerical tables. The cross-cultural and cross-disciplinary discussions on tables have opened up a number of very general, methodological questions on the structure, production, use and transmission of tables. 8 The intertwining of mathematical procedures and numerical tables that is so typical of Babylonian mathematical culture, is also found in some schools of medieval Arabic science and in some early modern European mathematicians. Further, some domains of both applied and pure mathematics are characterised by a strong interaction between procedures of computation and numerical tables, e.g., astronomy, ballistics, 18th century elliptic functions or number theory. An overview of the project can be found in Tournès (2011).

Much of the interest in algorithms is inspired directly or indirectly by the advent of the digital electronic computer, but surprisingly few historical studies exist on the particular interaction between algorithms, the computer and its users or programmers. It has often been argued that mathematics has or is changing since the advent of digital computing, but those historical studies have mostly focussed on the history of modelling and simulation (e.g. Galison, 1997or Dahan-Dalmedico, 2001), or on the mechanisation of proof (e.g. McKenzie, 1999). Concrete studies on how the computer has impacted on mathematics or on the use of algorithms are still scarce. It is still mostly a history written by its actors who write a short, often personal, mostly linear story of the development of an algorithm in a particular context.

Admittedly, the material context of the computer is often difficult to take into account in a reasonable and apprehensive way. Some of this know-how has not been explicitly recorded and much of it may quickly lead to over-detailed and technical stories. However, the co-evolution of algorithms and the computer is an intriguing field that invites historiographic exploration. For one of the first electronic digital computers, the ENIAC (1946), a number of such studies have appeared. Bullynck andDe Mol (2010) andDe Mol, Carlé andBullynck (2015) have shown the reciprocal dependencies that exist between the machine, its users and the mathematics used and transformed. It analyses questions as how speed and accuracy of an algorithm have to be balanced against simplicity of control and program expression, or how the architecture of a computer conditions the short and long term memory management of a procedure. A key idea is that the early computers functioned as a kind of meeting ground where people from various backgrounds came together and interacted and where specific styles of computing could develop. Haigh et al. (2014) bring a reconstruction of the production process and life cycle of the first Monte Carlo program that ran on the ENIAC rewired into stored-program mode.

Such studies can cover a longer period too, and track the influence of the rapid evolution of computer hardware and of programming languages and techniques on the choice, sequences and variations of algorithms within a specific discipline. Franck Varenne, 2007 andVarenne, 2008 has shown, for instance, that increasing speed and memory of computers, combined with the development of object-oriented programming, made it possible to study not only of the structure of plants, but also of the temporal unfolding of the structure itself through computer modelling. In a field that has gained in prominence due to cryptographic applications, computational number theory, Bullynck (2009) has described how a complex taxonomy of various factorisation algorithms slowly developed, not only echoing the size of the number to be factorised but also the variety of computer devices, from parallel to serial architecture, from hand-held calculators to supercomputers. In yet another field, that of chess programs, Ensmenger (2012) showed that over the years a rather brute-force minimax algorithm grew to be preferred over subtler, heuristic algorithms.

Conclusion.

It may be clear by now that the concept of algorithm has borne much fruit in recent historiography, but also that much work remains to be done. The book under review, Histoire d'algorithmes, still is a nice textbook to introduce the student to the fascinating field of the history of algorithms. This review essay has tried to supplement it by introducing a number of topics that have been explored in the recent historiography of algorithms and serves as an invitation to the reader to discover this thriving field. For those who want to dig deeper into a specific topic, the rather extensive list of references listed here may provide orientation and inspiration fur further reading.

Acknowledgements

This review essay was begun in collaboration with Jim Ritter while we were teaching a History of Algorithms course at Université Paris 8. Although I alone now feature as the author, this review owes much to this collaboration and to the inspiring discussions I had with Jim. All errors, interpretations and misconceptions are, of course, solely my responsibility.

 (Wu (1980, 43-44), cited in Hudeček (2012, 55)): "We set out the question and came up with a method of solution under inspiration from Ancient Chinese algebra. The reason is actually easy to understand; Chinese ancient mathematics was basically a mechanized mathematics.... The work of the present author on mechanization of mathematics is precisely a product of inspiration from these ideas and results, it is a direct continuation of our mathematics from the Nine Chapters up to the Song and Yuan periods." (This footnote is due to Jim Ritter.) 3 It may be remarked that the treatment of false position in the book under review (pp. 99-102) could be updated, fine-tuned and often improved by consulting this recent work on algorithms in non-Western traditions, viz. (Høyrup, 2002, pp. 59-80;101-107) and (Ritter, 2004, 179-184) for Babylonia; (Imhausen, 2003, pp. 39-42) and(Ritter, 2004, pp. 184-186) for Egypt; (Chemla and Shushun, 2004, pp. 549-599) for China.

4

Since 2009, an ANR-funded project "Pratiques déductives algorithmiques dans les mathématiques pré-algébriques" (ALGO) directed by F. Acerbi and B. Vitrac has brought together researchers on ancient mathematics to address these questions more systematically. See http://algo.hypotheses.org.gate3.inist.fr/a-propos.

5

Under the direction of A. Bernard at HASTEC a research group has been focussing on the role of series of problems in different (algorithmic) traditions ("Les séries de problèmes, un genre au croisement des cultures"). See http://problemata.hypotheses.org.gate3.inist.fr/. 6 E.g., during his work on his PhD, A. Heeffer (2006) has drawn up a database of problems and their copies and modifications in early algebra textbooks to study their migration and dissemination.

7

The list of early university courses on numerical analysis given in (Brezinski and Wuytack, 2001, pp. 4-20) confirms this. 8 See http://www.sphere.univ-paris-diderot.fr/spip.php?rubrique64.


N Akrour
email: akrour@latmos.ipsl.fr

F Mercier

C Mallet

L Barthes

Rainfall at Urban Hydrological Scales : Opportunist use of TV-SAT Receivers Network

Conclusion & perspectives

Filter term Gap with the observations • 𝑅 0 is the initial rainfall field discretized on a spatiotemporal grid. 𝑅 𝑛 is this initial field propagated by an advection numerical scheme from 𝑡 0 = 0 to 𝑡 𝑛 = 𝑛∆𝑡 • 𝑦 𝑛 contains the observations available at time 𝑛 (rain gages and KU-waves). • 𝐻 𝑛 is the operator projecting the field 𝑅 𝑛 on the observations space at time 𝑛 • 𝑅 𝑖𝑗 𝑛 is the mean rainfall around grid point (𝑖, 𝑗)

• 𝛼 is a coefficient to balance the weight of the 2 terms

A measurement system (Figure 1 ) allows estimating integrated rain attenuations along satellite-ground links (depending on the freezing level and on the relative position of the satellite from the receiver) with a high temporal resolution (15s) see [Barthes and Mallet (2013)]. A 4D-VAR assimilation algorithm (Figure 2) using an advection numerical scheme has been applied on 4 rain attenuation times series recorded by a single Ku band TV-SAT receiver during HyMeX experiment. The obtained rainfall fields [Mercier et al. (2015)] at resolution 10s and 0.5x0.5 km2 have been validated against co-located radar maps and rain gauges data.

• TV-SAT measurements networks can be used in assimilation algorithm to retrieve rainfall fields.

• A large number of links allows estimating original rain field including rain cell structure with a very accurate estimation • fine scale resolution retrieval (10 s and 500 x 500 m) • urban area all over the world (conditioned by the presence of TV SAT receivers) • TV-SAT measurements can be easily merged with other observations type (raingauges) in the assimilation algorithm.

The aim of this work is to rebuild rainfall fields at small scale (around 30 * 30𝑘𝑚 2 with 0,5 * 0,5𝑘𝑚 2 grid boxes) from measures of attenuation of electromagnetic waves. Theses waves come from geostationary TV satellites and are in the Ku band (around 12𝐺𝐻𝑧).

Opportunist measurement

Fig. 2 :

 2 Fig. 2 : Principle of the 4D-VAR data assimilation

Fig. 1 :

 1 Fig. 1 : Device measuring the attenuations of the microwaves of the KU band

Acknowlegement :

This work was supported by the French Programme National de Télédetection Spatiale (PNTS, http://www.insu.cnrs.fr/pnts), grant no. PNTS-2013-01, and by the IPSL (Institut Pierre-Simon Laplace)


François Léonard
email: leonard@ireq.ca

and Phase Unwrapping: Referencing the Phase at the Center of the Observation Window

Keywords: Continuous Signal Spectral Analysis, Tone Phase Estimation phase reference, phase unwrapping, DFT, FFT, IFFT, z-transform, spectral window and multi-tone parameter estimation, reconstruction, Cramér-Rao bound, Dirichlet kernel, high-resolution spectral estimation, high-accurancy spectral estimation

Introduction

In spectral analysis, the phase was, and still is, perceived as information that is more abstract and sophisticated than the amplitude and, therefore, a lesser degree of precision is accepted in its estimation. This view is broadened by the fact that, on the one hand, the phase is subjected to one or more  rotations under a spectral lobe and, on the other, the value of the phase is Hydro-Québec, 1999

Phase referenced to window center 2 less often used than the amplitude or frequency values. Despite the fact of some authors [START_REF] Rife | Single-tone Parameter Estimation from Discrete-time Observations[END_REF][START_REF] Ferrero | A New Approach to the Fourier Analysis of Periodic Signals for the Minimization of the Phase Errors[END_REF][START_REF] Offelli | The Influence of Windowing on the Accuracy of Multifrequency Signal Parameter Estimation[END_REF][START_REF] Schoukens | The Interpolated Fast Fourier Transform: A Comparative Study[END_REF] have decribed different ways to perform accurate phase estimation with phase referenced to the center of the observation window, many others do not used the latter's algorithms or similar algorithm with phase referenced to the window center. The knowledge about how to process the phase is not well established among scientific community and, indeed, for many the lack of precision in estimating the phase is a major impediment in its use.

The phase information is essential for many applications. Processing of the phase information over many channels is common in the field of mechanics. In mechanical vibration analysis and, more specifically, modal analysis [START_REF] Ewins | Modal Testing: Theory and Practice[END_REF], machinery diagnosis [START_REF] Lyon | Machinery Noise and Diagnostics[END_REF],

acoustic [START_REF] Bendat | Engineering applications of correlation and spectral analysis[END_REF] and sonar, phase information is vital for describing the phenomenon. It should be noted that modal analysis is largely used to experimentally validate finite-element modeled structures. In modal analysis, the phase is compared among different points of the structure for a given frequency mode. Similarly, in the field of electricity, simultaneous phase measurement over several channels is typically used in radar and power grid stability [START_REF] Missout | Measurement of Bus Voltage Angle Between Montreal and Sept-Iles[END_REF][START_REF]IEEE Standard for Synchrophasors for Power Systems[END_REF].

The phase is also of interest in phase modulation and demodulation in data communication, in which case phase modulation over a single channel is observed. In fact, regardless of the field of application, the value of the sinusoid phase on one channel is mainly (of use?) useful when:

-It is compared with another measurement performed simultaneously on another channel for the same given frequency;

-It is compared with another measurement performed elsewhere on the same channel for the same carrier or;

-The absolute phase value is request for a full reconstruction of the sinusoid found in the signal.

This last aspect concerning the reconstruction of the sinusoid is currently given very little consideration in the signal processing field. In fact, small proportions of scientists are currently able to reconstruct the sinusoid with an error corresponding to the noise embedded in the signal. In some cases, reconstruction of the sinusoids is required if one wishes to eliminate the sinusoids from the signal [START_REF] Danthez | Autocoherent spectrum : a usefull estimator for vibration analysis of rotating machinery. Accurate estimation and cancellation of pure tones[END_REF]. For instance, in detecting partial discharges in substation transformers, it is often useful to eliminate the radio transmission signal so that only partial discharge impulses remain [START_REF] Léonard | Referencing the Phase to the Center of the Spectral Window[END_REF][START_REF] Léonard | Pourquoi faut-il référencer la phase au milieu de la fenêtre spectrale[END_REF]. This time domain elimination of sinusoids give better results than subtraction of frequency components in frequency domain [START_REF] Wenzel | Noise suppression and data reduction for partial discharge measurements using orthogonal transformations[END_REF] since Hydro-Québec, 1999

Phase referenced to window center 3

we wipe not only the main spectral lobe but the whole spectral pattern corresponding to the sinusoid. In others cases, the sinusoid parameters are the requested information. For example, in the monitoring of hydraulic turbines [START_REF] Léonard | Hydro Equipment: Learning from Hydro-Québec's Monitoring Experience[END_REF], each vibration is stored in a database as a frequency, amplitude and phase. Through reconstruction, the orbit trajectory of the shaft in the guide bearing can then be synthesized at different dates and machine vibration behavior can be followed. This method used to compress, reconstruct and display the data is know as "holospectrum" [START_REF] Liangsheng | The Holospectrum: A New FFT Based Rotor Diagnostic Method[END_REF] in machinery vibration monitoring. The precise identification of a sinusoid in view of its later reconstruction thus finds many applications. This use of the phase estimation should become more common once everyone will be able to accurately extract the phase information.

We believe that there are still some gaps in the field of phase signal processing due to the simple paradigm where the phase is referenced at the beginning rather than in the center of the spectral window.

Effect of a phase reference shift on z-transform

To understand better the DFT phase behavior we must look at the z-phase distribution when few poles are present. For the next examples, we will used the finite-length sequence ztransform

X z   z k z k  x n   z n n0 N1  (1)
where the ratio z k z k shift the phase reference to the sample number k. Moreover, for enhance the amplitude graphic illustration, we have introduce a weighting function to the latter equation. We have plot the value

X z   z k z k  x n    z n n0 N1  z n 2 n0 N 1  (2)
instead of the value of 1. For our concern, this weighting function has no effect on phase distribution.

Hydro-Québec, 1999

Phase referenced to window center 4

Fig. 1 show the z-plane amplitude and phase for a single pole associated with a complex exponential located at f=+0.15 bin on 1.05 z-radius, for N= 20 time samples. The so-called pole circle is the virtual circle centered on the z-plane origin and passing through the pole.

With phase referenced to the first time sample, i.e. with k=0, we have phase circular slope inside the pole circle and a phase flatness outside the pole circle. On this figure and the followed one, each fringe correspond to a 2 radians phase step since the - value is black and + value is white. The phase increase in clockwise direction on a circle inside the pole circle illustrated on Fig. 1. To be more covenants for a connection between z-transform and a DFT of a real signal, we use a sinusoid instead a single complex exponential for the next developments. The Fig. 2a shown in the z-plane the two poles corresponding to a constant amplitude sinusoid at f=0.2055 bin. It should be note that the additivity theorem of poles in z-plane allowed complex amplitude addition, not phase addition. This explains the phase distribution of the double poles associated with a sinusoid. More specifically, for the example illustrated at Fig. 2b, the «petal» on the 180° z-plane direction appears more short than others one.

Hydro-Québec, 1999

Phase referenced to window center 5 The Figs.2b-2d shown the effect of phase reference on phase distribution in the z-plane. For the usual phase reference k=0, the z-phase increase clockwise with an N-1 radians/radian rate inside the pole circle. In the last ratio, the first "radian" is coming from the phase slope and the second is the z-plane polar angle. The "minus one" is generated from the  phase drift under pole. For the phase referenced to the window middle, k=N/2, the phase slope value (N-1)2 radians/radian is equally distributed between inside and outside pole circle with an opposite sign: we observe a anti-symmetry pattern lock on the pole circle. And lastly, for a phase referenced to the opposite window extremity, for k=N, the z-phase increase counterclockwise with an N-1 radians/radian rate outside the pole circle and have only a 2 shift inside the pole circle.

Hydro-Québec, 1999

Phase referenced to window center 6 We don't know if the z-plane phase distribution could be useful in some applications, but the DFT phase are request for many duties. The Figs. 3 and4 demonstrate numerically, respectively for k=0 and for k=N/2, that the phase value on the z-plane unit circle correspond to a zero padding DFT. Here, the small mismatches are coming from jitter on sampling position of the unit circle in a z-plane express as a 400 x 400 matrix: the coordinate values of circle are rounded to the nearest matrix indices. Observe how some small position discrepancies on z-plane, less than 1%, lead a visible phase curve mismatch on pole circle.

We must focus our attention on three important facts deduced from the Figs.1 to 4. First, if we are looking only one pole, the pattern reach his maximum complexity near the unit circle.

The radius position of this pattern are sensitive to any pole drift like it's occurring when a slight amplitude modulation decrease or increase a sinusoid amplitude in the observation window. This modulation has similar effect than do a circle cut with a different radius than the unit circle for a constant amplitude sinusoid. Moreover, in harmonic analysis with multiple quasi-constant tones, the corresponding poles are in the vicinity of the unit circle creating a complex phase pattern on this unit circle. Second, the phase shown very different results with k=N/2 than for k=0. This difference is amplified with unwrap phase (see next chapter for explanation about phase unwrapping). Under the pole main lobe this difference appear like a constant slope addition. Third, without the phase resolution enhance generated by zero padding on DFT, we have some "phase 2-skip" generated by a phase slope greather than  radians/spectral line. For example, if we cut the phase illustrated at Fig. 2b inside the unit circle, at radius ≈ 0.8 in the counterclockwise direction, with a chirp z-transform [START_REF] Rabiner | The Chirp z-Transform Algorithm[END_REF] (it's a DFT windowed by a exponential decreasing window), we obtain a phase shift of 2(N-1)/N radians between two consecutive lines: a phase unwrapping algorithm interprets this phase shift like a 2/N phase drift line to line. On unit circle, like the illustration on Fig. 3, we have few phase 2-skip occurring if we retain only the samples corresponding to the spectral lines. Conversely, no phase 2-skip could occurs in Fig. 4 since the maximum phase slope is less than  radians between two consecutive lines. However, with some poles a bit aside the unit circle or in the case of real field measurement, both phase references shows phase 2-skip without an appropriate zero padding to decrease line to line phase difference. Since the mean (average?) phase slope have the half value with k=N/2, the zero padding factor requested to avoid phase 2-skip is reduce too. Note that the zero padding factor is define like the ratio of DFT array length over original time array length. Finally, under an isolated pole main lobe, the phase differences line to line for k=0 is at the edge of phase 2-skip condition. On opposite, for k=N/2 the phase are in the best position, close to zero slope, to avoid phase 2skip.

In conclusion: the DFT cut the z-plane in the most phase sensitive region for a quasi-constant tone signal. The phase processing result is strongly function of the phase reference and phase 2-skip is common with phase unwrapping without zero padding on DFT. Phase 2-skip probability is reduced with phase referenced to middle time sample and can be minimize or avoid by zero padding technique.

Phase presentation and phase unwrapping in DFT frequency domain

Let us take the conventional image of the trigonometric circle to represent the phase value of a spectral line. This value usually originates from the arc tangent function "arctg{imaginary part, real part}," which yields the value of the angle in four quadrants. If we repeat this circle at each spectral line, we will obtain a cylinder on the frequency axis, as illustrated in Fig. 5. The frequency axis is located at the center of the cylinder and therefore is not in contact with the surface of the phase domain. It can be seen that a boundary line exists between the values  and -. In view of the difficulty of presenting this cylinder graphically, we will cut it along this boundary line. The choice is an arbitrary one because we could just as easily have taken the 0 radian line or any other angle value. Once it has been cut, the cylinder can be stood on end and displayed in a planar form, which is more convenient for screens and printers. Naturally, it takes a little imagination to be able to align the parts of the curves that have been cut in this "mapping." Fig. 5 presents a single example of a typical cut where the phase is expressed in the form of a curve as a function of the frequency and plotted on a cylinder centered on the frequency axis. In order to obtain a 2D presentation, all that is required is to cut the cylinder lengthwise, open it out and place it flat down. shortest path. The choice becomes more complex when two values are only about  apart and when we cannot tolerate a 2-skip in our application. Interpolation of the spectrum by adding zeros on either side of the time signal before applying the DFT makes it easier to select the shortest path. Instead of the latter "zero padding" solution, a numerical integration of the phase derivative can be chosen [START_REF] Tribolet | A New Phase Unwrapping Algorithm[END_REF]. The unwrapped phase can be displayed as shown in Fig. 5. Here, unwrapping allows the cut segments to be plotted in the correct trigonometric direction, i.e. clockwise or counterclockwise. The unwrapped phase can also be displayed without discontinuity on a scale covering several  radians [START_REF] Oppenheim | Digital Signal Processing[END_REF] (Fig. 11). In fact, if the boundary /- is cut in the same trigonometric direction several times, we accumulate as many rotations of 2 radians.

If we look back to the raw phase values illustrated at Figs. 3 and4, the vertical lines corresponding to some paths passing through /-boundary line do not really exist. Indeed, it is more suitable to show unwrapped phase values instead raw phase values. Accordingly, the unwrapped phase fitted in the ]-,] domain present a clear cut of the phase paths at - and  boundary lines without the latter ghost vertical lines.

Effect of time delay and time reference shift on the DFT phase

In the following, the spectral component nomenclature is used for the frequency object associated with a continuous-tone sinusoid. Where:

s nT    a 0  cos  0  nT   0   (3)
the sinusoid of amplitude a 0 , phase  0 and frequency f 0   0 2 using time t=nT and where T is the sampling period. For a rectangular window, the frequency object is the circular convolution over N spectral lines

S    a 0 exp  j 0       0    exp j 0       0   2        D T   (4)
between the sinusoid Fourier transform (FT) and the Dirichlet kernel [START_REF] Lanczos | Discourse on Fourier Series[END_REF] (see discussion in annex A)

D T   exp  j N  1 2 T     sin N 2 T     sin 1 2 T     , T   ,   , (5)
corresponding to the FT of the rectangular window defined as Hydro-Québec, 1999

Phase referenced to window center 10

w(n)  1.0 , n = 0,1,2,...,N -1 w(n)  0 , otherwhere , N even . (6)
A time delay of l samples consists in moving the signal's observation window as shown in Fig. 6a. This time delay only affects the left member of the convolution expressed in 4. It corresponds to moving the sinusoid phase with a value of l 0 T such that

S l,   a 0 2 exp  j  0  l 0 T         0   exp j  0  l 0 T         0      D T   (7)
is the component of the sinusoid observed with a delay of l samples.

On the other hand, a time shift of the DFT phase reference affects the both members of the convolution. To modify the phase reference, a "circular shift" is applied to the array that contains the samples transmitted to the DFT (Fig. 6b). The right or left part of the array is then truncated into k samples and joined to the other end of the array. This circular shift added a complex exponential product to the spectral components given at 4 like

S k,   exp jkT   a 0 exp  j 0       0    exp j 0       0   2        D T         (8)
with the resulting equation

S k,   a 0 2 exp  j 0   exp jkT   D    0  T   exp j 0   exp jkT   D    0  T     .(9)
Since the phase argument of the Dirichlet kernel

D k,T   exp  j N  1 2 T  jkT     sin N 2 T     sin 1 2 T     ,  T   ,   , (10)
modified by a "circular shift" of k samples is tilted by a linear phase shift kT, the two shifted kernels appearing in 9 are tilted by the same linear phase shift plus a constant phase shift  jk 0 T . In fact, we can write the sinusoid DFT 11) like a sum of two shifted and tilted kernels.

S k,   a 0 2 exp  j  0  k 0 T      D k,    0  T    exp j  0  k 0 T      D k,    0  T     (
Hydro-Québec, 1999

Phase referenced to window center 11 There is only one instance when a time delay and the time shift of the reference phase give the same result, i.e. eq.7 equivalent of eq.11, which is when the signal occurs periodically and the DFT exactly covers an integer number of periods. In the latter case, unless there is some non-periodic noise, only one sample per kernel associate with a spectral component is non-

nil, since D k,T   N    for NT  i  2 /  0 with i integers number.
If we discard this specific case resulting from a digital simulation and look at a measurement of a physical phenomenon, the time delay and reference time shift affect the phase in different ways.

More specifically, the phase relationship, which exists between the spectral lines that support a spectral component, remains practically unaffected with a time delay in the observation window. A shift of the observation window over time causes the positive-frequency component to be rotated in the opposite direction in relation to the negative-frequency component. An angular rotation of two Dirichlet kernels is then observed, with each kernel being associated with one of the two complex exponentials that form the sinusoid. If we look at the spectral lines of one of the two kernels, all of its lines undergo the same phase rotation. However, a time shift of the phase reference on k samples is equivalent to twisting the right end of the phase cylinder (Fig. 5) by kradians while keeping the left end set at 0 Hz. A time reference shift adds a linear phase shift between the spectral lines of a spectral component. As shown in Fig. 6 (right-hand side), in the first case a uniform vertical translation of the kernel phase occurs, while the second case involves a linear shift of the kernel phase.

Hydro-Québec, 1999

Phase referenced to window center 12

Main lobe phase plateau with phase referenced at the center of a spectral window

A N  1   2 sample shift of the phase reference provides the simplest kernel:

D k  N  1   2,T   sin N 2 T     sin 1 2 T     (12)
associated with a phase referenced at the center of a rectangular window. This kernel is "in plane" since they fit in the real plane and do not contain any imaginary part. Indeed, the phase values are 0 or .

The significance of the difference between 5 and 12 becomes evident when one looks at the distribution of the spectral line phase values related to a spectral component without zero padding. In Fig. 7, a plateau appears where  phase rotations are found and, inversely,  phase rotations appear where there is a plateau. If we excluded some numerical simulation cases and look at field measurements, the frequency of a sinusoid is located between two spectral lines. Finding the phase value presents a problem if the latter was referenced at the beginning of the window (k=0). In fact, interpolating the phase between two spectral lines is no easy matter when there is a phase rotation of  radians between these lines: any slight frequency estimation error will results in a greater phase deviation (Fig. 7). Consequently, referencing the phase at the center of the DFT sequence thus provides an advantage with the rectangular window. We will see that this advantage increases with the use of a spectral window such as cosine sum or a Gaussian shaped window. Before the demonstration about the presence of a large phase plateau under spectral lobe with some spectral windows, we will look at the effect of the applying window on constant tone. A window with a kernel W k,T   apply to the sinusoid described by 3 give a similar result, (13) than for a rectangular window (eq.11) and lead the same property: a circular shift in k samples introduces a linear distortion of the phase relationship between the spectral lines of the kernel, as illustrated in Fig. 2, plus a constant phase shift of  k 0 T on negative frequency kernel and of  k 0 T on positive frequency kernel.

S k,   a 0 2 exp  j  0  k 0 T      W k,    0  T    exp j  0  k 0 T     W k,    0  T     ,
To demonstrate existence of phase plateau, let us now look at the case of a spectral window (14) defined as a cosine sum like defined by Rife-Vincent [START_REF] Rife | Use the discrete Fourier transform in the measurement of frequencies and levels of tones[END_REF]. This definition includes the Hamming, Hanning, Blackman and Blackman-Harris window [START_REF] Harris | On the Use of Windows for Harmonic Analysis with the Discrete Fourier Transform[END_REF], but excludes Dolph-Chebychev and Kaiser-Bessel windows. For x  0, X  , the latter expression become

w x    c 0  c i  cos 2 X  i  x     i 1 I  , with x   X 2, X 2   , I  N / 2 and c i i 0 I   1.0,
w x    c 0  1   i c i cos 2 X i  x     i 1 I  (15)
Hydro-Québec, 1999

Phase referenced to window center 14

and match the DFT domain definition for the sampling process

x i  nT with n  0,...,N  1   (16)
where X  NT . Note that the latter window center is located at n=N/2 instead N  1   2 for the rectangular window since the window is symmetric over (with respect to/de part et d'autre) the maximum value of 1.0 at n=N/2. In Fourier analysis, the basic hypothesis is periodicity of the signal. If we repeat the latter window with a period X, we can write

W     1 2 1   i c i     2 X  i         2 X  i         i0 I  (17)
as the Fourier transform of 15. Apply a finite sampling correspond, in frequency domain, by the convolution of the window Fourier transform given at the Eq. 17 with the Dirichlet kernel given at the eq. 10 with k=1/2 to take account of the half sample shift of the window center (see annex A). If we add the k samples shift to insert the phase reference parameter, this sampling yield the window kernel

W k,T   exp jkT   1 2 1   i c i  D 1 2 ,T  2  i N      D 1 2 ,T  2 i N         i0 I  (18)
generated by a sum of Dirichlet kernels. This latter equation become

W k,T   exp jT k  N 2          1 2 c i  D N  1 2 ,T  2  i N      D N  1 2 ,T  2  i N         i 0 I  (19)
after some manipulations. Consequently, the window kernel W k,T   is an in-phase sum of "in plane" Dirichlet kernels give at eq.12 multiply by a complex exponential drive by an argument function of k  N / 2  . Indeed, since the c i values are positive, the summed Dirichlet kernels have their main lobes in phase inside the summation. For k=0, the complex exponential of 19 put an alternating signs line to line under a main spectral lobe like the example given at Fig. 8.a. On the opposite, when the phase is referenced to the window center with k  N  1   2 , we obtain a flat phase under a main spectral lobe (Fig. 8.b) since we have 1+2(I-1) in phase Dirichlet kernels joins side by side. For this latter case, a phase plateau is present over 1+2(I-1) spectral lines under the main lobe. For example, the Blackman-Harris 3-terms window shows a flat phase over 5 spectral lines.

Hydro-Québec, 1999

Phase referenced to window center 15

R I  t a) W(0, t) b) W(N/2, t)
Fig. 8: Referencing the phase to the beginning of a Rife-Vincent class window creates an alternating sign under the main kernel lobe (a), then instead, referencing the phase to the middle of the window leave a continuous phase value (b).

Similarly, for a gaussian spectral window defined by

w n    exp  1 2  n N 2      2       , for 0  n  N 2 (20)
and yielding the kernel [START_REF] Harris | On the Use of Windows for Harmonic Analysis with the Discrete Fourier Transform[END_REF] W T  

N 2 2  exp  1 2 1  T     2       , for   2.5 and T small, (21)
it is easier to visualize that the phase remains the same for T   . For instance, for =3.5, the phase is constant on more than 6 spectral lines.

We can generalize that, for symmetrical window like rectangular and bell-shaped spectral windows, referencing the phase at the center of the window simplifies the mathematical expression of the kernel by removing a complex exponential product and put a phase plateau under the main lobe. Many find that use of a low sidelobe spectral window increases the complexity of the form the phase will take because, when referenced at the beginning of the spectral window, the phase undergoes more than one rotation of 2 radians under the main lobe. As shown in the diagram at the top of Fig. 9, it becomes problematic to choose the frequency at which we have to take the phase value. On the other hand, referencing the phase at the center of the window allows us to obtain an almost horizontal flat section under the main lobe where it is easy to estimate the phase value.

The kernel associated with low sidelobe window increases both the widths of the lobe and the phase plateau horizontality. Factors that affect plateau horizontality include noise, leakage and modulation. The presence of noise under main lobe and leakage noise affect small amplitude components.

Leakage is mainly coming from sinusoids of neighboring frequencies. When we are not close to =0 or =/T, we can neglect the overlap of the sidelobe originating from the negative frequency components. Usually, the leakage can be well minimized with an appropriate window having a great sidelobe rejection. In many physical measurements, like done in mechanic field, the most horizontality disturbing factor is fluctuations in frequency and amplitude that the observed sinusoid may have. Sometimes, when a periodic fluctuation observed over less than few modulations beat, we don't have enough resolution to put the modulating components out the main lobe of the carrier.

Othertimes, the fluctuation is none-periodic creating a distortioned phase plateau like when we do a circle cut aside a pole in a z-plane.

Let us now look at a few practical aspects. For a DFT performed on N time samples, a scroll delay of N/2 samples is equivalent to twisting the cylinder by N/2 radians, i.e. a jump from a value of  at each spectral sample. From an algorithmic perspective, to move from one reference to another, from k=0 to k=N/2,  should be added to or subtracted from the phase of the uneven spectral lines, with the first line being line zero. Note that the correct phase value is that which originates from the interpolation of the value of the phase plateau at the sinusoid frequency. Further on we will present an algorithm which allows this frequency value to be found between two spectral lines. In the case of several sinusoids, a phase plateau is obtained under each of these, when we have no overlapping between two main lobes. The phase can then be recorded directly on the graph as shown in Fig. 10. However, for a DFT which includes a fairly large number of samples, it would be better to present the latter over an interval of ]-180°,180°] or ]-,] radians. In fact, the unwrapped DFT phase with over 100 samples often extends over an interval of several  radians. Consequently, for a DFT sequence with over 1000 samples, the phase plateaus are barely noticeable. Lastly, it should be remembered that the horizontality of the phase plateau is affected by a fluctuation in the amplitude or frequency. In the example in Fig. 7, the free vibration of the beam after a mechanical shock has a decreasing exponential amplitude envelope whereas the frequency increases by a fraction of the spectral line with a decrease in amplitude. This is why the phase plateau of the spectral component at 2770 Hz appears slightly inclined in this figure.

FFT high-accuracy tone-parameters estimation algorithm

A high-accuracy estimation algorithm provides the frequency, amplitude and phase values of a continuous-tone sine. For algorithms working in frequency domain, these values are Hydro-Québec, 1999

Phase referenced to window center 18 calculated using the spectral lines located under the main lobe associated to a spectral component in the FFT or DFT. The algorithm does not increase the resolving power, i.e. the distinction between the two spectral components which partly overlap. The high-accuracy estimation designation comes from the increase in accuracy of the sine parameter estimation. Without any high-accuracy estimation algorithm, the frequency of the largest-amplitude spectral line will be attributed to the sinusoid. There are several high-accuracy estimation algorithms, with the simplest ones being the "zero padding" algorithm and the calculation of the center of mass. Once the latter has been presented, we will then present the iterative algorithm that was developed. The "zero padding" technique consists in laying out an equal number of nil? zero? samples on either side of the time sequence before FFT processing when the phase is referenced to the window center (Fig. 11). Whenever, for phase referenced at the first time sample, the zeros are added at the end of the signal samples. For example, doubling the number of samples in such a manner allows the spectral component to be interpolated such that a new spectral line is found between each previous line. The calculation time and memory become exorbitant for a estimation better than 1% from the distance between two lines, i.e. 1% of 1/NT: the calculation time is then multiplied by 100 ln(100N)/ln(N) and 100 times more memory is required. Here, the spectrum detail is increased equally throughout the spectrum rather than only between two spectral lines at the top of the main lobe. A large part of the calculation effort and memory are then wasted. However, in practical applications, the zero padding is typically used to increase the number of samples for phase unwrapping and smooth graphic visual presentation of spectrum and spectrogram. The zero padding factors are usually set from 2 to 5 for the latter applications.

Hydro-Québec, 1999

Phase referenced to window center 19

Calculating the center of mass involves a simple and efficient estimate of the frequency between the two spectral lines. An FFT provides a sample of S(k,) where the value S i  S(N / 2,2 i NT ) is the complex amplitude of the spectral line i. Express in spectral line number, the frequency

f  i  S i S i * iP  S i S i *  , (22)
is resulting from a calculation of the center of mass over the distribution of main lobe energy of the lines defined by set P. Since the spectral amplitude is directly related to the sine amplitude, we can extract sine amplitude from S(k,). Like Hanrahan and Valkin [START_REF] Hanrahan | Fast method of estimation of signal parameters using the discrete Fourier transform[END_REF], for amplitude estimation we can used the amplitude sum

a 0  S(0,2 i NT) iP  w n   n0 N 1   2 N w 2 n   n 0 N 1  (23)
over main lobe. Another way is the use of Parceval relation for energy conservation where power of the sinusoid a 0 2 2 is equal to the sum of the power spectrum lines associated with the two sinusoid kernels. The amplitude estimator

a 0  4 S i S i * iP  N w 2 n   n 0 N 1  (24)
give a best results than 23 with less bias and less sensitive to noise, since the sum of square amplitudes put more weight on lines with high S/N ratio. Indeed, the last three equations have results not related to the phase reference. Finally, the phase estimation

 0  arg S int f      f  int f      arg S int f 1      arg S int f       (25)
is obtain from a linear interpolation with phase reference set on k=N/2, where int(f) is the integer part function of the frequency expressed in line number. The main source of error stems from truncating the summation. From a practical perspective, this sum extends over the few spectral lines to the top of the component. The weights of the lines which are not taken into account then causes an error in estimating the center of mass and involve a systematic underestimation of amplitude. When used rectangular window, this calculation is not recommended since it generates a systematic error of around 20% of the spectral interline spacing when a sum is done from the 8 spectral lines of greater amplitude. However, with a Hydro-Québec, 1999

Phase referenced to window center 20 high-sidelobe rejection window like the Blackman-Harris-4-term window, the center of mass is estimated at better than 0.6% of the interline spacing with only 5 spectral lines considered. With this last window, a sum on 8 spectral lines without any noise gives a relative accuracy in estimating the interline frequency of the order of a few ppm. However, when noise is present, the number of summed lines should be restricted to prevent the latter from affecting the calculation. The signal/noise ratio is at its maximum when only the maximum amplitude sample is retained for the sum in 22; the truncation error is then at its maximum since there is only one sample in the sum. A compromise must then be reached between the decrease in the truncation error and the contribution made by the noise when a spectral line is added. Such a compromise determines the algorithm's limitations. When a sinusoid with a slight modulation in amplitude or frequency is found, the calculation of the center of mass and the mean power of the component are sound, requires a minimal calculation effort, and delivers a result that is sometimes as valid as other more sophisticated algorithms which present the hypothesis of a continuous tone. In fact, when a tone that is not perfectly continuous is found, the interpolation approach cannot match a warped spectral main lobe.

To overcome the latter compromise between the truncation error and the noise contribution, we must take only one spectral line and take account of the relative weighting of this line over truncated lines. The solution calls for interpolation algorithm and has been published since mid 70's by Rife and Boorstyn [START_REF] Rife | Single-tone Parameter Estimation from Discrete-time Observations[END_REF]. One important thing about this paper it's his link to some other papers about the necessity of put the phase reference to the window center with IFFT. Rife and R. Boorstyn mentioned that the optimum value

t 0   N  1  / 2 
T "offers roughly a four-to-one reduction in the variance of the phase estimator over the natural choice of t 0  0 ". Using a different way, in a new approach for minimization of the phase error, Ferrero and Ottoboni find a similar conclusion [START_REF] Ferrero | A New Approach to the Fourier Analysis of Periodic Signals for the Minimization of the Phase Errors[END_REF]. About minimization of phase variance in the context of multifrequency signal parameter estimation with spectral windows, Offelli and Petri [START_REF] Offelli | The Influence of Windowing on the Accuracy of Multifrequency Signal Parameter Estimation[END_REF] wrote that the optimum value is n 0   1 2

 N with a "roughly ten-to-one variance reduction with respect to the usual choice n 0  0". Finally, in comparative study of IFFT algorithms, Schoukens, Pintelon and Van hamme ("h" minuscule sur hamme) demonstrated in their paper appendix [START_REF] Schoukens | The Interpolated Fast Fourier Transform: A Comparative Study[END_REF] that "the analytical formulation of rectangular window is considerably simplified if the time origin is chosen in the center of the rectangular window...". Moreover, these latter authors used this simplification in her maximum likelihood estimator (MLE).

Hydro-Québec, 1999

Phase referenced to window center 21 The IFFT algorithm presented (introduce?, dans le sens d'aborder partiellement) in this paper [START_REF] Léonard | Algorithme de haute-résolution spectrale: collé sur la limite ultime[END_REF] comprises three steps: a rough search, a coarse iterative search, and a final interpolation. Like many other IFFT algorithms, the interpolated frequency value is then used to interpolate the amplitude and phase. Consequently, a best frequency estimation leave/give? a best amplitude and phase estimate. The high-accuracy estimation algorithm uses the assumption that a continuous tone is present. This hypothesis limits to two the minimum number of spectral lines required for parameter estimation. Like many other algorithms, the rough search correspond to find(?) the indicia of the maximum amplitude spectral line. The coarse search consists of an iteration which maximizes the spectral correlation between the shape of the lobe of the component extracted from the signal and the shape of a synthesized lobe (Fig. 12). The correlation coefficient

 2  R l Q l  R l 2   Q l 2  , ( 23
)
between the magnitudes R l of the lobe of the measured signal and the magnitudes Q l of the synthesized lobe is calculated on few spectral lines of the largest magnitude. At each iteration, this correlation occurs for three separate frequencies: f i -df i , f i , f i +df i . The frequency with the strongest correlation amplitude becomes the central frequency f i+1 for the next correlation. The frequency range df i is reduced by half at each new iteration as illustrated in Hydro-Québec, 1999

Phase referenced to window center 22 Fig. 12. Consequently, the frequency uncertainty thus decreases by half at each iteration. At the beginning of iteration, the starting search range 2df 0 is set between 0.6 and 0.9 spectral lines and is centered on the line with the greatest amplitude. For instance, we can take df 0 =0.64 line to obtain an interline accuracy of 1% after six consecutive iterations.

Ultimately, however, i.e. after an infinite number of iterations, the maximum frequency excursion is 4df 0 of a spectral line. A crude frequency search with 22 can be used to substantially reduce the starting search range, i.e. the number of iterations required to obtain a given level of accuracy.

The success of this algorithm partly resides in the shape of the synthesized lobe. Ideally, the algorithm should be seen as an analog differential comparator where the input signals are processed in exactly the same way. This lobe may thus be synthesized using the FFT of a sinusoid with a frequency corresponding to f i -df i , f i or f i +df i and a phase of  i -d i ,  i or  i +d i . Thus, the iterative search is simultaneously performed on the frequency and phase since the phase has a slight influence on the lobe's magnitude profile. The same spectral window coefficients and FFT algorithm are then used on the same number of points.

Ultimately, the faults in centering the spectral window and the other systematic errors found in the processing of the signal are offset by the differential technique. Though it involves a prohibitive calculation time, the approach may be simplified by formulating two hypotheses:

-the magnitude of the main lobe is independent of the sinusoid phase;

-the shape of the main lobe only depends on the frequency interline position, without taking into account the value of the frequency expressed by the line number.

In fact, these hypotheses consist in disregarding the presence of an overlap between the positive and negative frequency components. Use of a spectral window with a substantial rejection of sidelobes allows these hypotheses to be met. Indeed, the relative error introduced with these hypotheses is then equivalent to the rejection value of the window's sidelobes. For field measurement, this error is often negligible in relation to the distortions of the lobe shape caused by fluctuations in amplitude and frequency of the tone.

As noted above, the hypothesis of a given magnitude which is not a function of the phase rules out the need for an iteration on the phase value. The iterative search is only done on the fractional part  of the frequency K+, where K corresponds to an integer number of the spectral line. If we replace the angular frequency  0 with 2(K+)/(NT) in ( 13) while disregarding the contribution of the negative frequencies, we obtain Hydro-Québec, 1999

Phase referenced to window center 23

S k,    0  a 0 2 exp jk 0 T  j 0   W k, T  2  N  2  N         (24)
such that

R l  S k,2  l / NT   a 0 2  W k, 2  l  2  2  N   (25)
presents the magnitude of the spectral lines under the main lobe with || small. In addition,

W k,2 N   W N / 2,2 N  , (26)
for a symmetrical window in relation to N/2: the imaginary part being null (nil?) in the kernel when the phase is referenced at N/2. Eq. 25 thus becomes

R l  S k,2  l / NT   a 0 2 W N / 2, 2  l  2  2  N   (27)
near the lobe summit since the contribution of the kernel found in the negative frequencies can be disregarded. Several approaches can be used to generate the shape of the lobe sampled by an FFT as a function of .

A numerical approach consists in realizing an FFT on a small number of sinusoid points of a frequency K'+ where the fractional part  corresponds to the desired position. K' is then a line number located far from the ends of the spectral band in order to minimize any overlapping of the positive and negative frequencies. Using this approximation, it is  the fractional part of the line position, which gives the shape of the lobe. To speed up the process, the shape of a spectral lobe between 0 and 0.5 spectral lines may also be tabulated [START_REF] Léonard | Algorithme de haute-résolution: version optimisée[END_REF]. Symmetrically, this table is used to locate the shape of a lobe from -1 to +1 spectral lines, i.e. throughout the search range. Lastly, the shape of the spectral lobe of Rife-Vincent windows class which consist of a cosine sum can be directly calculated using the mathematical expression stated in [START_REF] Lanczos | Discourse on Fourier Series[END_REF].

With the iterative process, each correlation is like a vertical cross-section of a mountain that gives the shape of the lobe as a function of the frequency . After 6 to 8 iterations (less than 1% of interline spacing), the surface progression is considered to be linear between the last two cross-sections. An interpolation leading to the final value may then occur.

Fig. 13 shows the last two lobe profile cross-sections realized by the iterative portion of the algorithm. The shape of the lobe to be identified is located somewhere between these two profiles. From the point of view of noise propagation, it is better to take into account the smallest possible number of lines at the top of the lobe. However, at least two spectral lines must be compared for the interpolation. Hence, interpolation used only two spectral lines for the fine search.

Q (f ) l 2 Q (f ) l 1 Q (f ) l+1 1 Q (f ) l+1 2 R l R l+1 f 2 f f 1 Iteration i+1 Iteration i Signal spectrum line l line l+1
Fig. 13: Linear interpolation of a lobe S i between two cross-sections found near the frequencies f 1 and f 2 .

Fig. 13 illustrates the case where the desired frequency f is located to the left of the spectral line of maximum amplitude Q l+1 (f 1 ) of the synthetic lobe with the strongest correlation and of frequency f 1 . Note that lines R l and R l+1 are not in the same scale as synthetic lines Q x (f x ).

The scale is not significant since only the relationship between the spectral lines of a given cross-section is used to calculate f. When the relationship R l /R l+1 and the relationships

Q l (f x )/Q l+1 (f x )
of the two adjacent synthetic lobes are known, one needs to determine the frequency f for which the interpolation is to be done. This frequency value is then used to interpolate the amplitude c Max which corresponds to the maximum-amplitude spectral line for a sinusoid with an amplitude of 1.0 located at the frequency f. The result will be the ratio S Max /c Max , i.e. the maximum of the spectrum over the synthetic maximum, which will determine the peak amplitude a 0 of the sinusoid described in (3).

The linearity hypothesis, where the spectral line amplitude ratio is a function of the difference in frequency, leads to the ratio

constant  f 2  f 1   Q l f 2   Q l 1 f 2    Q l f 1   Q l 1 f 1         . (20)
A simple rule of three then governs the interpolation of frequency Hydro-Québec, 1999

Phase referenced to window center 25

f  f 1  f 2  f 1   R l R l 1  Q l f 1   Q l1 f 1         Q l f 2   Q l1 f 2    Q l f 1   Q l 1 f 1         (21)
which is a function of the spectral line amplitude ratios. Note that the frequencies f, f 1 and f 2 may be line number units or expressed in hertz in the last equation. The interpolation of the spectral line of maximum amplitude is written

c Max  Q l 1 f 1   Q l 1 f 2   Q l 1 f 1      f  f 1 f 2  f 1       , at line number l+1, (22)
when Q 1+1 (f 1 ) corresponds to the amplitude of the maximum line of the synthetic lobe with the strongest correlation value. Note that the interpolation of c Max is different when the frequency is found to the right of the maximum line, in which case the following is used:

c Max  Q l f 1   Q l f 2   Q l f 1      f  f 1 f 2  f 1       , at line number l, (23)
since in this case Q 1 (f 1 ) corresponds to the line of maximum amplitude. The relationship R Max /c Max cannot be used directly without introducing a bias in the estimate of the amplitude. The fact that noise is found under the spectral lobe of the measurement causes the actual value of R Max to be overestimated. Statistically speaking, on a large number of measurements made using white noise that differs from one measurement to the next, it can be assumed that the energy of the noise under a spectral lobe is equivalent to half of the energy found under the straight line connecting the inflexion nodes R left and R right of the lobe (Fig. 14). The energy is adjusted to take into account the contribution of the noise such that

R Max 2  R X 2  R left 2  R right 2 4 with X  l,l  1   (24)
where R x corresponds to R i+1 or R i , depending on whether 22 or 23 is used. Lastly, the amplitude

a 0  R Max c Max N Q N S ( 25
)
associated with the tone is calculated while taking into account the number of samples N s of the signal's FFT and the number of samples N Q of the FFT which generates the synthetic lobe.

The phase interpolation

 0   l   l1   l   f  F l   F l 1  F l   , (26)
using  l and  l+1 , the phases below the spectral lines R 1 and R l+1 , is done using phases which do not cross the boundary line  /- (i.e. taking care to not mix -179° with +179°, in which case 360° is added to -179°). If the frequency f is expressed at the line number, then the frequencies F l and F l+1 will be in the same unit such that F l =l and F l+1 =l+1. It should be noted that the frequencies F l and F l+1 correspond to the frequencies of the spectral lines l and l+1, while frequencies f 1 and f 2 correspond to the cross-section frequencies related to the last two iterations. The interpolation described in [START_REF] Hanrahan | Estimation of Amplitudes, Frequency and Phase of Multi-tone Signal Components Using Low Sidelobe Time Windows[END_REF] thus takes place between two spectral lines while the interpolation of amplitude c Max takes place on a given spectral line between two iterations.

The test bench used to fine-tune the high-accuracy estimation algorithm is shown in Fig. 15.

It is used for two statistical calculations using the Monte-Carlo method. The first is used to calculate the variance on each estimated parameter, i.e.  f 2 ,   2 and  a 2 . To not introduce any bias in the calculation of the variance, a comparison is done of two estimations of a given sinusoid to which have been added different noise sequences prior to the FFT. The second statistical calculation provides the bias or accuracy using the mean spread between the estimated and input values. To increase the number of statistical tests, we have generated and Hydro-Québec, 1999

Phase referenced to window center 27 First of all, one must measure whether the Gaussian random-noise generator produces a noise that is sufficiently white. Note that the usual statistical generators "shuffle" the positions on the output sequence number for the random values in order to whiten the noise. Two random numbers are used for each random number requested. The first number determines which random number will be extracted from the array and the second random number replaces the first number extracted from the array. One must thus ensure that such a shuffler is found in the noise generator and we suggest that the whitening array be made larger than the FFT to be realized.

Hydro-Québec, 1999

Phase referenced to window center 28

We tested several versions similar to the algorithm presented above. For instance, we observed that although 24 eliminates any amplitude bias, it somewhat increases the dispersion. For the algorithm retained, no statistically significant bias can be observed on the frequency, amplitude or phase. With respect to the dispersion in response to a standard deviation  noise, the tested algorithm with the Blackman-Harris spectral window (-92dB)

allows us to move closer to the following bounds for a multiple-tone parameter estimation:

var f   9 2  2  a 0 2 T 2 N 2  1 PG (27)
for the variance in the estimation of the frequency in hertz (established empirically),

var a    2 PG ( 28
)
for the variance in the estimation of the amplitude, and

var      2 a 0 2 PG ( 29
)
for the variance in the estimation of the phase in radians, with the processing gain

PG  c i    2 c i 2  (30)
defined as the ratio of the coherent power gain over the noise power, with c i as the window time samples. Knowing that S N  a 0 2 2 2 , 27 and 29 can also be expressed as a function of the S/N ratio. Note that this bound is very similar to the Cramér-Rao bound (Eqs. 33-34).

The amplitude bound has been well established by Harris [START_REF] Harris | On the Use of Windows for Harmonic Analysis with the Discrete Fourier Transform[END_REF] and corresponds to the peakpicking of a spectral component whose frequency precisely matches a DFT spectral line. The phase bound is deduced from the amplitude bound. It is obvious that the amplitude and phase estimations come up against the same ultimate limit. In fact, a degree of uncertainty distributed isotropically in the complex plane, i.e. independent errors of the same magnitude in both axes, gives var     var a   a 0 2 with  in radians.

Different statistical estimates performed at 12,000 frequency, amplitude and phase give a variance that is twice that of the previous boundaries. The accuracy of the frequency and amplitude estimations is similar to that reported by others [START_REF] Ferrero | A New Approach to the Fourier Analysis of Periodic Signals for the Minimization of the Phase Errors[END_REF][START_REF] Offelli | The Influence of Windowing on the Accuracy of Multifrequency Signal Parameter Estimation[END_REF][START_REF] Schoukens | The Interpolated Fast Fourier Transform: A Comparative Study[END_REF][START_REF] Tabei | A Maximum Likelihood Estimator for Frequency and Decay Rate[END_REF][START_REF] Hanrahan | Estimation of Amplitudes, Frequency and Phase of Multi-tone Signal Components Using Low Sidelobe Time Windows[END_REF][START_REF] Rife | Multiple Tone Parameter Estimation from Discrete-Time Observations[END_REF][START_REF] Holm | Optimum FFT-Based Frequency Acquisition with Application to COSPAS-SARSAT[END_REF] and our bias estimator (Fig. 15) show no significative bias on frequency, amplitude or phase. On the other hand, we Hydro-Québec, 1999

Phase referenced to window center 29 believe that we are the only ones can use any kind of symmetrical windows. Also, the price of this iterative approach is the extensive amount of computation time.

Three factors contributed to the difference between the amplitude variance obtained and the amplitude lower bound. First, the processing gain 30 resulted from a calculation for the best case where =0, i.e. when there is an integer number of cycles in the observation window. The worst case corresponds to =0.5 and gives a processing gain of about 40% less.

Statistically, we obtain a mean value that is found between these two cases, i.e. an increase in the bound of about 20%. Second, the interpolation 22 of c Max involves use of the estimated frequency f which itself contains noise, thus substantially increasing the variance. Lastly, the adjustment for the amplitude bias 25 introduces an increase of 50% in the variance. These three values combined result in an increase of 85-95% of the actual value of the statistical bound given in 28 for an unbiased estimation of the amplitude of a sinusoid with a constant probability frequency density where   0.5, 0.5  . Seen from another perspective, it is as if the algorithm were using the information from about two spectral lines instead of only one, which ends up doubling the variance.

In our opinion, the best test consists in reconstructing the sinusoid using 3. On the one hand, this test involves doing a mix of the estimation errors in the actual frequency, amplitude and phase. It should be noted that these errors and theirs effect on reconstruction are not independent of each other: a parameter estimation error can be compensated or amplified by another error. On the other hand, combining these estimation errors gives a reconstruction error which is a function of time. In fact, it can be noted that the reconstruction error is not uniformly distributed throughout the time samples. Thanks to the algorithm's accuracy, we are able to reconstruct the sinusoid and subtract the sinusoid from the signal to the noise floor. The overall subtraction error is of similar magnitude as the noise and is at its smallest at the center of the time window.

Phase reference and Cramér-Rao bound

Like the bound expressed in Eqs. 27-29 for the Blackman-Harris window, the Cramér-Rao bound provides a lower limit for estimating the amplitude, frequency and phase error of a sinusoid when Gaussian white noise is found. However, this bound is the theoretical ultimate bound for a rectangular window and its demonstration begins in the time domain (see Rife and Boorstyn [START_REF] Rife | Single-tone Parameter Estimation from Discrete-time Observations[END_REF] for a basic explanation of the bound). The bound takes on different values according to whether certain parameters, such as the frequency, amplitude and phase, are Hydro-Québec, 1999

Phase referenced to window center 30 known in advance. For a noise with a standard deviation , a sinusoid with a peak amplitude a 0 , a sampling step T and a number N of samples, the variance in the frequency estimation will be

var f      2 4 2  a 0 2 T 2 k 2 N  2kP  Q   , (31)
with

P  n n 0 N 1   N N  1   2 , Q  n 2 n 0 N 1   N N  1   2N  1   6 and k  t R T , (32)
when the phase  is known and

var f   3 2  2  a 0 2 T 2 N N 2  1   ( 33
)
when it is not known, where t R corresponds to the point in time to which the estimated parameters are referenced. Similarly, the variance in the estimation of the phase, expressed in radians, is

var    f   2 a 0 2 N (34)
when the frequency is known and

var     12 2 k 2 N  2kP  Q   a 0 2 N 2 N 2  1   (35)
when the frequency is not known.

In almost all the applications of spectral estimation in a mechanical engineering context, none of the three parameters are known prior to being identified so that 33 and 35 predominate at the lower bounds. As pointed out by Rife and Boorstyn, the variance in the estimation phase expressed by 35 is at a minimum for

t R  N  1 2      T (36)
i.e., when the phase is referenced at the center of the observation window.

Hydro-Québec, 1999

Phase referenced to window center 31

We would like to introduce an important point here regarding the definition of the center of an observation window. On a signal with N-even samples, the center of the window is between two time samples. However, in a DFT, the window begin with sample #0 and ending at sample #N, excluding this latter. Hence, the center of the DFT is at N/2 rather than at (N-1)/2 (see annex A). Regarding the C-R bound, the demonstration of this bound begins in a temporal representation space, without the use of a DFT, which is why the center of the window corresponds to (N-1)/2 here.

We should, however, question Rife and Boorstyn's interpretation of 35; this variance combines both the estimate error in the phase parameter

   12 2 n  n c  n c  k   2 n0 N 1  a 0 2 N 2 N 2  1   . ( 39
)
After a few manipulations we obtain

var     12 2 n c  k   2  2 n c  k   n  n c   n0 N 1        a 0 2 N 2 N 2  1    12 2 n  n c   2 n0 N 1  a 0 2 N 2 N 2  1   (40)
where the second right-hand term corresponds to the variance in the estimate of the phase for k=n c , i.e. the variance expressed by 34.

Hydro-Québec, 1999

Phase referenced to window center 32

Since the sum

n  n c   n 0 N 1  is zero, 40 becomes var     12 2 n c  n R   2 a 0 2 N 2 N 2  1     2 a 0 2 N (41)
when it is expressed with a position referenced to the center of the observation window. This variance, expressed as var     var 2f

  T n R  n c     2  var    f , (42)
is equal to the sum of the variance in the error induced by the frequency estimation and the variance in the estimate of the phase parameter. Based on this last relationship, we can conclude that:

-Eq. 35 corresponds to 34 when k coincides with n c or when the frequency is known (i.e. var 2f   0).

-Eq. 35 does not yield the phase estimation error but rather the phase reconstruction error at position k; this error is a combination of both the estimation error in the phase parameter and the deviation induced by the frequency estimation error.

It is important to make a clear distinction between the phase parameter  =(n c T) and the reconstruction phase function (nT). The phase parameter is like the amplitude or frequency parameter, i.e. it is constant over the entire time window for an unmodulated tone. However, the phase parameter is related to a time reference; the frequency and phase parameters, on the other hand, have no time reference. The reason is that the phase parameter sets the time position for the entire sinusoid. We have come to the conclusion that the phase  of the tone must be referenced to the center of the window to eliminate the contribution of the frequency error. The phase function (nT) is generated using the phase and frequency parameters. The reconstruction error on the phase function (nT) is given by 42. This error is statistically minimal at the center of the time window and increases as one approaches the end of the time window.

Conclusion

Referencing a DFT phase in the center of an observation time window does not merely involve theoretical considerations but also has sizeable practical aspects. The phase plateau found under the spectral lobe substantially facilitates phase interpolation and reduces phase estimation error. From our interpretation of the Cramér-Rao bound and from our own experience with an high-accuracy estimation algorithm, like the experience related by different authors, referencing phase to the window center simplify mathematical expression and reduce phase standard deviation in multi-frequency signal parameter estimation. Moreover, put the phase reference at the window center reduce the probability of wrong path decision in the DFT phase unwrap process. Visually, z-plane and DFT phase graphic illustrations become more convenient: the phase referenced to the window center can be easily graphically read and interpreted.

Annex A

For the time window even sequence defined as

Fig. 1 :

 1 Fig. 1: Z-plane illustration of an increasing complex exponential of 3 cycles over 20 time samples data length (0.15 bin) located on radius 1.05, the radius on the pole circle.

Fig. 2 :

 2 Fig. 2: z-plane illustration of a sinusoid of 4.11 cycles over 20 time samples data length (0.2055 bin): a) z-amplitude, b) z-phase for k=0, c) z-phase for k=N/2, d) z-phase for k=N.

Fig. 3 :Fig. 4 :

 34 Fig. 3: Comparison of z-plane phase on unit circle and phase from zero padded DFT for phase referenced at the first time sample (k=0), and the corresponding unit circle unwrapped phase.

Fig. 5 :

 5 Fig. 5: 3D and 2D illustrations of the phase angle of a small portion of the spectrum

Fig. 6 :

 6 Fig. 6: Difference between a time delay on an observation window (a) and phase reference shift (b) of the DFT sequence for a constant tone at the  0 frequency.

Fig. 7 :

 7 Fig. 7: Typical illustration of the amplitude and phase of a spectral component of f 0 frequency estimated at the frequency f 0  , with the phase referenced at the beginning of the rectangular window (k=0) and at the center of the window (k=(N-1)/2).

Fig. 9 :

 9 Fig. 9: DFT phase with bellshaped windows, for a reference at the beginning of the window (k=0) and at the window's center (k=N/2).

Fig. 10 :

 10 Fig.10: Power spectrum and unwrapped phase of measured acceleration on a cantilever metal beam with some resonant vibration modes excited by a mechanical shock.

Fig. 11 :

 11 Fig. 11: Zero padding technique for phase referenced to the window center (a) and for phase referenced to the first time sample (b).

Fig. 12 :

 12 Fig. 12: Flow chart of the iterative search algorithm for the frequency.

Fig. 14 :

 14 Fig. 14: Inflexion nodes R left and R right on the spectral lobe.

Fig. 15 :

 15 Fig. 15: Test bench used to fine-tune and set the algorithm's final detection boundary.

   and the deviation induced by the frequency estimate error. In fact, looking once again at 35 and introducing the values relative to P and Q, introduce the position of the central sample in the window n c =(N-1)/2 such that var 

Fig. A. 1 :

 1 Fig. A.1: Position and corresponding span for a time windows symmetrical with respect to (N-1)/2 and respect to N/2, for an N time samples sequence.

Acknowledgements

The author thanks Georges Desrochers and Innocent Kamwa for the many consultations they offered as signal processing experts.