Georges-Pierre Bonneau
email: pierre.bonneau@imag.fr

Multiresolution Analysis with Non-Nested Spaces

Two m ultiresolution analyses (MRA) intended to be used in scienti c visualization, and that are both based on a non-nested set of approximating spaces, are presented. The notion of approximated r e nement is introduced to deal with non nested spaces. The rst MRA scheme, referred to as BLaC (Blending of Linear and Constant) wavelets is based on a one parameter family of wavelet bases that realizes a blend between the Haar and the linear wavelet bases. The approximated re nement is applied in the last part to build a second MRA scheme for data de ned on an arbitrary planar triangular mesh.

INTRODUCTION

The purpose of this paper is to present t wo m ultiresolution analyses (MRA) intended to be used in scienti c visualization, and that are both based on a non-nested set of approximating spaces. The need for non nested spaces arises from the fact that the required scaling functions do not ful ll any re nement equation. Therefore we i n troduce in the rst part the concept of approximated re nement equation, that allows to generalize the lter bank and exact reconstruction algorithms. The second part shows how t h i s concept enables to de ne a one parameter family of wavelet bases that realizes a blend between the Haar and the linear wavelet bases. These so-called BLaC (Blending of Linear and Constant) wavelets are used for the analysis of volume data de ned on uniform grids. The e ect of the blending parameter is shown on examples. In the last part of the paper, the approximated re nement is applied to built a MRA scheme for data de ned on an arbitrary planar triangular mesh. The ability to deal with arbitrary triangular meshes, without subdivision connectivity, can be achieved only through the use of non nested approximating spaces, as introduced in the rst part.

Multiresolution analysis with non nested approximating spaces

We assume that the reader is familiar with the basic concepts of multiresolution analysis. An introduction on this topic can be found in 4]. In this section, a general framework for constructing multiresolution analysis schemes based on a non nested sequence of approximating spaces is given.

Necessary conditions on the sequence of approximating spaces

The minimal requirements for this construction are the following: If L 2 ( ) is the set of all functions with nite energy de ned on a domain , the sequence V n of approximating subspaces of L 2 ( ) must satisfy:

V n is isomorph to a subspace of V n+1 (1)
V n is dense in L 2 ( ):

(2) Condition [START_REF] Bonneau | Wavelets: a multiresolution analysis with non-nested spaces[END_REF] means that the sequence of approximating spaces V n is \growing" in L 2 ( ). This condition replaces the nested condition V n V n+1 in the de nition of the classical multiresolution analysis schemes. If the spaces V n have nite dimensions, this condition may be written as dimV n dimV n+1 . Condition [START_REF] Eck | Multiresolution Analysis of Arbitrary Meshes[END_REF] implies that every function can be approximated with an arbitrary precision by functions in V n for some value n.

Approximated re nement equation

Let (' n j ) denote the scaling functions (i.e. (' n j ) j is a basis of V n ). If the nested condition V n V n+1 is satis ed, the link between two successive l e v els of resolution is provided by a so-called re nement equation, which expresses each scaling function ' n j as a linear combination of the ner scaling functions (' n+1 i ). In our setting, condition (1) implies that there exists functions ( 'n j ) j which form a basis of a space Ṽn , included in V n+1 . And the link between level n and n + 1 is now p r o vided through a two step process: First apply the isomorphism, i.e. replace the scaling functions (' n j ) j by t h e i n termediate functions ( 'n j ) j , and then apply the re nement equation which expresses 'n j as a linear combination of the ner scaling functions (' n+1 i ) i :

' n j $ 'n j = X i p n ij ' n+1 i : (3)
(3) is called the approximated r e nement equation because in the two next sections, where the framework is applied, 'n j is closed to ' n j .

Detail spaces and wavelet functions

The detail spaces W n are now de ned as a complementary of the space Ṽn in V n+1 :

V n+1 = Ṽn W n :
The wavelet functions ( n j ) j are basis functions of the space W n . T h us there exist coe cients (q n ij ) s u c h t h a t n j = X i q n ij ' n+1 i and coe cients (a n ij ), (b n ij ) s u c h that

' n+1 j = X i a n ij 'n i + X i b n ij n i :
From now, we assume that the approximating spaces V n have nite dimensions. Let us call P n , Q n , A n and B n the matrices with coe cients (p n ij ), (q n ij ), (a n ij ) a n d ( b n ij ). P Q are called synthesis matrices and A B are called analysis matrices. The columns of the matrices P and Q contain the coordinates of the intermediate functions 'n j and the wavelets n j in the basis of the ner scaling functions ' n+1 i . The rows of the matrices A and B contain the coordinates of the dual of the intermediate functions, and the dual of the wavelets, in the dual basis of the ner scaling functions. These matrices must satisfy the reconstruction condition P A+ QB = I : 2 A u t h o r p r e p r i n t Furthermore let G n+1 be the Gram-Schmidt matrix of the basis (' n+1 i ) (i.e. the matrix whose elements are < ' n+1 i ' n+1 j >).

Then the multiresolution analysis is called semi-orthogonal i P T G n+1 Q = 0 AG ;1 n+1 B T = 0 : (4) In this case P and Q can be calculated respectively from A and B by P = G ;1 n+1 A T (AG ;1 n+1 A T ) ;1 Q = G ;1 n+1 B T (BG ;1 n+1 B T ) ;1 :

(5) orthogonal i in addition to (4), one has P T G n+1 P = I Q T G n+1 Q = I and in this case P and Q are given by P = G ;1 n+1 A T Q = G ;1 n+1 B T : bi-orthogonal i BP = 0 AQ = 0 : One main di erence with the classical multiresolution analysis framework, is that the scaling functions (' n j ) are not the result of a subdivision scheme. They are not xed by the choice of the coe cients (p n ij ), and this is the reason why the Gram-Schmidt matrix G n+1 appears in the orthogonality equations.

Analysis and reconstruction algorithms

Given an approximation f n+1 = P i x n+1 i ' n+1 i of f in V n+1 , the next coarser approximation f n = P i x n i ' n i and the details g n = P i y n i n i are computed by (

x n i ) = A(x n+1 i ) (y n i ) = B(x n+1 i ) : (6)
The inverse transformation is given by (x n+1 i ) = P(x n i ) + Q(y n i ) :

(7) The analysis and reconstruction formulas (6) and (7) are the same as in the case when the nested condition V n V n+1 is satis ed, but they are related now to an approximation in two steps: 1 gives a geometric interpretation of one step of analysis and synthesis in the classical and new multiresolution framework, for the semi-orthogonal case.

Analysis: f n+1 approx. ;! e f n = P i x n i e ' n i isomorph. ! f n = P i x n i ' n i & g n Figure
f n+1 V n g n f n V n+1 V n V n fn ~fn+1 V n+1 g n f n f n f n+1 g n f n+1 f n g n f n Figure 1:
The classical and new multiresolution framework (semi-orthogonal case).

3

A u t h o r p r e p r i n t

Choice of the intermediate functions 'n j

A crucial point is the choice of the intermediate functions 'n j (i.e. the choice of the matrix P for the approximated re nement equation). In Section 3, this choice depends on the application, since 'n j is chosen to be close to ' n j , and so that the approximated re nement equation (3) blends between the Haar and linear re nement equations. In Section 4, the choice of 'n j is independent of the application, since it is chosen so that f n is always the best L 2 -approximation of f n+1 in V n . T h i s i s a c hieved i the smooth coe cients (x n i ) are the solution of the following linear system:

G n (x n i ) i = ( < ' n i X j x n+1 j ' n+1 j >) i
and this implies that A = G ;1 n (< ' n i ' n+1 j >) i j : (8) Since in Section 4, the multiresolution analysis veri es the semi-orthogonality conditions (4), the matrix P (and hence the functions 'n j ) are xed by [START_REF] Schr Oder | Spherical wavelets: E ciently representing functions on the sphere[END_REF] .

Therefore in the semi-orthogonal case, the choice (8) for the matrix A and ( 5) for the matrix P implies that f n is always the best approximation of f n+1 in V n . N e v ertheless, since the sequence of approximation spaces V n is not nested, this best approximation property does not hold anymore between f n and f n+k if k > 1. This implies theoretically that very coarse approximations could be far away from best approximations of the original data, and thereby inducing stability problems during the analysis. For the two following applications of sections 3 and 4, we didn't encount s u c h stability problems.

In conclusion, this section has introduced a theoretical framework, that can be used to handle practical situations where the scaling functions arising naturally from the application do not ful ll any re nement equation.

BLaC Wavelets

In this section, the previous results are applied to de ne a one-parameter family of wavelets that perform a blend between Haar wavelets and linear wavelets. The purpose of BLaC wavelets is to nd a compromise between the locality of the analysis (which is perfect for the Haar wavelets) and the regularity o f t h e approximations (which is better for the linear wavelets). Therefore BLaC stands for Blending of Linear and Constant. More details on the motivation for the de nition of BLaC wavelets can be found in 1].

In Section 3.1 the BLaC scaling functions and the approximated re nement equation, as de ned in the previous section, are introduced. In Section 3.2, the BLaC wavelets, assuming semi-orthogonality and minimal support, are computed. The BLaC scaling and wavelet functions, and the approximated re nement equations depend on a blending parameter 2 0 1], and are calculated so that for = 0, the MRA scheme is equal to the Haar scheme, and for = 1, to the linear scheme.

BLaC scaling functions and approximated re nement equation

The BLaC scaling function ' n i is built by dilation by a factor 2 n and translation of (i2 ;n ) of the function ' shown in gure 2: ' n i (x) = ' ; 2 n (x + i 2 ;n ) :

A u t h o r p r e p r i n t

Since we w ant our algorithms to work on data de ned on bounded domains, we select the indices ;1 : : : 2 n , and take the restriction to the interval 0 1] of the rst and last scaling functions, ' n ;1 and ' n 2 n . Figure 3 shows the scaling functions at the level 2 of resolution. The lack of symmetry between the rst and last scaling functions is due to the fact that there are 2 n ;1 Haar scaling functions, and 2 n linear scaling functions whose support intersects the interval 0 1]. From Figure 4 it is clear that ' n i can not be written as a linear combination of the ner scaling functions ' n+1 i : in the dashed area, ' n i has a slope, while all functions ' n+1 i are constant. Thus the new multiresolution analysis framework de ned in the previous section has to be applied. Figure 4: Need for non-nested approximating spaces. The intermediate functions 'n i must be chosen so that the approximated re nement equation ( 3) is a blending between the Haar and the linear re nement equations, and so that 'n i is close to ' n i . For Haar scaling functions, the exact re nement equation is ' n i = ' n+1 2i + ' n+1 2i+1 , a n d i n t h e l i n e a r c a s e i t i s ' n i = 1 2 ' n+1 2i + ' n+1 2i+1 + 1 2 ' n+1 2i+2 . Therefore we search for the following approximation of ' n i :

'n i = ' n+1 2i + ' n+1 2i+1 + ( 1 ; )' n+1 2i+2 k' n i ; 'n i k L 2 ( 0 1]) ! min : (9)
The constant depends only on the blending parameter . As expected, is equal to 1 for = 0 and is equal to 1=2 for = 1, leading respectively to the Haar and linear re nement equations. Figure 5 illustrates the approximated re nement equation.

A u t h o r p r e p r i n t

The matrix P is now xed as P = 2 6 6 6 6 6 6 6 6 4

BLaC wavelet functions

We impose two constraints on the BLaC wavelet functions ( n i ): they must be orthogonal to all intermediate functions ( 'n i ) (i.e. the semi-orthogonality condition (4) must hold) and they must have minimal support. n i is a linear combination of the scaling functions (' n+1 j ) j . Adding a function ' n+1 j in this linear combination gives one more degree of freedom to determine the wavelet, but also increases the support of the wavelet. The minimal degree of freedom necessary to ensure the semi-orthogonality is 5, i.e.:

n i = a' n+1 2i + b' n+1 2i+1 + c' n+1 2i+2 + d' n+1
2i+3 + e' n+1 2i+4 : Indeed, in this case the support of the function n i intersects only the support of four basis functions of Ṽn , and the semi-orthogonality condition reduces to: Figure 6: The BLaC Wavelets On the boundaries of the interval, the rst and last wavelets are linear combinations of only four scaling functions: and the semi-orthogonality condition P T G n+1 Q = 0 holds.

< n i 'n j >= 0 j = i ; 1 i i + 1 i + 2 :
n ;1 = b 0 ' n+1 ;1 + c 0 ' n+1 0 + d 0 ' n+1 1 + e 0 ' n+1 2 n 2 n = a 1 ' n+1 2 n+1 ;4 + b 1 ' n+1 2 n+1 ;3 + c 1 ' n+1 2 n+1 ;2 + d 1 ' n+1 2

Examples

Figure 7 shows the result of the partial reconstruction of a 2D (512 512) data set, with 1% of the detail coe cients, for three di erent v alues of . On the top left is the original data. On the top right w e used = 0 (i.e Haar analysis), on the bottom left = 0:5, and on the bottom right = 1 (i.e. linear analysis) was used.

Because linear wavelets are more regular than Haar wavelets, the result at the bottom right l o o k s smoother. On the other hand, the partial reconstruction using Haar wavelets at the top right is more sharp. The result of the reconstruction with = 0:5 at the bottom left clearly shows that BLaC wavelets enables to nd a compromise between the regularity and the locality of the analysis, or analogously between the sharpness and the smoothness.

Original data = 0 :0(Haar) = 0 :5 = 1 :0(Linear) A u t h o r p r e p r i n t

Multiresolution analysis on hierarchical triangulations

Large unstructured data sets originating from meteorological or geological measurements often need a special treatment because they are very expensive in storage and time. Although there exist a lot of di erent wavelet and multiresolution analysis constructions, they all are restricted to simple domains with subdivision connectivity. One exception is in 2], where the approach i s t o a p p r o ximate the given domain by a simpler one which has subdivision connectivity. For data de ned on uniform grids, domains with subdivision connectivity are constructed by 2-split, 4-split or 8-split of intervals, squares or cubes. Previous approaches on multiresolution analysis of data de ned on triangular meshes involve domains constructed by recursive 4-split of triangles [START_REF] Lounsbery | Multiresolution analysis for surfaces of arbitrary topological type[END_REF][START_REF] Stollnitz | \Wavelets for Computer Graphics, Theory and Applications[END_REF][START_REF] Schr Oder | Spherical wavelets: E ciently representing functions on the sphere[END_REF].

The mathematical reason for the need of the subdivision connectivity for MRA relies in the need of a nested set of approximation spaces. Indeed, the nested set condition [START_REF] Bonneau | Wavelets: a multiresolution analysis with non-nested spaces[END_REF], or analogously the re nement equation implies that the support of each scaling function at one level of resolution must be the union of the supports of some ner scaling functions. Since for triangular multiresolution, the support of the scaling functions are the triangles of the mesh, it follows that the nested set condition (1) implies the subdivision connectivity of the mesh over which the data is de ned. The development of Section 2 allows to construct multiresolution analysis on more general hierarchical triangulations, which don't have subdivision connectivity.

In Section 4.1. we describe the hierarchical triangulation structure that we c hose for our multiresolution analysis. In Section 4.2. we explain how the theoretical framework introduced in Section 2 can be applied to nd scaling and wavelet functions for such a hierarchical structure. The corresponding decomposition and reconstruction algorithms are introduced in Section 4.2. . Examples are shown in Section 4.3. .

The hierarchical triangulations and the approximation spaces

The MRA scheme de ned in the next section is dedicated to the analysis of piecewise constant functions de ned on an irregular triangular planar mesh. Such functions can be the result of a nite element m e s h calculus. At the nest level of resolution N, the domain is the disjoint union of the triangles T N k : = k2K(N) T N k and the nest approximation space V n is the set of functions de ned on and constant o ver the triangles T N k .

The triangulation of the domain at the level n is obtained by removing vertices from the triangulation at the level n+1 and re-triangulating the region equal to all triangles adjacent to the removed vertices. Such a construction de nes disjoint polygonal regions R n j in the domain , one for each set of adjacent removed vertices (see Figure 8 At the level of resolution n, the approximation space V n is the set of functions de ned on and constant on the triangles T n k .

No assumption is made on the order of suppression of the vertices, and on the algorithm for retriangulating the polygonal region. In particular this paper does not address the problem of choosing an order of suppression of the vertices. It assumes that this order is known. This paper solves the problem of de ning scaling and wavelet functions, together with smooth and detail coe cients, based on the knowledge of the hierarchy of triangulations.

In the current implementation the hierarchy is encoded with a graph structure. This graph has N levels, one for each level of resolution (i.e. one for each triangulation of the domain ). The nodes of the graph at the level n are the polygonal regions R n j . E a c h n o d e R n j contains two main lists of pointers: one for the child triangles T n+1 k , and one for the parent triangles T n k . F urthermore, a node R n j points onto a node R m j 0 at a coarser level m n, i f R m j 0 contains a child triangle which is a parent triangle of R n j (see Figure 9).

In the special case where the successive triangulations are constructed by 4 split, then each n o d e h a s only one parent triangle and four child triangles, and the graph structure reduces to a forest of triangle quadtrees (like the one used in 5].

Figure 2 :

 2 Figure 2: BLaC scaling function 4

Figure 3 :

 3 Figure 3: The BLaC scaling functions at level 2 of resolution.From Figure4it is clear that ' n i can not be written as a linear combination of the ner scaling functions ' n+1 i : in the dashed area, ' n i has a slope, while all functions ' n+1 i are constant. Thus the new multiresolution analysis framework de ned in the previous section has to be applied.

Figure 5 :

 5 Figure 5: BLaC re nement equation On the boundaries of the interval, the rst and last approximated re nement equations are given by ' n ;1 $ 'n ;1 = ' n+1

  is a linear homogeneous system of four equations in the ve u n k n o wns a b c d e. It is uniquely solved if we impose the L 2 -norm of the wavelet to be equal to 1: R j n i j 2 = 1. Note that the coe cients a b c d e are also independent of the resolution level n. Figure6shows one BLaC wavelet for various values of .

  n+1 ;1 : b 0 c 0 d 0 e 0 a 1 b 1 c 1 d 1 are calculated by i n verting the following linear homogeneous systems

Figure 7 :

 7 Figure 7: Partial reconstructions with BLaC wavelets. (1% of the coe cients) 7

Figure 8 :

 8 Figure 8: Polygonal regions de ned by v ertices removal. The triangle T n k are called the parent triangles of R n j , and the triangles T n+1 k are called child triangles of R n j .At the level of resolution n, the approximation space V n is the set of functions de ned on and constant on the triangles T n k .No assumption is made on the order of suppression of the vertices, and on the algorithm for retriangulating the polygonal region. In particular this paper does not address the problem of choosing an order of suppression of the vertices. It assumes that this order is known. This paper solves the problem of de ning scaling and wavelet functions, together with smooth and detail coe cients, based on the knowledge of the hierarchy of triangulations.In the current implementation the hierarchy is encoded with a graph structure. This graph has N levels, one for each level of resolution (i.e. one for each triangulation of the domain ). The nodes of the graph at the level n are the polygonal regions R n j . E a c h n o d e R n j contains two main lists of pointers: one for the child triangles T n+1 k , and one for the parent triangles T n k . F urthermore, a node R n j points onto a node R m j 0 at a coarser level m n, i f R m j 0 contains a child triangle which is a parent triangle of R n j (see Figure9).In the special case where the successive triangulations are constructed by 4 split, then each n o d e h a s only one parent triangle and four child triangles, and the graph structure reduces to a forest of triangle quadtrees (like the one used in 5].

3 Figure 9 :

 39 Figure 9: The hierarchical triangulations and the graph structure.

A u t h o r p r e p r i n t

Analysis and synthesis matrices

The scaling function ' n k is the function equal to 1 on the triangle T n k , and 0 otherwise:

Outside the polygonal regions R n j de ned in Section 4.1, the multiresolution analysis does nothing between level n and n+1: the scaling functions ' n k , ' n+1 k are the same, the smooth coe cients (x n k ), (x n+1 k ) also, and no detail coe cients need to be generated. It is therefore su cient to compute the analysis and synthesis matrices P Q A B locally to each region R n j . Let R n j = k2L(n+1 j)

T n k be such a region. A is rst computed, so that f n is always the best L 2 -approximation of f n+1 in V n , as explained in Section 2.5: A = G ;1 n < ' n i ' n+1 j >. G n is here a diagonal matrix whose i-th diagonal element is the area of the triangle T n i , a n d < ' n i ' n+1 j > is equal to the area of the intersection of the triangles T n i and T n+1 j . T h us A is the following (#L(n + 1 j ) #L(n j)) matrix:

area(T n i ) :

In the simple 4 split case, A is the row matrix ( 1

4 ). B is calculated so that the dual scaling and wavelet functions are orthogonal: AG ;1 n+1 B T = 0 (11) and so that the dual wavelet functions are orthonormal:

BG ;1 n+1 B T = I : This is done in two steps: rst x a submatrix of dimension #L(n j) #L(n j) o f B to be the identity and solve (11). Then apply a Gram-Schmidt algorithm to orthonormalize the rows of B. This choice for the analysis matrices A and B ensures that the approximation f n at the level n is the best L 2 approximation of f n+1 , the intermediate approximation fn at the level n is also the best L 2 approximation of f n+1 . Furthermore, this choice for the analysis matrices implies that the synthesis matrices P and Q can be calculated from A and B by

where G ;1 n+1 is here the diagonal matrix with i-th diagonal element equal to the area of the triangle T n+1 i .

Figure 10 shows the decomposition of a function f n+1 in a region R n j containing four child triangles T n+1 k and two p a r e n t triangles T n k . This region is created by the removal of one interior point. The decomposition is a two step process, as explained in Section 2.4: rst apply the analysis matrices to nd fn and the two w avelet components y n 0 n 0 , y n 1 n 1 . Then apply the isomorphism, i.e. replace the intermediate scaling functions 'n i by the scaling functions ' n i .

10

A u t h o r p r e p r i n t

Triangular decomposition and reconstruction algorithms

The graph structure described in Section 4.1. , and the analysis and synthesis matrices introduced in Section 4.2. enable to implement the three following algorithms.

Algorithm 1: Triangular wavelet decomposition. This algorithm starts with the values (x N i ) of a function f on the nest triangles T N k . Then for each level of resolution n and for each polygonal region R n j (as de ned in Section 4.1. ), the analysis matrices A and B introduced in Section 4.2. are used to compute the smooth coe cients (x n i ) and the detail coe cients (y n i ). These coe cients are stored in memory in place of the corresponding ner scaling coe cients (x n+1 i ). At the end remains only the scaling coe cients (x 0 i ) corresponding to the approximation of f at the lowest level of resolution, and the detail coe cients (y n i ) corresponding to all regions (R n j ). Algorithm 2: Reconstruction at a xed level n. This algorithm starts with the output of the wavelet decomposition algorithm. For each level of resolution k lower than n, and for each region R k j , the synthesis matrices P and Q are used to recover the scaling coe cients (x k+1 i ) from the coarser coe cients (x k i ) and the details (y k i ). Algorithm 3: Threshold reconstruction. This algorithm starts with the output of the wavelet decomposition algorithm, and with a threshold value . This algorithm select the subgraph of the graph structure introduced in Section 4.1. , which contains all regions R n j that have an associated detail coe cient whose absolute value is greater than the threshold value . Then the synthesis matrices P and Q are applied on all these selected regions to recover the scaling coe cients (x n+1 i ) from the coarser coe cients (x n i ) and detail coe cients (y n i ) t h a t have an absolute value greater than .

Examples

The hierarchical triangulations corresponding to the examples of Figure 11 were obtained by inserting randomly chosen vertices in the square domain = 0 512] 0 512], and constructing their Delaunay triangulation. The coarsest triangulation of the domain consists of two triangles covering it. The values of the test function at the nest level of resolution are the mean values of a 512 512 image on the nest triangles.

11

A u t h o r p r e p r i n t

The top picture of Figure 11 shows the original data, consisting of 48682 triangles together with their function value (a gray level between 0. and 256.).

The middle left picture of Figure 11 shows the result of the reconstruction of the data at a xed level (algorithm 2 of Section 4.3.), using approximatively 40% of the detail coe cients. The bottom left picture shows the corresponding triangulation.

The middle right picture shows the result of a threshold reconstruction (algorithm 3 of Section 4.3.) using approximatively the same number of detail coe cients. The corresponding triangulation (bottom right picture) illustrates that this algorithm selects more vertices in the areas of strong variations of the test function.

Original data set


Bernard Picinbono
email: bernard.picinbono@lss.supelec.fr

Measurements of Second-Order Properties of Point Processes

The second order statistical properties of point processes are described by the coincidence function which can be measured by a coincidence device but such measurements are long and complicated. We propose another method of measurement and we analyze its performances. The starting point is that the coincidence function can be deduced from the probability density functions of the life times (the distances between points) of the process. The idea is to transform the point process into a positive signal whose values are these distances. From an appropriate processing of this signal we deduce the coincidence function. For the validation of the method we use point processes for which the coincidence function is known. The agreement between theory and experiment is in general excellent. Finally the method is applied to measure the coincidence functions of some point processes for which no theoretical result is available.

Index Terms

Point processes, statistical measurements, signal processing, signal representation.

I. INTRODUCTION

Point processes (PP) play an important role in many areas of physics, statistics, and engineering sciences. For example in optical communications at low levels the only available information is the set of random instants at which photons are detected. Similarly the time instants at which telephone calls arrive at a switching center is a PP.

A PP is a random distribution of points in a space. If this space is the time axis we have a time PP and its random points P i are time instants t i sometimes called events. All the PPs considered below are time PPs. In what follows we consider only stationary PPs characterized by the fact that their statistical properties are invariant in time.

The complete statistical description of a PP is very complicated and in most applications one is obliged to use only first-or second-order statistical properties. The first-order description of a stationary PP is contained in its density λ which is the average number of points per unit of time. The transposition to PPs of the concepts of the correlation function requires a specific analysis. The first attempt to describe second-order properties of PPs was presented by Bartlett [START_REF] Bartlett | The spectral analysis of point processes[END_REF], [START_REF] Bartlett | An Introduction to Stochastic Processes[END_REF]. Since a PP has no correlation function, the idea was to introduce a similar function which was called the correlation density function. More recently it was indicated that this function also appears in the description of a PP by coincidence analysis [START_REF] Macchi | The coincidence approach to stochastic point processes[END_REF] and this is the reason to use the term coincidence function [START_REF] Picinbono | Coincidence analysis of point processes[END_REF] which contains all the second-order properties of a PP.

The coincidence function can be measured by a coincidence device, system widely used in Nuclear Physics. The measurements however are long and complicated and not adapted to various PPs. We propose another approach which does not start directly from the definition of the coincidence function but from its relation to the properties of the distances between successive points of the PP sometimes called life times [START_REF] Kryukov | Calculation of the correlation function and the spectral power density of random sampling[END_REF], [START_REF] Beutler | The spectral analysis of impulse processes[END_REF], [START_REF] Cox | Point Processes[END_REF]. Using the standard vocabulary of experimental nuclear physics, the device in order to measure the coincidence function is then no longer a coincidence circuit but a time to amplitude converter (TAC) which transforms distances between points into pulses of amplitude proportional to these distances.

These amplitudes constitute a discrete-time positive signal. By an appropriate statistical signal processing of these distances it it is possible to measure the coincidence function.

Coincidence function and coincidence measurements are at the basis of some applications in statistical optics and more lately in quantum communication [START_REF] Picinbono | Coincidence analysis of point processes[END_REF], [START_REF] Picinbono | Characterization of nonclassical optical fields by photodetection statistics[END_REF], [START_REF] Saleh | Photoelectron Statistics[END_REF]. For example the structure of the coinci-DRAFT February 13, 2018

dence function is a tool to decide whether or not an optical field is of classical or of non-classical nature [START_REF] Picinbono | Characterization of nonclassical optical fields by photodetection statistics[END_REF].

In Section II we present a short review of properties of coincidence functions and we describe the principle of their measurement from the analysis of life times. We discuss theoretically the precision of the method and its limitations. In Section III we test the method by using computer experiments. There are some PPs for which the coincidence function can be calculated theoretically. By making measurements on those PPs we can validate the behavior of the method and evaluate its precision. In general the experimental measurements are in excellent agreement with theoretical results. In the last section we apply the method in the case of some PPs for which the calculation of the coincidence function is almost impossible and we present and discuss some examples.

Point processes (PP) play an important role in many areas of physics, statistics, and engineering sciences. For example in optical communications at low levels the only available information is the set of random instants at which photons are detected. Similarly the time instants at which telephone calls arrive at a switching center is a PP.

A PP is a random distribution of points in a space. If this space is the time axis we have a time PP and its random points P i are time instants t i sometimes called events. All the PPs considered below are time PPs. In what follows we consider only stationary PPs characterized by the fact that their statistical properties are invariant in time.

The complete statistical description of a PP is very complicated and in most applications one is obliged to use only first-or second-order statistical properties. The first-order description of a stationary PP is contained in its density λ which is the average number of points per unit of time. The transposition to PPs of the concepts of the correlation function requires a specific analysis. The first attempt to describe second-order properties of PPs was presented by Bartlett [START_REF] Bartlett | The spectral analysis of point processes[END_REF], [START_REF] Bartlett | An Introduction to Stochastic Processes[END_REF]. Since a PP has no correlation function, the idea was to introduce a similar function which was called the correlation density function. More recently it was indicated that this function also appears in the description of a PP by coincidence analysis [START_REF] Macchi | The coincidence approach to stochastic point processes[END_REF] and this is the reason to use the term coincidence function [START_REF] Picinbono | Coincidence analysis of point processes[END_REF] which contains all the second-order properties of a PP.

The coincidence function can be measured by a coincidence device, system widely used in Nuclear Physics. The measurements however are long and complicated and not adapted to various PPs. We propose another approach which does not start directly from the definition of the coincidence function but from its relation to the properties of the distances between successive points of the PP sometimes called life times [START_REF] Kryukov | Calculation of the correlation function and the spectral power density of random sampling[END_REF], [START_REF] Beutler | The spectral analysis of impulse processes[END_REF], [START_REF] Cox | Point Processes[END_REF]. Using the standard vocabulary of experimental nuclear physics, the device in order to measure the coincidence function is then no longer a coincidence circuit but a time to amplitude converter February 13, 2018 DRAFT (TAC) which transforms distances between points into pulses of amplitude proportional to these distances.

These amplitudes constitute a discrete-time positive signal. By an appropriate statistical signal processing of these distances it it is possible to measure the coincidence function.

Coincidence function and coincidence measurements are at the basis of some applications in statistical optics and more lately in quantum communication [START_REF] Picinbono | Coincidence analysis of point processes[END_REF], [START_REF] Picinbono | Characterization of nonclassical optical fields by photodetection statistics[END_REF], [START_REF] Saleh | Photoelectron Statistics[END_REF]. For example the structure of the coincidence function is a tool to decide whether or not an optical field is of classical or of non-classical nature [START_REF] Picinbono | Characterization of nonclassical optical fields by photodetection statistics[END_REF].

In Section II we present a short review of properties of coincidence functions and we describe the principle of their measurement from the analysis of life times. We discuss theoretically the precision of the method and its limitations. In Section III we test the method by using computer experiments. There are some PPs for which the coincidence function can be calculated theoretically. By making measurements on those PPs we can validate the behavior of the method and evaluate its precision. In general the experimental measurements are in excellent agreement with theoretical results. In the last section we apply the method in the case of some PPs for which the calculation of the coincidence function is almost impossible and we present and discuss some examples.

II. DEFINITION AND MEASUREMENT OF THE COINCIDENCE FUNCTION

Let N (t) be the number of points of a PP in the interval [0, t[, where 0 is an arbitrary origin of time.

The increment dN (θ) = N (θ + dθ) -N (θ) is the number of points in the interval [θ, θ + dθ[. The second-order properties of a PP are characterized by its coincidence function c(t) defined by

E[dN (θ)dN (θ )] = c(θ -θ )dθdθ . (1)
The name comes from the fact that this function is generated by a coincident events. Indeed a coincidence at θ and θ is the event defined by [dN (θ) = 1] ∩ [dN (θ ) = 1]. For regular PPs, which are the only ones studied in this paper, the increment dN (θ) takes only the values 0 or 1 when dθ tends to 0.

As a consequence the expectation value appearing in (1) is the coincidence probability. The function b(t) = c(t)/λ, where λ is the density or the PP, is called the bunching function because it describes the bunching effect appearing in the PP obtained in photodetection [START_REF] Picinbono | Coincidence analysis of point processes[END_REF], [START_REF] Picinbono | Characterization of nonclassical optical fields by photodetection statistics[END_REF], [START_REF] Saleh | Photoelectron Statistics[END_REF]. It is also sometimes called the intensity function (see p. 69 of [START_REF] Cox | The Statitstical Anlysis of Series of Events[END_REF]).

Note that the coincidence function is an even function of t and that it tends in general to λ 2 when t tends to infinity. This property comes from the fact that for large values of θ -θ the increments dN (θ)

and dN (θ ) become in general uncorrelated and their mean value is λdθ. In the case of Poisson processes these increments are independent whatever the value of θ-θ , which implies that the coincidence function is constant and equal to λ 2 .

DRAFT February 13, 2018

The coincidence function satisfies a relation which will play a fundamental role in what follows. Let

x [n]
i be the random variable (RV) equal to the distance between a point t i of the PP and the nth point posterior to t i , or x

[n] i = t i+n -t i . This positive RV is sometimes called the life time of order n. Its distribution function F n (t) is the probability P

[x [n] i
≤ t] and the PDF f n (t) of this life time is the derivative of F n (t) with respect to t. For simplicity we assume that x

[n] i is a continuous RV or that F n (t) has a derivative for any t. Since x

[n] i > 0, then f n (t) = 0 when t < 0. It can be shown (see p. 69 of [START_REF] Cox | The Statitstical Anlysis of Series of Events[END_REF]) that the coincidence function satisfies

c(t) = λ ∞ n=1 f n (|t|).
(

) 2
The principle of the proof is presented in the Appendix. This equation appears in a rather different forms in [START_REF] Kryukov | Calculation of the correlation function and the spectral power density of random sampling[END_REF], [START_REF] Beutler | The spectral analysis of impulse processes[END_REF].

This relation is the basis for the estimation of the coincidence function. For this purpose we use a TAC system which transforms the sequence of random points t i of the PP into a sequence of values

x i = t i+1 -t i . This positive random signal x i describes completely the PP and any PP generates such a signal x i .

Suppose now that we have M observations of x i , 1 ≤ i ≤ M . The problem is to deduce from these observations an estimation of the coincidence function c(t). The measurement of λ appearing in ( 2) is obvious because the density is the inverse of the mean value of the distances x i . For the estimation of c(t) the first step is to replace the series (2) by a finite sum of S terms. As a consequence we arrive at the estimation of the truncated coincidence function c S (t) = λ S n=1 f n (|t|). The choice of the appropriate value of S will be discussed later.

For a given value of S such that M/S is an integer N , and for each k satisfying 1 ≤ k ≤ S, we deduce from the observation x i a set of S signals s

[k] i defined by s [k] i = x (k-1)N +i , 1 ≤ i ≤ N, (3)
and zero otherwise. From these signals we construct a set of S vectors X [k] with N components defined by

X [k] i = s [k] i + s [k] i+1 + ... + s [k] i+k-1 , 1 ≤ i ≤ N. (4)
The histogram of these components will clearly yield an estimation of the PDF f k (t). As a consequence the histogram of the components of the vector X = X [1] + X [2] + ... + X [S] yields an estimation of the truncated bunching function b S (t) = (1/λ)c S (t). More precisely, consider an interval [t, t + ∆T [ and let n k (t, ∆T ) be the number of samples X

[k]

i recorded in this interval. If ∆T is sufficiently small the mean value nk (t, ∆T ) of n k (t, ∆T ) is approximately N f k (t)∆T . Let n(t, ∆T ) the number of samples of the February 13, 2018 DRAFT components of the vector X recorded in the same interval. From the definition of X and c S (t) = λb S (t)

the mean value n(t, ∆T ) of n(t, ∆T ) is n(∆T ) = b S (t)∆T N = b S (t) M ∆T S . (5)
The histogram of X yields for each experiment a value of the recorded number n(t, ∆T ). If N is sufficiently large this number yields an estimation of n(t, ∆T ), and then of b S (t) or c S (t). This procedure necessarily introduces a statistical error analyzed below. Finally if we are interested only in the shape of the coincidence function, the histogram can be represented in arbitrary units and the factor M ∆T /S can be omitted.

Let us now discuss the influence of the various parameters appearing in this measurement. The two most important are the number N of samples recorded and the number S of terms used for the construction of the vector X or for the approximation of the series (2) by a sum of S terms. On a pure mathematical point of view, the approximation of a function like c(t) defined as a series of other functions requires the discussion of the uniform convergence of the series. There are PDFs for which this uniform convergence is guaranteed, but as the method is used in the case where the PDFs f n (t) of ( 2) are unknown it is impossible to introduce in advance more details on this point.

The number N determines the statistical precision of the histogram and the number S determines the precision of an approximation of a series by a finite sum. These numbers are not completely independent.

Indeed the total number of samples analyzed by the histogram that yields the estimation of c S (t) is

M = S × N .
It is this number which determines the duration of the measurement or the complexity of the computer analysis. It is often necessary to introduce an upper bound of this number and in this case we must choose N and S for a given value of M . This introduces a compromise analogous to the one existing between bias and variance in many statistical measurements, as for example spectral analysis.

This will be discussed more precisely in the next section.

There are two last parameters appearing in the construction of the histogram yielding the approximation of the coincidence function. The first is the range ∆R of values of t of the function c(t). This range must be limited for great values of t. Indeed it is clear that c S (t) tends to 0 when t → ∞ because it is a finite sum of PDFs while, as seen above, c(t) → λ 2 . The second is the number of bins of the histograms of values belonging to ∆R. Here also there is a well known compromise. Too many bins yields a better precision of the approximation of c S (t) but also decreases the number of samples belonging to each bin and therefore increases their variance. This compromise between bias and variance in estimation of the PDF of signals has been analyzed in numerous statistical signal processing papers.

DRAFT February 13, 2018

Finally it must be noticed that there is a systematic error in the k -1 last components of the vector

X [k] . In particular X [k] N = s [k] N =
x kN which is not a sum of k terms. This kind of error also appears in all the measurements of correlation functions of random signals. There are ways to correct this error, but this does not play any role in our methods because in our experiments N >> S.

III. PERFORMANCES OF THE METHOD

There are only a few PPs for which the coincidence function can be explicitly calculated. They will now be used to look at the agreement between theoretical results and those of computer experiments which yields an indication on the quality of the method.

A. Poisson Processes

The PDF of the life time of order n of a Poisson process is

f n (t) = λ (λt) n-1 (n -1)! exp(-λt). (6)
Inserting this expression in (2) yields c(t) = λ 2 , which illustrates the fact that a Poisson process has no memory. The error when approximating (2) by a finite sum calculated for λ = 1 is then S (t) = 1 -exp(-t) S n=1 t n /n!. By a simple calculation we see that we can approximate c(t), which is here equal to 1, by using 10 terms with an error smaller that 10 -2 for t < 4. With 15 terms this range becomes t < 7. Similarly by using 20 terms we can approximate c(t) with an error smaller than 10 -3 for t < 9. This value of S is the greatest used in the experiments described below.

Let us now present experimental results with Poisson processes. A stationary Poisson process is a renewal PP with an exponential PDF. Then the distances x i between successive points are independent and identically distributed (IID) random variables and their common PDF is given by [START_REF] Beutler | The spectral analysis of impulse processes[END_REF] with n = 1. It is easy to generate in a computer experiment such a sequence, which allows us to discuss the performance of the method in terms of the various parameters appearing in its construction. The experiments are realized with a Poisson process of density λ = 1. This implies that b S (t) = c S (t).

In Fig. 1

B. Erlang Processes

Erlang processes belong to a class of PPs in which the PDF of the life time of order 1 is

f 1 (t) = µ(µt) exp(-µt). (7)
The mean value of the life time is 2/µ and the density λ is then µ/2. According to [START_REF] Beutler | The spectral analysis of impulse processes[END_REF] this means that f 1 (t) is the PDF of the life time of order 2 of a Poisson process of density µ.

Since the PDF [START_REF] Cox | Point Processes[END_REF] does not define a specific PP there are various different Erlang processes. Their difference comes from the fact that, even if the PDFs of the RVs x i are the same, the other statistical properties are different. There are of course a great variety of Erlang processes, depending on the correlation between the x i s. We shall only consider the two simplest cases.

The first case is the renewal Erlang process. It is the renewal PP defined by the PDF [START_REF] Cox | Point Processes[END_REF]. It can be generated by the following procedure. Let u i and v i be two independent signals and suppose that they are sequences of IID RVs with the same exponential distribution. It is clear that

x i = u i + v i is a sequence
of IID RVs with the common PDF [START_REF] Cox | Point Processes[END_REF].

In order to calculate the coincidence or the bunching function b(t) = c(t)/λ we use the following procedure. The Laplace transform (LT) of ( 7) is

F (s) = [µ/(s + µ)] 2 .
The assumption that the PP is a renewal process yields that the LT of

f n (t) is F n (t) = [F (t)] n .
Then the LT of the bunching function is

B(s) = F (s)/[1 + F (s)] or B(s) = µ 2 (s + µ) 2 -µ 2 = µ 2 1 s - 1 s + 2µ . ( 8
)
As a result the bunching function is

b(t) = µ 2 [1 -exp(-2µt)] . (9)
The second Erlang process is defined as follows. Consider a Poisson process and let x i be the life time of order 1, or the distance between successive points. Let

y i = x i + x i+1 . It is clear that y i is the life DRAFT February 13, 2018
time of order 2 of the Poisson process and its PDF is given by [START_REF] Cox | Point Processes[END_REF]. The RVs y i however are no longer independent because, for example, x i appears in y i and y i-1 . Consider now the PP defined by the fact that y i is the distance between successive points. It is then an Erlang process but not a renewal process.

Changing the origin of time we can write y [1] = x 1 + x 2 and

y [n] = x 1 + 2(x 2 + x 3 + ... + x n ) + x n+1 , n > 1, (10)
where the x i s are IID RVS with the exponential distribution µ exp(-µt). We deduce that the LT of the PDF of y [n] , the life time of order n, is

F n (s) = µ s + µ 2 µ 2s + µ n-1 , n ≥ 1. (11)
It is easy to calculate the LT of the bunching function defined by

B(s) = ∞ n=1 F n (s) and its inverse LT is b(t) = µ 2 1 -(1 - t µ ) exp(-µt) , (12)
which must be compared to (9) obtained previously.

The experimental results for these two Erlang processes appear in Fig. 2. These figures use arbitrary units because we are only interested in the shape of the coincidence function in comparison with its theoretical value. Fig. 2.1 corresponds to the Erlang renewal process. The continuous curves are those deduced from [START_REF] Saleh | Photoelectron Statistics[END_REF] in Fig. 2.1 and (12) in Fig. 2.2 and the experimental points are obtained from an experiment with 10 6 samples and with µ = 1 and S = 10. We see that there is excellent agreement between theory and experiment. Even if the the difference between the two Erlang processes is rather small, it appears clearly.

All the results of this section show the good performance of the method of measurement of the coincidence function introduced in this paper.

C. Poisson Process with Input Dead Time

Dead time effects appear in almost all the PPs practically used. Dead time occur when two points of a PP are so close together that they cannot both be registered. This means that some points are erased.

If the value D of the dead time is very small compared to the mean distance between points, its effect can be neglected. On the other hand when the density µ of the process increases there is always a value of the density such that dead time must be taken into consideration.

The input dead time is characterized by the fact that each point t i of the PP generates an interval [t i , t i + D[ such that all the points of the process arriving in this interval are erased. It is clear that when the density µ increases the number of points erased increases and for very large density almost all the February 13, 2018 DRAFT points of the initial PP are erased. This corresponds to the classical congestion phenomenon. Note that the input dead time is sometimes called type II counter (see p. 101 of [START_REF] Cox | Point Processes[END_REF]).

The dead time effect transforms a given PP into another one. As in our approach a PP is defined by the life times, or the distance between successive points, the first problem is to calculate the life time after dead time in terms of the life time of the initial PP. The explicit expression of the result is almost impossible to obtain in closed form. A recursive algorithm, however, has been introduced and analyzed in [START_REF] Picinbono | Dead time effects in point processes[END_REF]. This algorithm is used for the estimation of the coincidence function when the initial PP is a

Poisson process. In this case the theoretical calculation of the coincidence function with input dead time is especially simple. For this it suffices to start from [START_REF] Bartlett | The spectral analysis of point processes[END_REF]. If θ -θ < D the expectation is zero because one of the two points is necessarily erased. Otherwise c(t) is constant and its value is λ 2 . Note that λ is the density of the PP with dead time. If the density of the initial Poisson process is µ a simple calculation yields that λ = µ exp(-µD). This expression explains the congestion phenomenon indicated above because λ tends to 0 when D tends to infinity. In the opposite case λ = µ if D = 0.

Experimental results are presented in Fig. 3. The density of the Poisson process is µ = 1 and two values of the dead time are used : 0.5 and 1. The number S of terms in the sum defining c(t) is 10, which ensures good precision of the results for t < 3. The fact that c(t) = 0 for t < D is well verified.

On the other hand the fact that c(t) is constant for t > D is also verified although statistical fluctuations of the estimation remain. These fluctuations are more important for D = 1 than for D = 0.5 because more points are erased in the second case than in the first. Finally, as in these experiments M = 3.6 10 6 , ∆T = 10 -2 , and S = 10, the application of ( 5) with b S (t) = µ exp(-µD) yields that the mean number of points n(∆T ) in each bin is for t > D equal to 2.183 or 1.324 for the cases 1 or 2 respectively. This is very well verified in the two histograms of Fig. 3.

IV. ANALYSIS OF SOME PARTICULAR POINT PROCESSES

In this section we shall analyze the second-order properties of some point processes interesting in various applications but for which theoretical analysis is difficult, which justifies the experimental approach.

A. Poisson Process with Output Dead Time

This dead time, also called type I counter (see p. 101 of [START_REF] Cox | Point Processes[END_REF]), is characterized by the fact that each point θ i of the Poisson process which is not erased generates an interval [θ i , θ i + D[ such that all the points t i of the PP in this interval are erased. This means that only the points not erased contribute to the DRAFT February 13, 2018 dead time. In this case when µ → ∞ the PP after dead time becomes a periodic process, or a sequence of points with equal distances D.

The calculation of the life time after output dead time in terms of the life time of the initial PP is very complicated. As previously, however, a recursive algorithm for this purpose has been introduced in [START_REF] Picinbono | Dead time effects in point processes[END_REF].

We shall use this algorithm for the estimation of the coincidence function after output dead time effect when the input is a Poisson process.

Let us first present how this coincidence function can be calculated in this case. It is easy to see that the PP after dead time is a renewal process entirely defined by the PDF of its life time. This PDF is

f (t) = u(t -D)µ exp[-µ(t -D)], (13)
where u(.) is the unit step function, D the value of the dead time and µ the density of the Poisson process. It results from this expression that the mean value of the life time after dead time is

D + 1/µ.
This quantity is the inverse of the density which then is λ = µ/[1 + µD]. We observe that when µ → ∞, λ → D, which was indicated above. Furthermore when D → ∞, λ → 0, which appears with any dead time effect.

As the PP analyzed is a renewal process, its PDF f n (t) is deduced from (13) by n convolutions and the result is

f n (t) = u(t -nD)µ [µ(t -nD)] n-1 (n -1)! exp[-µ(t -nD)]. (14)
In order to calculate the coincidence function we can use [START_REF] Bartlett | An Introduction to Stochastic Processes[END_REF]. There is no explicit analytic expression for c(t), but its numerical calculation is possible. Indeed because of the term u(t -nD) in (14) the series is always limited to a finite sum. As we must have t -nD > 0, the number of terms in this sum is the greatest integer smaller than t/D.

Numerical results of calculations are presented in Fig. 4 where some examples of the bunching function The experiments with output dead time are much more difficult. Indeed the algorithm yielding the life time after dead time presented in [START_REF] Picinbono | Dead time effects in point processes[END_REF] introduces much more complexity and calculation time. This limits February 13, 2018 DRAFT the possible values of M , and thus the statistical precision of the procedure. This effect is increased by the fact that, as in any experiment with dead time, many points are erased and the number of those playing a role in the estimation tends to 0 when D increases. These functions suggest the following comments. We note first that if b → 0, the PP tends toward a periodic PP, which means a periodic distribution of points with a period equal to 1. This is easily verified experimentally but has no practical interest. The value b = 0.1 is still rather small and a memory of the periodic structure remains. The different PDFs f n (t) appearing in the figure do not overlap, at least for t < 8 used in Fig. 6.1 . We observe clearly the rectangular and triangular functions f 1 (t) and f 2 (t). For large values of t not represented in the figure the Gaussian shapes overlap in such a way that their sum becomes constant and equal to the asymptotic value 1 of b(t). This, however, requires very great values of S and M which cannot be obtained in our experiments.

Experimental results appear in

C. Point Process with Uniform Distribution and Correlated Life Times

In this section we assume that the PDF of the life time is still rectangular, as in the previous one, but that the PP is no longer a renewal PP. This means that the distances between successive points are identically distributed (rectangular distribution) but no longer independent.

An algorithm for generating this kind of life time has been presented in [START_REF] Picinbono | Dead time effects in point processes[END_REF], and its principle is outlined in the Appendix. The normalized correlation function of the first order life time is exponential, or p |n| . If p = 0, we find again the case of the previous section.

Experimental results are presented in Fig. 7. The parameters of this figure are the same as in Fig. 6.2 and our objective is only to evaluate the effect of the correlation characterized by the parameter p. In all these figures f 1 (t) is the same, and this appears for values of t smaller than 1.333. We have then a good example of point processes with the same distribution of life time and entirely different second order properties characterized by the coincidence functions represented in Fig. 7.

V. CONCLUSION

The second order properties of a point process are entirely contained in its coincidence function. In this paper we have shown that this function can be measured by an appropriate processing of the distance between successive points obtained for example from a time-to-amplitude converter. The theoretical foundation of this processing is an expression that yields a relation between the coincidence function and the probability density function of the life time of the PP. For its practical realization we replace a series by a sum of a finite number of terms constructed from appropriate histograms of samples of the signal.

The method has been tested in the case of PPs where the coincidence function has an explicit theoretical expression as in the case of Poisson processes with or without dead time, and also of some Erlang PPs.

February 13, 2018 DRAFT

The experimental results are in excellent agreement with the theoretical calculations. Finally, when this method is used for the measurement of coincidence functions that cannot be obtained theoretically, the experimental results exhibit various interesting features.

This yields [START_REF] Bartlett | An Introduction to Stochastic Processes[END_REF].

The statistical signal used in Section IV C is constructed as follows. Let u k and v k be two independent IID random variables. Suppose furthermore that u k takes only the value 0 or 1 and let p be the probability that u k = 1. Let f (.) be the PDF of the v k . Consider now the signal x k = u k x k-1 + (1 -u k )v k . It is clear that if the PDF of x 1 is f (.), all the other x k have the same PDF. It can be shown that this is asymptotically verified whatever the PDF of x 1 . The x k s have then the same PDF but are correlated.

Their correlation function γ k satisfies the recursion γ k = pγ k-1 . This yields γ k = σ 2 p k , where σ 2 is the variance of the RVs v k . For p = 0, x k = v k , while for p = 1 all the RVs x i s are equal.

  we present experimental results obtained when processing M = 2.10 6 samples of the life times x i . The parameters of these histograms are S = 5, S = 10, S = 15, S = 20. The widths of the bins of the histograms are 10 -2 . This means that there are 2.10 3 values recorded in the interval [0, 20] of the figure. Three comments must be made.First we verify experimentally the fact that the histograms are constant in an interval increasing with S. In the last histogram obtained for S = 20 the methods yields excellent results for t < 12, which means a duration twelve times greater than the mean distance between points of the Poisson process. February 13, 2018 DRAFT Secondly, as the range of possible values increases with S, the statistical precision decreases because the number of samples recorded in each bin decreases. When comparing the first and the last histograms we observe an increase of the statistical fluctuations. The only way to decrease these fluctuations is to increase the total number M of samples analyzed. Thirdly the number of points recorded in each bin corresponds to the values deduced from (5). Applying this expression for t = 0 and noting that in our experiment c S (0) = c(0) = 1 because λ = 1, we obtain n(∆T ) = M ∆T /S. In this experiment M = 2.10 6 , ∆T = 10 -2 . The values of n(∆T ) corresponding to the four values of S used are then 4.10 4 , 2.10 3 , 1.33.10 3 , and 10 3 , which cleanly appears on the figure.

b

  (t) = (1/λ)c(t) are presented. The value of µ is 1 and the values of D are 1/4, 2/3, 1, 3/2, in such a way that the corresponding values of the density λ are 0.8, 0.6, 0.5, 0.4. We verify that b(t) → λ when t → ∞. In reality the values of t are limited to the interval 0 < t < 8 and for this interval the sum (2) is truncated to 12 terms. We see on this figure that when D = 1/4 the effect of the dead time on the coincidence function disappears when t > 2/3. For D = 1 this inequality becomes t > 4. Furthermore all these curves show a discontinuity of the derivative for the value t = 2D which is the time at which the second term f 2 (t) appears in the sum (2).

Fig. 5 . 5 and D = 1 . 5 . 3 .

 55153 Estimations of the coincidence function are presented for D = 0.The results are in a rather good agreement with those of Fig.4, but some statistical fluctuations remain. Their suppression would require more memory and computer time.B. Renewal Process with Uniform DistributionConsider a renewal PP in which the PDF f (t) of the distance between successive points is a rectangular function equal to 1/2b when t is in the interval [m -b, m + b] and zero otherwise. This kind of PDF appears in the description of jitter phenomena in which the distance between successive points is not strictly constant, which appears if b → 0. The mean value of the distance is clearly m and as a result the density λ of the corresponding PP is 1/m.The PDFs f n (t) can be calculated by successive convolutions. The calculation is possible but rapidly becomes very complicated and no explicit expression of the convolution is available. Two points however can be noted. First the PDF f 2 (t) of the second order life time is a triangular function equal to(1/4b 2 )[-| t-2m | +2b] for | t-2m |< 2,and zero otherwise. In particular f 2 (2m) = 1/2b. Secondly when n >> 1, f n (t) tends to have a Gaussian shape of the form N (nm, nσ 2 ), where σ 2 is the variance associated with the PDF f (t) and equal to b 2 /Experimental results are presented in Fig.6. The bunching function of the PP is estimated for the values m = 1. The values chosen for b are 1/10, 1/3, 1/2, and 1. The number of samples used is M = 3.10 6 , the number of terms in (2) is S = 9, and the width of the bins is ∆T = 10 -2 . As a consequence the coefficient M ∆T /S in (5) is (1/3)10 4 .

Fig. 6 .

 6 Fig. 6.2 corresponds to b = 1/3. This value has been chosen because it is the greatest value of b for

Fig. 7 . 1 corresponds

 71 to p = 0, the same situation as in Fig.6.2, and we find again the same form of coincidence function. When p increases this coincidence function is strongly modified and theoretical calculations are impossible because of the complexity of the model. But the effect of a correlation is perfectly clear in Figs.7.2, 7.3, and 7.4.

  VI. APPENDIXThe expectation[START_REF] Bartlett | The spectral analysis of point processes[END_REF] defining the coincidence function is equal to the probability P {[dN (θ) = 1] ∩ [dN (θ ) = 1]}. This probability can be writtenP {[dN (θ) = 1] ∩ [dN (θ ) = 1]} = P [dN (θ) = 1].P {[dN (θ ) = 1]|[dN (θ) = 1]}. (15)We deduce from the definition of the densityλ that P [dN (θ) = 1] = λdθ. The PDF f n (θ -θ ) is defined by f n (θ -θ )dθ = P {[dN (θ ) = 1]|A n (θ -θ) ∩ [dN (θ) = 1]},(16)where A n (θ -θ) is the event that there are n -1 points of the PP in[θ + dθ, θ [. It is clear that P {[dN (θ ) = 1]|[dN (θ) = 1]} = ∞ n=1 P {[dN (θ ) = 1]|A n (θ -θ) ∩ [dN (θ) = 1]}

  For this value of b the rectangular function f 1 (t) is equal to 3/2. Applying (5) with the value calculated above yields (1/2)10 4 = 5000, which clearly appears in Fig.6.2. The same can be said for f 2 (2) which has the same value. Finally for t → ∞ the bunching function tends to 1, which yields a number of samples in each bin of (1/3)10 4 , which appears in Figs.7.2, 7.3, and 7.4. Theoretical calculation of the bunching function corresponding to the values used in Figs.6.3 and 6.4 is almost impossible. This justifies the interest of experimental results displayed on these figures.

This value has been chosen because it is the greatest value of b for which f 1 (t) and f 2 (t) do not overlap. The corresponding limit point is t = 1+1/3, which clearly appears DRAFT February 13, 2018 on the figure.