Hani Abdeen

Stéphane Ducasse

Houari Sahraoui
email: sahraoui@iro.umontreal.ca

Ilham Alloui
email: ilham.alloui@univ-savoie.fr

Automatic Package Coupling and Cycle Minimization

Keywords: Reverse engineering, Re-engineering, Software modularization, Search algorithms

Object-oriented (OO) software is usually organized into subsystems using the concepts of package or module. Such modular structure helps applications to evolve when facing new requirements. However, studies show that as software evolves to meet requirements and environment changes, modularization quality degrades. To help maintainers improve the quality of software modularization we have designed and implemented a heuristic search-based approach for automatically optimizing inter-package connectivity (i.e., dependencies). In this paper, we present our approach and its underlying techniques and algorithm. We show through a case study how it enables maintainers to optimize OO package structure of source code. Our optimization approach is based on Simulated Annealing technique.

I. INTRODUCTION

In object-oriented languages such as Java, Smalltalk and C++, package structure allows people to organize their programs into subsystems. A well modularized system enables its evolution by supporting the replacement of its parts without impacting the complete system. A good organization of classes into identifiable and collaborating subsystems eases the understanding, maintenance, test and evolution of software systems [START_REF] Deremer | Programming in the large versus programming in the small[END_REF].

However code decays: as software evolves over time with the modification, addition and removal of new classes and dependencies, the modularization gradually drifts and looses quality [START_REF] Eick | Does code decay? assessing the evidence from change management data[END_REF]. A consequence is that some classes may not be placed in suitable packages [START_REF] Griswold | Automated assistance for program restructuring[END_REF]. To improve the quality of software modularization, optimizing the package structure and connectivity is required.

Software modularization is a graph partitioning problem [START_REF] Mitchell | On the evaluation of the bunch search-based software modularization algorithm[END_REF], [START_REF] Mitchell | Using interconnection style rules to infer software architecture relations[END_REF]. Since this last is known as a NP-hard problem [START_REF] Farrugia | Vertex-partitioning into fixed additive inducedhereditary properties is np-hard[END_REF], searching for good modularization by using deterministic procedures or exhaustive exploration of the search space is not feasible without additional heuristics [START_REF] Chapman | System design is an np-complete problem: Correspondence[END_REF], [START_REF] Mitchell | On the evaluation of the bunch search-based software modularization algorithm[END_REF]. We chose then an alternative approach based on heuristic search procedures to identify a good solution within a reasonable amount of computing time [START_REF] Ferland | Heuristic search methods for combinatorial programming problems[END_REF].

Heuristic search methods have already been successfully applied to the software modularization problem [START_REF] Clarke | Reformulating software engineering as a search problem[END_REF], [START_REF] Harman | The current state and future of search based software engineering[END_REF],

for example, to automatically decompose and modularize software systems. They are mainly based on clustering [START_REF] Abreu | Coupling and cohesion as modularization drivers: are we being over-persuaded?[END_REF], [START_REF] Bauer | Architecture-aware adaptive clustering of OO systems[END_REF], [START_REF] Lung | Program restructuring using clustering techniques[END_REF], [START_REF] Mancoridis | Bunch: A clustering tool for the recovery and maintenance of software system structures[END_REF], [START_REF] Mancoridis | Using automatic clustering to produce high-level system organizations of source code[END_REF], [START_REF] Serban | Restructuring software systems using clustering[END_REF], [START_REF] Mitchell | Using interconnection style rules to infer software architecture relations[END_REF] and evolutionary algorithms [START_REF] Doval | Automatic clustering of software systems using a genetic algorithm[END_REF], [START_REF] Harman | A new representation and crossover operator for search-based optimization of software modularization[END_REF], [START_REF] Liu | Using evolutionary algorithms to tackle large scale grouping problems[END_REF], [START_REF] Lutz | Evolving good hierarchical decompositions of complex systems[END_REF], [START_REF] Maini | Genetic algorithms for graph partitioning and incremental graph partitioning[END_REF], [START_REF] Mitchell | On the automatic modularization of software systems using the bunch tool[END_REF], [START_REF] Mitchell | On the evaluation of the bunch search-based software modularization algorithm[END_REF], [START_REF] Mitchell | Using interconnection style rules to infer software architecture relations[END_REF], [START_REF] Seng | Search-based improvement of subsystem decompositions[END_REF].

Few of these works address the problem of optimizing software modularization [START_REF] Harman | A new representation and crossover operator for search-based optimization of software modularization[END_REF], [START_REF] Mitchell | On the evaluation of the bunch search-based software modularization algorithm[END_REF], [START_REF] Mitchell | Using interconnection style rules to infer software architecture relations[END_REF], [START_REF] Seng | Search-based improvement of subsystem decompositions[END_REF], [START_REF] Tzerpos | The orphan adoption problem in architecture maintenance[END_REF]. Often existing approaches change (to various degrees) the existing package structure of an application. In such a case, it can be difficult for a software engineering to understand the resulting structure and to map it back to the situation he knows. Our approach is to support remodularisation of existing package structure by explicitly taking it into account and avoiding creating new package or related abstractions.

In this paper, we present an approach for automatically optimizing existing software modularizations by minimizing connectivity among packages, in particular cyclic-connectivity. The objective of the optimization process is inspired by well known package cohesion and coupling principles already discussed in [START_REF] Briand | A unified framework for cohesion measurement in object-oriented systems[END_REF], [START_REF] Fowler | Reducing coupling[END_REF], [START_REF] Martin | Agile Software Development. Principles, Patterns, and Practices[END_REF]. We limit ourselves to direct cyclic-connectivity and restrict our optimization actions to moving classes over existing packages. By definition, there is a direct cyclic-connectivity between two packages if they mutually depend on each other.

Our approach is based on Simulated Annealing [START_REF] Ferland | Heuristic search methods for combinatorial programming problems[END_REF], [START_REF] Kirkpatrick | Optimization by simulated annealing[END_REF], which is a neighborhood (local) search-based technique. Simulated Annealing is inspired by the annealing process in metallurgy [START_REF] Kirkpatrick | Optimization by simulated annealing[END_REF]. We chose this technique because, it suits well our problem, i.e., local optimization of an existing solution. Moreover, it has been shown to perform well in the context of automated OO class design improvement [START_REF] O'keeffe | Search-based software maintenance[END_REF], [START_REF] O'keeffe | Search-based refactoring for software maintenance[END_REF] and more generally, in the context of software clustering problems [START_REF] Mitchell | Using heuristic search techniques to extract design abstractions from source code[END_REF], [START_REF] Mitchell | On the evaluation of the bunch search-based software modularization algorithm[END_REF].

Contribution: we present an approach, using simulated annealing technique, for automatically reducing package coupling and cycles by only moving classes over packages while taking into account the existing class organization and package structure. In our approach, maintainers can define [START_REF] Abreu | Coupling and cohesion as modularization drivers: are we being over-persuaded?[END_REF] the maximal number of classes that can change their packages, [START_REF] Bauer | Architecture-aware adaptive clustering of OO systems[END_REF] the maximal number of classes that a package can contain, and (3) the classes that should not change their packages or/and the packages that should not be changed.

Section II defines the terminology we use and gives an overview about challenges in optimizing package structure. Section III presents the metrics used in evaluation functions of both modularization and package quality. We detail our optimization algorithm in Section IV. We validate our approach using real large applications and discuss the results in Section V. In Section VI we position our approach with related works, before concluding in Section VII.

II. BACKGROUND AND VOCABULARY

In this section, we introduce the terminology used in this paper and give an overview of some challenging issues we address.

A. Terminology

Before all, we define OO software Modularization M as a decomposition of the set of software classes M C into a set of packages M P . Fig. 1 shows two modularizations, both consists of 9 classes distributed over 3 packages: where Package represents an entity that only contains classes, where every class c belongs to only one package c p . We define package size p size by the number of its classes.

Every class c can be related to other classes, in consequence to their packages, through a set of Dependencies (c D ). This set consists of two subsets: Outgoing Dependencies (c Out.D ) and Incoming Dependencies (c Inc.D ). We denote a dependency that goes from a class c i to another one c j by the pair (c i , c j ). A dependency might be either method calls, class access, or class inheritance.

Every dependency is either internal if it is related to two classes belonging to the same package, or external if not. The set of dependencies related to a package p (p D ) represents the union of the set of internal dependencies (p Int.D ) and the set of external dependencies (p Ext.D ) relatively to p C , the set of p classes. The set p Ext.D consists of two subsets: dependencies that are either exiting p (p Ext.Out.D ) or that are pointing to p (p Ext.Inc.D ). The sets p Ext.Out.D and p Ext.Inc.D respectively relate p to its provider packages (p P ro.P ) and its client packages (p Cli.P ). M odularization 1 in Fig. 1 shows that there is one dependency exiting p1 and three dependencies pointing to p1: p1 To determine connectivity at the package level, we say that there is a Connection from a package y to another one z if there is n (n > 0) outgoing dependencies y pointing to z. M odularization 1 in Fig. 1 shows that there are two connections going from p2 to its provider packages: p2 Out.Con = [(p2, p1), (p2, p3)]; and there is one connection going to p2: p2 Inc.Con = [(p1, p2)].

Dependencies can form direct cyclic dependencies/connections between packages. If package y plays the role of client and provider for package z then we consider involved dependencies and connections between y and z as cyclic. We denote the set of cyclic dependencies related to a package p by p Cyc.D . The set p Cyc.D consists also of two subsets: cyclic-dependencies that are either exiting p (p Out.Cyc.D ) or pointing to p (p Inc.Cyc.D ). p Out.Cyc.D and p Inc.Cyc.D are dependencies causing cycles between packages (and not classes) in the context of the client-provider relation. Similarly we denote the set of cyclic connections related to p by p Cyc.Con : p Cyc.Con = p Out.Cyc.Con ∪ p Inc.Cyc.Con . M odularization 1 in Fig. 1 shows that p1 has one outgoing cyclic-dependency [(c 3 , c 6 )] and two incoming cyclic-dependencies [(c 7 , c 1 ), (c 7 , c 2 )]. Those cyclic-dependencies produce two cyclic-connections: p1 Cyc.Con = p2 Cyc.Con = [(p1, p2), (p2, p1)].

B. Challenges in Optimizing Package Structure

Optimizing package structure in existing large applications is a difficult problem because of the following reasons:

Large applications are usually very complex: they contain thousands of heavily inter-connected classes. Many of the dependencies are between classes belonging to different packages, which increases the inter-package connectivity. In such situation, the optimization problem is more difficult.

Classes usually are not well distributed over packages: in real applications, some packages contain large sets of classes and other packages contain few number of classes. As consequence, most application packages depend on those large packages. Furthermore, it is difficult to determine an ideal package size since it may depend on external factors such as the team structure, domain, or coding practice...

Optimizing some criteria may degrade others: minimizing inter-packages dependencies/connections may increase the number of noncyclic ones. e.g., Fig. 1 shows two modularizations, M odularization 1 and M odularization 2 . Both are composed of the same set of classes and packages. The difference between them is c 7 that is moved from p2 to p3 in M odularization 2 . In M odularization 2 there are 2 cyclic-dependencies [(c 7 , c 6 ), (c 6 , c 9 )] compared to 3 in M odularization 1 [(c 3 , c 6 ), (c 7 , c 1 ), (c 7 , c 2 )]. Thus moving c 7 has reduced the number of cyclic-dependencies. On the other hand, moving c 7 increases the number of inter-package dependencies. In M odularization 2 , there are 6 interpackage dependencies compared to 5 for M odularization 1 .

III. MODULARIZATION QUALITY

Our goal is to automatically optimize the decomposition of software system into packages so that the resulting organization of classes/packages, mainly, reduces connectivity and cyclic-connectivity between packages. This goal is inspired from well known quality principles already discussed in [START_REF] Briand | A unified framework for cohesion measurement in object-oriented systems[END_REF], [START_REF] Fowler | Reducing coupling[END_REF], [START_REF] Martin | Agile Software Development. Principles, Patterns, and Practices[END_REF] and in particular from the following principle: packages are desired to be loosely coupled and cohesive to a certain extent [START_REF] Fowler | Reducing coupling[END_REF]. In such a context, we need to define measures that evaluate package cohesion and coupling.

In addition, cyclic dependencies between packages are considered as an anti-pattern for package design [START_REF] Martin | Agile Software Development. Principles, Patterns, and Practices[END_REF]. In this section we define two suites of measures: the first is used when evaluating modularization quality and the second is used when evaluating modularity quality of single package. Note that all measures we define in this section take their value in the interval [0..1] where 1 is the optimal value and 0 is the worst value.

A. Measuring Modularization Quality

Inter-Package Dependencies: according to Common Closure Principle (CCP) [START_REF] Martin | Agile Software Development. Principles, Patterns, and Practices[END_REF], classes that change together should be grouped together. In such a context, optimizing modularization requires reducing the sum of inter-package dependencies (IP D =

CCQ(M) = 1 -IP D ICD
Inter-Package Connections: according to Common Reuse Principle (CRP) [START_REF] Martin | Agile Software Development. Principles, Patterns, and Practices[END_REF], classes that are reused together should be grouped together. In such a context, optimizing modularization requires reducing the sum of inter-package connections (IP C = [START_REF] Fowler | Reducing coupling[END_REF]. We define the measure CRQ to evaluate the Common Reuse Quality of a modularization M as follows:

|M P | i=1 |p i Out.Con |) [3],
CRQ(M) = 1 -IP C ICD
Inter-Package Cyclic-Dependencies: according to Acyclic Dependencies Principle (ADP) [START_REF] Martin | Agile Software Development. Principles, Patterns, and Practices[END_REF], dependencies between packages must not form cycles. In such a context, optimizing modularization requires reducing the sum of inter-package cyclic-dependencies (IP CD =

|M P | i=1 |p i Out.Cyc.D |).
We define the measure ADQ to measure the Acyclic Dependencies Quality of a modularization M as follows:

ADQ(M) = 1 -IP CD

ICD

Inter-Packages Cyclic-Connections: as for cyclic dependencies between packages, reducing cyclic connections between packages is required. For example, in M odularization 1 in Fig. 1, there are 3 cyclic dependencies [(c 3 , c 6 ), (c 7 , c 1 ), (c 7 , c 2 )] and 2 cyclic connections [(p1, p2), (p2, p1)]; moving c 7 to p 3 will reduce the number of cyclic-dependencies: in modularization 2 there are only 2 cyclic dependencies [(c 6 , c 9 ), (c 7 , c 6 )], but it remains 2 cyclic connections [(p2, p3), (p3, p2)]. We thus deduce that reducing inter-package cyclic dependencies does not necessarily reduce inter-package direct cyclic-connections (IP CC =

|M P | i=1 |p i Out.Cyc.Con |).
We define the measure ACQ to evaluate the Acyclic Connections Quality of a modularization M as follows:

ACQ(M) = 1 -IP CC ICD

B. Measuring Package Quality

In addition to measures presented in Section III-A, we define a set of measures that help us determine and quantify the quality of a single package within a given modularization. To normalize the value of those measures we use the number of dependencies related to the considered package

(|p D |) with |p D | > 0.
Package Cohesion: we relate package cohesion to the direct dependencies between its classes. In such a context, we consider that the cohesion of a package p is proportional to the number of internal dependencies within p (|p Int.D |). This is done according to the Common Closure Principle (CCP) [START_REF] Martin | Agile Software Development. Principles, Patterns, and Practices[END_REF]. We define the measure of package cohesion quality similarly to that in [START_REF] Abreu | Coupling and cohesion as modularization drivers: are we being over-persuaded?[END_REF] as follows:

CohesionQ(p) = |p Int.D | |p D |
Package Coupling: we relate package coupling to its efferent and afferent coupling (Ce, Ca) as defined by Martin in [START_REF] Martin | The tipping point: Stability and instability in OO design[END_REF]. Package Ce is the number of packages that this package depends upon (|p P ro.P |). Package Ca is the number of packages that depend upon this package (|p Cli.P |). According to the common reuse principle, we define the measure of package coupling quality using the number of package providers and clients as follows:

CouplingQ(p) = 1 -|p P ro.P ∪p Cli.P | |p D |
Package Cyclic-Dependencies: for automatically detecting packages that suffer from direct-cyclic dependencies we define a simple measure that evaluates the quality of package cyclic dependencies (CyclicDQ) using the number of package cyclic dependencies:

CyclicDQ(p) = 1 - |p Cyc.D | |p D |
Similarly we define another measure that evaluates package cyclic connections quality (CyclicDQ) using the number of package cyclic connections:

CyclicCQ(p) = 1 - |p Cyc.Con | |p D |

IV. OPTIMIZATION TECHNIQUE (METHODOLOGY)

To optimize package connectivity, we use an optimization procedure that starts with a given modularization and gradually modifies it, using small perturbations. At each step, the resulting modularization is evaluated to be possibly selected as an alternative modularization. The evaluation of modularization quality is based on measures defined in Section III-A. This section describes our optimization approach and algorithm.

A. Technique Overview

To address the problem of optimizing modularization, we use heuristic optimization technique based on simulated annealing algorithm [START_REF] Ferland | Heuristic search methods for combinatorial programming problems[END_REF], [START_REF] Kirkpatrick | Optimization by simulated annealing[END_REF]. Simulated annealing is an iterative procedure that belongs to the category of Neighborhood Search Techniques (N ST ). Algorithm 1 shows an overview of the optimization algorithm. The optimization process performs series of local searches with the global search parameter T current . T current represents in simulated annealing technique the current temperature of the annealing procedure which started with the value T start . A local search consists of num (num ≥ 1) searches of suboptimal solution. At each of them, a new modularization M trial is derived from a current one M current by applying to this later a modification. The derivation of M trial from M current is performed by the neighborhood function. Then, the algorithm evaluates the M trial and M current fitness using the fitness function F, where the bigger is the value of F(M), the better is modularization M: if M trial is better than M current then M trial becomes the M current ; then, if M current is better than the current best modularization M best , M current becomes the M best . At the end of each local search, the parameter T current decreases and another local search starts with the new value of T current . Decreasing T current is the responsibility of CoolingSchedule function. This latter is defined according to Keeffe etal. discussion [START_REF] O'keeffe | Search-based refactoring for software maintenance[END_REF] using a geometric cooling scheme: CoolingSchedule(T ) = 0.9975 * T . Local searches are repeated until reaching T stop (T current ≤ T stop ).

Algorithm 1 Optimization Algorithm

To circumvent the problem of local optima [START_REF] Ferland | Heuristic search methods for combinatorial programming problems[END_REF], a lessgood modularization can be accepted with some probability: a less-good modularization M trial can replace M current under some conditions AcceptanceCondition. Simulated annealing technique defines acceptance conditions in a way that the probability of accepting a less-good modularization decreases over time. We define AcceptanceCondition as follows: r > e -

B. Evaluating Modularization Quality (Fitness)

As for any search-based optimization problem, the definition of the fitness function represents a central concern as it guides the search. We define our fitness function as a combination of the measures defined in Section III-A. We define dependency quality (DQ) for a modularization M as the weighted average of Common Closure Quality (CCQ) and Acyclic Dependencies Quality (ADQ); and we define connection quality (CQ) for M as the weighted average of Common Reuse Quality (CRQ) and Acyclic Connections Quality (ACQ). To give higher intention to cyclic dependencies/connections between packages we define a factor of importance γ (γ = β α , β > α ≥ 1):

DQ(M) = α * CCQ(M)+β * ADQ(M) α+β CQ(M) = α * CRQ(M)+β * ACQ(M) α+β
Both functions DQ and CQ take their values in the interval [0..1] where 1 is the optimal value. The final fitness function is defined by the average of DQ and CQ:

F (M) = DQ(M)+CQ(M) 2
Our hypothesis is: optimizing F will reduce inter-package dependencies and connections, particularly cyclic ones. Furthermore, in addition to AcceptanceCondition for lessgood modularizations we defined in Section IV-A, the optimization process may accept a less-good resulting modularization only if the number of inter-package dependencies decreases (i.e., DQ increases). We expect such a decision facilitates the reduction of inter-package cyclic dependencies.

C. Modularization Constraints

In addition to the fitness function, our approach allows maintainers to define distinct constraints that should complete the evaluation process and guarantee maintainers' requirements. The rationale behind those constraints is to control the optimization process when optimizing a given modularization. e.g., putting a major partition of classes into one package can effectively reduce inter-package (cyclic-) dependencies/connections; such an approach is clearly not the best one to optimize software modularization. This section presents three constraints allowing control the optimization process. Section IV-D explains how the optimization process favors these constraints when deriving new modularizations.

1) Controlling package size: to avoid having very large and dominant packages, we introduce the following constraint: the size of every package (p size ) should always be smaller than a predefined number (size max ). We define size max for every package p relatively to its size in the original modularization p size V 0 : size max = δ + p size V 0 , δ ≥ 0. Maintainers can define δ according to the context of the concerned software system. We cannot determine upfront the good interval in which δ should be taken. In the scope of this paper, we define δ as the theoretical package size in M 0 which equals to the ratio:

|M0 C | |M0 P | .
It is worth to note that maintainers can define a different δ for each package: e.g., for a large package p, δ may be defined to 0; this way, p will never be larger than before.

2) Controlling modularization modification: maintainers should be able to define the limit of modifications that the optimization process can apply on the original modularization M 0 when it proceeds. In other words, the optimization process must take into account the maximal authorized distance (distance max ) between resulting modularizations and M 0 . In our context, for two modularizations that entail the same set of classes, we define the distance between them by the number of classes that changed packages. This way, distance max can be defined simply as the maximal number of classes that can change their packages.

3) Controlling modularization structure: moreover, we found that it is very helpful to allow maintainers decide whether some classes should not change their package and/or whether given packages should not be changed. We say that such classes/packages are frozen. This constraint is particularly helpful when maintainers know that a given package is well designed and should not be changed: e.g., if a small package p contains a couple of classes that extend classes from other package, p may be considered a well designed package, even if it is not cohesive.

D. Deriving New Modularization (Neighbor)

The neighborhood function (N ) is the second main concern of the optimization process. Defining N requires: (1) the definition of the set of modifications that N can use to derive new modularizations, (2) and the definition of a process that derives a new modularization from another one. This section presents our definition of N .

Since we search near optimal modularization by applying near minimal modification to the original modularization, we limit the set of modifications that N can use to only: moving a class c from its current package p source to another one p target . In this context, we say that c is the modification actor (c actor ). To minimize search-space we reduce the selectionspace of p target to the set of client and provider packages of c actor : p target ∈ (c actor P ro.P ∪ c actor Cli.P ).

We specify the derivation of a neighbor modularization of a modularization M by 4 sequential steps: (1) selecting p source , (2) selecting c actor , (3) selecting p target and then (4) moving c actor to p target . Selections in the first three steps are done arbitrary using a probability function. The probability function gives higher probability to the worst package into M P to be selected as p source , to the worst class into p source to be selected as c actor and to the nearest package to c actor to be selected as p target . The selection mechanism performs similarly to a roulette wheel selection, where each class/package is given a slice proportional to its probability to be selected and then we randomly take a position in the roulette and pick the corresponding class/package.

It is worth to note that packages and classes that are defined as frozen (Section IV-C3), do not belong to the selection spaces: a frozen package will never be a p source or p target , and a frozen class will never be a c actor .

The following subsections explain our definition of the probability of being selected as p source , c actor or p target . Note that in our definition of this probability we use the factor γ, as defined in the fitness function (Section IV-B), to pay more intention to cyclic dependencies/connections. 1) Selecting p source : the worst package in M P is, the highest probability to be selected it has. We relate package badness to the quality of its cohesion, coupling and of its external dependencies (i.e., the density of cyclic dependencies/connections related to the concerned package). We define the badness of package by using the measures: CohesionQ, CouplingQ, CyclicDQ and CyclicCQ (Section III-B): where we relate CohesionQ and CyclicDQ to package dependency quality (DQ). Also we relate CouplingQ and CyclicCQ to package connection quality (CQ). We define package quality functions, similarly to modularization quality functions defined in Section IV-B :

DQ(p) = α * CohesionQ(p)+β * CyclicDQ(p) α+β CQ(p) = α * CouplingQ(p)+β * CyclicCQ(p) α+β
Both functions DQ and CQ take their values in the interval [0..1] where 1 is the optimal value. Finally, we define package badness basing on the average of DQ and CQ:

Badness(p) = 1 -DQ(p)+CQ(p)
In addition to satisfy constraints discussed in Section IV-C, we define the probability of selecting a package p as p source by: ρ * Badness(p), where ρ is a factor takes its value in the interval [0..1]. It is the average of two sub-factors (ρ 1 , ρ 2 ):

• ρ 1 is based on p size : relatively to p size in the original modularization M 0 (p 0size ), a package whose size increased has a higher probability to be selected than a package whose size decreased. By doing so, we expect that the package size in resulting modularizations will be similar to that in the original modularization; • ρ 2 is based on the number of new classes into p: relatively to p 0 , a package that acquired the largest number of new classes (i.e., classes are not packaged in p 0 ) has the highest probability to be selected. By doing so, we favor moving classes that already changed their packages until they find their optimal package. 2) Selecting c actor : the worse a class in p source is, the highest probability to be selected it has. We relate class badness to the number of external dependencies related to the class (|c Ext.D |) and to the number of its external cyclicdependencies (related to its package p):

Badness(c) = α * |c Ext.D |+β * |c Ext.D ∩p Cyc.D | α+β
In addition, to satisfy the constraint of distance max (Section IV-C2), when the distance (d) between resulting modularizations and the original one increases, classes that have already changed their packages have higher probability to be c actor . In this context, we use the factor ρ = 1 -d distancemax . If ρ ≤ 0 then only classes which already changed their original packages can move. Only if 0 < ρ ≤ 1 then the optimization process can move more classes over packages but with a probability ρ. Thus we define the probability of selecting a class c as c actor as following:

0 ρ ≤ 0, not(isM oved(c)) ρ * Badness(c) ρ > 0, not(isM oved(c)) Badness(c) isM oved(c)
Where the predicate isM oved(c) is true if c is already moved from its original package.

3) Selecting p target : the nearest package to c actor is, the highest probability to be selected it has. We simply relate the nearness of a package p To satisfy the constraint on p size (size max defined in Section IV-C1), when package size increases its probability to be selected as p targert decreases. In this context, we use the factor ρ = 1-psize sizemax . If ρ ≤ 0 then the package size should not increase anymore. Only if 0 < ρ ≤ 1 then the package size can increase but with a probability ρ which decreases when p size increases. Thus we define the probability of selecting a package p as p target for a class c, as following:

0 ρ ≤ 0 ρ * N earness(p, c) 0 < ρ ≤ 1

V. EXPERIMENTS AND VALIDATION

To validate our optimization approach, we applied it to several software applications that differ in terms of: number of classes (|M C |), number of packages (|M P |), number of inter-class dependencies (ICD); number of interpackage dependencies (IP D), connections (IP C), cyclic dependencies (IP CD) and cyclic connections (IP CC). Table I shows information about the original modularization of those software applications.

Since the search process is not deterministic, we applied our algorithm 10 times for each software application and we calculated the average of modularization parameters cited in Table I. We used the parameters T start , T stop and num (Algorithm 1) with value 50, 1 and 30 respectively. On another hand, we weighted cyclic dependencies/connections to be three times more important than noncyclic dependencies/connections. We performed our experience twice: the first time, we did not used the constraint distance max . In the second time, we limited distance max to 5%, which means that only 5% of classes can change their original packages.

General optimization: Table III shows optimization results. In the resulting modularization for JEdit (JEdit 1 ), 10.2% of inter-package dependencies IP D, 23.3% of inter-package cyclic-dependencies IP CD, 24% of interpackage connections IP C and 37.2% of inter-package cyclicconnections IP CC have been removed. This significant improvement of inter-package connectivity was obtained by moving only 8.9% of the classes (d = 8.9%). Similarly for other case studies, the optimization process has improved original modularizations by moving a relatively small number of their classes. When limiting distance max to 5%, the algorithm obtained similar results.

Class distribution and package size: Table IV shows that some packages were empty in resulting modularizations -since their classes moved to other packages. For example, in ArgoU M L 1 , 25.4% of packages were empty. By inspecting packages in the original modularization we found that those empty packages are packages which have originally very small sizes (i.e., in average two or three classes) and have low quality for cohesion, coupling and/or cyclic dependencies. This conclusion were also true for the other case studies.

Fig. 2 shows an overview about package size and cohesion for the original modularization of ArgoU M L and for the resulting modularizations (Table III). We can see that empty packages in the resulting modularizations are packages whose sizes are small and whose cohesion is relatively worse. On another hand, Fig. 2 shows also that the size of some small packages, annotated by small packages, is increased in the resulting modularization. Dominant packages, annotated by dominant packages, is a main cause of bad distribution of classes: moving classes from dominant packages to small ones generally produces more dependencies and connections among packages. Maintainers can avoid moving classes from small packages to dominant ones by limiting the size max (Section IV-C1) of dominant packages to their original size: p sizemax = p size V 0 . This way, the dominant package size will never increase and the optimization process will search better modularizations by moving classes among/to smaller packages.

Fortunately, in the case of JEdit 1 , only 10.5% (2/19) of packages are empty (Table IV), where Table III shows that our optimization process has effectively optimize package connectivity in JEdit 1 . Now as a future work, we have to perform a deep manual validation since in presence of late-binding and frameworks, some small packages may extend larger ones and as such may have a real reason to exist. Note that defining such packages as f orzen (Section IV-C3) will keep those packages existing.

While some packages became empty, Table IV shows that the average package size ( |M C | |M P | ) for the resulting modularizations is really close to the average package size for the original ones Table I. Similarly, we can see that for the maximum package size (maxP size ). This shows that the optimization algorithm conserves the original system shape.

Package quality optimization: Table V shows that package quality average is also optimized: cohesion quality average (CohesionQ Avg ), coupling quality average (CouplingQ Avg ) and cyclic-dependency quality average (CyclicsDQ Avg ) for resulting packages are also almost all optimized, even if distance max is limited to only 5%. This can be seen also in Fig. 2.

In only one case (JEdit 1 and JEdit 2 ), the package coupling quality (CouplingQ) decreased with a very good improvement of CyclicsDQ and CohesionQ. We explain this by the fact that the optimization process gives more importance to inter-packages cyclic-dependencies. Indeed, CyclicDQ had a very bad value in the original modularization JEdit (Table II): the ratio of inter-package cyclic-dependencies ( IP CD ICD ) shows that 38.4% of inter-class dependencies form cyclic-dependencies between packages. Moreover there are 802 classes distributed over only 19 packages, so that the search space for generating new modularization is limited. Consistency of resulting modularizations: since our optimization approach uses random selection, different executions produce different modularizations. To evaluate the consistency of our optimization approach, we have applied it 10 times on each case study (Table I). As a result, each application has 10 modularizations [M 1 ..M 10 ]. Table VI shows the average distance between every pair (M i , M j ). For example, between resulting modularizations for JEdit, there are, in average, only 3% of classes that have not the same packages. For Jboss, only 5.6% of the classes have different packages in distinct resulting modularizations.

We mainly relate this very good consistency of resulting modularizations to the improvements we introduced to the neighbor function N Section IV-D (i.e., the probability function to being selected).

In conclusion, the obtained results are very convincing. For all the case studies, the new modularizations are clearly better than the original ones. Moreover, our optimization process produces very similar results.

VI. RELATED WORKS

Our work is mostly related to work on software modularization and decomposition [START_REF] Abreu | Coupling and cohesion as modularization drivers: are we being over-persuaded?[END_REF], [START_REF] Harman | A new representation and crossover operator for search-based optimization of software modularization[END_REF], [START_REF] Harman | Pareto optimal search based refactoring at the design level[END_REF], [START_REF] Mancoridis | Using automatic clustering to produce high-level system organizations of source code[END_REF], [START_REF] Mancoridis | Bunch: A clustering tool for the recovery and maintenance of software system structures[END_REF], [START_REF] Mitchell | Using heuristic search techniques to extract design abstractions from source code[END_REF], [START_REF] Mitchell | On the automatic modularization of software systems using the bunch tool[END_REF], [START_REF] Mitchell | On the evaluation of the bunch search-based software modularization algorithm[END_REF], [START_REF] Seng | Search-based improvement of subsystem decompositions[END_REF].

Mitchell etal. [START_REF] Mancoridis | Using automatic clustering to produce high-level system organizations of source code[END_REF], [START_REF] Mancoridis | Bunch: A clustering tool for the recovery and maintenance of software system structures[END_REF], [START_REF] Mitchell | Using heuristic search techniques to extract design abstractions from source code[END_REF] introduced a search-based approach based on hill-climbing clustering technique to cluster software modules (classes in our context). Their approach starts with an initial population of random modularizations. The clustering algorithm clusters each of the random modularization and selects the result with the largest quality as the suboptimal solution. Recently, they used Simulated Annealing technique to optimize resulting clusters [START_REF] Mitchell | Using heuristic search techniques to extract design abstractions from source code[END_REF], [START_REF] Mitchell | On the automatic modularization of software systems using the bunch tool[END_REF], [START_REF] Mitchell | On the evaluation of the bunch search-based software modularization algorithm[END_REF]. Their optimization approach creates new modularizations by moving randomly some classes (a block of classes) to new clusters. The goal of their approach is increasing cluster internal dependencies.

Harman etal. [START_REF] Harman | Pareto optimal search based refactoring at the design level[END_REF] introduces a non-exhaustive hill climbing approach to optimize and determine a sequence of class refactorings. Similarly to our approach, they also restricted their approach to only move methods (classes in our context) over existing classes (packages in our context). The goal of their approach is reducing the class coupling, basing on the Coupling Between Objects (CBO) metric [START_REF] Briand | A unified framework for cohesion measurement in object-oriented systems[END_REF]. To avoid having a very large classes, they also used the dispersion of methods over classes (the standard deviation of methods per class metric) as a factor to measure the quality of resulting class refactoring sequences.

Abreu etal. [START_REF] Abreu | Coupling and cohesion as modularization drivers: are we being over-persuaded?[END_REF] used hierarchical agglomerative clustering methods to decompose software classes into packages. Their clustering methods starts with a set of classes considering that each class is placed within a singleton cluster.

The goal of their approach is also increasing package internal dependencies (i.e., package cohesion). In addition to the package cohesion, they used the dispersion of classes over packages (i.e., package size dispersion) as a factor to measure the modularization quality.

Seng etal. [START_REF] Seng | Search-based improvement of subsystem decompositions[END_REF] and Harman etal. [START_REF] Harman | A new representation and crossover operator for search-based optimization of software modularization[END_REF] proposed genetic algorithms to partition software classes into subsystems (packages). Their algorithms start with an initial population of modularizations. These algorithms apply genetic operators on packages to modify current modularizations and/or create new modularizations into the population. The goal of both works is increasing package internal dependencies. Seng etal. consider also cyclic-dependencies between packages as anti-pattern for package design quality.

Our approach has several advantages compared to those works.

Considering explicitly the original modularization and controlling the optimization process: our approach tackles the problem of optimizing existing software modularizations rather than the problem of software re-modularization. Indeed, our optimization approach starts from one original modularization instead of an initial population of modularizations or a flat set of classes. Although we use an optimization technique similar to that in Mitchell etal. work, we restrict ours to moving classes over existing packages rather than creating new packages since we want to minimize the distance for a maintainer between the initial situation and the resulting one. Even if some prior works support the notion of importing a defined clustering [START_REF] Tzerpos | The orphan adoption problem in architecture maintenance[END_REF] and restrict modifications to only moving classes over existing packages [START_REF] Harman | Pareto optimal search based refactoring at the design level[END_REF], we did not found, in the software re-modularization literature, approaches that explicitly take into account the original modularization structure as we do. Our approach allows maintainers to specify a set of constraints (e.g., the number of classes that may change their package, the package maximal size and the packages/classes which should not be changed/moved) and therefore make them control the optimization process. Although those constraints are simple, they are very important and helpful for the automatic optimization of software modularization.

Searching better modularizations by doing near-minimal modifications: differently of those cited works, which search good modularizations without considering the distance between resulting modularizations and the original ones, we introduce a probability function that improves the derivation of neighbor modularizations by taking into account the distance between resulting modularizations and the original one, in addition to other constraints (Section IV-C) and package quality parameters (Section III-B). The great advantage of the probability function is finding better modularization by doing near-minimal modifications (Section V, the distance between resulting modularizations and the original one is very small).

Reducing package coupling and cycles in two levels (inter-package dependencies and connections): another advantage of our approach is that we use an evaluation function consisting of a combination of multiple measures. This allows us to have a much richer quality model than the approaches cited above which are mostly based on the unique goal of maximizing package internal dependencies. Although those approaches aim at reducing the global number of interpackage dependencies (e.g., the fitness measure MQ used by Mitchell etal. [START_REF] Mitchell | Using heuristic search techniques to extract design abstractions from source code[END_REF], [START_REF] Mitchell | On the automatic modularization of software systems using the bunch tool[END_REF], [START_REF] Mitchell | On the evaluation of the bunch search-based software modularization algorithm[END_REF]), they do not take into account the number of coupled packages (i.e., the number of inter-package connections). As consequence, they do not check whether package coupling is reduced or not along the optimization process. In addition, those cited prior works do not consider package cycles. Excepted Seng etal. [START_REF] Seng | Search-based improvement of subsystem decompositions[END_REF], they consider inter-package cyclic-dependencies, without taking in account inter-package cyclic-connections.

VII. CONCLUSION AND FUTURE WORK

In this paper, we addressed the problem of optimizing existing modularizations by reducing the connectivity, particularly the cyclic-connectivity, among packages. We proposed an optimization algorithm and a set of measures that our optimization process uses to automatically evaluate the quality of a modularization. When designing our optimization approach, we exploited several principles of package design quality to guide and to optimize the automatic derivation of new modularizations from an existing one. We limited the optimization process to only moving classes over existing packages. We also introduced constraints related to package size, to the number of classes that are allowed to change their packages and to the classes/packages that should not be moved/changed. The results obtained from 4 case studies on real large applications shown that our optimization algorithm has been able to reduce, significantly, package coupling and cycles, by moving a relatively small number of classes from their original packages. These results are important because the chosen applications have radically different original modularizations (in terms of number of classes/packages, inter-class/inter-package dependencies, etc.).

As future work, we intend to enhance our approach by: (1) supporting indirect cyclic-dependencies among packages, (2) taking into account visibility of classes and particular cases of classes (e.g., inner classes in Java). Indeed in this paper, we considered that classes are always public and then can change their packages, (3) setting up a real validation supported by proper statistics and qualitative analyses of resulting source code structures. We also envisage a deep comparison with the closest related work (e.g., Bunch).

Figure 1 .

 1 Figure 1. Example of two modularizations: different decompositions of the set of classes [c1..c9] into 3 packages [p1, p2, p3].

  Ext.Out.D = [(c 3 , c 6 )]; p1 Ext.Inc.D = [(c 7 , c 1 ), (c 7 , c 2 ), (c 8 , c 2 )]. It also shows that p1 has one provider package and two client packages: p1 P ro.P = [p2]; p1 Cli.P = [p2, p3].

  |M P | i=1 |p i Ext.Out.D |) [3], [11]. Since we do not change the dependencies between classes during our optimization process, we use the sum of inter-class dependencies (ICD = |M C | j=1 |c j Out.D |) as normalizer. We define the measure CCQ to evaluate the Common Closure Quality of a modularization M as follows:

  T current T start , r ∈ [0..1]. The value of r is generated randomly in the interval [0..1]. The function e -T current T start takes its value in the interval [0..1] ∀T current ≥ 0, and T current ≤ T start . It increases along the optimization process -since T current decreases. By doing so, the probability of accepting a less-good modularization decreases over time.

  to a class c to the number of dependencies that c has with p classes (|c D ∩ p D |) and to the number of dependencies related to c and present cyclic dependencies between c p and p (|c D ∩ p Cyc.D |): N earness(p, c) = α * |c D ∩p D |+β * |c D ∩p Cyc.D | α+β

Figure 2 .

 2 Figure 2. Package size and cohesion into ArgoUML original (dark gray) and resulting (light gray) modularizations. Packages have the same order in

  Tstop ≥ 1, Tstart > Tstop, num > 1 and M original Ensure: M best M best ← M original Mcurrent ← M

	M best ← Mcurrent
	end if
	else if AcceptanceCondition then
	-accepting a worse modularization-
	Mcurrent ← M trial
	end if
	end for
	-end of local search-
	Tcurrent ← CoolingSchedule(Tcurrent)
	end while
	-end of global search-
	Return M best .

Require: best

Tcurrent ← Tstart -starting global searchwhile Tcurrent > Tstop do -starting local searchfor i = 1 to num do -generating a new modularization and evaluating it-

M trial ← N eighborhood(Mcurrent) if F(M trial ) > F(Mcurrent) then Mcurrent ← M trial if F (Mcurrent) > F(M best ) then

Table I INFORMATION

 I ABOUT USED SOFTWARE APPLICATIONS.

	Original |MC |	ICD	|MP |	IP D	IP CD IP C IP CC maxPsize	|M C | |M P |
	JEdit	802	2683	19	1430	1032	110	26	173	42.2
	ArgoUML	1671	7432	76	5661	1406	517	63	156	22
	Jboss	3094	8859	455	7219	296	1898	41	80	6.8
	Azureus	4212	13945	380	10929	1319	2037	136	213	11

Table III OPTIMIZATIONS

 III ON INTER-PACKAGE CONNECTIVITY. THE TOP TABLE SHOWS THE PERCENT OF REDUCTION OF IPD, .. , IPCC (TABLE I) INTO RESULTING MODULARIZATIONS. THE BIGGEST NEGATIVE VALUE IS, THE BEST OPTIMIZATION IS. THE BOTTOM TABLE SHOWS THESE INFORMATION WHEN distancemax IS SPECIFIED AND LIMITED TO 5%.

Table V OPTIMIZATIONS

 V ON PACKAGE QUALITY. THE TOP TABLE SHOWS THE AVERAGE OPTIMIZATIONS ON PACKAGE QUALITY INTO RESULTING MODULARIZATIONS. VALUES ARE BASED ON TABLE II. THE BIGGEST POSITIVE VALUE IS, THE BEST OPTIMIZATION IS. THE BOTTOM TABLE SHOWS THESE INFORMATION WHEN distancemax IS SPECIFIED AND LIMITED TO 5%.

	Optimization1 CohesionQAvg CouplingQAvg CyclicsDQAvg
	JEdit1	+06.1%	-00.8%	+10.4%
	ArgoU M L1	+08.1%	+07.4%	+00.4%
	Jboss1	+08.5%	+11.6%	+01.8%
	Azureus1	+05.8%	+03.6%	+05.4%
	Optimization2 CohesionQAvg CouplingQAvg CyclicsDQAvg
	JEdit2	+05.4%	-02.2%	+05.2%
	ArgoU M L2	+06.0%	+07.4%	+00.1%
	Jboss2	+04.0%	+09.4%	+00.5%
	Azureus2	+03.9%	+04.6%	+02.2%

  TABLE SHOWS THE AVERAGE DISTANCE BETWEEN TEN RESULTING MODULARIZATIONS FOR EACH APPLICATION.

	Optimization1 DistanceAvg
	JEdit1	3%
	ArgoU M L1	4.1%
	Jboss1	5.6%
	Azureus1	4.9%

ACKNOWLEDGMENT

We gratefully acknowledge the financial support of the french ANR (National Research Agency) for the project "COOK: Réarchitecturisation des applications industrielles objets" (JC05 42872).

Optimization1

IP D IP CD IP C IP CC JEdit1(d = 8.9%) -10.2% -23.3% -24.0% -37.2% ArgoU M L1(d = 8.3%) -04.4% -09.0% -


Quoc Hung Ly

Amina Alaoui
email: alaoui@lami.enpc.fr

Silvano Erlicher

Laurent Baly

Towards a footwear design tool: influence of shoe midsole properties and ground stiffness on the impact force during running

Keywords: Running impact force, Dynamic model, Mid-soles dynamic properties, Ground stiffness, footwear design

. In particular, Liu and Nigg introduced at the lower level of their model, i.e. at the interface between the human body and the ground, a nonlinear element representing simultaneously the shoe mid-soles and the ground flexibility. The ground reaction force is modelled as the force supported by this nonlinear element, whose parameters are identified from several sets of experimental data. This approach proved to be robust and quite accurate. However, it does not explicitly take into account the shoe and the ground properties. It turns out to be impossible to study the influence of shoe materials on the impact force, for instance for footwear design purposes. In this paper, a modification of the Liu and Nigg's model is suggested, where the original nonlinear element is replaced with a bi-layered spring-damper-mass model: the first layer represents the shoe mid-sole and the second layer is associated with the ground.

Ground is modeled as an infinite elastic half-space. We have assumed a visco-elastic behaviour of the shoe material, so the damping of shoe material is taken into account. A methodology for the shoe soles characterization is proposed and used together with the proposed model. A parametric study is then conducted and the influence of the shoe properties on the impact force is quantified. Moreover, it is shown that impact forces are strongly affected by the ground stiffness, which should therefore be considered as an essential parameter in the footwear design.

Introduction

During running activity, the human body is exposed to repetitive impact forces, which may cause the development of lower back pain and running injuries for runners [START_REF] Clement | A survey of overuse running injuries[END_REF][START_REF] James | Injuries to runners[END_REF]. This is one of the reasons why the modelling and the reduction of impact forces are important axes of investigation in the field of biomechanical research. Other reasons are the need of improving muscle tuning and comfort during running. One can also observe that athletes adjust their kinematics and kinetics quickly when changing from one impact condition to another one.

The first studies on this subject dealt with the understanding of the intrinsic attenuation mechanisms in the human body. Indeed, two different categories can be distinguished in the technical literature: (i) the active shock absorption including joint positioning [START_REF] Bobert | Mechanical analysis of the landing phase in heel toe running[END_REF][START_REF] Denoth | Load on the locomotor system and modelling[END_REF] and muscle activity [START_REF] Christina | Effect of localized muscle fatigue on vertical ground reaction forces and ankle joint motion during running[END_REF] and (ii) the passive shock absorption including synovial fluid, bone, heel pad and articular cartilage. A different mechanism has also been identified, involving the knee flexion amplitude during running: it has been proven that it also plays a dominant role in the shock absorption [START_REF] Bobert | Mechanical analysis of the landing phase in heel toe running[END_REF][START_REF] Denoth | Load on the locomotor system and modelling[END_REF][START_REF] Lafortune | Dominant role of interface over knee angle for cushioning impact loading and regulating initial leg stiffness[END_REF].

The effectiveness of these shock absorption mechanisms is bounded by biomechanical or morphological constraints, and it is very difficult to modify them in order to increase cushioning during running. Hence, a further increase of the global attenuation should be obtained by rather acting on the shoe properties (midsoles mainly) and properly accounting for ground stiffness.

Concerning the midsole influence, some experimental studies concluded that there is a strong correlation between the attenuation of impact force peaks and the decrease of shoe sole stiffness [START_REF] Lafortune | Dominant role of interface over knee angle for cushioning impact loading and regulating initial leg stiffness[END_REF][START_REF] Aerts | Deformation characteristics of the heel region of the shod foot during a simulated heel strike: the effect of varying midsoles hardness[END_REF]. Other authors showed that shoe sole stiffness has no or small influence on the impact force peaks [START_REF] Clarke | The effects of shoe cushioning upon ground reaction forces in running[END_REF][START_REF] Nigg | External Force Measurements with Sport Shoes and Playing Surfaces[END_REF][START_REF] Snel | Shock Absorbing Characteristics of Running Shoes during Actual Running[END_REF][START_REF] Devita | The effects of Time on Selected Ground Reaction Force Parameters[END_REF]. The different experimental protocols (fixed leg or heel or fixed foot-shoe combination) can partially explain these discrepancies. However, the main reason of this lack of consensus may result from the way the notion of shoe sole stiffness is defined in each analysis. Repetitive impacts during running are a dynamic phenomenon. Hence, a pertinent description of the running impact forces cannot be exclusively based on the knowledge of shoe-sole stiffness, but the damping should be also accounted for. Let us quote, for example, the work of [START_REF] Cole | Loading of the joints of the extremities during the impact phase in running[END_REF] who described the impact force as a nonlinear function of the impact rate taking into account the material properties of heel pad, shoe and contact surface. The same definition was also used by [START_REF] Nigg | The effect of the muscle stiffness and damping on simulated impact force peaks during running[END_REF][START_REF] Liu | A mechanical model to determine the influence of masses and mass distribution on the impact force during running[END_REF] to model the impact loading during running.

In the literature, the impact forces are studied by both experimental investigations and developing mechanical models of the body-shoe-ground system. Simple spring-damper-mass models have been successfully developed by [START_REF] Alexander | Mechanism in animal movement[END_REF][START_REF] Cavagna | The determinants of step frequency in running, trotting, and hoping in man and other vertebrates[END_REF][START_REF] Farley | Leg stiffness and strike frequency in human running[END_REF][START_REF] Ferris | Interaction of leg stiffness and surface stiffness during human hoping[END_REF]Ito et al., 1983;[START_REF] Mcmahon | The mechanics of running: how does stiffness couple with speed[END_REF]Kim et al., 1994;Nigg and Anton, 1994). However, all these models neglected non rigid masses of the body, like muscles and other soft tissues. Indeed, human body is constituted by both rigid parts (bones) and soft tissues linked each other through elastic and viscous connections. [START_REF] Liu | A mechanical model to determine the influence of masses and mass distribution on the impact force during running[END_REF] and [START_REF] Nigg | The effect of the muscle stiffness and damping on simulated impact force peaks during running[END_REF] proposed a model of the human body for describing running taking into account both rigid and wobbling masses. These authors used their model to study the influence of masses, mass ratios, stiffness constants and damping coefficients on the simulated impact force. A few years later, Zadpoor and co-workers [START_REF] Zadpoor | A mechanical model to determine the influence of masses and mass distribution on the impact force during running -a discussion[END_REF][START_REF] Zadpoor | A model-based parametric study of impact force during running[END_REF] brought some corrections to the parameters (impact velocities, coefficients of the function of the ground reaction force) used by Liu and

Nigg.

The objective of our study is to develop a tool for the design of shoe soles for running. The work is focused on modelling the impact of the heel on the ground and takes into account the different mechanical and biomechanical parameters involved during this impact For the biomechanics of the human body, we started from the work of Liu & Nigg and the modifications suggested by Zadpoor & Nikooyan (LNZN model). However, in this model, characterization of the shoe soles is obtained using the human pendulum test and thus doesn't correspond to an intrinsic property (neither material nor structural) of soles. Furthermore the LNZN model does not take explicitly into account the stiffness of soil.

To create this design tool, we have coupled the LNZN model with a bi-layered subsystem representing the mass, stiffness and damping of the shoe sole and the ground stiffness allowing us to integrate visco-elastic properties of soles and estimates of soil stiffness. The mechanical testing protocol of soles properties is precisely described and the ground stiffness is evaluated by a FE calculation.

We then set the dynamic equations of this new model that we have programmed into the Matlab code to study the influence of the various parameters involved.

Methods

In this section, two models are presented: (i) the [START_REF] Liu | A mechanical model to determine the influence of masses and mass distribution on the impact force during running[END_REF] model with the modifications suggested by [START_REF] Zadpoor | A model-based parametric study of impact force during running[END_REF] (LNZN model) and (ii) the new model proposed in this paper. The main novelty introduced here consists in adding to the LNZN model a bi-layered spring-damper-mass subsystem, associated with midsoles and ground. As it will be seen hereafter, this modification allows a clear interpretation of the role of the ground and of shoe sole properties on the running impact force.

LNZN model

The original model of Liu and Nigg is represented by the scheme of Figure 1a. The rigid and wobbling masses of the system are described in the caption of the Figure, and the parameters are listed in Table 1. The corresponding equations of motion read:

                           , . 4 . 3 4 4 3 ) 5 4 ( 4 .. 4 4 , . 4 . 3 4 . 3 . 1 1 4 3 ) 5 4 ( 3 2 3 3 1 1 3 .. 3 3 , . 2 . 1 2 3 2 3 2 1 2 2 .. 2 2 , . 2 . 1 2 . 3 . 1 1 2 1 2 3 1 1 1 .. 1 1 x x c x x k k g m x m x x c x x c x x k k x x k x x k g m x m x x c x x k x x k g m x m x x c x x c x x k x x k g F g m x m                                    (1)
where g is the gravity acceleration. The vertical contact force F g is defined as follows:

    0 0 0 1 1 . 1 1 1              x F x x x c x a A F g e d b C g (2)
It is a nonlinear viscoelastic force, which represents simultaneously the material properties of the heel pad, shoe and contact surface. The parameters a, b, c, d and e were assumed to be shoe-specific, since the heel pad and the ground were unchanged during the tests. Two types of shoe-foot models (soft and hard) were used in the simulations, with parameter values determined by pendulum impact tests (cf. Table 2). The constant A c was taken to be 2, assuming that the average contact area during ground contact was twice the contact area between the impact pendulum and the foot in the experiments of [START_REF] Aerts | Deformation characteristics of the heel region of the shod foot during a simulated heel strike: the effect of varying midsoles hardness[END_REF]. The same vertical velocities (0.6m/s) at touchdown were chosen for all body masses. [START_REF] Zadpoor | A model-based parametric study of impact force during running[END_REF] showed that some parameters of Liu and Nigg's model had to be modified in order to give a better fitting of the experimental results of the impact force. These authors proposed the parameter values and the touchdown velocities listed in Table 3. In summary, The LNZN model is defined by Eqs. ( 1)-( 2) and Table 3.

The new model

Starting from the LNZN model, a new mechanical model is proposed here, characterized by bi-layered subsystem representing the mass, stiffness and damping of the shoe and the ground stiffness (Figure 1b). As a result, it is possible to distinguish the ground reaction force F g from the contact force between the foot and the shoe sole, indicated by F s .

The theoretical aspects of the new model were first published in [START_REF] Ly | Influence of shoe midsoles dynamic properties and ground stiffness on the impact force during running[END_REF]. In this first version, we considered linear viscoelastic behaviour for shoe soles. The model proposed in this study includes a new method of characterization and identification of the viscoelastic behaviour of shoe soles taking into account geometric nonlinearities. Furthermore, this new model presents also a method of estimation of ground stiffness from the literature data.

Estimation of ground stiffness:

Ground is modeled as an infinite elastic half-space. A unit load F is applied over a circular area of 40 mm radius, corresponding to the average size of the contact area of sole/ground, in the form of uniform pressure.

The ground stiffness is then calculated as:

max max 1 u u F k s  
where u max represents the maximum displacement of the soil under load. Calculations were made using ABAQUS FE software.

Using literature data on elastic modulus for soils, we computed estimations of the ground stiffness using two approaches. For the study, we then considered six different soils whose properties are summarized in Table 4. Details are given in Appendix A.

Viscoelastic properties of shoe soles and model dynamics equations:

The shoe sole is modelled as a unique structure submitted to the impact of the heel during running. The load is thus considered, according to ASTM F 1976, as an impact load applied though a spherical indenter on the rear part of the sole (Figure 2). We have assumed a visco-elastic behaviour of the shoe material, where the geometric nonlinearity, due to contact between the slightly convex shape of the rear of the sole and the spherical shape of the impactor, is taken into account for calculating the soles elastic stiffness k 6 following a Hertz type contact law:

1 0 6 6 * *           b R u u k F
where R 0 is the radius of the indenter According to the scheme of Figure 1b, the expression of the ground reaction force reads

5 x k F s g   (3)
whereas the force between the foot and the shoe sole becomes:

                                     . 5 . 1 6 5 * 5 1 * 6 1 0 1 6 x x c R x x x x k s F b (4)
This approach makes it possible to insert into the model the shoe sole properties (m 5 , k 6 , b 6 , c 6 ), as well as the ground stiffness (k s ). Hence, the influence of these quantities on the impact force F s exerted on the foot can be easily evaluated (see equation 4). The equations of motion for the new system read:

                                           , . 5 . 1 6 5 * 5 1 6 5 .. 5 5 , . 4 . 3 4 4 3 ) 5 4 ( 4 .. 4 4 , . 3 . 4 4 . 3 . 1 1 3 4 ) 5 4 ( 3 2 3 3 1 1 3 .. 3 3 , . 2 . 1 2 2 3 3 2 1 2 2 .. 2 2 , . 1 . 5 6 . 1 . 2 2 . 1 . 3 1 5 * 5 1 6 1 2 2 1 3 1 1 .. 1 1 5 1 6 0 1 1 6 0 1 x x c x k R x x x x k g m x m x x c x x k k g m x m x x c x x c x x k k x x k x x k g m x m x x c x x k x x k g m x m x x c x x c x x c R x x x x k x x k x x k g m x m s b b                                                                 (5)
A MATLAB code has been developed to solve this equation system. Let us remind that these simulations do not include take off-phase running. Therefore, the vertical movement is simulated only for the first 250 ms after touchdown. The parameters and initial conditions are the same as those of the LNZN model for all masses representing the human body (Table 3). Concerning the shoe (m 5 ), the assumption that the foot is perfectly stuck on the shoe at the initial time is done. Hence, the shoe touchdown velocity is taken equal to that of the lower rigid body mass (m 1 ).

Determination of viscoelastic properties of shoe soles

The horizontal speed of the race determines the frequency of impacts, which then determines the frequency of the mechanical testing protocol for the characterization of viscoelastic shoe insoles. In this study, the viscoelastic characteristics are obtained for a frequency of 1.5 Hz.

As for the simplified model of the human body (LNZN model), the shoe sole is modelled by a mass (m 5 )spring (k 6 ) -damper (c 6 ) system. The exponent (b 6 ) is arising from the geometric nonlinearity, due to contact between the shoe sole and the indenter.

For this study, we have tested 4 running shoes. They all have the same shape and they differ only by the hardness of the midsoles material: 40 Asker C, 50 Asker C, 60 Asker C and 65 Asker C. The tests were carried out using a ZWICK ROELL HC25 hydraulic machine, and the load applied via a cylindrical indenter with a spherical tip (Figure 2). An oscillating sinusoidal force is applied on the rear area of the sole, corresponding to the supporting surface of the heel at initial impact. From the stabilized cycle, we determine the parameters k 6 (f) and b 6 (f) relating to the elastic behaviour of the sole, by interpolating the whole cycle as a power law. The value of energy dissipation during a loading cycle is then used to calculate the value of c 6 (f). An example of the testing results and the evaluation of the shoe sole mechanical parameters is shown on Figure 3.

Results

Mechanical properties of shoe soles

The results presented hereafter are related to the 4 selected shoes and correspond to a load amplitude of 1200N (F max = 1400N, F min = 200N) and a frequency of 1.5 Hz. The values of the parameters k 6 (f), b 6 (f) and c 6 (f) are computed as described in the previous section. The results are gathered in Table 5.

Reaction forces

Having now the values for the new parameters of our model, we can compute, as for the original model, the ground reaction force F g and the velocities but also the impact force F s exerted on the foot and thus the difference F g -F s corresponding to specific influence of the shoe. Results of simulated ground reaction forces are presented on Figures 4 and5, corresponding to two different soil stiffness values. On Figures 6 and7, we have plotted the difference F g -F s , that can be seen as the intrinsic cushioning of the shoe sole, for the same soils.

Loading rate and passive force peak level

The evolution of the ground reaction force allows us to determine the passive peak (first force peak level) and to calculate the loading rate (between 10% and 90% of the passive peak). As we can notice on Figures 8a and8b, the evolution of the passive peak versus the ground stiffness presents a minimum while the loading rate increases in an almost linear way. These evolutions are summarized on Figures 9a and9b for all studied shoes.

Discussion

The main purpose of this paper is a proper modelling of the impact force during running. However, the model also provides other information, like the acceleration of the lower mass. Moreover, the presence of a layer associated with the shoe-sole, allows the deformation of the shoe-sole material to be computed.

Influence of the ground stiffness k s

The usual models for evaluating the impact force do not explicitly account for the ground stiffness effect. However, the simple spring element introduced here (stiffness k s ) shows the great importance of the ground on the impact force. From Figures 9a and9b, one can see that the passive force peak and the loading rate, as well as the time evolution of the reaction force, are deeply affected by the ground stiffness. In particular, for all shoes, low value of ground stiffness can cancel out the first impact force peak. This result is important for shoe design, in order to better adapt shoe cushioning according to the running surface. Another conclusion, highlighted through Figures 6, 7 and 9c, is that the shoe cushioning is not really meaningful for soft soils, as the intrinsic gain ((Fg-Fs) max) is less than 2%, while it represents almost 10% for hard surface. Runners adjust leg stiffness for their first step on a new surface [START_REF] Ferris | Runners adjust leg stiffness for their first step on a new running surface[END_REF]. However, supposing that the leg stiffness does not change during running, it has been shown in this paper that the variations of the impact forces due to the modifications of the ground stiffness are significant and that the ground stiffness should be considered as an essential parameter in shoe cushioning design criteria.

Influence of shoe sole

As we can notice on Figures 4 and5, shoe characteristics have very small influence on loading rate but contribute significantly to the first force peak level, mainly for stiff soils (fig. 7.). Besides, for a ground of stiffness higher than 350 kN/m, we observe a good correlation between the attenuation of impact force peaks and the decrease of shoe sole stiffness (fig. 9). This result confirms works led by [START_REF] Lafortune | Dominant role of interface over knee angle for cushioning impact loading and regulating initial leg stiffness[END_REF] or [START_REF] Aerts | Deformation characteristics of the heel region of the shod foot during a simulated heel strike: the effect of varying midsoles hardness[END_REF]. While for ground stiffness lower than 350 kN/m, it seems that shoe sole stiffness has almost no influence on the impact force peaks. We find here the conclusions of [START_REF] Clarke | The effects of shoe cushioning upon ground reaction forces in running[END_REF], [START_REF] Nigg | External Force Measurements with Sport Shoes and Playing Surfaces[END_REF] or [START_REF] Snel | Shock Absorbing Characteristics of Running Shoes during Actual Running[END_REF]. The ground stiffness could thus be one reason in this lack of consensus concerning the correlation between shoe sole stiffness and the attenuation of impact force peaks.

Acceleration of lower rigid mass (m 1 )

Tibial acceleration is often used as criterion for evaluating the shoe-sole cushioning. Typical values of tibial accelerations are given in Table 6 [START_REF] Lafortune | Three-dimensional acceleration of the tibia during walking and running[END_REF]. Measured acceleration magnitudes are affected by many factors such as location and fixation of sensors, running velocity, running level and also by wearing conditions. Hence, experimental results in literature cannot be easily compared. Moreover, the "physical" location of the mass m 1 of the model corresponds only approximately to tibia.

Figure 10 illustrates and compares the time evolution of the acceleration (mass m 1 ), with k s = 880 kN/m and a load frequency of 1.5 Hz. Table 6 shows the minimum and maximum mass m 1 acceleration predicted by the model as a function of the ground stiffness, for a stride frequency of 1.5Hz.

It is difficult to compare quantitatively these values of acceleration with those of the literature, because we have no information about the subjects involved in the tests of the authors quoted on Table 7. Also, the ground stiffness is not estimated and shoes cannot be compared because they were not characterized in the same way.

However, the values of acceleration are of the same order than those in the literature.

Conclusions

A mechanical model has been developed, able to predict the time evolution of the vertical impact running force. It provides results consistent with those reported in the literature. In addition, it suggests a new approach in the modelling of the mechanical interaction between the human body and the ground through the shoes: the reaction force is no longer represented by a nonlinear function depending on coefficients without a clear physical meaning, but is defined explicitly accounting for the ground and shoe-soles material properties.

A first result obtained by this model concerns the quantitative evaluation of the ground stiffness influence on the impact running forces. Moreover, it has been proposed an approach for the characterization of the shoe-soles mechanical properties accounting for the influence of geometrical non-linearity due to the contact. This methodology has been applied to a sample of 4 shoes, but it can be used for every type of shoe material and thickness. Hence, it constitutes a very effective footwear design tool, which can be used independently or together with the proposed model. Table 1. The parameters (masses, spring constants, damping coefficients) of the system Table 2. The parameters of vertical reaction force.

List of Tables

Table 3. Parameters of the original ground reaction force model, the modified ground reaction force model and the modified touchdown velocities. Table 4. Properties of the studied soils Table 5. Mechanical properties of studied shoes at 1.5 Hz Table 6. Values of peak positive axial acceleration measured for shod running, after [START_REF] Lafortune | Three-dimensional acceleration of the tibia during walking and running[END_REF] Table 7. Values of peak acceleration (mass m 1 ) as a function of the stiffness of the running surface Figure 1 (a) A simplified spring-damper-mass model for the human body during a running step: lower body rigid mass (m 1 ) and wobbling mass (m 2 ); upper body rigid mass (m 3 ) and wobbling mass (m 4 ); compressive spring (k 1 ) and damper (c 1 ) that connect the upper and lower body rigid masses; spring (k 3 ) and spring-damper unit (k 2 , c 2 ) connecting the lower wobbling mass to the upper and lower body rigid masses; spring (k 5 ) and spring-damper unit (k 4 , c 4 ) connecting the upper wobbling mass to the upper rigid mass. (b) A new spring-mass-damper model: shoe rigid mass (m5); spring-damper unit (k6, b6, c6) connecting the lower rigid mass to the shoe rigid mass; compressive spring (ks) representing the ground stiffness Figure 2 The indenter used for loading the shoe sole and the characterization of mechanical parameters. Surface Acceleration (g) [START_REF] Clarke | The effects of shoe cushioning upon ground reaction forces in running[END_REF] 3.8 Treadmill 6.9-7.41 Clarke et al. (1985a) 3.8 Treadmill 9.55 Clarke et al. (1985b) 3.35-5.36 Treadmill 6.2-10.7 Nigg et al. (1974) Not reported Various 6.0-11 Valiant (1990) 3.83 Treadmill 7.1 [START_REF] Lafortune | Three-dimensional acceleration of the tibia during walking and running[END_REF] 3.5-4.7 Treadmill 2.9-5

Figure 2 The indenter used for loading the shoe sole and the characterization of mechanical parameters.

Figure 3

 3 Figure 3Example of testing results for the determination of shoe mechanical parameters.

Figure 4

 4 Figure 4 Ground reaction force Fg for soil S3 (Ground stiffness: 210 kN/m)

Figure 6

 6 Figure 6 Shoe cushioning (Fg -Fs) for soil S3 (Ground stiffness: 210kN/m).

Figure 7

 7 Figure 7 Shoe cushioning (Fg -Fs) for soil S6 (Ground stiffness: 880 kN/m).

  Figure 8 (a) Influence of ground stiffness on the simulated vertical ground reaction force for hard shoes; (b) Influence of ground stiffness on the simulated vertical ground reaction force for soft shoes.

  Figure 9 (a) Evolution of the passive force peak according to the ground stiffness; (b) Evolution of the loading rate according to the ground stiffness; (c) Evolution of the maximum shoe cushioning (Fg-Fs) according to the ground stiffness; (d) Evolution of the acceleration of m 1 according to the ground stiffness.

Figure 10

 10 Figure 10Acceleration of the mass m 1 (g) for soil S6 (Ground stiffness: 880 kN/m); load frequency 1.5 Hz.

  Figure 1 (a) A simplified spring-damper-mass model for the human body during a running step: lower body rigid mass (m 1 ) and wobbling mass (m 2 ); upper body rigid mass (m 3 ) and wobbling mass (m 4 ); compressive spring (k 1 ) and damper (c 1 ) that connect the upper and lower body rigid masses; spring (k 3 ) and spring-damper unit (k 2 , c 2 ) connecting the lower wobbling mass to the upper and lower body rigid masses; spring (k 5 ) and spring-damper unit (k 4 , c 4 ) connecting the upper wobbling mass to the upper rigid mass. (b) A new spring-mass-damper model: shoe rigid mass (m5); spring-damper unit (k 6 , b 6 , c 6 ) connecting the lower rigid mass to the shoe rigid mass; compressive spring (ks) representing the ground stiffness

Figure 3 Figure 5

 35 Figure 3Example of testing results for the determination of shoe mechanical parameters.

FgFigure 6 FgFigure 7

 67 Figure 6Shoe cushioning (F g -F s ) for soil S3 (Ground stiffness: 210 kN/m).

First

  Figure 9 (a) Evolution of the passive force peak according to the ground stiffness; (b) Evolution of the loading rate according to the ground stiffness; (c) Evolution of the maximum shoe cushioning (F g -F s ) according to the ground stiffness; (d) Evolution of the acceleration of m 1 according to the ground stiffness.

  Figure 10Acceleration of the mass m 1 [g] for soil S6 (ground stiffness k s = 880 kN/m); load frequency 1.5 Hz.

Table 2 .

 2 The parameters of vertical reaction force and touchdown velocities

		a	b	c	d	e
	Soft shoe	1.0*10 6	1.56	2.0*10 4	0.73	1.0
	Hard shoe	1.0*10 6	1.38	2.0*10 4	0.75	1.0
		m 1	m 2	m 3	m 4
	Velocities 0.6 m/s	0.6 m/s	0.6 m/s	0.6 m/s

Table 3 .

 3 The modified parameters of vertical reaction force and the modified touchdown velocities

		a	b	c	d	e
	Soft shoe	0.6*10 6	1.56	2.0*10 4	0.73	1.0
	Hard shoe	0.6*10 6	1.38	2.0*10 4	0.75	1.0
		m 1	m 2	m 3	m 4
	Velocities 0.96 m/s	0.96 m/s	2.0 m/s	2.0 m/s

Table 4 .

 4 Properties of the studied soils

	Reference			Material	Min Young Modulus (MPa)	Max Young Modulus (MPa)	Mean Young Modulus (MPa)	Poisson ratio	Ground stiffness k s (kN/m)
	SETRA			BBTM			5500	0.3
	Schosser			S1: multilayer	Sand+ clay + gravel (cf. appendix A.)	99
	Schosser		S2= S1 + 25mm BBTM			143
	Bridge Software Institute		S3: Soft soil	3.59	7.18	5.39	0.4	210
	Bridge Software Institute	S4= S3 + 25mm BBTM			359
	Bridge Software Institute		S5: Medium soil	7.18	14.36	10.77	0.45	429
	Bridge Software Institute		S6: Hard soil	14.36	28.73	21.55	0.5	880
	Table 5. Mechanical properties of studied shoes at 1.5 Hz
	Shoe	k 6 (N/m)		b 6	c 6 (Ns/m)
	40 Asker C	2.22E+05	2.1679	9.14E+02
	50 Asker C	2.18E+05	1.9519	1.37E+03
	60 Asker C	3.62E+05	2.3018	1.87E+03
	65 Asker C	4.03E+05	2.4346	2.17E+03

Table 6 .

 6 Values of peak positive axial acceleration measured for shod running, after[START_REF] Lafortune | Three-dimensional acceleration of the tibia during walking and running[END_REF].

Table 7 .

 7 Values of peak acceleration (mass m 1 ) as a function of the stiffness of the running surface

	Ground stiffness	Acceleration
	(kN.m -1 )	(g)
	99	2.34 -4.03
	143	2.91 -5.13
	210	3.6 -6.54
	359	4.51 -8.53
	429	4.97 -9.56
	880	6.31 -13.04


Pascale Bréhonnet
email: pascale.brehonnet@univ-brest.fr

N Tanguy

P Vilbé

L-C Calvez

Autocorrelation-Invariant Discrete-time Functions and Associated Orthogonal Sequences

Keywords: Autocorrelation, Orthogonal functions, Discrete-time functions

The paper studies frequency domain characterization and properties of autocorrelation discrete time functions. Such functions may be useful in the design and synthesis of signaling waveforms with prescribed autocorrelation functions.

Introduction

It has long been known that the orthonormal set of causal discrete time functions known as the discrete Laguerre functions [START_REF] Gottlieb | Concerning some polynomials orthogonal on a finite or enumerable set of points[END_REF][START_REF] King | Digital Laguerre filters[END_REF] an Autocorrelation-Invariant Discrete (AID) time function. A similar curiosity which holds for continuous-time Laguerre functions has led Gonsalves [START_REF] Gonsalves | Autocorrelation-Invariant Functions and an Associated Orthogonal Set of Functions[END_REF] to investigate an associated orthogonal set of Autocorrelation-Invariant Continuous time functions. AI functions find their use in the design and synthesis of signaling waveforms. For example, if certain specifications on the autocorrelation function of a discrete-time signal can be met by an AID function

[ ] k f 0 for 0 ≥ k , then [ ] n f 0
is immediately an acceptable discrete-time waveform and so are the other members of the orthogonal set derived from [ ] n f 0 which have the same autocorrelation function. Orthonormal functions have been used for reducing the complexity of an echo canceler [START_REF] Davidson | Reduced Complexity Echo Cancellation Using Orthonormal Functions[END_REF]. Some applications of Laguerre functions can be found in the references given in [START_REF] Tanguy | Online Optimization of the Time Scale in Adaptative Laguerre-Based Filters[END_REF]. The purpose of this paper is to investigate the existence and properties of a class AID C of real causal finite energy discrete-time functions [ ]

n f
with rational z-transform

( ) z F and autocorrelation function [ ] [ ] k f k ff = φ .
(1)

Frequency domain characterization

Considering the z-transform of (1) yields

( ) ( ) [ ] 0 1 1 f z F z F z F z F - ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ + = ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ . ( 2
)
Without loss of generality, the rational z-transform ( )

z F
can be written as

( ) ( ) z R z F - = 1 (3)
where ( )

z R
is an irreducible rational function

( ) ( ) ( ) z A z P z R = ˆ (4)
with ( ) ( )

A P deg deg ≤ . On account of the finite energy assumption, ( ) z R
has all its poles inside the unit circle in the z-plane. Substituting (3) into (2) yields

( ) [ ] 0 1 1 f z R z R - = ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ . (5)
Assuming [ ]

1 0 = f implies ( ) 0 = z R and then ( ) 1 = z F
, a trivial solution of (2) which presents no particular interest. Therefore, we will subsequently assume [ ] 1 0 ≠ f which implies ( ) ( )

A P deg deg = , because in case ( ) ( ) A P deg deg <
, the initial value theorem for the ztransform applied to (3) yields [ ] 1 0 = f . In view of (4), (5) may be written

( ) [ ] ( ) ( ) ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ - = ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ z A z A f z P z P 1 0 1 1 .
Multiplying both sides by m z , ( ) ( )

P A m deg deg = =
, and denoting by ( ) ( )

z A z z A m R 1 = the reciprocal polynomial, we obtain ( ) ( ) [ ] ( ) ( ) ( ) z A z A f z P z P R R 0 1 - = (6)
Because ( )

z P and ( ) z A are relatively prime ( ( ) z R
is irreducible by assumption), for a given stable polynomial ( ) z A with non-zero constant term the only acceptable solution is

( ) ( ) z A z P R λ = (7)
where λ denotes a real constant to be determined. Notice that if the constant term in ( )

z A is zero, there is no ( ) z P
for which (6) holds. Therefore, in the following, we will assume that the constant term in the polynomial ( )

z A is non-zero. Substituting (7) into (6) yields [ ] 0 1 2 f - = λ where [ ] 0 f
depends on λ . Indeed, in view of ( 3) and ( 4)

( ) ( ) ( ) z A z A z F R λ - = 1 .
Applying the initial value theorem for the z-transform of a causal function yields

[ ] m f λα - = 1 0 with ( ) ( ) ( ) ( ) 0 0 lim R R z m A A z A z A = = ∞ → α (8) hence m λα λ = 2 and the only acceptable solution is m α λ = . Therefore ( ) ( ) ( ) z A z A z F R m α - = 1 . ( 9
)
With the polynomial

( ) z A m 1 - defined as ( ) ( ) ( ) z A z A z zA R m m α - = - 1 (9) can be rewritten as ( ) ( ) ( ) z A z zA z F m 1 - = (10)
It is worth noting that

( ) z A m 1 -
is the first polynomial generated when testing stability of ( )

z A
with the widely used Aström-Raible table [START_REF] Aström | Introduction to stochastic control theory[END_REF][START_REF] Raible | A simplification of Jury's tabular form[END_REF]. Thus, every stable polynomial ( )

z A with non-zero constant term readily yields a function AID C f ∈
for which (1) holds. This function is given by its z-transform (10).

An orthogonal set

[ ] { } n f i can be generated by passing [ ] [ ] n f n f 0 =
through identical all-pass filters in Laguerre fashion, each filter having transfer function ( ) ( ) ( )

z A z A z T R =
. Some of the interesting properties of these functions are summarized in the following section.

AID-time functions properties

In the following, the cross-correlation function of two real finite-energy discrete-time functions [ ]

n f i and [ ] n f j , denoted by [ ] k ij φ , is defined as [ ] [ ] [ ] ∑ +∞ -∞ = + = n j i ij k n f n f k φ If the real causal finite-energy discrete-time functions [ ] n f i have the z-transforms ( ) ( ) ( ) ( ) ( ) i R m i z A z A z A z zA z F ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ = -1 , K , 2 , 1 , 0 = i where ( ) z A is a stable polynomial of degree 1 ≥ m
with non-zero constant term, ( )

z A R
is the reciprocal polynomial and

( ) z A m 1 - is the first polynomial, of degree 1 - m ,
n f 0 is [ ] [ ] k f k 0 00 = φ hence AID C f ∈ 0 .
This has been proved in section 2. Because an autocorrelation function is maximum at the origin, the following holds

[ ] [ ] 2 0 0 1 0 m f k f α - = ≤
k f k ii 0 = φ On account of ( ) ( ) ( ) ( ) z A z A z A z A R R = 1 1
the following holds for the energy spectra

( ) ( ) ( ) ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ = ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ = Φ z F z F z F z F z i i ii 1 1 0 0 which in turn yields [ ] [ ] k k ii 00 φ φ = .
As a corollary notice that all the [ ]

n f i have the same energy [ ] [ ] [ ] 2 0 0 2 2 1 0 0 m ii n i i f n f f α φ - = = = = ∑ ∞ = .
Figures 1 and2 illustrate the auto and cross-correlation for a set of three functions in the case where 0 f is a decaying exponential ( ( )

5 . 0 - = z z A
). The theoretical infinite sequences have been truncated to 20 terms.

n f i and [ ] n f j is [ ] [ ] [ ] k f k f k i j m i j ij 1 - - - - = α φ , i j >
which depends only on the difference i jand is zero for 0

≤ k .
This can be demonstrated as follows. The cross-energy spectrum of [ ]

n f i and [ ] n f j is ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) z F z F z A z A z A z A z A z A z A z A z A z A z A z A z F z F z i j m i j i j R m i j R R m i j R R m R m j i ij 1 1 1 1 1 1 - - - - - - - - = ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ - ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ - = ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ - ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ - = ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ = Φ α α α α α
and hence, after z-transform inversion, the above expression of the cross-correlation function results.

Because [ ] k ij φ is a linear combination of two causal functions, it equals zero for 0 < k . Letting 0 = k [ ] [ ] [ ] 0 0 0 1 - - - - = i j m i j ij f f α φ , i j >
The initial value theorem for the z-transform leads to

[ ] ( ) ( ) i m m i z i z F f α α 2 1 lim 0 - = = ∞ → from which [ ] 0 0 = ij φ follows (for i j > ). Corollary: the [ ] { } n f i
form an orthogonal set.

The following orthogonality relation holds

[ ] [ ] ∑ +∞ = = + 0 0 n j i k n f n f , j i > , 0 ≥ k .
The infinite sum on the left hand side is

n f i for 2 , 1 , 0 = i have the same autocorrelation. The cross-correlation of [ ] [ ] n f n f i 0 = and [ ] [ ] n f n f j 1 = on one hand and that of [ ] [ ] n f n f i 1 = and [ ] [ ] n f n f j 2 =
on the other hand are the same. This is due to property 3, cross-correlation only depends on the difference

1 = -i j . [ ] [ ] [ ] k f k f k i j m i j ij 1 - - - - = α φ , i j >
which depends only on the difference i jand is zero for

0 ≤ k . 4°) If [ ] n f 0 has a zero at N n = , then [ ] n f i and [ ] N n f i + are orthogonal. Indeed, from 2° [ ] [ ] [ ] ∑ +∞ = = + 0 0 n i i k f k n f n f , 0 > k Letting N k = , orthogonality of [ ] n f i and [ ] N n f i + follows from [ ] 0 0 = N f . Corollary: the zero crossings of AID C f ∈ 0 cannot be periodic
The proof is by contradiction. The integers ,... ,..., , 2

1 q N N N are zero crossings of [ ] n f 0 if for all 1 ≥ q a) [ ] 0 0 = q N f and [ ] [ ] 0 1 1 0 0 < + - q q N f N f b) [ ] n f 0 has constant sign in the interval [ [ 1 ,
0 N and in any of the intervals ] [

1 , + q q N N .
Assume that the zero crossings are given by the infinite sequence ,...

n f i with [ ] n f j is [ ] [ ] [ ] [ ] n f n f n f n f j i m j i j i 1 + + + - = * α
which depends only on the sum j i + .

Starting with the z-transform ( ) ( )

z F z F j i
, the proof is similar to that of the cross-correlation property in 3°.

6°) Consider a second function

AID C g ∈ 0 and let [ ] n g j have the z-transform ( ) ( ) ( ) ( ) ( ) j R R r j z B z B z B z B z G ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ - = β 1 with ( ) z B
a stable polynomial of degree 1 ≥ r with non-zero constant term, ( ) ( )

0 0 R r B B = β , and let [ ] n c denote the discrete convolution [ ] [ ] [ ] [ ] n g n f n c n c j i ij * = = .
Then, the autocorrelation of [ ] n c is given by the nice relationship This can be demonstrated as follows. Because ( ) ( ) ( )

[ ] [ ] [ ] k g k f k cc 0 0 * = φ (
z G z F z C j i = the energy spectrum ( ) ( ) ( ) z C z C z cc 1 = Φ may be rearranged as ( ) ( ) ( ) ( ) ( ) ( ) ( ) [ ] ( ) ( ) [ ] [ ] { } [ ] { } n g Z n f Z z G z G z F z F z G z G z F z F z j j i i cc 0 0 0 0 0 0 1 1 1 1 = = = Φ
Z-transform inversion leads to the expression (11) and the proof is achieved. The property can be easily generalized for more than two functions. Since the autocorrelation function of is strictly positive : two functions in AID C cannot be orthogonal..

Conclusion

In this paper, the existence and interesting properties of Autocorrelation-Invariant Discrete (AID) time functions with prescribed spectrum have been established. Associated sets of orthogonal sequences with identical spectrum are readily derived. The derived properties may find their use in the design and implementation of communication systems.

  can be generated by passing [ ] ( ) [ ] step function), through a cascade of iidentical all-pass discrete-time filters having transfer function i th member of the set. Then each member has the same autocorrelation function since they have identical energy spectra. Furthermore, this autocorrelation [ ] which is, in this sense, invariant under the operation of autocorrelation. For this reason we call [ ] n l 0

  Figure 1 -Time functions 2 1 0 , , f f f

  Figure 2 -auto and cross-correlation functions of 2 1 0 , , f f f

  respectively recognized as the autocorrelation functions of [ ]

  is independent of i , j , all the [ ] n c ij have the same energy i.e.

  product of two arbitrary functions 0 f and 0 g in AID C

  generated via Aström-Raible table, then the following properties are verified:

1°) The autocorrelation function of [ ]