Maroua Mehri
email: maroua.mehri@univ-lr.fr

Petra Gomez-Krämer

Pierre Héroux
email: pierre.heroux@univ-rouen.fr

Rémy Mullot
email: remy.mullot@univ-lr.fr

Old document image segmentation using the autocorrelation function and multiresolution analysis

Keywords: Segmentation, autocorrelation, directional rose, multiresolution, consensus clustering

Recent progress in the digitization of heterogeneous collections of ancient documents has rekindled new challenges in information retrieval in digital libraries and document layout analysis. Therefore, in order to control the quality of historical document image digitization and to meet the need of a characterization of their content using intermediate level metadata (between image and document structure), we propose a fast automatic layout segmentation of old document images based on five descriptors. Those descriptors, based on the autocorrelation function, are obtained by multiresolution analysis and used afterwards in a specific clustering method. The method proposed in this article has the advantage that it is performed without any hypothesis on the document structure, either about the document model (physical structure), or the typographical parameters (logical structure). It is also parameter-free since it automatically adapts to the image content. In this paper, firstly, we detail our proposal to characterize the content of old documents by extracting the autocorrelation features in the different areas of a page and at several resolutions. Then, we show that is possible to automatically find the homogeneous regions defined by similar indices of autocorrelation without knowledge about the number of clusters using adapted hierarchical ascendant classification and consensus clustering approaches. To assess our method, we apply our algorithm on 316 old document images, which encompass six centuries (1200-1900) of French history, in order to demonstrate the performance of our proposal in terms of segmentation and characterization of heterogeneous corpus content. Moreover, we define a new evaluation metric, the homogeneity measure, which aims at evaluating the segmentation and characterization accuracy of our methodology. We find a 85% of mean homogeneity accuracy. Those results help to represent a document by a hierarchy of layout structure and content, and to define one or more signatures for each page, on the basis of a hierarchical representation of homogeneous blocks and their topology.

I. INTRODUCTION

A document has a structure which provides additional information. Without document structure, it would be difficult to index and retrieve correctly the information contained in the document. In this context, document structure analysis remains a fundamental step and crucial stage in any indexing and retrieval system such as optical character recognition (OCR) and graphic recognition modules. Several researches have been carried out to characterize the document layout with the result of structuring it into three different levels: the physical layout, the intermediate functional structure, and the logical structure analysis [START_REF] Mao | Document structure analysis algorithms: a literature survey[END_REF]. Firstly, the physical level defines both the typography and document organization. The typography sets the style of different kinds of informative regions (fonts, colors, lines, frames, etc.) and the form layout (line spacing, alignment, etc.). The document organization describes the layout of all the visual elements (characters, words, lines, blocks, columns, and non-text regions) which compose it and the topographical relationships between those elements (hierarchy, inclusion, neighborhood position). Secondly, the intermediate functional structure is a physical interpreted one which is adapted to the recognition of the logical structure. Finally, the logical level concerns the interpretation of different parts that compose the document and specifies the logical relationship between them [START_REF] Mullot | Les documents écrits: De la numérisation à l'indexation par le contenu[END_REF].

In this work, we focus only on the first segmentation level i.e. the physical level. Several works have been presented on document image segmentation and characterization tools. There are two main kinds of page segmentation and characterization methods: the first one includes data-driven, model-driven, and hybrid methods and the second one is based on texture analysis.

A. Data-driven, model-driven, and hybrid methods

Data-driven, model-driven, and hybrid methods are devoted to contemporary document recognition and are significantly widespread in the literature because those methods are based on a strong a priori knowledge such as the repetitiveness of document structure in a corpus. This family of document structure analysis and page segmentation methods, mixing a priori knowledge and image analysis, can be classified in three main categories: data-driven, model-driven, and hybrid methods [START_REF] Mao | Document structure analysis algorithms: a literature survey[END_REF], [START_REF] Mullot | Les documents écrits: De la numérisation à l'indexation par le contenu[END_REF].

The first class of methods, known as data-driven segmentation, does not include (or little) knowledge of the model. Those approaches are based on low-level data mining of pixels (color, position, etc.). They rely on the study of the space between pixels and the grouping of pixels in order to segment the different elements of the page content into blocks. For example, the Run Length Smearing Algorithm (RLSA) [START_REF] Wong | Document Analysis System[END_REF] studies the spaces between black pixels in order to link together neighboring black areas. There are certain limitations of data-driven segmentation methods: Firstly, they are based on the definition of complex criteria and rules. Secondly, those methods are sensitive to noise and not robust to slanted texts. Thus, the data-driven approaches are suitable for documents whose areas are clearly demarcated and rectangular. Furthermore, the pertinence of this segmentation category depends on the particular layout and idiosyncrasies of the documents.

The second class of methods (model-driven methods) is guided by a model of the document. Often used for well-defined and invariant structured documents, those approaches are based on strong a priori knowledge to guide the segmentation and recognition. For example, the XY-CUT algorithm [START_REF] Khedekar | Text -Image Separation in Devanagari Documents[END_REF] consists in computing the horizontal and vertical projection profiles (corresponding to the sum of the pixels along the horizontal axis and the vertical axis) of the whole document image and in iteratively splitting them into smaller ranges until a condition about hollow projections (interline spaces) has been satisfied. This requires the definition of criteria for cutting (and possibly fusion). Although the modeldriven approaches are generally faster, they are not welladapted to complex layout documents.

Finally, the hybrid methods are often non-sequential and combine data-driven and model-driven algorithms. For example, the hybrid method [START_REF] Pavlidis | Page Segmentation and Classification[END_REF] uses a split-andmerge strategy. However, by combining tools from datadriven and model-driven approaches, hybrid methods can deal with a wide variety of documents and cope with complicated page segmentation problems, but many parameters and thresholds must be adjusted.

B. Segmentation methods based on texture analysis

To handle the drawbacks of the methods described above, new alternatives for document structure analysis based on texture have been developed, which ensure the segmentation and the characterization of information contained in a document. Those methods are pixelbased and do not require neither a document model nor a priori information relative to the semantic and physical characteristics of the document class. Thus, texture-based segmentation methods characterize generally complex documents aiming at segmenting the different elements of the page content into homogeneous blocks. A texture is defined as a spatially uniform distribution of local gray-value variations. In the literature [START_REF] Allier | Texture Feature Characterization for Logical Pre-labeling[END_REF], texture is defined as a suitable measure for the analysis of the block contents of the physical layout. Texture analysis methods have been used in image segmentation in order to extract textural characteristics. Texture-based segmentation methods can be classified into four categories: methods based on probabilistic models, and geometric, frequential, and statistical methods.

Markov random fields [START_REF] Nicolas | Handwritten document segmentation using hidden Markov random fields[END_REF] and fractals [START_REF] Ferrell | Application of fractal encoding techniques for image segmentation[END_REF] are both the most commonly used tools based on probabilistic models. This category of texture-based segmentation methods is complex to implement. There are many difficulties in the learning phase and a long computation time is required. Markov random fields are perfectly adapted to documents with high variability in terms of the layout and the quality of the scanned document, which yields good performances in handwritten documents. However, Markov random fields are not robust since the learning phase is only valid for one type of document at a time. Fractal dimensions compute measures of texture roughness and repeatability of a pattern. Fractals are considered as a useful tool for image segmentation when the image characteristics tend to be predictable and repetitive and in which the objects to segment tend to be irregular or different from the background.

The second class of texture analysis based segmentation methods is known as geometric methods. Those methods are used to describe intricate patterns to retrieve and to characterize the notion of a texton. Among the classics of geometric methods, momentbased texture segmentation [START_REF] Tuceryan | Moment Based Texture Segmentation[END_REF] is one of the well-known methods. Anyway, moment-based texture segmentation is not sufficient to discriminate all types of texture and the algorithm needs a non-linear transformation of the images. Another geometric method [START_REF] Varshney | Block-segmentation and Classification of Grayscale Postal Images[END_REF] aims at extracting three classes: text, background and design from postal images based on six features derived from wavelet transforms. Even if the proposed algorithm has a good recognition rate, one of the features must be adjusted manually and the efficiency (computation time) of the algorithm is limited.

Frequential methods like Gabor filters, Fourier transforms, and wavelets are widely used in indexing and segmentation of natural images. For instance, the unsupervised texture segmentation algorithm [START_REF] Jain | Text Segmentation Using Gabor Filters for Automatic Document Processing[END_REF] is used to segment an input image into regions of homogeneous texture based on a bank of Gabor filters. Gabor filters have the advantage of reducing the computational complexity and are suitable for document texture analysis. One of the limitations of such an algorithm based on a fixed set of Gabor filters is that many parameters must be fixed [START_REF] Raju | Text Localization and Extraction from Complex Color Images[END_REF]. Another frequential approach [START_REF] Sabharwal | Indexing image databases using wavelet and discrete Fourier transform[END_REF] combines the wavelet and the Fourier transform to index image databases. Although their algorithm is faster and more robust than the separate use of the discrete Fourier and wavelet transforms, the computation time is directly dependent on the level of wavelet decomposition. The method described in Ref. [START_REF] Qiao | Document image segmentation using Gabor wavelet and kernel-based methods[END_REF], combines kernel-based methods and a Gabor wavelet to segment document images scanned from popular newspapers and journals. According to the authors, the analysis of multiresolution and multiorientation properties of an image is ensured, but the effectiveness and computational complexity of the algorithm is no longer preserved and a proper postprocessing is needed to improve the segmentation result.

The last class is statistical methods that have the advantage of being simple to implement and their effectiveness is proved. The GLCM (Grey Level Cooccurrence Matrix) [START_REF] Haralick | Textural features for image classification[END_REF] is one of the classics of texture analysis based statistical segmentation methods. By computing some indices on the GLCM [START_REF] Hall-Beyer | GLCM Texture: A Tutorial[END_REF], the texture regularity and repetitiveness are characterized. The authors of [START_REF] Caponetti | Document page segmentation using neuro-fuzzy approach[END_REF] propose a document page segmentation method using a neuro-fuzzy methodology. Another approach is the extraction algorithm of texture features [START_REF] Journet | Document image characterization using a multiresolution analysis of the texture: application to old documents[END_REF] that is devoted to the analysis of documents. The computed texture features are based on frequencies and the autocorrelation function. This method gives good information on the principal orientations and periodicities of the texture allowing to characterize the content of images without any assumption on the image structure or properties. Although their results are promising, their algorithm is computationally expensive because it is carried out for each pixel and the size of the analysis window is a critical parameter that is difficult to determine. The authors of [START_REF] Uttama | Segmentation and Retrieval of Ancient Graphic Documents[END_REF] introduce a novel lettrine segmentation method based a combination of different texture analysis approaches: GLCM [START_REF] Haralick | Textural features for image classification[END_REF], autocorrelation function [START_REF] Journet | Document image characterization using a multiresolution analysis of the texture: application to old documents[END_REF], etc.

The authors of [START_REF] Lemaitre | Multiresolution Cooperation Improves Document Structure Recognition[END_REF] claim that multiresolution in document image analysis and pyramid methods [START_REF] Contassot-Vivier | Multiresolution approach for image processing[END_REF] in image processing allow to perceive texture at different scales and provide rich information like the characteristics of gray level distribution. There are two ways to implement multiresolution algorithms: The first solution consists in fixing the window size and changing the size of the image. The second solution is to keep the size of the original image and simply vary the size of the analysis window.

C. Scope

This work is a part of the DIGIDOC project (Document Image diGitisation with Interactive DescriptiOn Capability) 1 . Generally, the DIGIDOC project aims at simplifying and improving the archiving, processing, comparison and indexing of digitized old document images. Specifically, the objective of the DIGIDOC project is to integrate a module in scanners that will provide in addition to the scanned image a set of descriptors computed on it. Those descriptors dedicated to the acquisition, storage, analysis, and indexing of the scanned documents, will adapt the quality of the scanned document with respect to its content and to the subsequent use of the document image. This would help ensuring better interaction with scanners and offers new tools for document analysis.

The work presented in this paper is a part of our goal to provide a similarity measure between pages by defining one or more signatures, e.g. a graph of homogeneous regions, for each digitized page. Hence, a set of metadata characterizing the physical structure of pages in terms of homogeneous areas and topological relationships has to be proposed. In this context, we present a method to distinguish similar regions of the analyzed document without expressing any hypothesis about its physical structure and its logical structure. According to the literature, texture analysis based segmentation methods are the most appropriate choice if no assumption on document structure (model) or the typographical parameters (font size) should be made. Therefore, we present a fast and automatic layout segmentation to determine homogeneous areas in old document images based on texture analysis. Our proposal consists firstly in characterizing the content of old documents by extracting the autocorrelation features in the different areas of a page and at several resolutions. Then, we show that is possible to find automatically the homogeneous regions defined by similar indices of autocorrelation without knowledge about the number of clusters using adapted hierarchical ascendant classification (HAC) and consensus clustering (CC) approaches. This paper is organized as follows. We detail in Section II our method for the segmentation of old documents and the characterization of their content by extracting five autocorrelation features in the different areas of the digitized page and at four resolutions. To assess our method, we present in Section III, the pertinence of our proposition with segmentation experiments which show homogeneous regions defined by similar indices of autocorrelation and we quantify the relevance of our layout analysis method with a new evaluation metric, the homogeneity accuracy. Our conclusion and future work are presented in Section IV.

II. SEGMENTATION USING THE AUTOCORRELATION FUNCTION AND MULTIRESOLUTION

We propose a fast automatic layout segmentation of old document images based on non-parametric tools: the autocorrelation function and multiresolution analysis. Our method focuses on the extraction of texture information without binarizing the image. We avoid a binarization of the image because it causes a loss of information specifically textural information. The pertinence of the segmentation experiments of old documents [START_REF] Journet | Document image characterization using a multiresolution analysis of the texture: application to old documents[END_REF] that are based on the autocorrelation function leads us to work with autocorrelation features. The authors use the autocorrelation results to construct a rose of directions [START_REF] Bres | Contributions à la quantification des critères de transparence et d'anisotropie par une approche globale: application au contrôle de qualité de matériaux composites[END_REF]. In order to obtain the orientation information, three features have been derived from the directonal rose: its main orientation, the intensity of the autocorrelation function for the main orientation, and the standard deviation of the intensities of the directional rose [START_REF] Journet | Document image characterization using a multiresolution analysis of the texture: application to old documents[END_REF]. In order to extract texture information, we chose to compute five descriptors based on the autocorrelation function carried out in the different areas of a gray-scale page and at different resolutions. The proposed approach is pixel-based and independent of the model or type of document image. It does not require any a priori knowledge on the content or font styles of the document. Thus, our method is adapted to all kinds of document. The proposed approach is illustrated in Figure 1. It is composed of two main stages which are detailed below: Firstly, the computation of autocorrelation features for each page which are then used in an unsupervised clustering approach (block 1 on Figure 1) in order to determine and characterize the homogeneous regions in the document (block 3 on Figure 1). Secondly, we estimate the true number of clusters of the homogeneous regions defined by similar autocorrelation indices by performing the consensus clustering method applied to a number of pixels chosen randomly from a few pages of a book (block 2 on Figure 1). The clustering approach is more significant when we compute it on an entire book instead of processing each page of a book individually as our goal is to compare and index the content of digitized book.

Autocorrelation features + Gray level

A. Autocorrelation feature computation

The first stage of our proposal (see Figure 1) is to compute the autocorrelation features. We propose a feature vector composed of five texture indices all extracted from the autocorrelation function and based on a nonparametric tool, the directional rose [START_REF] Bres | Contributions à la quantification des critères de transparence et d'anisotropie par une approche globale: application au contrôle de qualité de matériaux composites[END_REF]. The extraction of those five texture attributes helps us to characterize the content of the digitized document. The goal of our feature extraction process is to propose a set of metadata characterizing the physical structure of pages in terms of homogeneous areas and topological relationships. The texture features are computed at various sizes of analysis windows in order to adopt a multiresolution approach. We decide to keep the size of the original image and varying the size of the analysis window. The sliding window is shifted horizontally and vertically scanning the whole image. In order to avoid side effects, a quick and easy way to compute texture features on the whole image, is border replication. The extraction of autocorrelation descriptors per block is carried out at four different sizes of sliding windows: (16 × 16), (32 × 32), (64 × 64), and (128 × 128). Therefore, we obtain 20 numeric values (5 texture indices × 4 sliding window sizes) for each selected pixel from the digitized document image (see Figure 2). We will detail in the following the five texture indices.

The textual regions in the document are considered as textured areas while non-text contents in the document, such as blank spaces, graphics, and noise, are considered as regions with different textures. Therefore, we use the directional rose derived from the autocorrelation function to characterize homogeneous regions in a digitized document. The autocorrelation function

= F F T -1 ([F F T [I(x, y)] F F T * [I(x, y)]]) (1)
where I(x + α, y + β) is the translation of the analysis window of an image I(x, y) by α and β pixels along the horizontal and vertical axes respectively, defined on the plane Ω. F F T , (.) * , and (.) -1 denote respectively the Fast Fourier transform, the complex conjugate, and the inverse transform.

From the autocorrelation function, the directional rose can be deduced, that has been initially proposed by [START_REF] Bres | Contributions à la quantification des critères de transparence et d'anisotropie par une approche globale: application au contrôle de qualité de matériaux composites[END_REF]. The rose of directions is a polar diagram based on the analysis of the autocorrelation results. In order to identify the main orientation of the analyzed image, the rose of directions is computed for each orientation by summing up the different values of the autocorrelation function (see equation ( 1)):

R I (x,y) (Θ i ) = Di R I(α,β) (x,y) (2)
where Θ i ∈ [0, 180] is the selected orientation of the set of the possible orientations D i , which is represented by a straight line passing through (x, y) and the angle Θ i . In order to select only the relative variations of all contributions for each direction, the authors of [START_REF] Journet | Document image characterization using a multiresolution analysis of the texture: application to old documents[END_REF] present a normalization of the rose of directions. The definition of the relative sum R I (x,y) (Θ i ) is:

R I (x,y) (Θ i ) = R I (x,y) (Θ i ) -R I min R I max -R I min (3)
with R I max = R I min and where R I min and R I max represent respectively the minimum and maximum value of R I (x,y) (Θ i ), computed both on the analysis window of an image I(x, y).

In order to illustrate the relevance of discriminating textual regions from graphical ones in the analyzed document and to determine the main orientation of a texture, we present in Figure 3 the results of the directional rose for four different textures. It can be observed that the shape of the rose is different for each type of texture. For textual regions such as (c), the shape of the rose depends on the orientation of the text and the main information. The horizontal orientation (0 • and 180 • ) is clearly identifiable in (g). For the drawing (d), the directional rose (h) is deformed. A study of the variance of the shape of the rose applied to old grayscale documents has shown that the variety of textures does not define a homogeneous model of the rose of directions [START_REF] Journet | Document image characterization using a multiresolution analysis of the texture: application to old documents[END_REF]. The computation of the rose helps us to extract significant and relevant indices. The authors of [START_REF] Journet | Document image characterization using a multiresolution analysis of the texture: application to old documents[END_REF] define three texture features linked to the orientation information in order to analyze the digitized document and to describe their contents. The first texture feature F (1) (x,y) corresponds to the main angle of the rose of directions extracted from its maximal intensity. It is normalized by the deviation from the horizontal angle in order to avoid handling circular data. It is given by:

F (1) (x,y) = 180 -argmax Θi∈[0,180] (R I (x,y) (Θ i )) (4)
Likewise the second texture feature F

(x,y) corresponds to the intensity of the autocorrelation function for the main orientation (equation ( 4)), which is computed on the non-normalized value of the autocorrelation function (equation ( 2)). This feature evaluates the level of anisotropy of the analysis window of an image I(x, y) since the directional rose ensures the association of gray levels of pixels in a specific direction. It is computed as:

F (2) (x,y) = R I (x,y) ( argmax Θi∈[0,180] (R I (x,y) (Θ i ))) (5)
The third texture index F

(x,y) characterizes the overall shape of the rose. F 1, it means that the main orientation is significantly more prevalent than the other orientations. Otherwise, if the variance of the rose intensities is high, it signifies that the rose is deformed and a large number of orientations are present in different proportions (graphic blocks). Hence, the third texture attribute is:

F (3) (x,y) = σ 2 (R I (x,y) (Θ i )) (6)
where

Θ i ∈ [0, 180] \argmax(R I (x,y) (Θ i ))
and σ represents the standard deviation estimator. 2 where µ, σ and n are respectively the mean value, standard deviation and 179 orientation values.

σ 2 = 1 n-1 n i=1 (R I (x,y) (Θ i )) 2 -n n-1 (µ)
In addition to the three texture features that are related to the orientation information of the autocorrelation function, we introduce two other texture attributes also in relationship with the autocorrelation function: the mean stroke width and height [START_REF] Ouji | Chromatic / Achromatic Separation in Noisy Document Images[END_REF] of an image. In contrast to the initial work which computes the features in horizontal and vertical direction [START_REF] Ouji | Chromatic / Achromatic Separation in Noisy Document Images[END_REF], we propose estimating the mean stroke width and height accurately along the axis of the main angle of the directional rose in order to indicate the order of magnitude of the main strokes thickness.

The next texture index corresponds to the estimation of the mean stroke width F (4) (x,y) . It is computed from a derivative of the autocorrelation function along the axis of the directional rose's main angle Θ (see equation ( 4)) if Θ ∈ [START_REF] Varshney | Block-segmentation and Classification of Grayscale Postal Images[END_REF]80], otherwise the mean stroke width is estimated along the horizontal axis. If the growth rate of the sequence defined in the equation [START_REF] Nicolas | Handwritten document segmentation using hidden Markov random fields[END_REF] becomes lower than 10%, we estimate the mean stroke width, otherwise we continue to compute the sequence (equation ( 7)) until we reach the horizontal borders of the analyzed sliding window. Thus, we define F (4) (x,y) as:

F (4) (x,y) = Θ∈[10,80] I(x, y) -T Θ (α,0) (I( y |tan(Θ)| , y)) (7)
where T Θ (α,0) (I(., .)) is the translation of the analysis window of an image I by α pixels along the axis of the main angle Θ of the directional rose and Θ = F (1) (x,y) . The computation of the last texture attribute is similar to that of the fourth texture index

F (4) (x,y) . F (5) (x,y)
is an estimation of the mean stroke height computed along the axis of the directional rose's main angle Θ if Θ ∈ [START_REF] Varshney | Block-segmentation and Classification of Grayscale Postal Images[END_REF]80], otherwise the mean stroke height is estimated along the vertical axis. If the growth rate of the sequence defined in the equation ( 8) becomes lower than 10%, we estimate the mean stroke height, otherwise we continue to compute the sequence (equation( 8)) until we reach the vertical borders of the analyzed sliding window. F

B. Unsupervised clustering approach

After the computation of the autocorrelation features of all pages, we propose a non-parametric unsupervised clustering method (block 1 on Figure 1) which aims at determining the homogeneous regions. The homogenous regions are assumed to have similar autocorrelation indices. They are obtained without any knowledge on the number of clusters using hierarchical ascendant classification (HAC) and consensus clustering (CC). Firstly, we perform HAC in order to discriminate the foreground cluster which is considered as the most representative and significant pixels. But, the great variability of the old document content, like the various fonts and different types of drawing (ornaments, illuminations, drop caps, initial letters, frames, stamps, etc.), remains a difficulty in correctly predicting the optimal number of clusters. Thus, we propose to perform CC applied to a number of pixels chosen randomly from a few pages of a book in order to estimate the true number of clusters. Finally, we apply a third clustering step using HAC in order to classify all foreground pixels of the document in homogenous regions.

First of all, we perform a segmentation step with the help of non-supervised techniques which aims at extracting two clusters. One represents the information of the foreground (noise, text fields, drawings, etc.) and the other represents the background. The authors of [START_REF] Nguyen | Stroke feature extraction for lettrine indexing[END_REF] apply HAC on stroke features in order to classify the strokes of initial letters with interesting classification results. We chose to perform an adapted HAC using the Ward criterion [START_REF] Ward | Hierarchical Grouping to Optimize an Objective Function[END_REF], applied only on the average gray level of the analyzed sliding window and setting the maximum number of clusters to two. In the following, we consider only the foreground cluster in order to deal with representative and significant pixels and to limit the amount of data.

After the extraction of the foreground pixels, a second clustering step is performed using a non-supervised technique. The objective is to determine the exact number of clusters and to find the homogeneous regions defined by similar autocorrelation indices. Typically, conventional unsupervised clustering techniques like partitioning and hierarchical methods [START_REF] Duda | Pattern Classification. 2nd Edition[END_REF] can not determine the ideal number of clusters. However, using the CC method [START_REF] Simpson | Merged consensus clustering to assess and improve class discovery with microarray data[END_REF], the ideal number of clusters can be obtained. The CC consists in calculating a consensus matrix obtained by iterating multiple runs of clustering algorithms with random and re-sampled clustering options. The consensus matrix analyzes the consistency of clustering results from five different clustering algorithms: agglomerative hierarchical clustering (agnes) [START_REF] Kaufman | Finding Groups in Data: An Introduction to Cluster Analysis[END_REF], divisive analysis clustering (diana) [START_REF] Kaufman | Finding Groups in Data: An Introduction to Cluster Analysis[END_REF], partitioning around representative objects (pam) [START_REF] Kaufman | Finding Groups in Data: An Introduction to Cluster Analysis[END_REF], k-means clustering (k-means) [START_REF] Macqueen | Some Methods for Classification and Analysis of Multivariate Observations[END_REF], and hierarchical cluster analysis (hac) [START_REF] Lance | A general theory of classificatory sorting strategies 1. Hierarchical systems[END_REF]. The authors of [START_REF] Simpson | Merged consensus clustering to assess and improve class discovery with microarray data[END_REF] show that hierarchical clustering methods are highly sensitive to outliers while partitioning ones are relatively insensitive. Therefore, they propose to use a new merged consensus clustering approach by applying a weighted averaging of the clustering results to estimate the true number of clusters. The clustering approach is more significant when we compute it on an entire book instead of processing each page of a book individually as our goal is to compare and index the content of a digitized book. Thus, in order to estimate the true number of clusters k, we use the consensus merge method [START_REF] Simpson | Merged consensus clustering to assess and improve class discovery with microarray data[END_REF] by determining k for five different clustering methods (agnes, diana, pam, kmeans, hac) only for a set of randomly selected pixels of few pages of a book. By weighting the different clustering methods, we mitigate extremes in consensus values that can be created by the sensitivity of some algorithms. Finally, the optimal number of clusters corresponds to the largest change in area under the cumulative density curve ∆k for the merge consensus matrix. An example of ∆k is shown in Figure 4. The optimal number of clusters k opt = 2 is estimated by finding the peak in the ∆k values of the merge curve.

Once the number of clusters is known, we perform a third clustering step in order to classify all foreground pixels of the document in the homogenous regions. HAC is applied using the optimal number of clusters k opt and the Ward criterion. Figure 5 illustrates the final result of the three clustering steps where two clusters were obtained. From (b) and (c), we can see that the page has been segmented into graphic regions (blue), which correspond to an ornament and a drop cap, and textual regions (red). (d) represents the proposed ground-truth (graphic regions (blue), textual regions (font1-red, font2-green)).

III. EXPERIMENTAL RESULTS

To assess our method, we present in this section the results of our algorithm on heterogeneous content. Moreover, we define a new evaluation metric, the homogeneity measure, which aims at evaluating the segmentation and characterization accuracy of our method.

A. Evaluation Corpus

In the context of a collaboration with the BnF within the DIGIDOC project, we had access to Gallica digital library2 for old documents. The heterogeneity of pages, the database size and the degradation of certain documents of this corpus are examples representing specific issues and scientific challenges. This corpus includes many particularities of old documents, for instance a great variability of the page layout: complicated layout, random alignment, specific fonts, presence of embellishments (ornaments, illuminations, drop caps, initial letters, frames, etc.), variations in spacing between the characters, words, lines paragraphs and margins, and the superimposition of informations layers (stamps, handwritten notes, noise, back-to-front interference). For a first evaluation of our approach, we have selected 316 pages from 13 books of two categories: 7 printed monographs and 6 manuscripts that encompass six centuries (1200-1900) of French history. For each category, we have decided to select three types of page content: 110 pages containing only two fonts, 100 pages containing graphics and single font texts, and 106 pages containing graphics and text with two different fonts.

B. Method evaluation

Our method reaches very satisfying results when comparing visually the segmentation results (see Figure 5). Indeed, this method of assessing the effectiveness of a segmentation method is inherently a subjective evaluation and we need to assess and evaluate the effectiveness using an appropriate quantitative metric. The lack of appropriate quantitative measures for segmentation quality and the difficulty of defining criteria for specific application-dependent segmentation are the shortcomings that limit researchers in an objective unsupervised evaluation of their results.

Evaluation of segmentation and region classification requires a ground-truth which is performed using the Ground-truthing Editor (GEDI) 3 , a public domain document image annotation tool, providing us spatial boundaries of regions with labels. GEDI is considered as an annotation toolkit, which allows us to define manually bounding boxes drawn around each selected zone in order to define in page regions with different labels.

Several methods have been presented to measure the performance of segmentation methods [START_REF] Wontaek | Performance Evaluation Tools for Zone Segmentation and Classification (PETS)[END_REF]. For instance, the authors of [START_REF] Ge | New benchmark for image segmentation evaluation[END_REF] evaluate their segmentation 3 http://gedigroundtruth.sourceforge.net/ using Jaccard coefficient. However, this coefficient is not suitable to assess the accuracy of our method because our goal is not an accurate pixel-based segmentation, but we are interested in finding the homogeneous regions defined by similar indices of autocorrelation. Thus, by defining our ground-truth, we define manually rectangular regions drawn around each selected zone and precise different labels when regions with different fonts (see Figure 5-(d)). Furthermore, each rectangular region is characterized by its location in the page, its height, its width, and a label. Then, we define an evaluation metric based on an homogeneity measure, which aims at evaluating the accuracy of our methodology in terms of matching regions. This metric is based on spatial overlaps of the ground-truth rectangle and the segmentation result. It takes into account the label correspondence between the ground-truth and result regions, based on the set of pixels they contain. Our evaluation accuracy is defined by the following equation: .., g j , ..., g m } are respectively the sets of result blocks and rectangular regions of the ground-truth. L B = {l B1 , l B2 , ..., l Bi , ..., l Bn } corresponds to a set of labels obtained with our clustering methodology.

H(B, G) = 1 |G| j max 1≤k≤kopt (|b i , (b i ∈ g j ) ∧ (l Bi = k)|) |{b i ∈ g j } | (
The results of our homogeneity measure (see equation ( 9)) are presented in Table I mean homogeneity accuracy. The overall results are quite satisfying especially for the manuscript document category which contains textual (one and two fonts) and non-textual regions. The mean homogeneity accuracy is 94% for the manuscript document category (one font and graphics). One assumption can be that the manuscript documents contain drawing regions that are more compact and homogeneous than the printed document ones. By comparing the average of homogeneity measure for different document categories, we observe that a greater H(B, G) for pages containing graphics and single font texts. This yields that the extracted autocorrelation features are able to distinguish textual regions from graphical ones. We conclude that the autocorrelation descriptors are suitable to distinguish textual regions from graphical ones of the analyzed document without expressing any hypothesis either about its physical structure or its logical structure.

IV. CONCLUSIONS AND FUTURE WORK

We proposed in this paper a method for the segmentation and characterization of old document images without any a priori knowledge. The segmentation is based on a feature vector that is composed of texture indices, all based on the autocorrelation function and the directional rose. The texture features are automatically extracted in the different areas of a page and at several resolutions. The robustness of the extracted features is used in a non-parametric unsupervised clustering method which aims at determining the homogeneous regions which are assumed to be defined by similar indices of autocorrelation. Moreover, the number of clusters does not need to be known in advance as it is determined automatically.

The proposed method has been evaluated with promising results on 316 pages of old documents. We conclude that the used autocorrelation features allow a good discrimination of the foreground layers of the document, particularly of two classes: text and graphics. Results show that the autocorrelation descriptors are suitable to distinguish textual regions from graphical ones of the analyzed document.

The first aspect of future work will be to improve the classification by integrating a rejection mechanism taking into account a compactness criterion, e.g. the Mahalanobis distance, to classify all the pixels of a digitized document. Furthermore, we will elaborate frequential and statistical texture features in order to refine the segmentation and ensure the distinction between different text fonts and various graphic types. Besides, our clustering method performs classification on all pages of the document which can be very costly in case of large documents e.g. books. Thus, incremental clustering methods will be elaborated in order to reduce computation times. This work gives reason to several perspectives. For instance, by representing a document with a hierarchical layout structure, a signature can be defined for each page. This page signature can be used in similarity measure for pages for the categorization of pages and retrieval.

Fig. 1 .

 1 Fig. 1. Presentation of our automatic layout segmentation of old document images.

Fig. 2 .

 2 Fig. 2. Feature extraction.

Fig. 3 .

 3 Fig. 3. Examples of directional roses. {(a),(b),(c),(d)} are the original images and {(e),(f),(g),(h)} are respectively their roses of directions.

  y) is the variance of the intensities of the rose, except for the orientation of the maximal intensity. If F

  y) = Θ∈[10,80] I(x, y) -T Θ (0,β) (I(x, x * |tan(Θ)|)) (8) where T Θ (0,β) (I(., .)) is the translation of the analysis window of an image I by β pixels along the axis of the main angle Θ of the directional rose and Θ = F (1) (x,y) .

Fig. 4 .

 4 Fig. 4. Consensus clustering: Plot of ∆k changes in area under the cumulative density curve for the consensus matrix for each clustering experiment against number of clusters k.

Fig. 5 .

 5 Fig. 5. Example of segmentation result: (a) original grayscale image, (b) cluster representing the graphics, (c) cluster representing the text, and (d) ground-truth.

  9) where |.| is the number of pixels in the given block. B = {b 1 , b 2 , ..., b i , ..., b n } and G = {g 1 , g 2 , .

  . We find a 85% of Document category Document content Number of pages µ(H) σ(H) Max(H) Min(H)

		One font and graphics	50	0,94	0,03	0,99	0,83
	Manuscript	Two fonts and graphics	56	0,84	0,05	0,97	0,64
	H(B,G)	Only two fonts One font and graphics	50 50	0,87 0,84	0,05 0,14	0,96 1,00	0,72 0,60
	Printed	Two fonts and graphics	50	0,80	0,05	0,92	0,62
		Only two fonts	60	0,80	0,10	0,98	0,49

TABLE I .

 I HOMOGENEITY METRIC. µ(H) AND σ(H) ARE RESPECTIVELY THE MEAN VALUE AND STANDARD DEVIATION VALUES OF THE HOMOGENEITY MEASURE.

The DIGIDOC project is funded by the ANR (French National Research Agency), referenced under ANR-10-CORD-0020. For more details, see http://www.agence-nationale-recherche. fr/en/anr-funded-project/?tx lwmsuivibilan pi2[CODE] =ANR-10-CORD-0020

http://gallica.bnf.fr

V. ACKNOWLEDGMENTS

The authors would like to thank Geneviève CRON of the BnF for providing access to the Gallica digital library.


A Gracia

S Azzopardi

E Woirgard

Thermal performance evaluation of SiC power devices packaging

Devices based on wide-band gap semiconductors such as SiC, GaN allow high power densities, size reduction, high integration and elevated operating temperatures. In this study, we present solutions for high temperature power packages based on finite element simulations and experimental approach. Analytical methods are used to investigate the best choice among the different selected materials for the substrates, the solder joints and the baseplate. Various power assemblies using SiC diodes have been fabricated. Their thermal performances are evaluated through the measurements of the thermal impedance and the thermal resistance and also finite elements simulations. The correlation between the simulation results and the experimental ones are discussed.

Introduction

Silicon carbide (SiC) devices are promising components for high performance power system applications, because SiC has a band gap of 3eV which is about three times higher than conventionally used Si. The features of SiC offer the advantage of high temperature operation at more than 150 ºC which is higher than the maximum operating temperature of Si devices [START_REF] Actitty | High voltage, high temperature power electronic utility applications[END_REF]. A SiC device has been already studied with an operating temperature of 400 ºC [START_REF] Daimon | Operation of Schottky-barrier Field-effect Transistors of 3C-SiC up to 400°C[END_REF][START_REF] Chen | Demonstration of 500℃ AC Amplifier Based on SiC MESFET and Ceramic Packaging[END_REF]. High temperature operation makes it possible to eliminate the need of a cooling system, resulting in a small system size [START_REF] Lang | A Novel Chip Joint Method for High Temperature Operated SiC Power Modules[END_REF]. A high-density power system is expected to be realized by developing a high temperature module. Then to use SiC devices operating at high temperature, special care will be needed for power package. In this paper, we propose to analyse their thermal performances by experimental and simulation methods.

High temperature device packaging concepts

The selection of materials for a high temperature package (figure 1) is based on a combination of material properties: conductivity, Coefficient of Thermal Expansion (CTE), elasticity, and cost. CTE is one the most important property: then materials are selected for their high Thermal Conductivity (TC) and low CTE.

Figure 1 -Structure of power electronic module [START_REF] Dupont | Contribution a l'étude de la durée de vie des assemblages de puissance dans des environnements haute température et avec des cycles thermiques de grande amplitude[END_REF] We use a selection of materials suitable (

Experimental approach

The measurement of thermal resistance takes place in two phases using a thermal analyzer ANATECH PHASE 10 (figure 2). The first one is the calibration of the Temperature Sensitive Parameter (TSP) allowing the measurement of temperature at the level of the voltage V d , and the second is the measure of the thermal resistance. This will be explained more in detail hereafter.

Figure 2 -Test bench of thermal resistance measurement

We used the test vehicle of the whole packaging to measure the thermal resistance. The test vehicles are made by the silicon carbide diode, the die attach, the ceramic substrate with different metallization (table 1), the solder joints (SnAg 3 ) between the substrate and the baseplate, and the baseplate (CuC 40 ). The silicon carbide diode is wired by 6 aluminium wires (figure 3).

Figure 3 -Test vehicle of the whole package

The

The calibration of the Temperature Sensitive Parameter (TSP)

The thermal resistance measurement of power semiconductor is based on the calibration of the TSP. In the case of silicon carbide diode, there is a linear relationship between the on-state voltage drop across a P-N junction and temperature (Equation 1) [START_REF] Sulima | Contribution à la modélisation analytique tridimensionnelle de l'auto-échauffement dans les transistors bipolaires à hétérojonctionde type Si/SiGe[END_REF].

(1)

With, V d : voltage across the P-N junction (V) I d : direct current through the P-N junction (A) I s : saturation current (A) K : Boltzman constant (eV.K -1 ) T : local temperature (K) q : charge (e)

When the calibration is finished, we know the exact equation of V d =f(T), and it is possible to obtain at any time the junction temperature from the value of the voltage V d.

The SiC diode is connected to a current generator, a constant current "I 0 " of 1mA which can avoid the selfheating is injected into the diode. The test vehicles are placed in a dielectric liquid bath GALDEN D02, and a thermocouple is placed nearly by the TV (figure 4).

Figure 4 -Schematic description of the calibration of the TSP [7]

Firstly, the temperature of the dielectric liquid is increased gradually by the heater to 90°C. Once this temperature is reached, the heating of the bath is stopped. Then, the step of cooling down is carried out by a fan to come down to room temperature (25°C). During this period which may last several hours, several measurements are obtained at every 5°C step of temperature from 90°C down to 40°C and the voltage across the junction. We obtain the results of the calibration of TSP (figure 5) We can observe that the voltage V d decreases when the junction temperature (T j ) increases. The relationship between the junction temperature (T j ) and the voltage (V d ) can be expressed mathematically by the equation 2, the slope line is -1,9 mV/°C for SiC device:

] I I ln[ q KT V s d d  I 0 0,
(

With, T j : junction temperature (°C) m : slope of the line (V.°C -1 ) V d : voltage across the P-N junction (V) V 0 : voltage of the source coordinate (V)

Thermal resistance measurement

Once the calibration of TSP is finished, we can begin measuring the thermal resistance, which consists in two phases: one is the phase of heating the SiC diode; the other is the measurement of the thermal resistance (figure 6).

Figure 6 -Schematic descriptions of two phases to measure thermal resistance [7]

It is important to consider the safe operating area of the component, a power of about 9 W is applied by a generator to SiC diode. The test vehicle is connected by an electronic card interface to measure the thermal resistance of whole package. A thermocouple is placed between the controlled cold plate (T ref =25°C) and the interlayer double sided adhesive. The thermal resistance of the assembly is the sum of all thermal resistances of each material (figure 7).

Figure 7 -schematic descriptions of different thermal resistances for a test vehicle (heat transfer by conduction)

We can define: -R th-JC as the junction to case thermal resistance. It is defined between the junction temperature T j (maximum temperature in the chip) and temperature T c (bottom surface temperature of the baseplate).

-R th-INT as the thermal resistance of the adhesive layer. The test vehicle is fixed on the heat spreader with an interlayer double-sided adhesive.

-R th-HS as the thermal resistance between the heat spreader and the air.

The table below presents the technical characteristics of an interlayer double-sided adhesive:

Thermal conductivity : 1 W/m.K Thickness : 0.15 mm Thermal resistance : 0.17 °C/W * Operating temperature : -45°C / +170°C Isolation voltage :

3 kV *we take into account the substrate surface. Table 3 -Specifications for the thermal interface doublesided adhesive [START_REF]Radiosprares components[END_REF]

Experimental results

Thermal resistance

We obtain the results of thermal resistance for the matrix (Table 1) of test vehicles, and have been summarized in the table 4: The thermal resistance results (table 4) illustrate that TV2 model (with die attach solder AuGe 12 , substrate ceramic Si 3 N 4 Cu metallization) has a low thermal resistance (1.82°C/W), it might have a good performance assembly. TV6 model (with die attach solder SnPb 95 , substrate ceramic AlN Al metallization) shows the heat dissipation is less effective compared to others. So it appears that the solder AuGe 12 seems to be the best choice in terms of thermal performance for configurations of identical metallic substrates.

TVs Power (W) T j (°C) T ref (°C) V j (V) R th (°C/W)

Numerical simulations approach

The general design strategy is to choose a set of materials that appears to match by looking at the material properties (table 5), and then develop a package which drives these properties to the best use.

Material Using ANSYS finite elements software, a 3D thermal simulation of power assembly has been built. For reducing simulation time, only a quarter of the power assembly is simulated. Three-dimensional solid elements (linear brick and tetrahedron) were used for meshing the geometry (figure 11).The interfaces being supposed perfect. The boundary condition is imposed on the bottom surface of the assembly (baseplate) with a constant temperature 25°C, this temperature is the reference temperature T ref .

A power of about 9 W is injected into SiC diode.

Thermal simulations results

The junction temperature T j is the maximum temperature in the die, localized at centre of upper part of the SiC device (figure 12).

The heat flux balance across assemblies and the thermal resistances were evaluated from the thermal simulations. The thermal resistances were calculated from the chip maximal temperature and the reference temperature, measured when the steady state is reached (Equation 4).

Correlation between the simulation results and the experimental results

Figure 14 shows the comparison between experimental results and simulation ones for different models of thermal resistance assembly.

Figure 14 -Comparison between experimental results and simulation ones for different models of thermal resistance assembly

The two curves show that they have the same trend for the different models. But the experimental results have higher values than the simulation results due to the thermal resistance of double-sided adhesive thermal interface not include in the simulation. Additionally, the interfaces of assembly are not perfect, it also increases the thermal resistance values.

Conclusion

This study is to develop high temperature packaging technologies based on available materials, using ANSYS finite elements software to propose design rules and carry out the thermal simulation. We have taken into account the different selections of materials with the most critical design parameters to define the model. Thermal resistances are calculated by thermal simulation and have been measured by experimental approach in order to verify the simulation results and thermal performances. Thermal impedance method is used to determine the failure of the assembly after aging test showing the sensibility of test vehicle to thermal fatigue and the weakness of the substrate to baseplate report process.

Figure 5 -

 5 Figure 5 -Calibration results of TSP (voltage-temperature relationship)

Figure 11 -

 11 Figure 11 -Meshed geometric model assembly

Figure 12 -Figure 13 -

 1213 Figure 12 -Temperature distributions in Kelvin in 3D model assembly Thermal simulation results (R th ) with different selections of materials (table 1) is summarized in the figure 13.

Table 1 -

 1

	SiC die Die attach Substrate Metallization
	TV1 Diode	SnPb 95	Si 3 N 4	Cu
	TV2 Diode	AuGe 12	Si 3 N 4	Cu
	TV3 Diode	AuGe 12	AlN	Al
	TV4 Diode	AuSi 3	Si 3 N 4	Cu
	TV5 Diode	AuSi 3	AlN	Al
	TV6 Diode	SnPb 95	AlN	Al

table 1) for high temperature package assembly, including die attach (Au 88 Ge 12 , Au 97 Si 3, Sn 5 Pb 95 ) , substrate (Si 3 N 4 , AlN), wire bonding (Al, Au) and baseplate (CuC, Cu-diamond).Various materials combinations are proposed in order to try to find and optimize the performance and the reliability of the assembly. From a thermal point of view, the thermal resistance (R th ) of power package can be considered as one of the key parameters for high temperature device packaging concepts. Matrix of test vehicles with SiC diodes

Table 2 -

 2 geometrical parameters of the test vehicle are indicated in the table below: The geometric parameters of the test vehicle

		Material	Surface (mm 2 )	Thickness (μm)
	wire bond	Al	125μm diameter	-
	SiC chip	SiC	3.64*3.64	400
	die attach	AuGe 12 /AuSi 3	3.8*3.8	50
		/SnPb 95
	substrate	Si 3 N 4	30*30	320
		AlN	30*30	1000
	metallization	Cu	-	300
		Al	-	400
	solder joint	SnAg 3	30*30	100
	baseplate	CuC 40	40*40	4000

Table 4 -

 4 Experimental results of thermal resistance for different models of assemblies

	d V		m		T	j		0 V

Table 5 -

 5

Thermal properties of the elements

[9][START_REF] Coppola | Survey on Hign-Temperature Packaging Materials for SiC-Based Power Electronics Modules[END_REF]


Parametrization of matrix-valued lossless functions based on boundary interpolation

Ralf Peeters and Martine Olivi and Bernard Hanzon

Abstract-This paper is concerned with parametrization issues for rational lossless matrix valued functions. In the same vein as previous works, interpolation theory with metric constraints is used to ensure the lossless property. We consider here boundary interpolation and provide a new parametrization of balanced canonical forms in which the parameters are angular derivatives. We finally investigate the possibility to parametrize orthogonal wavelets with vanishing moments using these results.

I. INTRODUCTION

Lossless systems and their transfer functions play a central role in system theory mainly, but not only, due to the Douglas-Shapiro-Shields factorization. Lossless matrixvalued functions generalize, in some sense, the notion of denominator to the matrix case. As such they are involved in many representations of stable systems. They are also present in many applications, as the scattering matrix of a resonant filter or the polyphase matrix of a filter bank, for example. Their parametrization is an important issue for many purposes going from optimization, model reduction to physical parameters recovery. In the same vein as previous works [START_REF] Hanzon | Balanced Parametrizations of Stable SISO All-Pass Systems in Discrete-Time[END_REF], [START_REF] Hanzon | Balanced realizations of discrete-time stable all-pass systems and the tangential Schur algorithm[END_REF], we aim to stress the connections between two complementary descriptions: canonical forms on one hand and Schur type algorithms on the other hand. Thereby, we get parametrizations which combine the computational interest of canonical forms and the conceptual interest of Schur analysis.

In [START_REF] Hanzon | Balanced realizations of discrete-time stable all-pass systems and the tangential Schur algorithm[END_REF], by making appropriate choices in the Schur algorithm, construction of balanced canonical forms is achieved in a recursive way with unitary matrix multiplications. Each step of the recursion involves an interpolation condition at a point located within the analyticity domain. However, for several reasons, it may prove helpful to allow for interpolation points on the boundary of the analyticity domain (the circle in discrete-time and the imaginary axis in continuoustime). In [START_REF] Olivi | Lossless scalar functions: boundary interpolation, Schur algorithm and Ober's canonical form[END_REF], the tridiagonal canonical form of Ober in continuous-time was obtained from a recursion involving interpolation conditions at infinity. Vanishing moments, diagonal Markov parameters, can also be interpreted in term of boundary interpolation conditions.

In this paper, we extend the work of [START_REF] Hanzon | Balanced realizations of discrete-time stable all-pass systems and the tangential Schur algorithm[END_REF] on discrete-time lossless systems to allow for interpolation conditions on the R. Peeters

II. PARAMETRIZATION OF LOSSLESS MATRICES FROM

THE TANGENTIAL SCHUR ALGORITHM.

In this section, we briefly recall the results of [START_REF] Hanzon | Balanced realizations of discrete-time stable all-pass systems and the tangential Schur algorithm[END_REF]. The tangential Schur algorithm recursively handles an interpolation condition of the form

F (1/ w)u = v, (1)
in which w ∈ D, u and v are m-vectors such that u = 1 and v < 1.

The solutions of this interpolation problem (see e.g. [START_REF] Ball | Interpolation of rational matrix functions[END_REF]) are given by means of a linear fractional transformation (LFT) of the form

T Θ (G) = (Θ 22 G + Θ 21 )(Θ 12 G + Θ 11 ) -1 , (2)
where

Θ(z) = Θ 11 (z) Θ 12 (z) Θ 21 (z) Θ 22 (z)
is a rational J-inner function of size 2m × 2m. More precisely, Theorem 2.1: All the rational lossless functions F (z) satisfying the interpolation condition (1) can be represented as a linear fractional transformation (LFT)

F = T Θw,u,v (G) (3)
of a lossless function G(z). The symbol Θ w,u,v is the Jinner function given, in terms of the interpolation data, by

Θ w,u,v (z) =     I 2m + b w (z) b w (ν) -1 u v u v * J 1 -v 2     H, (4)
where b w (z) is the Blaschke factor

b w (z) = z -w 1 -wz , (5)
and J the matrix

J = I m 0 0 -I m . (6)
The unit complex number ν and the constant J-unitary matrix H (HJH * = J) are arbitrary and should be specified. Moreover, we have that

deg F = deg G + 1. Proof: see [4]
The tangential Schur algorithm yields from a lossless function F (z) of McMillan degree n a finite sequence of interpolation vectors associated with -a sequence w = (w n , w n-1 , . . . , w 1 ) of points of the unit disk -a sequence u = (u n , u n-1 , . . . , u 1 ) of interpolation direction vectors It works as follows: let

F n = F . Then, for k = n, n -1, . . . put v k = F k (1/ wk )u k ; • if v k < 1, let F k-1 = T -1 Θw k ,u k ,v k F k , (7)
in which

Θ w k ,u k ,v k (z) is defined by (4). • if v k has norm 1, then stop.
A local coordinate map ϕ(w, u) can be associated with the sequences w and u by

ϕ(w, u) : F (z) ∈ V(w, u) → (v 1 , v 2 , . . . , v k , F 0 ),
where the domain V(w, u) of the map consists with functions for which

v k < 1 for k = n, n -1 . . . , 1.
The freedom in the choice of ν and H has been used in the literature to serve various purposes. In [START_REF] Hanzon | Balanced realizations of discrete-time stable all-pass systems and the tangential Schur algorithm[END_REF] it is used to get a simple computation of balanced canonical forms. We introduce the Halmos extension H(E) of a contractive matrix E

H(E) = (I -EE * ) -1/2 E(I -E * E) -1/2 E * (I -EE * ) -1/2 (I -E * E) -1/2 = (I -EE * ) -1/2 (I -EE * ) -1/2 E (I -E * E) -1/2 E * (I -E * E) -1/2 . ( 8
D C B A = V 0 0 I n-1   D 0 C 0 1 0 B 0 A   U * 0 0 I n-1 , (10)
in which U and V are the unitary matrices

U = I m -(1 + w κ(w, v))uu * α(w, v) u α(w, v) u * wκ(w, v) , (11)
V = I m -(1 -κ(w, v)) vv * v 2 α(w, v) v -α w,v v * κ(w, v) (12)
with

α(w, v) = 1 -|w| 2 1 -|w| 2 v 2 , ( 13
) κ(w, v) = 1 -v 2 1 -|w| 2 v 2 . ( 14
)
Proof: see [START_REF] Hanzon | Balanced realizations of discrete-time stable all-pass systems and the tangential Schur algorithm[END_REF] Note that formula (9) is obtained from (4) by choosing

ν = 1 + w 1 + w , H = H(uv * ) -1 H( wuv * ).

III. PARAMETRIZATION OF LOSSLESS MATRICES FROM BOUNDARY DATA

Boundary interpolation of rational Schur functions is studied in [START_REF] Ball | Interpolation of rational matrix functions[END_REF]. We give here a specific treatment for lossless matrix functions. We are thus interested in rational lossless functions F that satisfy an interpolation condition at a given point ξ on the unit circle T. As it turns out, also the first derivative of F at ξ plays a role.

A. The angular derivative

Note that at a point ξ ∈ T it holds that F (ξ) is unitary, so that F (ξ) -1 = F (ξ) * . This gives the identity

ξF (ξ) * F ′ (ξ) = ξF ′ (ξ) * F (ξ),
which proves that the expression ξF (ξ) * F ′ (ξ) is actually Hermitian.

This also shows that a lossless function F and its derivative F ′ have a particular relationship at a given point ξ on the unit circle. More precisely, we have the following result.

Lemma 3.1: Let F (z) be a m × m rational lossless function, ξ ∈ T a point on the unit circle, and u ∈ C m a given vector. Then the scalar expression

ρ = -ξu * F (ξ) * F ′ (ξ)u, (15)
is nonnegative real: ρ ≥ 0.

Proof: To show this, we may use the Potapov decomposition [START_REF] Potapov | The multiplicative structure of J-contractive matrix functions[END_REF] of

F F = U 0 B w1,u1 B w2,u2 . . . B wn,un (16)
where w 1 , w 2 , . . . , w n belong to the open unit disk, u 1 , u 2 , . . . , u n are unit vectors,

B w,u = I + b -1 w -1 uu * , b w (z) -1 = 1 -wz z -w , (17)
are elementary (degree one) Potapov factors, and U 0 is a constant unitary matrix. First note that F ′ is computed as:

F ′ = U 0 n k=1 k-1 l=1 B w l ,u l (b -1 w k ) ′ u k u * k n l=k+1 B w l ,u l ,
in which we have that

(b -1 w k ) ′ (z) = - 1 -|w k | 2 (z -w k ) 2 .
At the point ξ ∈ T we shall introduce the notation

U k = B w k ,u k (ξ).
Note that b w k (ξ) has modulus 1, and that U k is unitary. We now have:

F (ξ) = U 0 U 1 . . . U n , and
F ′ (ξ) = - n k=1 1 -|w k | 2 (ξ -w k ) 2 k-1 l=0 U l u k u * k n l=k+1 U l ,
Therefore, using the relation

U * k u k u * k = b -1 w k (ξ)u k u * k = 1 -ξw k ξ -wk u k u * k ,
we get:

-ξF (ξ) * F ′ (ξ) = n k=1 1 -|w k | 2 |ξ -w k | 2 Û * k u k u * k Ûk , with Ûk = n l=k+1 U l .
Obviously this expression is positive semi-definite, as it is a positively weighted sum of positive semi-definite matrix products. Therefore the scalar quantity ρ = -ξu * F (ξ) * F ′ (ξ)u is nonnegative real. Note that from the available expression for ξF ′ (ξ)F (ξ) * we can quickly deduce the conditions under which ρ becomes equal to zero. This only happens when each of the terms in the weighted sum becomes equal to zero. Note that the weights are all strictly positive. Hence: ρ = 0 occurs if and only if

u * n u = 0, u * n-1 U n u = 0, u * n-2 U n-1 U n u = 0, . . . . . . . . . u * 1 U 2 . . . U n-1 U n u = 0.
But because of the structure of the matrix U n , the condition u * n u = 0 implies that U n u = u, whence the condition u * n-1 U n u = 0 becomes equivalent to the condition u * n-1 u = 0. Continuing this argument it follows that ρ = 0 if and only if

u * n u = u * n-1 u = • • • = u * 1 u = 0.
This allows us to draw the following conclusions. (i) In the case n = 0 no conditions apply and ρ = 0. Indeed, the inner function F (z) = U 0 is constant and its derivative F ′ (z) is the zero matrix. (ii) In the scalar case (with n > 0) the fact that u n = 1 implies that ρ > 0 if the vector u is nonzero. (See also [START_REF] Ball | Interpolation of rational matrix functions[END_REF][Lemma 21.1.1.]) (iii) In the multivariable case (with n > 0) each of the conditions u k * u = 0 generically does not hold, so that generically ρ > 0. If u is to be chosen from a basis for the space C p , then at least one such choice exists which achieves u k * u = 0. (iv) In the multivariable case (with n > 0) if the vectors u 1 , u 2 , . . . , u n happen to span the space C p (which requires n ≥ p), then ρ > 0 for every nonzero vector u . Otherwise, (and certainly for 0 < n < p) a nontrivial subspace of vectors u does exist for which ρ = 0. (v) The conditions under which ρ = 0 does or does not hold, are completely independent of the value of the interpolation point ξ. Therefore, a suitable choice of u which achieves strict positivity of ρ can be made independently of the choice of an interpolation point ξ .

B. Boundary interpolation of lossless functions.

A well-posed boundary interpolation problem includes the value of the angular derivative and can be stated as follows. Find all the lossless functions F (z) such that

F (ξ)u = v v * F ′ (ξ)u = -ξρ (18)
where ξ belongs to the unit circle, u and v are unit vectors and ρ > 0 is the angular derivative. Any solution of the (general angular derivative) interpolation problem (18) can be represented as a linear fractional transformation of a lossless function G such that G(ξ)u = v and the degree of G is equal to the degree of F minus one. The J-inner symbol associated with the LFT is, up to a J-unitary constant factor

Φ ξ,u,v,ρ = I 2m - 1 2ρ c ξ u v u v * J, (19)
in which c ξ is the Caratheodory function

c ξ (z) = ξ + z ξ -z . (20)
Note that an elementary J-inner factor which is analytic in the open unit disk must be either of the form (4) or of the form (19) if its pole is on the disk [START_REF] Dym | J-contractive matrix functions, reproducing kernel spaces and interpolation[END_REF]. Theorem 3.1: Part I -increasing the degree Let ξ ∈ T, u ∈ C m , u = 1, v ∈ C m with v = 1, and ρ ∈ R + . Let Φ ξ,u,v,ρ be the elementary J-inner factor given by ( 19). Let G be a m × m rational lossless function of McMillan degree n such that G(ξ)u = v. Let F be defined by

F = T Φ ξ,u,v,ρ (G). Then F is a m × m rational lossless function of McMillan degree n + 1 such that F (ξ)u = v, v * F ′ (ξ)u = -ξρ.
Part II -decreasing the degree Let F be a m × m rational lossless function of McMillan degree n + 1 ≥ 1. Let ξ ∈ T and let u ∈ C m with u = 1 be such that ρ := -ξu * F (ξ) * F ′ (ξ)u is strictly positive. Let v := F (ξ)u (which satisfies v = 1), and let Φ ξ,u,v,ρ be the elementary J-inner factor (19). Let G be defined by

G = T -1 Φ ξ,u,v,ρ (F ). Then G is a m × m rational lossless function of McMillan degree n such that G(ξ)u = v.
Proof: Part I. The linear fractional transformation F = T Φ ξ,u,v,ρ (G) can be worked out as follows:

F = G + 1 2ρ c ξ v(v * G -u * ) I m + 1 2ρ c ξ u(v * G -u * ) -1 = G + 1 2ρ c ξ v(v * G -u * ) I m - c ξ 2ρ + c ξ (v * Gu -1) u(v * G -u * ) yielding F (z) = G(z) - (z + ξ)(G(z)u -v)(v * G(z) -u * ) 2ρ(ξ -z) + (z + ξ)(v * G(z)u -1)
Therefore we have that

F (z)u = v + 2ρ(ξ -z)(G(z)u -v) 2ρ(ξ -z) + (z + ξ)(v * G(z)u -1)
.

At the point z = ξ it follows that F (ξ)u = v, which shows the first interpolation condition to be satisfied. For the second interpolation condition we have to perform differentiation. At the point z = ξ we have:

F ′ (ξ)u = - ρ(G(ξ)u -v) ξ(v * G(ξ)u -1) . (21)
(Note that division by zero does not occur, because of the requirement G(ξ)u = v). Upon premultiplication by v * it is obtained that:

v * F ′ (ξ)u = -ξρ,
which shows the second interpolation condition to be satisfied.

Part II.

The linear fractional transformation

F = T Φ -1 ξ,u,v,ρ
(Q) can be worked out as under Part I, upon noting that algebraically the inverse of the J-inner matrix is obtained as

Φ -1 ξ,u,v,ρ = Φ ξ,u,v,-ρ
i.e., by reversing the sign of ρ. Therefore:

G(z) = F (z) + (z + ξ)(F (z)u -v)(v * F (z) -u * ) 2ρ(ξ -z) -(z + ξ)(v * F (z)u -1)
From the definition v = F (ξ)u it is clear that pole-zero cancellation takes place at z = ξ. In fact, a double pole-zero cancellation takes place at that point. To see this, we proceed as follows. As under Part I it holds that

G(z)u = v + 2ρ(ξ -z)(F (z)u -v) 2ρ(ξ -z) -(z + ξ)(v * F (z)u -1)
.

which can be rewritten as:

G(z)u = v + 2ρ F (z)u-v z-ξ 2ρ+(z+ξ) v * F (z)u-1 z-ξ z-ξ
.

From the definitions F (ξ)u = v and ρ = -ξv * F ′ (ξ)u, we now have that:

G(z)u = v + 2ρF ′ (ξ)u -2ξv * F ′ (ξ)u+(z+ξ) v * (F (z)-F (ξ))u z-ξ z-ξ = v + 2ρF ′ (ξ)u v * F ′ (ξ) + (z + ξ) F (z)-F (ξ)-F ′ (ξ)(z-ξ) (z-ξ) 2 u .
At the point z = ξ it follows that :

G(ξ)u = v + 2ρF ′ (ξ)u -ξρ + ξv * F ′′ (ξ)u .
Upon multiplication by v, this yields

v * G(ξ)u = 1 + 2 ξρ 2 ξρ -ξv * F ′′ (ξ)u
where the second term is nonzero (as it has a nonzero numerator), so that it follows that G(ξ)u = v.

C. Computing balanced realizations

In [START_REF] Peeters | Linear fractional transformations and balanced realizations of discrete time stable allpass systems[END_REF], a unified framework is presented in which linear fractional transformations on transfer functions are represented by corresponding linear fractional transformations on state-space realization matrices.

D Ĉ B Â = T ∆ (R e ),
where

R e =   D 0 C 0 I m 0 B 0 A   is an extended realization matrix for G and ∆ =         D 11 C 1 0 D 12 0 0 B 1 A 0 B 2 0 0 0 0 I n 0 0 0 D 21 C 2 0 D 22 0 0 0 0 0 0 I m 0 0 0 0 0 0 I n         Proof: see [6, th. 8]
This theorem can be applied directly to the situation at hand. A (J-balanced) state-space realization (A, B, C, D) of Φ ξ,u,v,ρ is given by

A = ξ B = 1 √ ρ u v * J C = -ξ √ ρ u v D = I 2m -1 √ 2ρ u v u v * J. Theorem 3.3: A balanced state-space realization ( Â, B, Ĉ, D) for F = T Θ (G) of McMillan degree n + 1 is computed from a given balanced realization (A, B, C, D) of G by D Ĉ B Â =   D 0 C 0 -ξ 0 B 0 A   +   Du -v 2 √ ρξ Bu     D * v -u 2 √ ρ C * v   * 1 -v * Du + 2ρ (
v k = F k (ξ k )u k , ρ k = -ξ k v * k F ′ (ξ k )u k . If ρ k > 0, then put F k-1 = T -1 Φ ξ k ,u k ,v k ,ρ k (F k ), (23)
otherwise stop. Recall that generically, ρ k > 0, as explained in section III-A. The sequences (v 1 , v 2 , . . . , v n ) and (ρ 1 , ρ 2 , . . . ρ n ) can be used as parameters, the function F (z) is then computed by reversing the algorithm. However, all the sequences of v i 's are not allowed since the conditions

v k = F k-1 (ξ k )u k , k = 1, . . . , n,
must be satisfied.

IV. BOUNDARY INTERPOLATION AS A LIMIT OF CLASSICAL INTERPOLATION

We now prove that, for an appropriate choice of the Junitary constant factor, the linear fractional transformation for the boundary problem arises through a suitable limiting process from a sequence of LFT's for interior problems. With this choice, (22) recovers the form described in [START_REF] Hanzon | Balanced realizations of discrete-time stable all-pass systems and the tangential Schur algorithm[END_REF], that is a product of unitary matrices. In this section, we normalize the Blaschke product as

ζ w (z) = - |w| w b w (z). ( 24
)
and the J-inner function (4) as

Θ w,u,v (z) = I 2p + (ζ w (z) -1) 1 -v 2 u v u v * J. (25)
It is obtained by choosing ν = -w |w| and H = I 2m in (4). Lemma 4.1: The J-inner matrix Θ w,u,v given by ( 25) is related to the J-inner matrix Θw,u,v through the relation

Θw,u,v (z) = Θ w,u,v (z)H 1 -|w| 1 -|w| v 2 uv * I p + (-w |w| -1)uu * 0 0 I p ,
where H(E) denotes the Halmos extension ( 8) of a contractive matrix E.

Proof: The two J-inner functions may only differ by a right J-unitary constant matrix H

Θw,u,v (z) = Θ w,u,v (z)H, which satisfies H = Θw,u,v -w |w| . Note that b w -w |w| = -w
|w| and that for every contractive matrix E and unitary matrix P , we have

P 0 0 I H(E) = H(P E) P 0 0 I
We thus have

H = H(uv * ) I p + (-w |w| -1)uu * 0 0 I p H( wuv * ) = H(uv * )H(-|w|uv * ) I p + (-w |w| -1)uu * 0 0 I p It is easily verified that H(uv * )H(-|w|uv * ) = H 1 -|w| 1 -|w| v 2 uv *
Proposition 2: Let Θ w,u,v be the J-inner function defined by (25) and F be a rational lossless function; let ξ ∈ T, and

v = F (ξ)u, ρ = -ξv * F ′ (ξ)u. ( 26
)
Assume ρ > 0 and let Φ ξ,u,v,ρ be the J-inner function defined by (19). We then have

lim w→ξ 1 -F (1/ w)u 2 1 -|w| 2 = ρ ( 27
)
and for z = ξ,

lim w→ξ Θ w,u,F (1/ w)u (z) = Φ ξ,u,v,ρ (z) (28) lim w→ξ Θw,u,F (1/ w)u (z) = Φ ξ,u,v,ρ (z) Ĥ, (29)
where

Ĥ = H uv * 1 + 2ρ I m -(1 + ξ)uu * 0 0 I m Proof: Write the Taylor series of F (z) about ξ ∈ T F (z) = F (ξ) + (z -ξ)F ′ (ξ) + o(|z -ξ|).
Using (26) we get

1 -F (z)u 2 = 1 -u * F (z) * F (z)u = -2Re [(z -ξ)v * F ′ (ξ)u] + o(|z -ξ| 2 ),
and since lim

w→ξ 1 -|w| Re(ξ/w -1) = 1
we get (27). Note that

ζ w (z) -1 = (|w| -1) |w|z + w w -|w| 2 z , so that lim w→ξ ζ w (z) -1 1 -F (1/ w)u 2 = - 1 2ρ c ξ (z)
which achieves the proof of (28). Writing

1 -|w| v 2 = 1 2 [(1 -|w|)(1 + v 2 ) + (1 + |w|)(1 -v 2 ),
we get

lim w→ξ 1 -|w| v 2 1 -|w| = 1 + 2ρ,
which gives (29).

D Ĉ B Â = V 0 0 I n   D 0 C 0 1 0 B 0 A   Û * 0 0 I n .
(30) where Û and V are unitary (p+1)×(p+1) complex matrices depending on ξ, u, v, ρ as follows

Û =   1 1+ρ u I p -(1 + ξ ρ 1+ρ )uu * ξ ρ 1+ρ 1 1+ρ u *   , (31) V =   1 1+ρ v I m -(1 - ρ 1+ρ )uu * ρ 1+ρ - 1 1+ρ v *   . (32)
Proof: As previously, the limits of ( 13) and ( 14) are easily computed

lim w→ξ α(w, v) = 1 1 + ρ , lim w→ξ κ(w, v) = ρ 1 + ρ ,
and the matrices U and V in Theorem 1 converge to Û and V . Formula (30) can also be checked directly using Theorem 3.2.

V. PARAMETRIZATION OF ORTHOGONAL WAVELETS WITH VANISHING MOMENTS

We start from the theory of orthogonal wavelets as discussed in [START_REF] Strang | Wavelets and filter banks Wellesley Cambridge press[END_REF]. However, we do not require the wavelets to have a compact support. We consider a multi-resolution structure for L 2 (R), consisting of a nested sequence of linear subspaces V k , k ∈ Z satisfying the properties (i) nesting:

∀ k ∈ Z, V k ⊂ V k+1 (ii) empty intersection: ∩ k∈Z V k = 0 (iii) completness: ∪ k∈Z V k = L 2 (R) (iv) shift invariance: ∀ f ∈ L 2 (R), ∀ k, m ∈ Z, f (t) ∈ V k ⇒ f (t-m) ∈ V k (v) scale invariance: ∀ f ∈ L 2 (R), ∀ k ∈ Z, f (t) ∈ V k ⇒ f (2t) ∈ V k+1
The spaces V k are called approximation spaces. The detail spaces W k are the orthogonal complement of V k in the enveloping spaces V k+1 . We shall hypothesize that there exists a function φ(z) (resp. ψ(z)), called the scaling function (resp. wavelet function), which generates a shift invariant orthonormal basis {φ(t -m), m ∈ Z} (resp. {ψ(t -m), m ∈ Z}) for the approximation space V 0 (resp. W 0 ). It follows from the properties of the multi-resolution structure that there exists sequences of coefficients C k and D k , k ∈ Z, which allows to represent the scaling and the wavelet functions as a linear combination of the orthogonal basis functions of V 1 :

φ(t) = √ 2 k∈Z C k φ(2t + k), ψ(t) = √ 2 k∈Z D k φ(2t + k), (33)
with

C k = √ 2 φ(t)φ(2t + k) dt, D k = √ 2 ψ(t)φ(2t + k) dt.
Orthonormality of the basis for the spaces V k and W k allows us to derive important properties of the coefficient sequences (C k ) and (D k ). The reverse approach, which consists in starting from such coefficients and looking for the scaling and wavelet functions, is a highly non-trivial task that has received quite some attention in the literature. When a signal s(t) is to be decomposed into an approximation signal a(t) ∈ V k and a detail signal d(t) ∈ W k satisfying s(t) = a(t) + d(t), this can be performed directly from the coefficient sequences of these signals with respect to the various orthonormal bases. Let

s(t) = √ 2 k∈Z s k φ(2t -k) (34)
a(t) = k∈Z a k φ(t -k) (35)
d(t) = k∈Z d k φ(t -k). (36)
If we introduce the filters

C(z) = k∈Z C k z -k , D(z) = k∈Z D k z -k , (37)
it follows that

a k = m∈Z C 2k-m s m , d k = m∈Z D 2k-m s m .
The filters C(z) and D(z) are first used to filter the sequence (s k ), after which down-sampling is applied to the outcome. An equivalent procedure, by first splitting the sequence (s k )

into its even and its odd phase, can be interpreted in terms of the associated polyphase filter H(z). It is given by

H(z) = C e (z) C o (z) D e (z) D o (z) (38)
where

Ce(z) = P k∈Z C 2k z -k De(z) = P k∈Z D 2k z -k Co(z) = P k∈Z C 2k+1 z -k Do(z) = P k∈Z D 2k+1 z -k (39)
Filtering the input sequence u k = s 2k s 2k+1 by H(z) produces the output sequence

y k = a k d k .
The orthogonality conditions on the filter coefficients are in fact equivalent to the condition that the polyphase matrix H(z) is lossless. This provides us with a very useful and concise description of the class of filter banks leading to orthogonal wavelets, which allows us to parametrize this class directly with orthonormality built in.

A wavelet ψ(t) is said to have a vanishing moment of order p, if it holds that

t p ψ(t) dt = 0. ( 40
)
If a wavelet has m vanishing moments, for order p = 0, 1, . . . , m -1, then the scaling and wavelets functions have continuous derivative up to order m -1. These are very desirable properties in view of practical applications. A vanishing moment of order 0 is a standard admissibility condition for all wavelets. One famous class of orthogonal wavelets exhibiting vanishing moments, is the class of Daubechies wavelets. This class is obtained by restricting to FIR filters C(z) and D(z) of finite order 2n -1, which allows to built in n vanishing moments. We shall see that vanishing moment conditions can be reformulated in terms of boundary interpolation conditions for the lossless polyphase matrix function H(z).

We introduce the notation

v p = t p φ(t) dt, w p = t p ψ(t) dt.
Note that orthonormality of the basis impose that v 0 = 1.

Multiplying the equations (33) by t p and integrating, we get the following relations

v m w m = m p=0 m p v m-p C (p) (1) D (p) (1) , (41)
where m p is the binomial coefficient and C(z) and D(z) the filters (37) and (.) (p) denotes the pth derivative.

To derive interpolation conditions on the polyphase matrix H(z) we shall use the equality

C(z) D(z) = H(z 2 ) 1 z -1

A. One vanishing moment

For m = 0, (41) yields

1 w0 v0 = H(1) 1 √ 2 1 √ 2
.

Since H(z) is lossless, it takes orthogonal values at every point on the unit circle, and in particular H(1) is orthogonal. It is thus norm preserving and we must have w 0 = 0. Therefore, a vanishing moment of order 0 is automatically built. It also follows that

H(1) = 1 √ 2 1 √ 2 -1 √ 2 1 √ 2
, up to a change of sign in the second row.

If we impose the McMillan degree of H(z) to be 1, then only a single solution exists: H(z) = H(1). It is the polyphase matrix associated with the Haar wavelet.

To proceed with vanishing moments of higher order, we focus our attention on the lossless matrix

H(z) = H(z)H(1) T ,
which satisfies the property H(1) = I 2 .

B. Two vanishing moments

For m = 1 and imposing w 1 = 0, (41) yields, after some computations, the following interpolation condition for the first derivative of H(z)

H(1) (1) 1 0 = 1 4 2v 1 + 1 -1 = -ρ 1 -1 4 . (42)
From section III-B, we have a description of all the lossless functions satisfying boundary interpolation conditions of the form (18). The lossless matrix H satisfies such a condition

with ξ = 1, u = v = 1 0 and ρ 1 = 1 4 β 1 = -1 4 (2v 1 + 1)
whihc is a positive number if and only if v 1 < 1/2. Then β 1 appears as a free parameter and we may use Theorem 3.1 to describe the set of solutions in terms of a linear fractional transformation of an arbitrary lossless function G. The condition for a vanishing moment of order 1 happens to be 0 1 H(1) (1)

1 0 = - 1 4 . (43)
Using the symmetry of H(1) (1), this condition translates into an equivalent condition on the structure of the matrix G(1). Proposition 3: Any 2 × 2 lossless function H satisfying H(1) = I 2 together with (43) can be represented as H = T Φ1 (G), where Φ 1 is given by (19) with ξ = 1, ρ = ρ 1 and u

= v = 1 0 T and G(z) is a 2 × 2 lossless function satisfying G(1) =   1-β 2 1 1+β 2 1 -2β1 1+β 2 1 -2β1 1+β 2 1 - 1-β 2 1 1+β 2 1   , β 1 = 4ρ 1 (44)
Proof: The value of G( 1) is obtained form (43) using (21) and the fact that G(1) is orthogonal.

If we consider the case where the McMillan degree of H equals 1, and there are two vanishing moments built in, we then arrive at an explicit parametrization

H(z) = I 2 - 2(z -1) z(1 + 4β 1 + β 2 1 ) -(1 -4β 1 + β 2 1 ) β 1 1 β 1 1 T .
This shows that the pole of H is located at z =

1-4β1+β 2 1 1+4β1+β 2 1
, which for β 1 > 0 has a real value in the interval ] -1, 1[. We arrive at a lossless FIR system if and only if this pole is located at z = 0. This implies that β 1 = 2 ± √ 3. The first solution corresponds to the polyphase matrix associated with the Daubechies db2 wavelet, the second solution provides an alternative with reversed filter coefficients.

C. Three vanishing moments

We define the lossless matrix G by G(z) = G(z)G [START_REF] Ball | Interpolation of rational matrix functions[END_REF] where G(1) is given by (44). For m = 2 and imposing w 2 = 0, (41) yields the following condition for a vanishing moment of order 2 1 0 H(2) (1) 0 1 = 1 4

β 1 2 + 1 . ( 45
)
This condition translates into the equivalent condition on G .

1 0 G(2) (1) 1 -β 1 = β 1 4 1 -β 2 1 1 -β 2 1 . (46
Proof: The proof relies on technical computations with Mathematica.

If we consider the case where the McMillan degree of H equals 2, and there are three vanishing moments built in, we then arrive at an explicit parametrization in terms of ρ 1 and ρ 2 . Imposing the lossless function to be FIR, and working out the expressions with Mathematica, we find two pairs of values (ρ 1 , ρ 2 ) one of which gives the Dauberchies db3 wavelet.

VI. CONCLUSION AND FUTURE WORKS

In this paper we have got parametrization of (balanced realizations) of lossless functions in which the parameters are boundary interpolation value and angular derivatives. This is interesting form an application viewpoint since it seems that a number of meaningful quantities (vanishing moments, Markov parameters, etc.) can be interpreted in terms boundary interpolation. In this connection, we have got explicit parametrizations of 2 × 2 polyphase matrices of arbitrary order n with 4 vanishing moments built in (the results have been presented up to 3), in terms of angular derivative (positive) parameters. However, the conditions are cleverly handled in an unusual recursive fashion that we still do not completely understand. In our opinion, restriction to the FIR class would be an interesting target for further research.

These results could also be used to derive parametrizations for Schur functions. Indeed, a Schur function can be viewed as a sub-block of a lossless function of higher dimension. A parametrization of Schur functions could thus be deduced from that of lossless functions. The use of boundary interpolation could be proved necessary in this approach to handle functions which are not strictly Schur.

  is with the Department of Knowledge Engineering, Universiteit Maastrich, PO Box 616, 6200 MD Maastricht, The Netherlands ralf.peeters@math.unimaas.nl M. Olivi is with INRIA Sophia-Antipolis, BP 93, 06902 Sophia-Antipolis Cedex, France Martine.Olivi@sophia.inria.fr B. Hanzon is with the School of Mathematical Sciences, University College Cork, Cork, Ireland b.hanzon@ucc.ie unit circle.

) Proposition 1 :

 1 Let Θw,u,v = H(uv * ) I m + (b w (z) -1)uu * 0 0 I m H( wuv * ).

( 9 )

 9 A balanced realization ( Ã, B, C, D) of G = T Θw,u,v (G) can be computed from a balanced realization ( Ã, B, C, D) of G by the formula

Theorem 3 . 2 :

 32 Let G be a m × m rational matrix function which is proper and let (A, B, C, D) be a minimal statespace realization of G. Let Θ = Θ 11 Θ 12 Θ 21 Θ 22 be a 2m × 2m rational matrix function of McMillan degree d (block-partitioned into blocks of size m) and such that Θ(1/z) is proper. Let (A, B, C, D) = A, B 1 B 2 , C 1 C 2 , D 11 D 12 D 21 D 22 be a block partitioned minimal state-space realization of Θ(1/z). Then a state-space realization ( Â, B, Ĉ, D) for F = T Θ (G) is provided by

  22) D. Parametrization issue. Theorem 3.1 is the basis for a Schur type recursive algorithm and for a parametrization of lossless functions of degree n in terms of a sequence of boundary interpolation data. Let σ = (ξ k , u k ) k=1,...,n , |ξ k | = 1, u k = 1, be a sequence of boundary interpolation data. Given a lossless function F (z) of McMillan degree n, the algorithm handles the sequence of boundary interpolation conditions as follows: F n = F and for k = n, n -1, . . . 2 let

Theorem 4 . 1 :

 41 Let (A, B, C, D) be a balanced realization of a lossless function G. A balanced realization ( Â, B, Ĉ, D) of Ĝ = T Φξ,u,v,ρ (G) is given by

) Proposition 4 :F ( 1

 41 Any 2×2 lossless function H(z) satisfying the three vanishing moment conditions H(1) = I 2 , (43) and (45) can be represented as H = T Φ1 (T Φ2 (F )K 1 ) where Φ i is the J-inner matrix (19) with ξ = 1, ρ = ρ i and u = v = 1 0T , K 1 is given by (44) and F is a 2 × 2 lossless function satisfying