Jakob Puchinger
email: jakobp@csse.unimelb.edu.au

Peter J Stuckey

Mark Wallace

Sebastian Brand
email: sbrand@csse.unimelb.edu.au

From High-Level Model to Branch-and-Price Solution in G12

The G12 project is developing a software environment for stating and solving combinatorial problems by mapping a high-level model of the problem to an efficient combination of solving methods. Model annotations are used to control this process. In this paper we explain the mapping to branch-and-price solving. G12 supports the selection of specialised sub-problem solvers, the aggregation of identical subproblems, automatic disaggregation when required by search, and the use of specialised branching rules. We demonstrate the benefits of the G12 framework on three examples: a trucking problem, cutting stock, and two-dimensional bin packing.

Introduction

Combinatorial optimisation problems are easy to state, but hard to solve, and they arise in a huge variety of applications. Branch-and-price is one of many powerful methods for solving them. This paper describes how Dantzig-Wolfe decomposition, column generation and branch-and-price are integrated into the hybrid optimisation platform G12 [START_REF] Stuckey | The G12 project: Mapping solver independent models to efficient solutions[END_REF]. The G12 project is developing a software environment for stating and solving combinatorial problems by mapping a highlevel model of the problem to an efficient combination of solving methods. We call such a combination of methods a hybrid algorithm. Because there is no method for choosing the best way to solve a given problem, we believe the (human) problem solver must be able to experiment with different hybrid algorithms. To meet this purpose the G12 project is developing user-controlled mappings from a high level model to different solving methods. These mappings must satisfy three conflicting objectives. They must be efficient, enabling the human problem solver to tightly control the behaviour of the algorithm if necessary for performance; flexible, allowing plug-and-play between different sub-algorithms; easy-to-use and easy-to-change for efficient experimentation with alternative hybrid algorithms.

The mapping to branch-and-price presented in this paper is designed to meet all three objectives (in reverse order):

-The user can select branch-and-price and control its behaviour by annotating a high-level model of the problem. -The generated algorithm can use a separate solver for the subproblem. The user can control the decomposition and select the subproblem solver by further annotations. -Inefficiencies arising as a result of identical subproblems are avoided by aggregating them, but the user is still enabled to express search control in terms of variables in the original model. The system also supports specialised branching rules, allowing fine-grained control of search where necessary.

The G12 platform consists of three major components, the modelling language Zinc [START_REF] Garcia De La Banda | The modelling language Zinc[END_REF], the model transformation language Cadmium [START_REF] Duck | ACD term rewriting[END_REF], and several internal and external solvers written and/or interfaced using the general-purpose programming language Mercury [START_REF] Somogyi | The execution algorithm of Mercury, an efficient purely declarative logic programming language[END_REF]. All solvers and solver instances are specified in terms of their specific capabilities, i.e. the type of problems they can solve, the type of information they can return, and how they solve a problem.

On the Mercury level these specifications are described using type classes. Basic solvers such as a Finite Domain Constraint Programming (FD) solver or a Linear Programming (LP) solver can be used as building blocks for other solvers. The column generation and branch-and-price modules use and implement such solver type classes. This system of pluggable components allows us to quickly design new hybrid algorithms and to combine existing solvers in innovative ways.

Trucking problem. Consider the following trucking problem, inspired by [START_REF] Boland | A column generation approach to delivery planning over time with inhomogeneous service providers and service interval constraints[END_REF]. We are given N trucks; each truck has a cost and an amount of material it can load. We are further given T time periods; in each time period a given demand of material has to be shipped. Each truck also has constraints on usage: in each consecutive k time periods it must be used at least l and at most u times. The Zinc model of the problem follows: At each time period we need to choose which trucks to use in order to ship enough material and satisfy the usage limits. The sum set(s, f ) function returns Σ e∈s f (e), while the sequence([y 1 , . . . , y n ], l, u, k) constrains the sum of each subsequence of length k,

y i + • • • + y i-k-1 , 1 ≤ i ≤ n -k + 1
to be between l and u inclusive. As it stands this model is directly executable in an FD solver that supports set variables. There exist specialised propagators for sum set and sequence. In Zinc we can control the search by adding an annotation on the solve item, for example:

solve :: set_search(x, "first_fail", "indomain", "complete") minimize sum(p in Periods)(sum_set(x[p], Loads) >= Demand[p]);

which indicates we label the set variables with smallest domain first (first fail) by first trying to exclude an unknown element of the set and then including it (indomain) in a complete search.

Dantzig-Wolfe Decomposition and Column Generation

Dantzig-Wolfe decomposition is a standard way to decompose an integer programming model into a master problem and one or several subproblems [START_REF] Dantzig | Decomposition principle for linear programs[END_REF][START_REF] Desaulniers | Column Generation. GERAD 25th Anniversary Series[END_REF]. The bound on the objective resulting from the LP relaxation of the decomposed model is usually stronger than that of the original formulation (if the subproblem does not have the integrality property). This can result in a smaller search space in LP-based branch-and-bound algorithms.

The Original Problem has the form

OP: minimise k∈K c k x k subject to k∈K A k j x k ≥ b j ∀j = 1 . . . M x k ∈ D k k ∈ K.
The D k are finite sets of vectors in Z N k + implicitly defined by additional constraints. We view the elements of D k to be indexed using an index set P k ; that is, we have D k = {d k p | p ∈ P k }. We can then alternatively write

D k = {e k ∈ R N k | e k = p∈P k d k p λ k p , p∈P k λ k p = 1; λ k p ∈ {0, 1} ∀p ∈ P k }.
Substituting the x k by the λ k p in OP, we obtain the Master Problem:

MP: minimise k∈K p∈P k c k d k p λ k p subject to k∈K p∈P k A k j d k p λ k p ≥ b j ∀j = 1 . . . M (1)
p∈P k λ k p = 1 k ∈ K (2)
λ k p ∈ {0, 1} ∀p ∈ P k , k ∈ K.
Dantzig-Wolfe decomposition typically results in a Master Problem with many variables. To deal with a possibly exponential number of variables, delayed column generation [START_REF] Desaulniers | Column Generation. GERAD 25th Anniversary Series[END_REF] is used. Starting from a restricted LP-relaxation of the original problem, the Restricted Master Problem (RMP), variables (columns) are lazily included in order to find an optimal solution.

The simplex algorithm for solving linear programs proceeds from one basic feasible solution to the next one, always in direction of a potential improvement of the objective function. This is achieved by adding a variable with profitable reduced cost to the basis and by removing some other variable from it. Reduced costs can be seen as an optimistic estimate of the amount of improvement achieved by a unit increase of their corresponding variable. This is the crucial property of the simplex algorithm exploited in column generation. For every D k , a subproblem is solved to determine such variables. In case of a minimisation problem, the objective is to find feasible columns d k with negative reduced cost:

(c k -πA k )d k -µ k
where π are the dual variable values corresponding to the constraints (1) and µ k is the dual value of the kth convexity constraint (2). We do not need to find a column with maximum profit; adding a "good" column is sufficient.

Solving with G12

The G12 system allows one to take a model written in Zinc, transform it to various underlying solvers using Cadmium, and then execute it. We can use standard or user-defined Cadmium transformations. Mappings from Zinc to FD or LP models are available [START_REF] Brand | Flexible, rule-based constraint model linearisation[END_REF]. To control these transformations the user can annotate the model. The trucking problem example illustrates the use of an annotation to define search for an FD solver.

At the solver programming language (Mercury) level, G12 defines interfaces to solvers such as an FD solver, a continuous interval constraint solver, and linear programming solvers using type classes. Various implementations of these interfaces are provided, e.g. for LP/MIP solvers such as CPLEX, COIN-OR/OSI, and others. The Dantzig-Wolfe decomposition column generation, default branch-and-bound, and branch-and-price solvers heavily rely on the LP solver interfaces. These interfaces provide standard predicates for variable creation, constraint posting, setting an objective function, and LP and MIP optimisation.

The advantage of this architecture is that we can easily plug different LP solvers into modules such as column generation and branch-and-bound.

Dantzig-Wolfe Decomposition and Column Generation in Zinc

In order to use Dantzig-Wolfe decomposition and column generation on a highlevel model in G12, we need to annotate the model to explain: what parts define the sub-problems, which solver is to be used for each subproblem, and which solver is to be used for the master problem.

For instance, the trucking problem example can be annotated as follows: which exposes which variables x will be used in column generation. For each Period a subproblem is defined in terms of its constraints and solver. Note that we could have used a more specialised solver here since the subproblem is a knapsack problem. Finally, the solver for the master problem and the search specification, branch-and-bound selecting the most fractional variable first and performing a standard split, are attached to the solve item.

We then perform a Dantzig-Wolfe decomposition on the model, separating original, master, and subproblem variables, as well as adding constraints linking those variables: The so called master variables are place-holders representing the implicit sums of the λ variables p∈P k d k p λ k p as introduced in MP. Note that the search is still expressed in terms of the original problem variables.

Since column generation is to be used, the transformation must linearise the master constraints and objective function. The subproblem solver could use the original set representation of the variables, but for this example it too requires linearisation of the sub-problem constraints.

We can linearise the master and subproblem constraints giving linear definitions for the sum set and sequence globals. This can be done in Zinc as:

function var int:sum_set(var set of $T:s, array[$T] of int: cost) = sum(e in index_set(cost))(cost[e] * bool2int(e in s)); predicate sequence(array[int] of var int:y, int:l, int:u, int:k) = forall(i in min(index_set(y)) .. max(index_set(y)) -k + 1)( let { var int: s = sum(j in i .. i + k -1)( y[j] ) } in s >= l /\ s <= u);
where index set returns the set of indices of its array argument, and bool2int coerces a Boolean to 0..1. Finally, we transform the array of set variables x to a two-dimensional array of 0..1 variables such that x[p,t] = 1 if t in x[p].

Implementation

Column generation works by first transforming the original model to the form demanded by the solvers. Then it builds the subproblems and attaches them to the requested solvers. These solvers must support optimisation with a linear objective function, and preferably support it in an incremental way.

Then the restricted master problem is defined and attached to a solver that supports delayed column generation: currently LP solvers, although we are working on adding a hybrid volume algorithm/LP solver [START_REF] Barahona | The volume algorithm: producing primal solutions with a subgradient method[END_REF][START_REF] Anbil | Column generation and the airline crew pairing problem[END_REF].

The G12 Dantzig-Wolfe decomposition and column generation solver interface implements most of the standard functionality of the G12 LP solver interface. From the outside it looks mostly like a standard LP solver set up with the original problem using the original (linearised) variables. The mapping between the original variables and the master problem variables is straight-forward; we simply set

x k = p∈P k d k p λ k p .
The main difference lies in the initialisation of the column generation module. First the subproblem solver instances have to be added, then the variables to be decomposed are created, and finally the master problem constraints are posted.

Similarly to the simplex algorithm, column generation requires an initial feasible solution. If it is not provided by the user, we introduce artificial variables in order to determine it automatically. At the end of this first phase the artificial variables are removed from the problem [START_REF] Vanderbeck | Branching in branch-and-price: a generic scheme[END_REF].

Since the column generation algorithm alone only solves the LP-relaxed version of the problem, we have to branch in order to guarantee integrality of the variables. The default branch-and-bound module is a simple, standard linear programming based branch-and-bound algorithm branching on the original model variables, which does not affect the subproblem structure [START_REF] Villeneuve | On compact formulations for integer programs solved by column generation[END_REF].

The additional branching constraints could of course render the RMP infeasible. But, since we are usually not dealing with the complete master problem, additional columns could restore feasibility of the RMP. Such columns are obtained by solving a problem very similar to the pricing problem [START_REF] Gunluk | A branch-and-price algorithm and new test problems for spectrum auctions[END_REF].

The availability of the original variables in the column generation solver is the key to being able to use this solver in further hybrids. We can use it with an arbitrary search strategy on the original variables, or for example in combination with an FD solver, by communicating bounds on the original variables.

The Trucking Problem

We solved several different instances of our trucking example showing the advantages of using DW-decomposition. Table 1 shows results on five different instances displaying the number of search nodes and the time required for solving the model using an FD solver, using a linearised model with branch-and-bound (LP-BB), and using DW-decomposition and column generation (DW). For the trucking example the DW-decomposition is so strong that it yielded the optimal integral solution in the root node without a need to branch; so instead of nodes we show the number of columns generated for the DW-decomposed problem. We also display the value of the LP-relaxation at the root node for the linear models.

For this problem we used our own branch-and-bound module using CPLEX as LP solver as well as IP subproblem solver. In general, any kind of LP solver (with G12 interfaces) can be used as master solver, and also any kind of subproblem solver is supported.

Identical Subproblems

Dantzig-Wolfe decomposition often results in highly symmetrical models because of structurally identical subproblems, i.e. the objective coefficients, the master problem constraints and the subproblem constraints are identical. A typical example for such a model is the cutting stock problem [START_REF] Kantorovich | Mathematical methods of organizing and planning production[END_REF][START_REF] Gilmore | A linear programming approach to the cuttingstock problem (part I)[END_REF].

Aggregating Identical Subproblems

Solving problems with identical subproblems by the pure Dantzig-Wolfe approach can be quite inefficient. This issue is usually overcome by aggregating the identical subproblems. The set K of subproblem indices is partitioned into sets K s by grouping the indices of identical subproblems; s ranges over some S. We turn

A s j d s p λ s p ≥ b j ∀j = 1 . . . M p∈P s λ s p = |K s | s ∈ S λ s p ≤ |K s |, λ s p ∈ Z + ∀p ∈ P s , s ∈ S.

Automatic Disaggregation when Branching on Original Variables

The direct mapping between the original variables and the newly introduced variables is not obvious anymore. In the aggregated case we have

x k = p∈P s λ s p d s p /|K s |.
Unfortunately, this usually leads to highly fractional values for the original variables, even if the λ s p variables take integer values. We therefore first decompose the λ s p values into (non-aggregated) λ k p values preserving integrality as much as possible, and then we use the mapping for the non-aggregated case.

In order to allow branching on the original variables we have to disaggregate the problem as required by the branching. The column generation module allows one to post any kind of linear constraint on the original problem variables without affecting the subproblem structure. Each aggregated subproblem appearing in these constraints is automatically disaggregated and considered by the column generation iterations in the subsequent nodes. Given K identical subproblems, if a constraint is posted involving an original variable belonging to the kth subproblem, this subproblem becomes different to the others and is disaggregated (while the remaining K -1 subproblems are kept aggregated). In order to implement this complex behaviour, the column generation module maintains a mapping between the original variables and their associated subproblems. It also tracks the aggregation status of all the subproblems by keeping a list of active subproblems. The disaggregations are rolled back upon backtracking.

The Cutting Stock Problem

In the cutting stock problem, we are given N items with associated lengths and demands. We are further given stock pieces with length L and an upper bound K on the number of required stock pieces for satisfying the demand (a trivial upper bound is the sum over all the demands). The following Zinc model corresponds to the formulation by Kantorovich [START_REF] Kantorovich | Mathematical methods of organizing and planning production[END_REF]:

] | k in 1..K]); constraint forall(i in 1..N)(sum([items[k, i] | k in 1..K]) >= i_demand[i]); constraint forall( k in 1..K)( sum(i in 1..N)(items[k,i] * i_length[i]) <= pieces[k] * L :: colgen_subproblem_constraint(k, "knapsack"));
The original model is a linear program. The annotations for the column generation variables and subproblems are as before. But this time we introduce a new annotation colgen symmetric which annotates a type. This indicates that the model is symmetric in this dimension and the resulting column generation should aggregate in this dimension. A Cadmium transformation can then be used to create an aggregated version of the variables and constraints as follows: The colgen link constraints associate the aggregated master and subproblem variables with the original problem variables. The m pieces and m items variables are place-holders representing the implicit sums of aggregated λ variables p∈P s d s p λ s p as introduced in the AMP. The s pieces and s items variables are the actual subproblem variables. This model is similar to the wellknown column generation formulation first described by Gilmore and Gomory [START_REF] Gilmore | A linear programming approach to the cuttingstock problem (part I)[END_REF], although that does not retain the original variables. Note that it is conceivable to use Cadmium to detect symmetries and automatically add colgen symmetric annotations.

In the following experiment we evaluate possible differences when using the aggregated and the non-aggregated DW-decomposition. The results are shown in Table 2. We display in percent how often an optimal solution, a feasible solution, or no solution was found. We further give average objective values and number of explored nodes where at least a feasible solution was found. Average run-times over all the instances are also shown. The maximum run-time per instance was 5 minutes. We used CPLEX as LP solver and a specialised dynamic programming algorithm implemented in Mercury for solving the knapsack subproblems. The CPLEX MIP solver was used as primal heuristic to solve the restricted master problem to integrality at every 100th node with a time limit of 10 seconds, as specified using the colgen_ph annotation. The instances used are randomly generated using CUTGEN1 [START_REF] Gau | CUTGEN1: a problem generator for the standard onedimensional cutting stock problem[END_REF]. Instances of Classes 1-12 have stock length L = 1000; each class consists of 10 instances.

For almost all classes, aggregating identical subproblems presents an advantage in the number of solved instances, solution quality and solving time.

Specialised Branching Rules

In order to overcome symmetry issues, specialised branching rules for specific problem types were developed; see e.g. [START_REF] Barnhart | Branch-and-price: Column generation for solving huge integer programs[END_REF]. They usually require changes to the subproblems during the branch-and-bound process. G12 enables users to implement such specialised branching rules, changing the structure of the subproblems, but preserving aggregations. The column generation module allows one to ask for fractional columns of the DW-decomposed model. It returns their values as well as their entries in the constraint matrix of the master problem. Using this information the user can define specialised branching rules by introducing constraint branches on subproblem variables. In the master problem these constraint branches can be enforced by setting forbidden columns to zero in their respective branch. The column generation module provides a predicate by which the user can specify a list of column patterns that have to be set to zero. In our current system the specialised branching rules are implemented in Mercury. We are working on extensions to Zinc so that users will be able to specify such rules at the modelling level.

The Two-Dimensional Bin Packing Problem

In order to demonstrate the effectiveness of specialised branching rules we implemented a simple, well-known rule for the two-dimensional bin packing problem. It is similar to the one described in [START_REF] Puchinger | Models and algorithms for three-stage twodimensional bin packing[END_REF], which is based on a well known branching rule for set partitioning [START_REF] Ryan | An integer programming approach to scheduling[END_REF]. The solution space is divided by branching on whether two different items are in the same bin. We always choose the two highest items appearing in a pattern whose corresponding column generation master variable λ has an LP solution value closest to 0.5.

In the two-dimensional bin packing problem (2DBPP), we are given N rectangular items of given height and width. These items have to be placed on (or cut out) of bins of height H and width W, of which there are at most K. The variant we consider here does not allow items to be rotated; only level packings are allowed. Each bin can be divided into several levels, and each level contains the items beside each other [START_REF] Lodi | Models and bounds for two-dimensional level packing problems[END_REF]. For ease of modelling, we assume that the items are sorted by non-increasing heights. The formulation in Zinc is as follows: The bp annotation to the solve item tells the system to use the branch-and-price algorithm choosing the most fractional master variable and using the specialised branching rule.

Table 3 displays the results of applying standard branching on the original variables or using the specialised branching rule. We tested these approaches on the set of 500 instances described in [START_REF] Lodi | Models and bounds for two-dimensional level packing problems[END_REF]. They are divided in 10 classes of 50 instances each, with item numbers ranging from 20 to 100 in each class. While many instances could be solved to optimality in the root node, our specialised branching rules did reach optimal solutions more often in the given limited runtime.

Related Work and Conclusion

The practical usefulness of column generation and branch-and-price has been well-established over the last 20 years [START_REF] Desaulniers | Column Generation. GERAD 25th Anniversary Series[END_REF][START_REF] Barnhart | Branch-and-price: Column generation for solving huge integer programs[END_REF]. More recently it has emerged that column generation provides an ideal method for combining approaches, such as constraint programming, local search, and integer/linear programming. Columns can be generated by constraint programming or application-specific algorithms, while the master problem is handled using branch-and-price [START_REF] Junker | A framework for constraint programming based column generation[END_REF][START_REF] Yunes | A hybrid approach for solving large scale crew scheduling problems[END_REF][START_REF] Rousseau | Solving VRPTWs with constraint programming based column generation[END_REF][START_REF] Puchinger | Models and algorithms for three-stage twodimensional bin packing[END_REF].

For systems such as G12 that support hybrid algorithms, Dantzig-Wolfe decomposition, column generation and branch-and-price provide an elegant way for the different solving techniques to be combined. However, the specification of this form of hybrid is quite complex, as it requires adaptation of simplexbased approaches to support the lazy generation of columns. Thus systems such as ABACUS [START_REF] Jünger | The ABACUS system for branch-and-cut-and-price algorithms in integer programming and combinatorial optimization[END_REF], MINTO [START_REF] Nemhauser | MINTO, a Mixed INTeger Optimizer[END_REF], OPL script [START_REF] Van Hentenryck | OPL Script: Composing and controlling models[END_REF], MAESTRO [START_REF] Chabrier | Génération de Colonnes et de Coupes utilisant des sous-problèmes de plus court chemin[END_REF], COIN/BCP [START_REF] Ralphs | COIN/BCP user's manual[END_REF], and SCIP [START_REF] Achterberg | SCIP -a framework to integrate constraint and mixed integer programming[END_REF] offer facilities to support the implementation of branch-and-price on top of generic integer/linear programming packages. However, these systems still require the user to understand the technical details of branch-and-price: the purpose was to support algorithm implementation rather than problem modelling.

Certainly column generation is technical, but for people trying to solve combinatorial problems the most important requirement is to be able to try out an algorithm, or more generally a hybrid algorithm, quickly and easily without rewriting the problem specification. The first attempt to provide a column generation library was in ECL i PS e [START_REF] Eremin | Using Dual Values to Integrate Row and Column Generation into Constraint Logic Programming[END_REF]. This system introduced the idea of an aggregate variable appearing in the master problem to represent a set of values returned as columns from multiple solutions to identical subproblems. However this library assumes a fixed set of variables in each subproblem, and precludes search choices which break some of the subproblem symmetries. In order to achieve tight control over branch-and-price, sophisticated ECL i PS e users have required special adaptations of the column generation library in order to be able to work directly with low level primitives [START_REF] Papadakos | Integrated airline scheduling[END_REF].

The facility to annotate the same Zinc model in two different ways, as in the examples above, and thus have the problem solved by the FD solver, or by column generation according to the annotation, is completely novel. Moreover the facility to perform search on user variables and have any symmetries which are dynamically broken during search still correctly and efficiently handled automatically by the column generation solver is also new. Thirdly the facility to define specialised search still using the mapping provided by the library provides the full flexibility needed by the expert user.

The G12 scheme is to add annotations to a conceptual problem model, and thus turn it into a design model that maps to a specific algorithm. Annotating a constraint, occurring in the conceptual model, with the (name of the) solver that will handle it, is a simple example of this scheme.

Column generation is an interesting challenge because it does not naturally fit into the above scheme. Certainly we view the column generation module as a solver in the normal way (as discussed in Section 3). However annotating a constraint with the column generation solver is not enough: the solver needs to know which subproblem the constraint belongs to, the master problem or the subproblem. Moreover there is not one column generation solver: the master problem might be sent to one underlying solver and the subproblem to another. Finally branch-and-price search is closely connected with the column generation solver, and annotations to control the search can be crucial to the performance of the algorithm.

Each requirement has been satisfied in Zinc by having a sufficiently expressive annotation language. For example an annotation with a compound term (colgen_subproblem_constraint(p, "mip")) was used to specify the subproblem solver in Section 3.1, and the search was specified by multiple annotations.

The next particular challenge of column generation is that the variables (and constraints) used in the conceptual model of the problem are quite different from those needed in the design model. Our column generation module automates this mapping using G12's Cadmium mapping language. To ensure the annotations are still meaningful with respect to the new variables, the annotations have to be transformed by Cadmium in the same way. Moreover the search control as illustrated in Section 4 must be mapped to search steps expressed in terms of the design model variables.

The greatest design and implementation challenge was to have these still work, fully automatically, when handling symmetry by generating aggregated variables (used when solving the subproblem) and dynamically disaggregating some of them during search. Indeed, each symmetry-breaking search step causes the design model to be updated so as to operate on a new set of variables.

The design model is expressed in terms of a simplified version of Zinc, illustrated in Section 4.3. The specification of our language for expressing design models is still fluid, and so currently the translation to the Mercury codewhich is very similar, but uses different syntax -is by hand.

One interesting challenge arising out of this work is how to automatically detect identical subproblems. This is a completely novel form of automated symmetry detection, which is of significant practical value, as the results in Table 2 reveal.

We plan to implement the search annotation transformations necessary to enable specialised branching schemes to be expressed in Zinc. We also plan to build in an implementation of the generic branching scheme described in [START_REF] Vanderbeck | Branching in branch-and-price: a generic scheme[END_REF]. We further intend to address issues related to adding multiple columns and column pool management.

Finally we envisage to explore the use of the column generation module for solving a subproblem within a larger problem -thus supporting, for example, a combination of row and column generation.

  array[Periods] of var set of Trucks: x :: colgen_var; constraint forall(p in Periods)( sum_set(x[p], Load) >= Demand[p] :: colgen_subproblem_constraint(p, "mip")); solve :: colgen_solver("lp") :: lp_bb(x, most_frac, std_split) minimize sum(p in Periods)(sum_set(x[p], Cost));

  type Pieces = 1..K :: colgen_symmetric; int: N; type Items = 1..N; int: L; array[Items] of int: i_length; array[Items] of int: i_demand; array[Pieces] of var 0..1: pieces :: colgen_var; array[Pieces, Items] of var int: items :: colgen_var; solve :: lp_bb([pieces, items], most_frac, std_split) :: colgen_ph(100, 10) :: colgen_solver("lp") minimize sum([ pieces[k

  ] of int: ItemWidth; int: H; array[Items] of int: ItemHeight; array[Bins] of var 0..1: bin :: colgen_var; array[Bins, Items] of var 0..1: item :: colgen_var; solve :: bp([bin, item], most_frac_master, special_split) :: colgen_ph("mip", 100, 10) :: colgen_solver("lp") minimize sum(k in Bins)(bin[k]); constraint forall(j in Items)(sum(k in Bins)(item[k, j]) >= 1); constraint forall(k in Bins)( is_feasible_packing(bin[k], [item[k, j] | j in Items]) :: colgen_subproblem_constraint(k, "mip")); set of tuple(Items, Items): Idx = {(i, j) | i, j in Items where j >= i}; predicate is_feasible_packing(var 0..1: l_bin, array[Items] of var 0..1: l_item) = let { array[Idx] of var 0..1: x } in forall (i in Items)( sum(j in i..N)(ItemWidth[j] * x[i, j]) <= W * x[i, i]) /\ sum(i in Items)(ItemHeight[i] * x[i, i]) <= l_bin * H /\ forall(j in Items)(l_item[j] = sum(i in 1..j)(x[i, j]));

Table 1 .

 1 The trucking example: finite domain model versus linearised branch-andbound versus DW-decomposition

	Instance	FD	LP-BB			DW
	Trucks Periods Nodes Time Nodes LP opt. Time Columns LP/IP opt. Time
	4	6	4655 0.80s 3282 177.0	0.55s	19	220.0	0.18s
	4	6	5860 0.85s 1992 177.0	0.47s	12	210.0	0.16s
	4	6	4607 0.77s 3102 177.0	0.55s	20	224.0	0.18s
	4	8	39848 5.04s 25646 267.0	2.64s	24	324.0	0.18s
	6	7	2361926 215.90s 194000 244.8 18.75s	18	287.0	0.18s

  ≤ λ s p ≤ |K s | and p∈P s λ s p = |K s |. The Master Problem MP becomes the Aggregated Master Problem:

	k∈K s	p∈P k d k p λ k p	into	p∈P s d s p λ s p
	where λ s p are integer variables satisfying 0 AMP: minimise s∈S p∈P s c s d s p λ s p
	subject to	s∈S	p∈P s

Table 2 .

 2 Results for cutting-stock with a maximum run-time of 5 min.

	Class Items			No Aggregation				Aggregation
			Opt % Feas % No %	Obj Nodes Time [s] Opt % Feas % No%	Obj Nodes Time[s]
	Class1	10	30	70	0	12.70 3325.80 210.40 30	70	0	12.60 3596.60 209.95
	Class2	10	70	10	20	118.75	125 100.89 90	10	0	112.90 283.80 59.36
	Class3	20	30	0	70	23.33 766.67 242.52 20	80	0	24.50 823.10 250.05
	Class4	20	0	0	100	n.a.	n.a. 298.63 10	30	60 222.50	400 268.17
	Class5	10	100	0	0	49.50 75.20	6.07 100	0	0	49.50	0	0.32
	Class6	10	80	10	10	518.56 38.89	68.39 100	0	0	494.90 143.40 21.84
	Class7	20	70	20	10	90.22 212.22 105.18 90	10	0	90 225.90 50.00
	Class8	20	60	0	40	947.83 16.67 184.24 90	10	0	893.50 40.60 30.51
	Class9	10	100	0	0	64	20	2.04 100	0	0	64	50	1.79
	Class10 10	80	10	10	657.67 43.78	70.08 90	10	0	639.70	169 39.27
	Class11 20	70	10	20	117.75 104.75	95.10 80	20	0	115.50 253.60 60.15
	Class12 20	70	10	20 1182.25 10.25 154.79 80	20	0 1146.90 120.60 50.06
	Average	63.33 11.67 25 330.74 457.60 128.19 73.33 21.67 5	327.46 514.61 86.79

Table 3 .

 3 Results for two-dimensional bin packing with a maximum run-time of 5 min.

	Class			Std. Branching				Sp. Branching
		Opt % Feas % No % Obj Nodes Time [s] Opt % Feas % No% Obj Nodes Time[s]
	Class1	68	22	10 19.49 45.87 109.90 90	8	2 39.90 41.14 53.54
	Class2	26	0	74	1.31	0 223.24 30	2	68 64.19 6.19 203.08
	Class3	70	10	20 13.05	10 116.37 84	8	8 13.85 11.87 82.90
	Class4	26	0	74	1.31	0 228.76 26	0	74	1.31	0 228.74
	Class5	84	6	10 17.40 8.89	69.65 90	2	8 17.61 3.93 53.13
	Class6	24	0	76	1.08	0 228.03 24	0	76	1.08	0 227.97
	Class7	76	16	8 16.30 33.70	80.52 88	10	2 16.78 57.18 57.52
	Class8	78	10	12 15.73 14.77	89.48 84	6	10 15.98 13.60 77.04
	Class9	96	4	0 42.62 6.72	13.94 100	0	0 42.60 0.32	2.17
	Class10	48	4	48	7.46 18.15 155.95 52	0	48	7.46 7.54 149.39
	Average 59.6	7.2	33.2 17.95 17.58 131.58 66.8	3.6	29.6 23.65 18.38 113.55

Acknowledgements

We would like to thank the members of the G12 team at NICTA VRL for helpful discussions and implementation work.

NICTA is funded by the Australian Government as represented by the Department of Broadband, Communications and the Digital Economy and the Australian Research Council through the ICT Centre of Excellence program.


Bruno Torrésani

An Overview of Wavelet Analysis and Time-Frequency Analysis (a minicourse)

We describe several aspects of wavelet analysis and more general methods of time-frequency analysis, emphasizing applications to signal analysis and processing problems.

Introduction

Time-Frequency analysis has been introduced from the need of developing a mathematical version of musical notation. Classical analysis is generally based upon a representation of functions, either the usual one (hereafter termed the "time representation") or the Fourier representation ("frequency representation"): formally f (ω) = f (t)e -iωt dt .

It is known that the frequency representation yields a lot of difficult mathematical problems (see e.g. [START_REF] Zygmund | Trigonometric Series[END_REF]). However, the frequency representation has also been criticized for its lack of physical significance. From a physical point of view, it does not make sense to think of a function or a signal as being a superposition of functions which do not possess any time localization properties. Let us take the example of a musical signal. Our ear interprets it as a series of notes, or "atoms of sound" appearing at given times, with a finite duration and a given height (the fundamental frequency). If the signal contains a given note once, say A, a Fourier representation of the signal will exhibit a peak at the corresponding frequency, without any indication of location or duration. We now quote J. Ville:

"...the representation is mathematically correct because the phases of the tones close to A have managed to suppress it by interference phenomena before it is heard, and to enforce it, again by interferences, when it is heard... However this is a deformation of reality: when the A is not heard, it is simply because it has not been played yet..."

This motivates the search for a mathematical analogue of the musical notation: a way of representing functions in terms of joint time-frequency variables.

1 email address: Bruno.Torresani@sophia.inria.fr 1 This programme of time-frequency representations has been achieved in various ways and various contexts. The first contributions seem to be due to Ville on one hand, who proposed to use the Wigner distribution as a "timefrequency density", and Gabor on the other hand, whose approach was based upon decompositions into "time-frequency atoms", generated as time and frequency translates of a Gaussian function. Since then, Ville's approach has led to the theory of quadratic time-frequency representations, very popular in the signal analysis community. Gabor-type expansions are also very popular, and have more recently enjoyed a renewed interest, since the development of wavelet theory in the early eighties.

A. Grossmann and J. Morlet realized in 1983 that any square-integrable function may be expanded into wavelets, namely translates and dilates of a single function (the mother wavelet), provided that the latter possess some mild oscillation properties. This simple result, and the remark that the corresponding expansions are in fact very close to Gabor's expansions, showed that wavelet expansions had an important potential as a tool for signal analysis. Grossmann and his collaborators also provided a beautiful interpretation of wavelet and Gabor expansions in terms of square-integrable group representations, emphasizing the importance of symmetry groups in the construction.

At about the same time, motivated by image processing problems, and the construction by Y. Meyer of an orthonormal basis of L 2 (R) consisting in a discrete system of wavelets (see [START_REF] Meyer | Ondelettes et opérateurs,I: Ondelettes; II: Opérateurs de Calderón-Zygmund[END_REF] for example), S. Mallat [START_REF] Mallat | A theory for multiresolution signal decomposition: the wavelet representation[END_REF] developed the theory of Multiresolution analysis, providing a sort of algorithmic framework for generating wavelet bases. This discovery had in fact a tremendous impact in various fields, as it realized a synthesis between several apparently disconnected fields: image analysis and computer vision, mathematical analysis, signal coding,... A major achievement was the construction by I. Daubechies [START_REF] Daubechies | Orthonormal bases of compactly supported wavelets[END_REF] of orthonormal bases of compactly supported wavelets. The construction was generalized in many respects since then: construction of new compactly supported wavelet orthonormal bases (see [15] for a general account), biorthogonal bases [START_REF] Cohen | Biorthogonal bases of compactly supported wavelets[END_REF], wavelet bases on intervals and open sets, generalized subband decompositions [START_REF] Coifman | Remarques sur l'analyse de Fourier à fenêtre[END_REF][START_REF] Wickerhauser | Adapted Wavelet Analysis, from Theory to Software[END_REF] (wavelet packets, local trigonometric bases,...), and many other developments.

The wavelet theory has also evolved very much along lines motivated by applications. Besides signal and image processing applications, which were in some sense the first motivations for wavelet-based methods, new applications have gradually emerged. Without being exhaustive again, let us quote applications to numerical analysis [START_REF] Beylkin | Fast wavelet transforms and numerical algorithms I[END_REF], dynamical systems [START_REF] Bacry | Singularity spectrum of fractal signals from wavelet transform: exact results[END_REF] and turbulence, quantum physics, and more recently statistics. All these subjects brought additional contributions to the general theory of wavelets.

The present paper is a written version of lectures given by the author at the Dubna workshop on "self-similar systems", in July 1998. The goal of the lectures was to present a general introduction to the wavelet theory and some more recent developments, with a necessarily "coarse sampling" of the subject. Obviously, the theory is by now too large to fit into a single survey paper. Therefore, we have chosen to limit our presentation to specific aspects, mostly related to signal processing problems. We start our presentation with the introduction of wavelet transforms in the context of time-frequency analysis, in parallel with the so-called "quadratic time-frequency transforms", and the continuous Gabor transform. After describing some classical time-frequency problems, we illustrate the 2 Time-frequency analysis; continuous transforms

The whole story of wavelet analysis started with the introduction of the continuous wavelet transform in a seminal article by A. Grossmann and J. Morlet [START_REF] Grossmann | Decomposition of Hardy functions into square integrable wavelets of constant shape[END_REF]. Many earlier contributions should in fact be quoted, both in the signal processing and mathematical literature, but it was only in that paper that wavelet transform was considered as such.

Like the continuous Gabor transform (also called short-time Fourier transform), the continuous wavelet transform is a prototype of the so-called linear time frequency representations, which provide decompositions of functions as superpositions of elementary waveforms. Those linear time-frequency representations are often compared with the quadratic time-frequency representations, a family of tools which are considered more powerful in some respects, but generally lack of robustness in complex situations. We start our analysis by a short account of the quadratic time-frequency representations.

Quadratic time-frequency transforms

Definitions and general properties The first instance of time-frequency transform is the ambiguity function, familiar to radar specialists (see e.g. [START_REF]Radar and Sonar, Part I[END_REF]). The (cross) ambiguity function A f,g (b, ω) of a pair of functions f (t) and g(t) may be introduced as the solution of the following problem: let g(t) be a fixed reference signal, and let f (t) be an observation, assumed to be a noisy version of a time and frequency shifted copy of g(t):

f (t) = Ag(t -τ )e iξt + n(t) ,
where A is a constant and n(t) is a random perturbation, supposed to be a white noise for the sake of simplicity. Let us assume that g ∈ L 2 (R). If the random component is absent, the problem may be formulated as follows: find a family of continuous linear forms L b,ω such that for all τ, ξ), L b,ω (g τ,ξ ) attains its maximum for (b, ω)

= (τ, ξ). If n(t) is present, L b,ω (f ) = AL b,ω (g τ,ξ ) + L b,ω (n).
According to Riesz's representation theorem, there exists a function

ϕ b,ω (t) ∈ L 2 (R) such that for all f ∈ L 2 (R), L b,ω (f ) = f, ϕ b,ω . Then, denoting by E {X} the expectation of a random variable X, E |L b,ω (n)| 2 = ||ϕ b,ω || 2 ,
and the optimal ϕ b,ω is defined as the one which maximizes the Signal to Noise Ratio (SNR)

ρ 2 = | g τ,ξ , ϕ τ,ξ | 2 /E |L b,ω (n)| 2 .
The solution ϕ b,ω = Kg b,ω (where K = 0 is a constant) is an immediate consequence of the Cauchy-Schwarz inequality. Therefore, the parameters τ and ξ may be obtained by maximizing with respect to b, ω the square-modulus of

f (t)g(t -b)e -iωt dt .
The latter is essentially a time-frequency cross-correlation of f (t) and g(t), obtained by considering scalar products of f (t) with time and frequency shifted copies of g(t) of the form e iωt g(t -b) (with b, ω ∈ R). Therefore, it measures how close f (t) is to time-frequency shifted copies of g(t).

The definition (including the case of random time series) is given below.

DEFINITION 2.1 1. Let f ∈ L 2 (R).
Its ambiguity function is defined by

A f (τ, ξ) = f (t + τ /2)f (t -τ /2)e -iξt dt . (1)
2. Let {X t , t ∈ R} a second order random time series. Then its ambiguity function is defined by

A X (τ, ξ) = E X t+τ /2 X t-τ /2 e iξt dt . (2)
It is easily seen 4 . More generally, it may be shown that A ∈ L p (R 2 ) for all p ∈ [1, ∞] (bounds for the corresponding L p (R 2 ) norms have been derived by E. Lieb [START_REF] Lieb | Integral bounds for radar ambiguity functions and Wigner distributions[END_REF]). More generally, the definition of ambiguity function may be extended to distributions, and it may be shown that the "ambiguity distribution" of a distribution φ

that if f ∈ L 2 (R), then A f is a bounded function (with ||A f || ∞ ≤ ||f || 2 ). In addition, a direct calculation shows that if f ∈ L 2 (R), then A f ∈ L 2 (R 2 ), and that ||A f || 2 2 = 2π||f ||
∈ S ′ (R) is a distribution A φ ∈ S ′ (R 2 ).
The ambiguity and cross-ambiguity functions possess several interesting properties. Among them, the so-called Moyal's formula is extremely important:

PROPOSITION 2.1 Let f, f ′ , g, g ′ ∈ L 2 (R). Then A f,g , A f ′ ,g ′ ∈ L 2 (R 2 ), and A f,g , A f ′ ,g ′ = 2π f, f ′ g, g ′ . (3)
The proof of the proposition follows from a simple calculation.

The non-deterministic version may be given a similar interpretation. Its properties depend on the properties of the covariance operator C of the process, defined by its matrix elements: for all f, g ∈ D(R), Cf, g = E X, g X, f .

For example, if C extends to a Hilbert-Schmidt operator, which we denote by

C ∈ L 2 , then A ∈ L 2 (R 2 ).
A convenient way of expressing the ambiguity function makes use of the translation operator T b , defined by T b f (t) = f (t -b) and the modulation operator E ω defined by

E ω f (t) = e iωt f (t): A f,g (τ, ξ) = T -τ 2 E -ξ 2 f, T τ 2 E ξ 2
g . This shows in particular that the ambiguity function A f (τ, ξ) only provide estimates for the spreading of the analyzed object in the joint time-frequency plane, but not on its localization in that space. Such an analysis is done by the Wigner-Ville (WV) distribution, defined by

DEFINITION 2.2 1. Let f ∈ L 2 (R). Its Wigner-Ville distribution is the func- tion of two real variables W f defined by W f (b, ω) = f (b + τ /2)f (b -τ /2)e -iωτ dτ (5)
2. Let {X t , t ∈ R} a second order time series. Then its Wigner-Ville distribution (or Wigner spectrum) is defined by

E X (b, ω) = E X b+τ /2 X b-τ /2 e -iωτ dτ (6)
REMARK 2.1 It is readily seen that the Wigner-Ville function and the ambiguity function are related via a symplectic Fourier transform

A(τ, ξ) = 1 2π W(b, ω)e -i(ξb-ωτ ) db dω , (7)
W(b, ω) = 1 2π A(τ, ξ)e i(ξb-ωτ ) dτ dξ .

The same holds true in the non-deterministic context. Therefore, the Wigner function is square-integrable as soon as the ambiguity function is so. Therefore,

W f ∈ L 2 (R 2 ) as soon as f ∈ L 2 (R). In fact, Lieb's estimates show that when f ∈ L 2 (R), W f ∈ L p (R 2 ) for all p ∈ [1, ∞]. REMARK 2.2
The Wigner-Ville distribution may also be defined when f is a distribution. Indeed, if f ∈ S ′ (R), it may be shown (see [START_REF] Folland | Harmonic Analysis of Phase Space[END_REF] for example) that W f ∈ S ′ (R 2 ). In fact, the Wigner-Ville distribution is nothing but the Weyl symbol of the operator of orthogonal projection onto f , which may be written:

P f g(t) = 1 2π||f || 2 W f t + b 2
, ω e iω(t-b) g(b)dbdω .

Properties of the Wigner-Ville distribution

The WV distributions enjoy remarkable properties. The first one reflects the orthogonality relations of the Ambiguity functions, and go under the name of Moyal's formulae: let

f, f ′ , g, g ′ ∈ L 2 (R). Then W f,g , W f ′ ,g ′ ∈ L 2 (R 2 ) and W f,g , W f ′ ,g ′ = 2π f, f ′ g ′ , g . (9)
The orthogonality relations of the Wigner-Ville coefficients are a direct consequence of the corresponding relations for ambiguity functions, and the symplectic Fourier transform formulas [START_REF] Beylkin | Fast wavelet transforms and numerical algorithms I[END_REF] and [START_REF]Radar and Sonar, Part I[END_REF].

A main property of the WV distribution is its covariance with respect to a certain number of simple transformations. Namely:

1. Translations: if g(t) = f (t -τ ), W g (b, ω) = W f (b -τ, ω). 2. Modulations: if g(t) = e iλt f (t), then W g (b, ω) = W f (b, ω -λ). 3. Rescalings: if g(t) = 1 √ a f t a , then W g (b, ω) = W f b a , aω .
4. Time-frequency rotations: Let θ be such that cos θ = 0, and set

g θ (t) = 1 √ cos θ 1 2π e i tan θ(t 2 +ξ 2 )/2 e itξ/ cos θ f (ξ)dξ . Then W g θ (b, ω) = W f (b cos θ + ω sin θ, -b sin θ + ω cos θ) .
More generally, the Wigner transform is covariant under a general group of transformations, called the metaplectic group. We refer to [START_REF] Folland | Harmonic Analysis of Phase Space[END_REF] for more details.

Examples

Time-frequency localization The examples that best illustrate the optimal localization properties of the WV distribution are the pure oscillations, i.e. the distributionof the form f (t) = e iλt . The WV transform of such an f (t) has to be defined as a two-dimensional distribution, and one easily shows that

W f (b, ω) = 2πδ(ω -λ).
More generally, such optimal localization properties are preserved by the simple transformations alluded to in the previous section. While the effect of translations, modulations and rescalings are easy to understand, let us stress the role played by time-frequency rotations. It is easily shown that the socalled linear chirps, i.e. the distributions with a linearly time-varying frequency e iλt+αt 2 /2 may be obtained by appropriate translation, modulation, rescaling and time-frequency rotation of a pure oscillation. Therefore, its WV distribution inherits the perfect localization properties from those of the pure oscillations, and one obtains W f (b, ω) = 2πδ(ω -(λ + αb)). An example of such a behavior is provided in Fig. 1, which represents the WV distribution of a linear chirp.

Interferences As a quadratic functional of the function f (t), the WV distribution yields interference terms. Namely, let

f 1 , f 2 ∈ L 2 (R), and let f = f 1 +f 2 .
Then one immediately sees that

W f (b, ω) = W f1 (b, ω) + W f2 (b, ω) + 2ℜW f1,f2 (b, ω)
Even in the case where both W f1 (b, ω) and W f2 (b, ω) are sharply localized in the (b, ω) plane, the "cross term" 2ℜW f1,f2 (b, ω) introduces an extra component in the WV transform of f , which degrades the resolution. An example is given in Fig. 2, where we display the WV distribution of the sum of three pulses.

The treatment of such interference terms has received a considerable attention in the signal processing literature during the past 10 years. Interferences terms are generally attenuated by appropriate smoothings of the WV distribution. Examples are provided by the following classes of time-frequency distributions:

1. The Cohen's class: let Π ∈ L 1 (R 2 ) be a fixed function. Associate with it the following time-frequency distribution

ρ f (b, ω) = Π(τ, ξ)W(b -τ, ω -ξ)dτ dξ . (10)
ρ f (b, ω) is expected to be smoother than the Wigner-Ville distribution, and in particular to contain much less important interference terms. Notice that since it is obtained from a WV distribution via a two-dimensional convolution (a translation invariant operation), ρ f (b, ω) inherits from W(b, ω) the time-frequency translations covariance.

2. The affine class: Given a function of two variables Π(b, ω), introduce

ρ f (b, a) = Π(τ, ξ)W b -τ a , aξ dτ dξ . (11)
ρ f (b, a) is a time-scale representation of f . With the same arguments as before, it may be shown that it inherits from the WV distributions its properties of covariance with respect to translations and rescalings.

The need for linear time-frequency transform The quadratic time-frequency distributions we just described have many interesting properties. We refer to [START_REF] Flandrin | Temps-Fréquence. Traité des Nouvelles Technologies[END_REF] for a detailed analysis of such properties in a signal processing perspective. However, they are often difficult to interpret in practical situations. For that reason, simpler alternatives such as linear time-frequency transforms are often prefered. We describe below the simplest two examples of such transforms. We first focus on the case of continuous transforms, and postpone the description of the discretization problem to a subsequent section.

Continuous Gabor transform (CGT)

Definitions The simplest localized version of Fourier analysis is provided by windowed Fourier transform, whose main idea is to localize the signal first by multiplying it by a smooth and localized window, and then perform a Fourier transform. More precisely, the construction goes as follows. Start from a function g ∈ L 2 (R), such that ||g|| = 0, and associate with it the following family of Gaborlets

g (b,ω) (t) = e iω(t-b) g(t -b) . ( 12
)
DEFINITION 2.3 Let g ∈ L 2 (R) be a window. The continuous Gabor transform of a finite-energy signal f ∈ L 2 (R) is defined by the integral transform G f (b, ω) = f, g (b,ω) = f (t) g(t -b)e -iω(t-b) dt . (13)
Gaborlets yield decomposition formulas for functions in L 2 (R), as follows.

THEOREM 2.1 Let g ∈ L 2 (R) be a non trivial window (i.e. ||g|| = 0.) Then every f ∈ L 2 (R) admits the decomposition

f (t) = 1 2π||g|| 2 ∞ -∞ ∞ -∞ G f (b, ω)g (b,ω) (t)dbdω , (14)
where equality holds in the weak L 2 (R) sense.

In other words, the mapping

L 2 (R) ∋ f ֒→ 1 ||g|| √ 2π G f ∈ L 2 (R 2
) is an isometry between L 2 (R) and L 2 (R 2 ), and the inverse mapping is provided by the adjoint mapping.

The case of random time series DEFINITION 2.4 Let {X t , t ∈ R} be a mean zero second order random time series, and let g ∈ L 2 (R) be a window. The CGT of X t is the random timefrequency series defined by

G X (b, ω) = X, g (b,ω) (15)
By definition, the covariance operator C of the time series is defined by its matrix elements Cf, g = E { X, g f, X }. Obviously,

E G f (b, ω)G f (b ′ , ω ′ ) = Cg (b ′ ,ω ′ ) , g (b,ω) (16)
The case where the time series {X t , t ∈ R} under consideration is (second order) stationary is particularly interesting. By definition, the covariance operator is in such a case a convolution operator, and one immediately sees that

E G f (b, ω)G f (b ′ , ω) = 1 2π e iη(b-b ′ ) E(η) |ĝ(η -ω)| 2 dη ,
where E(η) stands for the spectral density of the time series. However, the main interest of the CGT lies in its potential for handling non stationary situations (in the stationary case, the covariance operator is a convolution operator, which is perfectly handled by Fourier methods). Particularly interesting is the case of the so-called locally stationary time series, which are basically random time series whose covariance operator is "almost diagonal" in an appropriate Gabor representation. Such a situation has been discussed by various authors in various contexts (see e.g. [START_REF] Kozek | Spectral Estimation in Non-Stationary Environments[END_REF][START_REF] Mallat | Adaptive covariance estimation of locally stationary processes[END_REF][START_REF] Carmona | Practical Time-Frequency Analysis: Gabor and wavelet transforms, with an implementation in S[END_REF]33]).

examples The CGT and similar tools have been quite popular in the speech processing literature, because of its capability of handling the so-called locally harmonic signals, namely signals which may be modeled in the form

f (t) = K k=1 A k (t)e iφ h (t) , (17)
where the local amplitudes a k (t) and the local frequencies ω k (t) = φ ′ k (t) are assumed to be slowly varying. Such signals are called locally harmonic when, in addition, the local frequencies are close to be integer multiples of a fundamental frequency (the pitch frequency: ω k (t) ≈ kω 1 (t). The CGT of such signals have been studied by several authors. For example we have the following estimates

LEMMA 2.1 Let g ∈ L 2 (R) be a window such that |tg(t)|dt < ∞ and t 2 |g(t)|dt < ∞. Let f (t) = A(t)e iφ(t)
be such that A, φ ∈ C 2 (R), and let G f (b, ω) denote the continuous Gabor coefficients of f with respect to the window g. Then

G f (b, ω) = A(b)e iφ(b) ĝ(φ ′ (b) -ω) [1 + R(b, ω)] , (18)
and the remainder R(b, ω) is bounded as

|R(b, ω)| ≤ K 1 |A ′ (b)| |A(b)| + K 2 sup u |A ′′ (u)| |A(b)| + K 3 sup u |φ ′′ (u)| |A(b)| .
A proof of a statement close to that one may be found in [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF]. See also [START_REF] Delprat | Asymptotic wavelet and Gabor analysis: extraction of instantaneous frequencies[END_REF][START_REF] Ryan | Introduction to Continuous Wavelet Analysis[END_REF][START_REF] Carmona | Practical Time-Frequency Analysis: Gabor and wavelet transforms, with an implementation in S[END_REF] for similar approaches.

The meaning of such statements is the following: as soon as the amplitudes A(t) and frequencies φ ′ (t) are slowly varying enough, one may write

G f (b, ω) ≈ A(b)e iφ(b) ĝ(φ ′ (b) -ω) .
Assume now that g(t) is a smooth function, located near the origin t = 0 (typically a Gaussian function.) The latter expression shows that for each value of the time variable b, the modulus of the CGT attains its maximum on a curve (the so-called ridge) of equation ω = φ ′ (b), i.e. describing the instantaneous frequency of the function f (t).

We refer to [START_REF] Carmona | Practical Time-Frequency Analysis: Gabor and wavelet transforms, with an implementation in S[END_REF] for a large number of numerical illustrations of such a situation. We limit our illustrations to an example of speech signal, in which several such components are present. Thanks to the linearity of the CGT, the CGT of a signal of the type ( 17) is of the form

G f (b, ω) ≈ K k=1 A k (b)e iφ k (b) ĝ(φ ′ k (b) -ω) ,
and as soon as g is such that the various φ ′ k (b) -φ ′ ℓ (b) for ℓ = k are small enough, the CGT of such a signal localizes itself near K different ridges. An example of such a behavior is given in Fig. 3, where we display 625 milliseconds of speech signal: /How are you ?/ (top) and the modulus of the CGT. Observe the localisation near the ridges.

Continuous Wavelet Transform (CWT)

An alternative to Gabor transform was proposed more recently by Grossmann and Morlet [START_REF] Grossmann | Decomposition of Hardy functions into square integrable wavelets of constant shape[END_REF]. The main idea was to improve the time resolution of Gabor transform, by changing the rule for generating the "basis functions". This may be done by replacing the modulation operation used to generate Gaborlets by a scaling operation.

Wavelet transform Let ψ ∈ L 1 (R) ∩ L 2 (R) be a fixed function (in fact, it is sufficient to assume ψ ∈ L 1 (R), but for convenience also assume that ψ ∈ L 2 (R). This extra assumption ensures the boundedness of the wavelet transform.) From now on it will be called the analyzing wavelet. It is also sometimes called the mother wavelet of the analysis. The corresponding family of wavelets is the family {ψ (b,a) ; b ∈ R, a ∈ R * + } of shifted and scaled copies of ψ defined as follows. If b ∈ R and a ∈ R * + we set:

ψ (b,a) (t) = 1 a ψ t -b a , t ∈ R . (19)
The wavelet ψ (b,a) can be viewed as a copy of the original wavelet ψ rescaled by a and centered around the "time" b. Given an analyzing wavelet ψ, the associated continuous wavelet transform is defined as follows

DEFINITION 2.5 Let ψ ∈ L 1 (R) ∩ L 2 (R)
be an analyzing wavelet. The continuous wavelet transform (CWT for short) of a finite-energy signal f (t) is defined by the integral:

T f (b, a) = f, ψ (b,a) = 1 a f (t)ψ t -b a dt . (20)
Like Gaborlets, wavelets may form complete sets of functions in L 2 (R), and we have in particular

THEOREM 2.2 Let ψ ∈ L 1 (R) ∩ L 2 (R)
, be such that the number c ψ defined by:

c ψ = ∞ 0 | ψ(aξ)| 2 da a ( 21
)
is finite, nonzero and independent of ξ ∈ R. Then every f ∈ L 2 (R) admits the decomposition

f (t) = 1 c ψ ∞ -∞ ∞ 0 T f (b, a)ψ (b,a) (t) da a db , (22)
where the convergence holds in the strong L 2 (R) sense.

In particular, we also have "energy conservation":

if f ∈ L 2 (R), then T f ∈ L 2 (R × R * + , db da a ), and
||T f || 2 = c ψ ||f || 2 .
Notice that in [START_REF] Duflo | On the regular representation of a nonunimodular locally compact group[END_REF], the constant c ψ can only depend on the sign of ξ ∈ R, therefore assuming independence wrt ξ is a simple symmetry assumption. The fact that 0 < c ψ < ∞ implies that ψ(0) = 0, so that the wavelet ψ(t) has to oscillate enough to be of zero mean.

The CWT of random time series may be defined in a way similar to the CGT.

Examples The CWT has a behavior similar to that of the CGT in many respects. The main difference lies in the fact that wavelets are extremely precise at small scales (where they lose frequency resolution), and more frequency localized at large scales (where time resolution is lost). A main application of this fact is the analysis of regularity (see below). A visualization of this effect may be found in Fig. 3 (bottom), where we display the modulus of the CWT of the speech signal shown at the top of the figure. As may be seen, at large scales, the wavelets have a sufficient frequency resolution to analyze carefully the first harmonics (namely, the pitch frequency and the first harmonic). For smaller scales, frequency resolution is lost, and the same wavelet is unable to "separate" several harmonic components. This results in interferences between the harmonic components, which yield the oscillations of the modulus in the b direction that appear on the image. We refer to [START_REF] Carmona | Practical Time-Frequency Analysis: Gabor and wavelet transforms, with an implementation in S[END_REF] for a more detailed analysis of such applications.

Linear transforms as approximations

Interestingly enough, there is a strong connection between the Wigner-Ville distributions and the linear decompositions we just reviewed (or more precisely their squared modulus). We express this in the following

PROPOSITION 2.2 Let f ∈ L 2 (R)
, and let T f , G f and E f denote respectively its CWT, CGT and WV transforms. Then the following two properties are true:

1. The squared modulus |G f (b, ω)| 2 of the CGT is a smoothed version of the WV distribution W f of f : |G f (b, ω)| 2 = W f (τ, ξ)W g (τ -b, ξ -ω)dτ dξ (23)
2. The squared modulus |T f (b, a)| 2 of the CWT is an affine smoothing of the WV distribution W f of f :

|T f (b, a)| 2 = W f (τ, ξ)W ψ τ -b a , aξ dτ dξ (24)
These results follow from an easy calculation. See e.g. [START_REF] Flandrin | Temps-Fréquence. Traité des Nouvelles Technologies[END_REF] for more details.

Littlewood-Paley transform

We close this section with a short account of a simplified version of wavelets decomposition, known for a long time under the name of Littlewood-Paley decomposition. Littlewood-Paley decomposition also provide a first step towards the discretization of the continuous wavelet transforms, by a discretization of the scale variable only. Let ψ ∈ L 1 (R) ∩ L 2 (R) be such that 0 < c ψ < ∞, and let us assume that φ is such that φ(aξ)

2 = ∞ a ψ(uξ) 2 du u . ( 25
)
Assuming an obvious definition for φ (b,a) (t), we set

S f (b, a) = f, φ (b,a) (26)
for all f ∈ L 2 (R). Then equation ( 22) may be replaced by

f (t) = 1 c ψ lim a→0 S f (b, a)φ (b,a) (t)db = 1 c ψ S f (b, a)φ (b,a) (t)db + ∞ -∞ a 0 T f (b, u)ψ (b,u) (t) du u db .
(27) The function φ is called a (bilinear) scaling function associated with the wavelet ψ. We refer to [START_REF] Duval-Destin | Continuous wavelet decompositions: multiresolution and contrast analysis[END_REF] for more details.

It is important to emphasize the difference between wavelet coefficients T f and scaling function coefficients S f . While the coefficients S f (b, a) provide an approximation of the function f at scale a, T f (b, a) give details at scale a, i.e., essentially differences between two consecutive approximations.

In fact, the integral with respect to the scale variable in the continuous wavelet decompositions may easily be replaced with a discrete sum, by the following trick. Let Ψ(t) be such that

Ψ(ω) 2 = 1 1/2 ψ(aω) 2 da a (28)
The, setting for all f ∈ L 2 (R)

Σ j (f )(t) = 1 c ψ S f (b, 2 -j )φ (b,2 -j ) (t)db (29) D f (b, a) = f, Ψ (b,a) (30)
∆ j (f )(t) = 1 c ψ D f (b, 2 -j )Ψ (b,2 -j ) (t)db (31)
we obtain the following decomposition: for all f ∈ L 2 (R),

f = Σ j0 (f ) + ∞ j=j0+1 ∆ j (f ) = ∞ j=-∞ ∆ j (f ) (32)
Such decompositions are called Littlewood-Paley decompositions, and play a great role in several contexts, for example in the study of partial differential equations. They are also extremely useful in the context of signal analysis, for they provide a version of wavelet analysis in which the translation invariance properties of the continuous transform are preserved. Littlewood-Paley decompositions may also be thought of as a first step towards multiresolution decompositions, which will be described in Section 4

Wavelets and the characterization of regularity

The wavelet transform is often compared to a mathematical microscope, because of its capability of zooming in. In fact, another crucial property of wavelets is their blindness to certain "regular" behaviors. Indeed, we have seen already that a wavelet has necessarily zero integral. Hence, two functions which differ only by a constant have the same wavelet coefficients: wavelets are blind to constants. If necessary, one may make them blind to higher order polynomial behavior, by imposing vanishing moments conditions, of the type

t m ψ(t)dt = 0 , m = 0, 1, . . . M .
This is best illustrated by the study of Hölder regularity properties of functions, which we sketch below.

Global regularity properties

We start with an example of characterization of some global regularity properties. Let us first recall some definitions. For 0 < α < 1, let C α denote the Banach space of functions f satisfying

|f (t) -f (s)| ≤ C|t -s| α , t, s ∈ R, (33)
for some constant 0 < C < ∞. Such functions are said to be Hölder continuous of order α. Then we have

THEOREM 3.1 Let ψ ∈ L 1 (R) be such that c α = |t| α |ψ(t)| dt < ∞ , and sup |δ|<1 |ψ ′ (u + δ)| du < ∞ . ( 34
)
Then f ∈ C α if and only if

|T f (b, a)| ≤ C ′ a α , a ∈ R, (35)
for some finite positive constant C ′ independent of b.

We will not give a complete proof of this result here (see, e.g., [START_REF] Frazier | Littlewood-Paley Theory and the Classification of Function Spaces[END_REF]), but in order to illustrate the role of the vanishing moments, let us sketch the proof of the first half. Let ψ be a wavelet satisfying the above properties, and let f ∈ C α . Then since ψ(t)dt = 0, we have

T f (b, a) = 1 a [f (t) -f (b)] ψ t-b a dt, and
|T f (b, a)| ≤ C a |t -b| α ψ t -b a dt ≤ Cc α a α .
which proves the first part of the theorem. The converse is proved in a similar way and makes use of the second assumption on the wavelet.

REMARK 3.1 If one associates to any f ∈ C α the infimum f α of the constants C such that (33) holds, it may be proved that • α provides C α with a norm, which makes it a Banach space (notice that functions in this Banach space, are only defined modulo additive constants). Then Theorem 3.1 expresses the fact that for a suitably chosen wavelet, the number inf a -α |T f (b, a)| defines an equivalent norm on C α . This is one of the simplest examples of the characterization of function spaces by mean of the wavelet transform. Elaborating on such arguments leads to the characterization of wider classes of functional spaces (see, e.g., [START_REF] Meyer | Ondelettes et opérateurs,I: Ondelettes; II: Opérateurs de Calderón-Zygmund[END_REF] or [START_REF] Frazier | Littlewood-Paley Theory and the Classification of Function Spaces[END_REF] for a review).

The theorem above only addresses the case of Hölder spaces of exponent 0 < α < 1. More generally, let us assume that n < α < n + 1. Then the function f belongs to the Hölder space C α if for any t there exists a polynomial P n of order n and a constant C such that for |h| small enough,

|f (t) -P n (h)| ≤ C|h| α . ( 36
)
Equivalently, f is C α if its derivative of order n is C α-n . To adapt to this new situation, extra assumptions on the wavelet are needed. To see this, let us consider the following toy example. Assume that f has M -1 continuous derivatives in a neighborhood U of a point t = t 0 . Then for every t ∈ U we may write

f (t) = f (t 0 ) + (t -t 0 )f ′ (t 0 ) + • • • + 1 (M -1)! (t -t 0 ) M -1 f (M -1) (t 0 ) + r(t) .
If the wavelet ψ (assumed to have compact support for the sake of the present argument) has M vanishing moments, then for sufficiently small values of the scale parameter a (i.e., such that supp(ψ (b,a) ) ⊂ U for some b), we have

T f (b, a) = f, ψ (b,a) = r, ψ (b,a) .
Then we can say that, at least in some sense, the wavelet does not "see" the regular behavior of f (t) near t = t 0 and "focuses" on the potentially singular part r(t). If in addition we assume that f (t) is C M near t = t 0 , then we have r(t) = (t -t 0 ) M ρ(t), with |ρ(t)| ≤ sup |f (M ) |/M !, and we immediately conclude that T f (b, a) ∼ O(a M ) as a → 0 .

More generally, Theorem 3.1 is still valid if one assumes in addition that the wavelet ψ has n vanishing moments (and a corresponding generalization of condition [START_REF] Hogan | Extensions of the Heisenberg group by dilations and frames[END_REF]). The proofs go along the same lines, except that instead of using the zero integral property of ψ to subtract f (b) from f (t), one uses all the vanishing moment assumption to subtract P n (t -b) from f (t) inside the integral.

Elaborating on such arguments leads to the characterization of Hölder spaces of functions via wavelet transform (see, e.g., [START_REF] Meyer | Ondelettes et opérateurs,I: Ondelettes; II: Opérateurs de Calderón-Zygmund[END_REF][START_REF] Frazier | Littlewood-Paley Theory and the Classification of Function Spaces[END_REF]) and the whole body of work on the characterization of singularities of a given signal [START_REF] Mallat | Singularities detection and processing with wavelets[END_REF]. More elaborate arguments also yield characterizations of more general functional classes, such as Sobolev, Besov and Triebel classes.

Local regularity

The previous results are global in the sense that they require uniform regularity throughout the real line. Another classical example involves local regularity. For a given point t we say that the function f is locally Hölder continuous of order α at t if it satisfies

|f (t) -f (s)| ≤ C|t -s| α ( 37
)
for some constant C and for |t -s| small enough. Then we have the "local" counterpart of the previous result:

THEOREM 3.2 Let ψ ∈ L 1 (
R) be a wavelet with M + 1 vanishing moments such that for any m, there is a constant C m such that ψ (k) (t) ≤ Cm 1+|t| m . 1. Let the function f be locally Hölder continuous of order α at t 0 . Then

|T f (b, a)| ≤ C ′ (a α + |b -t 0 | α ) ( 38
)
for some constant C ′ .

2. Conversely, Let α < M be a non-integer number, and assume that there exist numbers α ′ < α and A such that for all a, b

|T (b, a)| ≤ Aa α 1 + t 0 -b a α ′ . ( 39
)
Then f is locally Hölder continuous of order α at t 0 .

The proof of the theorem uses the ingredients we have developed above. Let us start with the first part, and consider T (b, a)

= 1 a (f (t) -P (t)) ψ t-b a dt. |T (b, a)| ≤ A a |t -t 0 | α ψ t -b a dt ≤ 2 α A a (|t -b| α + |b -t 0 | α ) ψ t -b a dt ≤ 2 α A a α |u| α |ψ(u)|du + |t 0 -b| α ||ψ|| 1 ,
which proves the first assertion. For the converse assertion, one first needs to exhibit the polynomial P (t). Let us consider the Littlewood-Paley blocks

∆ j (t) = ∞ -∞ 2 j+1 2 j T (b, a)ψ (b,a) (t) da a db
According to the assumptions made on ψ, each block may be bounded as

|∆ j (t)| ≤ K2 jα ∞ -∞ 1 + 2 α ′ |u| α ′ + 2 α ′ t-t0 2 j α ′ 1 + |u| m du .
Taking m > α ′ + 1 yields the following estimate:

LEMMA 3.1 |∆ j (t)| ≤ K2 jα 1 + t -t 0 2 j α ′ . ( 40
)
In the same way, we obtain the estimates LEMMA 3.2 For all k smaller than the integer part of α,

∆ (k) j (t) ≤ K2 j(α-k) 1 + t -t 0 2 j α ′ . ( 41
)
Let us now consider the truncated Taylor expansion

P t0 (t) = ⌊α⌋ k=0   ∞ j=-∞ ∆ (k) j (t 0 )   (t -t 0 ) k k!
The behavior of this sum as j → -∞ is easily controlled thanks to Lemma 3.2.

In addition, the fact that

|T (b, a)| ≤ ||f || ||ψ||/ √ a implies |∆ (k) j (t 0 )| ≤ K2 -j(k+ 1 2 )
, which takes care of the limit j → ∞. Then, P t0 (t) is well defined. Consider now f (t) -P t0 (t). Let J be such that 2 J-1 ≤ |t -t 0 | < 2 J . It follows from Taylor's formula that

∞ j=J   ∆ j (t)- ⌊α⌋ k=0 ∆ (k) j (t 0 ) (t -t 0 ) k k!   ≤ |t -t 0 | ⌊α⌋+1 (⌊α⌋ + 1)! sup s∈[t,t0] ∆ (⌊α⌋+1) j (s) ≤ K|t-t 0 | α
In addition, the sum from -∞ to J -1 is bounded by

K J-1 j=-∞   2 jα 1 + t -t 0 2 j α ′ + ⌊α⌋ k=0 2 j(α-k) |t -t 0 | k k!   ≤ K   2 Jα + 2 J(α-α ′ ) |t -t 0 | α ′ + ⌊α⌋ k=0 2 J(α-k) |t -t 0 | k k!   ≤ K   2 α |t -t 0 | α + 2 α-α ′ |t -t 0 | α + ⌊α⌋ k=0 2 α-k k! |t -t 0 | α   ≤ K|t -t 0 | α ,
where K is a constant, which may change from line to line. This sketches the proof. More details may be found for example in [START_REF] Jaffard | Wavelet Methods for Pointwise Regularity and Local Oscillations of Functions[END_REF] and [START_REF] Mallat | Singularities detection and processing with wavelets[END_REF] which contain reviews especially tailored to mathematicians and signal analysts, respectively. REMARK 3.2 It is instructive to examine the meaning of conditions of the type [START_REF] Kozek | Spectral Estimation in Non-Stationary Environments[END_REF]. Let us consider the cone in the (b, a) half-plane defined by the condition |b -t 0 | < a. Within this cone we have |T f (b, a)| = O(a α ) as a → 0. Outside the cone, the behavior is governed by the distance of b to the point t 0 . These two behaviors are generally different, and have to be studied independently. However, it is shown in [START_REF] Mallat | Singularities detection and processing with wavelets[END_REF] that non-oscillating singularities may be characterized by the behavior of their wavelet transform within the cone.

REMARK 3.3
It is also shown in [START_REF] Mallat | Singularities detection and processing with wavelets[END_REF] that rapidly oscillating singularities cannot be characterized by the behavior of their wavelet transform in the cone. A typical example is given by the function f (t) = sin 1/t whose instantaneous frequency tends to infinity as x → 0. The wavelet transform modulus is maximum on a curve of equation b = Ka 2 for some constant K depending only on the wavelet, and this curve is not in the cone. In such a case, the oscillations have to be analyzed carefully (see, e.g., [START_REF] Jaffard | Wavelet Methods for Pointwise Regularity and Local Oscillations of Functions[END_REF].)

REMARK 3.4
We close this discussion with a more "practical" remark. Most often, data are only available in the form of discrete signals, i.e. only for discrete values of the variable (with a fixed resolution given by the sampling frequency). Therefore, the notions of singularities and Hölder exponents are strictly speaking meaningless (as is the limit a → 0.) Nevertheless, one can say that the behavior of the wavelet coefficients across scales provides a good way of describing the regularity of functions whose samples coincide with the observations up to a given resolution.

Multiresolution wavelet analysis

Let us now turn to an approach which is completely different, namely the multiresolution approach. This approach is based upon a general methodology (the multiresolution analysis, which automatically generates orthonormal bases of wavelets. Remarkably enough, such an approach turned out to be extremely close to image processing methods (subband filtering methods), which led to fast pyramidal algorithms for wavelet decompositions. We shall not go into the details of the theory, which received a considerable attention in the recent years, and refer to [15,[START_REF] Vetterli | Wavelets and SubBand Coding[END_REF][START_REF] Wickerhauser | Adapted Wavelet Analysis, from Theory to Software[END_REF] for details. We rather sketch the construction and give a few consequences and examples.

Multiresolution analysis and associated wavelet basis

Let us start by introducing the notion of multiresolution analysis.

DEFINITION 4.1 A multiresolution analysis (MRA) of L 2 (R) is a collection of nested closed subspaces V j ⊂ L 2 (R) • • • ⊂ V -2 ⊂ V -1 ⊂ V 0 ⊂ V 1 ⊂ V 2 . . . (42)
such that the following properties hold:

1. ∪V j = L 2 (R) and ∩V j = {0}. 2. If f ∈ V 0 , then f ( • -k) ∈ V 0 for all k ∈ Z; f ∈ V j if and only if f ( • /2) ∈ V j-1 .
3. There exists a function χ ∈ V 0 such that the collection of the integer translates χ(

• -k) for k ∈ Z is a Riesz basis of V 0 .
The function χ is called a scaling function, and the Riesz basis may be orthonormalized by Gram's procedure by setting

φ(ω) = χ(ω) k | χ(ω + 2πk)| 2 . ( 43
)
COROLLARY 4.1 With the same notation as before, the collection {φ( • -k); k ∈ Z} is an orthonormal basis of V 0 .

The function φ(t) is also called a scaling function. An immediate consequence of the definition of φ in ( 43) is that the Fourier transform satisfies

k | φ(ω + 2πk)| 2 = 1 . (44)
From now on, we shall essentially work with the scaling function φ associated with orthonormal basis. The inclusion of the V j spaces shows that φ(t) (as well as χ(t)) may be expressed as a linear combination of the functions 2 1/2 φ(2t-k), which form an orthonormal basis of V 1 . The consequence is the so-called twoscale difference equation (or refinement equation)

φ(t) = √ 2 k h k φ(2t + k) . (45)
Let now the 2π-periodic function H(ω) be defined by

H(ω) = 1 √ 2 k h k e ikω , (46)
Clearly, we have φ(2ω) = H(ω) φ(ω). Let W j denote the orthogonal complement of V j in V j+1 . The cornerstone of the multiresolution theory is the existence of a function ψ ∈ W 0 such that the collection {ψ( • -k), k ∈ Z} is an orthonormal basis of W 0 . More precisely, if we define the 2π-periodic function G(ω) by

G(ω) = e iω H(ω + π), (47)
and if we denote by 2 -1/2 g k its Fourier coefficients, i.e., if we set

G(ω) = 1 √ 2 k g k e ikω , g k = -(-1) k h 1-k , k ∈ Z , (48)
then we have

THEOREM 4.1 If ψ(t) is defined by ψ(t) = √ 2 k g k φ(2t + k), t ∈ R, (49)
then the collection {ψ( • -k); k ∈ Z} is an orthonormal basis of W 0 .

See, e.g., [START_REF] Meyer | Ondelettes et opérateurs,I: Ondelettes; II: Opérateurs de Calderón-Zygmund[END_REF] for a proof. Notice that ψ constructed in this way is in fact an analyzing wavelet in the sense of the continuous wavelet transform introduced earlier. Formula (49) is a second example of refinement equation. Set

ψ jk (t) = 2 j/2 ψ 2 j t -k , φ jk (t) = 2 j/2 φ 2 j t -k . (50)
By definition, it follows that for a given j ∈ Z, the collection {φ jk ; k ∈ Z} (resp. {ψ jk ; k ∈ Z}) is an orthonormal basis of V j (resp. W j ), and we have COROLLARY 4.2 With the same notation as before, the collection of wavelets {ψ jk ; j, k ∈ Z} is an orthonormal basis of L 2 (R).

Quadrature mirror filters

The remarkable property of H(ω) and G(ω) is that they satisfy the so-called Quadrature Mirror Filter (QMF) condition, which is a direct consequence of equations ( 44) and ( 45):

LEMMA 4.1
The periodic functions H(ω) and G(ω) defined in [START_REF] Meyer | Ondelettes et opérateurs,I: Ondelettes; II: Opérateurs de Calderón-Zygmund[END_REF] and (47)

satisfy |H(ω)| 2 + |G(ω)| 2 = 1, ω ∈ R. (51)
Let now f ∈ L 2 (R) be any square-integrable function, and set

s jk = f, φ jk , and d jk = f, ψ jk (52)
Then, an immediate consequence of the refinement equations and Lemma 4.1 is

LEMMA 4.2
1. The coefficients {d jk,j,j∈Z } and {s jk,j,j∈Z } may be computed recursively via the following pyramidal algorithm

s j-1 k = ℓ h 2k-ℓ s jℓ , and d j-1 k = ℓ g 2k-ℓ s jℓ . ( 53
)
2. The finer resolution coefficients may be obtained from the coarser ones by

s j+1 ℓ = k (h 2k-ℓ s jℓ + g 2k-ℓ d jℓ ) . (54)
These two equations are of fundamental practical importance, for they automatically yield fast algorithms for discrete wavelet transform and inverse transform.

Let us also mention that such pyramidal algorithms have been generalized to yield new orthonormal bases of L 2 (R), known as wavelet packet bases. The interested reader may refer to [START_REF] Wickerhauser | Adapted Wavelet Analysis, from Theory to Software[END_REF] for a detailed account of such generalizations.

Discrete Littlewood-Paley transform

Let us finally quote for the record a variant of discrete wavelet transform, which is essentially a discrete version of the Littlewood-Paley transform we described above. The goal of such a transform is to restore a kind of (discrete) translation invariance. Let us consider a multiresolution analysis, with scaling function φ(t) and wavelet ψ(t), and introduce the following family of functions, defined by

ψ j k (t) = 2 j/2 ψ 2 j (t -k) , φ j k (t) = 2 j/2 φ 2 j (t -k) . (55)
This family forms an overcomplete set. For any f ∈ L 2 (R), define

D j k = f, ψ j k , and S j k = f, φ j k . (56)
An immediate calculation yields the following pyramidal algorithm for computing such coefficients:

S j-1 k = ℓ h k-2 j-1 ℓ S j ℓ , and D j-1 k = ℓ g k-2 j-1 ℓ S j ℓ . (57)
These coefficients yield an alternate wavelet representation for the function f (t), suitable for problems in which translation invariance is an important issue. See for example [START_REF] Mallat | Singularities detection and processing with wavelets[END_REF] for some applications in the image processing context.

Implementing Symmetries

As we have seen already, the continuous (linear or quadratic) time-frequency transforms described at the beginning of this paper possess important built-in symmetry properties. Such covariance properties, which are crucial for practical purposes (in particular for signal processing applications) are in fact intrinsic to the transforms, and even characterize them. This fact is conveniently described by means of group-theoretical methods. Indeed, the set of simple transformations used to generate the wavelets from a single one in general inherits the structure of a group G (as is the case for instance for translations, modulations or dilations). The deep connection between the usual wavelet decompositions, coherent states theory and the theory of square-integrable group representations was emphasized by H. Moscovici and A. Verona [START_REF] Moscovici | Coherent states representations of nilpotent Lie groups[END_REF][START_REF] Moscovici | Coherent states and square-integrable representations[END_REF] and A. Grossmann, J. Morlet and T. Paul independently [START_REF] Grossmann | Transforms associated with square integrable group representations I[END_REF]. As a result, all such wavelets are generated from a representation of a group of simple transformations G, in such a way that the representation is unitarily equivalent to a subrepresentation of the regular representation, i.e. a representation of the group G onto L 2 (G).

Square-integrable group representations

The connexion between time-frequency representation theorems and squareintegrable group representations was realized by A. Grossmann, J. Morlet and T. Paul in [START_REF] Grossmann | Transforms associated with square integrable group representations I[END_REF]. We briefly sketch the construction. DEFINITION 5.1 Let G be a separable locally compact Lie group, and let π be a unitary strongly continuous representation of G on the Hilbert space H. π is said to be square-integrable (or to belong to the discrete series of G) if π is irreducible, and if there exists at least a vector v ∈ H such that

0 < G | π(g)v, v | 2 dµ(g) < ∞ . ( 58
)
Such a vector is said to be admissible.

Square-integrable group representations have been extensively studied in the literature, in particular for compact groups, locally compact unimodular groups and non-unimodular locally compact groups. The results may be summarized in the following theorem, due to Duflo and Moore: THEOREM 5.1 Let π be a square-integrable strongly continuous unitary representation of the locally compact group G on H. Then there exists a positive self-adjoint operator C such that for any admissible vectors v 1 , v 2 ∈ H and for any

u 1 , u 2 ∈ H G u 1 , π(g).v 1 π(g).v 2 , u 2 dµ(g) = C 1/2 v 2 , C 1/2 v 1 u 1 , u 2 (59)
Moreover, the set of admissible vectors coincides with the domain of C.

Let λ : G → U (L 2 (G)) be the left-regular representation of G: if f ∈ L 2 (G), (λ(h)f ) (g) = f (h -1 g) . (60)
Theorem 5.1 shows that a representation π of G is square integrable if and only if is is unitarily equivalent to a subrepresentation of the left-regular representation λ. The corresponding intertwiners can be realized as follows. If v is an admissible vector in H, and v ′ ∈ H, introduce the Schur coefficients, i.e the matrix coefficients of elements of G:

c v,v ′ (g) = v ′ , π(g)v , g ∈ G (61)
Let T be the left transform [START_REF] Grossmann | Transforms associated with square integrable group representations I[END_REF], i.e. the map defined by

T : u ∈ H → T u = c v,u (•) ∈ L 2 (G) (62)
T realizes the intertwinning between π and λ:

T • π = λ • T (63)
The idea of Grossmann, Morlet and Paul was to use (59) and (62) for the analysis of functions, in the case where H is a function space. This was the starting point of many applications, especially in a signal analysis context. The left transform T is used to obtain another representation of functions, and (63) expresses the covariance of the transform. Notice that the continuity assumption is fulfilled by construction.

The classical examples

Gabor functions We start with the case of the so-called canonical coherent states, generated from the n-dimensional Weyl-Heisenberg group

G W H = R 2n × S 1 , with group operation (q, p, ϕ) • (q ′ , p ′ , ϕ ′ ) = (q + q ′ , p + p ′ , ϕ + ϕ ′ + p.q ′ [mod 2π]) (64)
It follows from the Stone-Von-Neumann theorem that any irreducible unitary representations of G W H is unitarily equivalent to one of the following form (see [START_REF] Mackey | Theory of unitary group representations[END_REF][START_REF] Schempp | Harmonic analysis on the Heisenberg nilpotent Lie group[END_REF]).

If f ∈ L 2 (R n ) π(q, p, ϕ) • f (x) = e iµ(ϕ-p.(x-q)) f (x -q)
for some µ ∈ Z * . Let us focus on the one-dimensional situation for the sake of simplicity, and assume µ = 1. The representation is square-integrable, and the corresponding left transform is the Gabor transform described in Definition 2.3, up to a phase factor e iϕ . More precisely, the Gaborlets read g (b,ω) = π(b, ω, 0)g .

This can be interpreted as a square-integrable projective representation of the abelian group R 2n . Notice that in such a case, one does not have full covariance by time-frequency shifts, precisely because of this phase factor (twisted covariance.)

One-dimensional wavelets Let G af f = R × R * + denote
(notice the slight change in the normalization.) This representation is squareintegrable, and the corresponding left transform is nothing but the affine wavelet transform described in Definition 2.5. The covariance equation (63) simply expresses that the affine wavelet transform of a dilated and shifted copy of a function f (t) coincides with a dilated and shifted version of the wavelet transform of f (t). Finally, the admissibility of a vector ψ ∈ L 2 (R) reduces to [START_REF] Duflo | On the regular representation of a nonunimodular locally compact group[END_REF].

In a signal analysis context, the translation parameter is interpreted as a position (or time) parameter, and the scale parameter a as a frequency parameter (more precisely the inverse of a frequency parameter). One is then led to a time-frequency representation.

N -dimensional wavelets To generalize the previous construction to the ndimensional case, one faces an irreducibility problem. The n-dimensional affine group G af f = R n × R * + does not act irreducibly on L 2 (R n ). To restore irreducibility, R. Murenzi [START_REF] Murenzi | Ondelettes multidimensionnelles et application à l'analyse d'images[END_REF] proposed to consider the semi-direct product of the affine group by SO(n). The resulting group, denoted by IG(n) then yields wavelets of the form

ψ (b,a,r) (x) = a -n/2 ψ r -1 x -b a , (66)
and the following resolution of the identity: for all f ∈ L 2 (R n ) f = C ψ f, ψ (b,a,r) ψ (b,a,r) db a n da a dθ .

(67)

Generalizations

Several extensions and generalizations have been developed. A first way of generalization concerns the study of more general classes of groups possessing square-integrable representations. Let us mention for example the works of Bernier and Taylor [START_REF] Bernier | Wavelets and square integrable representations[END_REF], Führ [START_REF] Führ | Wavelet frames and admissibilityin higher dimensions[END_REF] and Aniello et al [START_REF] Aniello | Wavelet transforms and discrete frames associated to semidirect products[END_REF] on semi-direct products of the form H × R n , where H is a closed subgroup of GL(n, R), with group law (h, b)(h ′ , b ′ ) = (hh ′ , b + hb ′ ). The above mentioned authors have in particular given necessary and sufficient conditions on the group H which ensure that the representation under consideration is square integrable. The application of the general theory to the case of discrete groups has also been studied. We refer to [START_REF] Flornes | Wavelets on Discrete Fields[END_REF][START_REF] Antoine | An algebraic approach to discrete dilations[END_REF] for more details.

In more general situations, i.e. for more general groups, the general theory does not apply, since the representation under consideration is not squareintegrable. However, a representation may become square-integrable when restricted to an appropriate quotient group. The general theory has been developed by Antoine, Ali and Gazeau, and may be found in [START_REF] Ali | Coherent states and their generalizations: an overview[END_REF], where several particular cases are also studied.

Conclusion

We have described in this paper a (limited) number of aspects of wavelet analysis and time-frequency analysis, with emphasis on some signal processing related problems. Our goal was to give an idea of the broadness of the theory, as well as the huge range of applications and potential applications. The interested reader may refer to several excellent textbooks available in the literature (up to now, more than a hundred monographs have been published on the subject).

Before closing this review, let us stress that there are still many problems, connected to wavelet theory and applications, which are still to be solved. For the record, let us mention a few of them.

First, the connection between continuous and discrete wavelet systems is not completely understood. We have seen the algebraic and geometrical origin of continuous wavelet decompositions. The multiresolution approach seems to be also extremely constrained by algebraic arguments, which should be developed further (see [START_REF] Antoine | An algebraic approach to discrete dilations[END_REF]).

Also, a very important problem (from the practical point of view) is that of adaptive decomposition. We have seen that the "time-frequency toolbox" offers several different ways of decomposing a signal, which raises the problem of selecting the most appropriate one. This point has been studied by several authors (see for example [START_REF] Wickerhauser | Adapted Wavelet Analysis, from Theory to Software[END_REF]), and should be one of the main research areas in the near future.

  the one-dimensional affine group, or "ax + b" group with group operation (b, a) • (b ′ , a ′ ) = (b + ab ′ , aa ′ ) It may be proved G af f has only two inequivalent irreducible unitary representations on H 2 ± (R) = {f ∈ L 2 (R), f (±ω) = 0 ∀ω ≥ 0}, of the form [π(b, a) • f ] (

FIGURES

  FIGURES

Figure 1 :

 1 Figure 1: Example of a WV transform: the case of a linear chirp (coded with gray levels).

Figure 2 :

 2 Figure 2: Interferences with the WV transform: the case of the sum of three pulses.

Figure 3 :

 3 Figure 3: Gabor (middle) and wavelet (bottom) transforms of 625ms of speech signal (top).