Marc Louis

Maurice François
email: francois@lmt.ens-cachan.fr

M L M François Lmt

Marc François

Lmt

A plasticity model with yield surface distortion for non proportional loading

Keywords:

In order to enhance the modeling of metallic materials behavior in non proportional loadings, a modification of the classical elastic-plastic models including distortion of the yield surface is proposed. The new yield criterion uses the same norm as in the classical von Mises based criteria, and a "distorted stress" S d replacing the usual stress deviator S. The obtained yield surface is then "egg-shaped" similar to those experimentally observed and depends on only one new material parameter. The theory is built in such a way as to recover the classical one for proportional loading. An identification procedure is proposed to obtain the material parameters. Simulations and experiments are compared for a 2024 T4 aluminum alloy for both proportional and nonproportional tension-torsion loading paths.

INTRODUCTION

Most of the phenomenological plasticity models for metals use von Mises based yield criteria. This implies that the yield surface is a hypersphere whose center is the back stress X in the 5-dimensional space of stress deviator tensors. For proportional loadings, the exact shape of the yield surface has no importance since all the deviators governing plasticity remain colinear tensors. In contrast, for non proportional loadings, the stress path may affect any point of the yield surface, and this shape becomes very important as its normal gives the direction of the plastic flow through normality [START_REF] Hill | The mathematical theory of plasticity[END_REF]).

Phillips and co-workers [START_REF] Phillips | Experimental investigation of corners in the yield surface[END_REF]; [START_REF] Phillips | The effect of the loading path on the yield surface at elevated temperatures[END_REF]) carried out the first experiments on pure aluminum and obtained very complex shapes for yield surfaces: for each straight loading path, a corner forms in the front end and a flat zone forms in the rear. The model developed by [START_REF] Kurtyka | A geometric description of distortionnal plastic hardening of deviatoric metals[END_REF][START_REF] Kurtyka | Evolution equations for distortionnal plastic hardening[END_REF] is able to model such complex distortion with geometric considerations and material constants.

Many authors [START_REF] Bui | Ecrouissage des métaux[END_REF]; Shiratori et al. [1976aShiratori et al. [ , 1976b]]; [START_REF] Winstone | Influence of the prestress on the yield surface of the cast nickel superalloy MAR-M002 at elevated temperature[END_REF]; [START_REF] Rousset | Surface seuil de plasticité: détermination automatique et modélisation[END_REF]; [START_REF] Rousset | Sur la déformation des surfaces seuil en palasticité cyclique[END_REF]; [START_REF] Cheng | Experimental determination of strain-induced anisotropy during nonproportional straining of an Al/Mg alloy at room temperature[END_REF]; [START_REF] Wu | On the experimental determination of yield surfaces and some results on annealed 304 stainless steel[END_REF]; [START_REF] Khan | An experimental study on subsequent yield surface after finite shear prestraining[END_REF]; [START_REF] Boucher | Experimental studies of yield surfaces of aluminium alloy and low carbon steel under complex biaxial loading[END_REF]; [START_REF] Ishikawa | Subsequent yield surface probed from its current center[END_REF]) have also measured or modeled [START_REF] Gupta | An alternative formulation for interpolating experimental yield surfaces[END_REF]) distortion of yield surfaces. However, they have observed more "egg-shaped" surfaces (see Figs. 1 and2). This simple shape will be our first assumption and will be more precisely defined later. The second hypothesis made here is that the egg-axis is the backstress X, in other words that the yield surface remains invariant respect to a rotation (in the deviatoric space) around X. [START_REF] Dahan | Evolution model of the elastic domain by any loading[END_REF] introduce a new internal variable to describe the distortion with (S-X) as the eggaxis (S being the stress deviator). The choice between these two directions cannot be made on the basis of available testing and the chosen direction seem to fit the experiments correctly. A third hypothesis is that the distortion from the sphere to the egg is proportional to the ratio of the norm of the backstress X to the new kinematic hardening limit X l . This can be seen as a natural first step in modeling and leads to classical theory while X l is set at an infinite value.

The entire model is developed within the framework of irreversible thermodynamics [START_REF] Haplhen | Sur les matériaux standarts généralisés[END_REF]; [START_REF] Lemaitre | Mechanics of Solids Materials[END_REF]) involving the Helmotz free energy , the yield function f and the dissipation pseudo-potential F. The formulation allows an easy implementation in computerized finite element codes for a low calculation cost. An implicit Newton scheme for stress driven transformations is described and has been programmed.

CLASSICAL ELASTIC-PLASTIC MODEL

For details regarding the theory we refer to Lemaitre and [START_REF] Lemaitre | Mechanics of Solids Materials[END_REF]. Variables are the Cauchy's stress s associated with the elastic (reversible) strain e e ; the opposite of the stress deviator -S associated with the plastic (irreversible) strain e p ; the backstress X associated with the plastic strain a (kinematic hardening) and the isotropic hardening R associated with the plastic stain r (isotropic hardening). Constants are the stiffness tensor K and the plastic flow constant C. H(r) is the integral of the isotropic hardening function h(r). The thermodynamic forces A k =(s, -S, X, R) are obtained by differentiation of the following Helmoltz free energy  with respect to the state variables V k =(e e , e p , a, r).

(e e , a, r) = 1 2 e e :K:e e + 1 2 C a:a + H(r) (1.1)

s = K:e e (1.2) X = C a (1.3) R = H'(r) = h(r) (1.4)
The yield surface is described by the yield function f based on the von Mises normative expression. In this paper, all classical corrective terms in 3/2 or 2/3 (whose role is to simplify the equations obtained in tension-compression) have been removed as it has been preferred to simplify tensorial equations used in computerized calculation than some manual calculation used only for identifications. The yield stress  y is defined then to be 2/3 of the classical value identified in tension-compression. This leads to description of the yield criterion as an hypersphere of radius R+ y .

f(S, X, R) = || S -X || -R - y (1.5)
If the plastic flow condition [START_REF] François | Comportement mécanique des matériaux[END_REF]) is satisfied, the following equations governs the plastic flow.

(f = 0) & ( ∂f ∂S : dS > 0 )  (plastic flow)

(1.6)

The dissipation pseudo-potential F differs from the yield function f (Eqn 1.5) by a quadratic term in X which provides the non-linear kinematic hardening [START_REF] Armstrong | A mathematical representation of the multiaxial Baushinger effect[END_REF]) weighted by the corresponding material constant .

F(S, X, R) = f(S, X, R)+  2C X:X (1.7)
The state variable evolution dV k is given by derivation of the dissipation pseudopotential with respect to the associated thermodynamic forces A k and the plastic multiplier d.

dV k = - ∂F ∂A k d (1.8)
We deduce from Eqn (1.8) Hill's normality rule [1950] when considering the thermodynamic force -S (Eqn 1.9). The following flow rules take into account Eqn (1.7) in order to refer only to the yield function f instead of F.

de p = ∂f ∂S d (1.9) da = - ∂f ∂X d -  C X d (1.10) dr = - ∂f ∂R d (1.11)
The value of the plastic multiplier d is given by the consistency equation df=0.

d =       C ∂f ∂X : ∂f ∂X +  X: ∂f ∂X + h'(h -1 (R))       ∂f ∂R 2 -1 ∂f ∂S : dS (1.12)
The following paper presents a modification of this theory involving only the yield function (Eqn 1.5). Others equations remain unaltered.

EXPERIMENTS AND OBSERVATIONS ON THE YIELD SURFACES

The experimental results used in this paper have been reported by [START_REF] Rousset | Surface seuil de plasticité: détermination automatique et modélisation[END_REF]. These involved 2024-T4 aluminum alloy using tubular thin-walled specimens allowing two dimensional loading in tension and torsion. The tensile stress deviator is either classically written in the canonic base as (e 1 e 1 ) or in the tensorial base of the deviators presented below as s 1 e 1 . The shear deviatoric stress is either classically written in the canonic base as (e 1 e 2 +e 2 e 1 ) or in the tensorial base of the deviators presented below as s 2 e 2 . Correspondence between  and s 1 ,  and s 2 is given Eqn (2.3). We first sketch this tensorial base of deviators different from Ilyushin's one [START_REF] Tanaka | A nonproportionality parameter and a cyclic viscoplastic constitutive model taking into account amplitude dependances and memory effects of isotropic hardening[END_REF]).

Tensorial base

Most of the measurements of a material testing are related to tensorial values, the stress S and the plastic strain e p in the case of plasticity. In our case of two dimensional loadings, one can plot (in 2 dimensions) two curves related to variables with different physical dimension (i.e. stress to strain) and two others related to variables of the same tensor (stress to stress; plastic strain to plastic strain).

The first way to plot such curves is the use of the components related to the canonic base of the tensors, in our example  11 to  12 for stresses,  p 11 to  p 12 for strains. This leads to an affinity in the representation of the yield surface: the spherical von Mises criterion appears as an ellipse and, in the general case, the affinity can be confused with the distortion.

The second way is to use classical correction factors: (;  3 ) for stresses and ( p ;  p / 3 ) for plastic strains (with = 11 ; = 12 ;  p = p 11 ;  p =2 p 12 ). These corrective terms lead to a circular representation of the von Mises criterion but are defined only for the tension or torsion loadings; for example, the superposition of a (non acting in plasticity) hydrostatic pressure leads to a modification of these corrective terms.

The third way recalled here is to use an orthonormal tensorial base [START_REF] Mehrabadi | Eigentensors of linear anisotropic elastic materials[END_REF]). The five tensors e 1 …e 5 are deviatoric: e 1 is relative to a tension along the e 1 axis, e 2 …e 5 are pure shears, and e 6 is spherical. Their structure is given below with respect to the canonical base.

e 1 = 1 6         2 0 0 0 -1 0 0 0 -1 e 2 = 1 2         0 1 0 1 0 0 0 0 0 … e 5 = 1 2         0 0 0 0 1 0 0 0 -1 e 6 = 1 3         1 0 0 0 1 0 0 0 1 e I : e J =  IJ (2.1)
As e 1 …e 5 is a base of the deviatoric subspace and e 6 is the hydrostatic base, we have the following properties (H means the hydrostatic part of the stress s):

s = S + H ; S =  i=1:5
(s:e i ) e i ; H = (s:e 6 ) e 6 (2.2)

The presented testings have been made in tension-torsion whose classical notations for stresses are  for  11 and  for  12 . The corresponding projections in the proposed base are noted s 1 =s:e 1 and s 2 =s:e 2 with the following correspondence:

s 1 = 2 3  ; s 2 = 2  (2.3)
In the same way, classical notation for plastic strains are  p for  p 11 and  p for 2 p 12 . The corresponding projections in the proposed base are noted e 1 =e p :e 1 and e 2 =e p :e 2 with the following correspondence:

e 1 = 3 2  p ; e 2 =  p 2 (2.4)
Use of a tensorial base leads to consider the plot of two components of a unique tensor as the cross section of the tensorial space by a well chosen hyperplane [START_REF] Rychlewski | On Hookes's Law[END_REF]). All the properties valid in the space of deviators will remain in the 2dimensional representation. For example, the von Mises based criterion used in the classical model will be represented as a circle whose radius is equal to R+ y , independent of the nature of the loading; an angle defined in the tensorial space (for exaple  defined later in Eqn (3.4)) is measurable on a plot such as Figs.

(1 and 2).

Combined with the simplified criteria as in Eqn (1.5) and later in Eqn (3.1) in which the corrective terms in 3/2 are suppressed these bases lead to many other simplifications. For example, the slope of the stress-strain curve (Fig. 4) at the beginning of plasticity is exactly C for any loading (not only the tensile one); the cumulative plastic strain p defined as   de p :de p , is exactly the length of the curve in the plastic strain to plastic strain graph.

Experimental results

The first experiment (Figs. 1 and4) is a monotonic tensile one, stopped at four different steps (O, A, I, B) for the measurement of yield surfaces at the following stresses =(0, 280, 320, 350 MPa) or s 1 =(0, 229, 261, 286 MPa) (Eqn 2.3). Only the steps (0, A, B) have been plotted for clarity.

The second experiment (Fig. 2) is a non proportional tension-torsion loading stopped at four steps (0, 0 1 , I, A, B) with the respective stresses =(0, 275, 275, 275, 275) MPa and =(0, 0, 75, 120, -130) MPa or s 1 =(0, 225, 225, 225, 225) MPa and s 2 =(0, 0, 106, 170, -184) MPa (Eqn 2.3). Only the steps (0, A, B) have been plotted for clarity.

At each step, a "straight loading" (a straight line in the deviatoric stress space) was prescribed from the approximate center of the yield surface (approximating X) until the detection of a non linearity in the stress to strain curve; the corresponding offset (strain) was set at 5.10 -5 . For each state M, the loading stress S M is plotted by a dark circle; each i th experimental point of the yield surface S iM for the state M is plotted as an open circle (Figs. 1 and2). In order to minimize the effect of the accumulated plasticity involved in the detection of the yield surface, the loading path was designed such as two successive yield stresses be as far as possible in the stress space (i.e. describing a star for five measurements).

It can be easily seen that the classical models, in which the yield surface is described by a von Mises based criterion, leading to a circular representation in Figs.

(1 and 2), deviate significantly from the shape determined by experiments.

NEW PLASTICITY MODEL FOR ELASTIC-PLASTIC METALS

In this section we propose a new yield criterion from which we deduce evolution laws. Its construction is illustrated in Fig. (3). The distorted stress S d is defined so that such as all the deviatoric stresses S for which the yield function f is null belong to a "hyper-egg" (first hypothesis presented in the introduction). We decompose the deviatoric stress S into its part S x collinear to X and it's orthogonal part S o according to:

S = S x + S o ; S x = S:X || X || 2 X (3.2)
We define the distortion to be parabolic (for simplicity) with respect to ||S o || (can be seen in real size on Fig. 3).

S d = S + S o :S o 2 X l (R +  y ) X (3.3)
The loci of S d for f=0 correspond to a von Mises' hypersphere whereas the yield surface, i.e. the loci of S for f=0, correspond to an egg-shaped surface curve for any section by an hyperplane containing X (Fig. 3) and corresponds to a circle of radius (R+ y ) for any section by a hyperplane orthogonal to X. The egg's symmetry axis is the backstress X (second hypothesis). The distortion is linearly proportional to the ratio of the norm of X to the new constant introduced, the "kinematic hardening limit" X l (third hypothesis). This leads naturally to a virgin state (in which the back stress is null) satisfying the von Mises criterion and the existence of a limit state for which the distortion is complete. The distortion is maximum when ||X||=X l . At this step, the elastic domain has a locally flat zone (at stress a in Fig. 3), but remains convex as for all states. Note that the use of classical non linear kinematic hardening rules described in section (1) lead automatically to such limit value (C/) for ||X||. The backstress X represents, as for the classical reference model, the middle of the elastic domain in the case of proportional loading. Another choice would be to write the model in order to define X as the center of the maximum diameter of the egg but this leads to an influence of the distortion even on proportional loadings, making identification more difficult [START_REF] François | A simple model for the determination of yield surfaces in non proportional loadings[END_REF]).

In proportional loadings, as tensorial variables (S, X, e p ) remain colinear, the othogonal part S o remains null, S d remains always equal to S and the classical model described in section (1) is recovered.

In the new model, the normal n of the yield surface is not collinear to S-X (as is the case for the classical model) when the loading is non proportional (when S does not remain collinear to X); this will affect behavior of the model during non proportional loadings involving plasticity, as in two-dimensional ratchetting for instance.

Another point of view can be given while using as parameter the following angle  between (S d -X) and X represented in real size in Fig. (3) which allows us to describe every possible state for the material with a reduced set of variables:

cos() = (S d -X):X ||S d -X|| ||X|| (3.4)
Let us suppose a state defined by (S, X, R). The tensors S and X define a hyperplane in the 5-dimensional space of deviators with an orthonormal basis x, o defined as:

x = X ||X|| ; o = S o ||S o || (3.5)
The yield surface in this plane can be described respect to the angle  in the base (x, o).

The state is then fully described by X (the norm of X), R and (with 0≤≤). Let us consider now a case involving plasticity, i.e. f=0 (Eqn 1.6). The following equations give a simple geometrical interpretation of the yield surface and show clearly the link to classical theory (X l  leads to circular yield surfaces).

S d -X = (R+ y ) (cos() x + sin() o) (3.6) S -X = (R+ y )             cos() -sin 2 () X 2X l x + sin() o (3.7)
Introducing the normal of the yield surface relative to the distorted stress n d (Fig. 3), the following gradients can be obtained with respect to the thermodynamic forces (-S, X, R):

n d = ∂f ∂S d = S d -X ||S d -X|| (3.8) ∂f ∂S = n d + n d :X X l (R+ y ) S o ∂f ∂X = -       1 - S o :S o 2X l (R+ y ) n d -       n d :S x X l (R+ y ) S o ∂f ∂R = -       S o :S o 2X l (R+ y ) 2 n d :X -1 (3.9)
It can be again easily verified that the form of the classical model can be found either when S o is null (uniaxial proportional case) or when X l tends to infinity. All these equations can be rewritten with the angle  (Eqn 3.4) in the (x, o) base.

n d = cos() x + sin() o (3.10) ∂f ∂S = cos() x + sin()       1+cos() X X l o ∂f ∂X = -       1- (R+ y )sin 2 () 2X l cos() x -       1 + R+ y 2X l ( ) 3cos 2 ()-1 + X cos() 2X l       2-sin 2 () R+ y X l sin() o ∂f ∂R = - sin 2 () cos() X 2X l -1 (3.11)
In the classical model, dr is always equal to dp, the increment of accumulated plastic strain p defined as dp=||de p ||. Fig. (5) shows the ratio dp/dr with respect to the angle  in the case ||X||=X l for which it's minima/maxima are the farthest from 1 but dp/dr always remains between 0.811 and 1.190 and is equal to 1 for proportional loadings (=0 or =) and for =/2. This indicates that we can keep the constants of the isotropic hardening function h(r) unchanged from the ones identified for the classical model and one can say that the isotropic hardening is still driven by the accumulated plastic strain.

We have not yet proven analytically the positiveness of the dissipation dD: dD = S:de p -X:da -R dr

dD =       S: ∂f ∂S + X: ∂f ∂X +  C X:X + R ∂f ∂R d (3.12)
We use here a numerical study. The evolution equation (Eqns 1.9, 1.10 and 1.11) and the expression of gradients (Eqn 3.11) show that all the involved tensors are in the (x, o) plane even if dS is out of plane since only its projection is used in the expression of d (Eqn. 1.12). The plastic multiplier d is calculated while setting ∂f/∂S:dS to an arbitrary positive value (Eqn.1.6). The minimum of dD in the (x, o) plane is numerically determined with respect to a reduced set of variables [, X, R] and material constants [ y , X l , C, ] (and the 2 constants [k, m] involved in the isotropic hardening power law h(r)=k r 1/m used here). By normalizing by  y all the parameters proportional to a stress (i.e. X, R,  y , X l , C and k), we consider 8 parameters. For each of them, two values are considered, namely the minimum and maximum values:

(0 ≤  ≤ ), (0 ≤ X < C/), (0 ≤ R < L y ), (C/ ≤ X l < L y ), (0 < C < L y ), (0 ≤  < L), (0 ≤ k < L y ), (0 < m < L),
where L is a "large" value (10 5 ). To avoid reaching local minima, we carried out the minimization with 2 8 initial values corresponding to the combinations for the 8 parameters with respect to their minimum / maximum values. Within these bounds, the minimum of the dissipation has been found positive. Therefore no new restriction on the choice of material constants is found to fulfill the second principle of thermodynamics.

An implicit scheme of Newton-Raphson type has been developed in order to model a stress path on a representative volume element. Calculation is then simple as the dissipation equations are stress driven and as elastic strain and Young's modulus are not involved. The Figs.

(1 and 2) show in plain lines the simulated yield surfaces for each state, obtained with the material constants identified in section 4.

IDENTIFICATION

The material constants to be identified are the yield stress  y , the non linear kinematic hardening parameters (C, ), the constant X l associated to the distortional effect and the isotropic hardening function h(r). The optimization procedure is first made using yield surface measurements (in the stress space) and then using the stress-strain curve. The experiment used for identification is the monotonic tensile test (Figs. 1 and4). We denote S iM the i th experimentally measured stress belonging to the yield surface for the state M. The yield stress  y has been identified as the average norm of measured yield surface stresses S iO on virgin material (O denotes the virgin state, and the brackets denotes the average on the term i):

 y = < || S iO || > (4.1)
Let us introduce the following distance d M from the measured stresses S iM to any theoretical yield surface defined by X and R at state M as follows:

d M (X, R; X l ,  y ) = < | f(S iM , X, R; X l ,  y ) | > (4.2)
This definition is close (equal if there is no distortion) to the average of real distances (in the stress space) from S iM to the yield surface defined by (X, R; X l ,  y ). As the considered tensile testing (used for identification) imposes that X is collinear to the stress deviator e 1 (Eqn 2.1), only the norm X=||X|| remains unknown. For each step M, the measured yield stresses S iM allow us to determine numerically the three values (X M , R M , X lM ) which define the closest theoretical yield surface to experimental data (the closest egg) as follows:

(X M , R M , X lM ) = arg     min (X, R, X l ) ( ) d M (Xe 1 , R; X l ,  y ) (4.
3)

The kinematic hardening limit X l is chosen as the average of the X lM values obtained for each step. As the length of the egg-shaped yield surface along the e 1 axis is equal to 2( y +R M ) in tension, the R M can be plotted versus the experimentally obtained accumulated plastic strain p M =||e pM || in our monotonic tension. This allows us to fit the isotropic hardening rule R=h(r) (as r=p in the considered monotonic case, see section 3).

In this paper we have used a classical power law, R=kr (1/m) which needs at least three steps to be identified. The remaining C and  constants are identified with the stress versus plastic strain curve (Fig. 4); C is the slope of the curve at the beginning of plasticity and  is defined by the saturation value for X, i.e. C/. One can remark that the "natural value" for the kinematic hardening limit introduced, X l is X l =C/. In the example, X l =130 MPa and C/=114.3 MPa, this means that the complete distortion (locally flat back) cannot be reached. The values obtained for 2024-T4 aluminum alloy are given for a power isotropic hardening law R=kr (1/m) in Table (1).

COMPARISON OF THE SIMULATIONS WITH THE EXPERIMENTS

Let us compare theoretical (solid lines) to experimental (circles) results for the monotonic tensile test described in section 3. The constants used (Table 1) have been obtained from the identification process described above. Concerning the stress versus plastic strain curve (Fig. 4), the first step A is not accurately fit, due to the simple hardening rules used (the distortion modification is not involved in this curve for monotonic loading). The average distance (Eqn 4.2) between simulation and measurements is 5.99 MPa for the four steps. The same distance, computed with the same material constants but with X l set to an infinity (classical model) is 11.18 MPa. In this case, the yield surface could be represented (on Fig. 1) by a circle coinciding to the egg-shape of the presented model at the stresses collinear to X (see Fig. 3). Some points on the back side (opposite to the corner) of yield surfaces are not very well fitted by the new model. In fact, these points are much more difficult to detect precisely: the stress strain curve presents, at the beginning of plasticity, a very progressive curvature in the unloading direction instead of the (easy to detect) sharp transition in the loading direction. Some experimental yield stresses do not comply with the symmetry with respect to the e 1 tensorial direction (horizontal axis on Fig. 1) that is expected. This may be due to an initial anisotropy of the material or some other experimental problem.

The second experiment, a non proportional tension-torsion loading, is described in section 2. The theoretical curves are obtained with the values obtained by the identification process described earlier on the previous monotonic tensile testing (Table 1). The average distance (Eqn 4.2) from measured yield surfaces to simulated ones is 11.25 MPa. Considering all the same constants identical but X l infinite (classical model) the distance is 13.46 MPa. The corresponding yield surfaces have circular representation on Fig.

(3), they are not drawn for clarity. The observation made before about back side points for tensile testing remains true, but the error is considerably bigger than in the monotonic case; this leads to the supposition that the isotropic hardening rule, as a monotonic function of r (or p in the classic case, see section 3) overestimates the size of the elasticity domain in the case of non proportional loading. Finally, one can observe that the experimental points seem to be distant from the theoretical curve proportionally to the intensity of the shear stress. This is demonstrated also for the initial (von Mises type) yield surface. The aluminum may obey to a criterion distinguishing the tension's deviator from other shears, e.g. Tresca's criterion. But these criteria do not respect [START_REF] Ilyushin | On the relation between stresses and small strains in the mechanics of continua[END_REF] postulate thus prohibiting the very useful intrinsic writing of plasticity in deviatoric subspace.

CONCLUSION

The presented theory allows a simple modeling of the distortion of yield surfaces within the classical thermodynamical framework; it has many possible applications such as metal forming, fatigue, forming limit curves, etc, and constitutes a generalization of the classical theory. It can be easily adapted to more sophisticated models including, for example, the models improved for ratchetting [START_REF] Chaboche | Modeling of ratchetting: evaluation of various approaches[END_REF]), with a reasonable augmentation of calculation time. Some variants can be envisaged: it is possible to link the distortion effect to a variable other than X, for example one of the X i for multiple kinematic hardening models, or to the plastic strain e p . Some two-dimensional ratchetting experiments are planned as the disortion effect takes on a great importance in this type of loading (mainly through the direction of the normal of the yield surface).

Fig. 1 .

 1 Fig. 1. Yield surfaces for monotonic loading.

Fig. 2 .

 2 Fig. 2. Yield surfaces for non proportional tension-torsion testing.

Fig. 3 .

 3 Fig. 3. Schematic of the yield surface.In order to keep the very simple formulation given by the von Mises' norm, we introduce a new "distorted stress" S d which replaces the stress deviator S involved in the classical expression. One can remark that S d is not a new thermodynamic force, but only a function of existing variables. The other terms involved in the yield function f are the back stress X, the isotropic hardening function R and the yield stress  y which is also the radius of the initial yield surface.

Fig. 4 .

 4 Fig. 4. Stress to plastic strain for tensile testing.

Fig. 5 .

 5 Fig. 5. Ratio dp/dr as a function of .

  Fig. (1) allows the comparison of yield surfaces in the tensiontorsion [ e 1 e 2 ] hyperplane of deviatoric stresses (according to tensorial base, Eqn 2.1).

Table 1 .

 1 Material constants identified for the 2024-T4 aluminium alloy.

	page 13


Peter L Søndergaard

An Efficient Algorithm for the Discrete Gabor Transform using full length Windows

 to the Gabor analysis/synthesis operator. The factorization provides a fast method for computing the discrete Gabor transform (DGT) and several algorithms associated with it. The factorization algorithm should be used when the involved window and signal have the same length. An optimized implementation of the algorithm is freely available for download.

I. INTRODUCTION

The finite, discrete Gabor transform (DGT) of a signal f of length L is given by c (m, n, w) = L-1 l=0 f (l, w)g (l -an)e -2πiml/M .

(

) 1
Here g is a window (filter prototype) that localizes the signal in time and in frequency. The DGT is equivalent to a Fourier modulated filter bank with M channels and decimation in time a, [START_REF] Bölcskei | Equivalence of DFT filter banks and Gabor expansions[END_REF].

Efficient computation of a DGT can be done by several methods: If the window g has short support (consists of relatively few filter taps), a filter bank based approach can be used. We shall instead focus on the case when g and f are equally long. The main advantage of the algorithm presented is its ease of use: The running time is guaranteed to be small even for long windows. This allows for the practical use of non-compactly supported windows like the Gaussian and its tight and dual windows without truncating them.

In the case when the window and signal have the same length, a factorization of the frame operator matrix was found by Zibulski and Zeevi in [START_REF] Zeevi | Oversampling in the Gabor scheme[END_REF]. The method was initially developed in the L 2 (R) setting, and was adapted for the finite, discrete setting by Bastiaans and Geilen in [START_REF] Bastiaans | On the discrete Gabor transform and the discrete Zak transform[END_REF]. They extended it to also cover the analysis/synthesis operator. A simple, but not so efficient, method was developed for the Gabor analysis/synthesis operator by Prinz in [START_REF] Prinz | Calculating the dual Gabor window for general sampling sets[END_REF]. Strohmer [START_REF] Strohmer | Numerical algorithms for discrete Gabor expansions[END_REF] improved the method and obtained the lowest known computational complexity for computing the Gabor frame operator. This paper extends Strohmer's method to also cover the Gabor analysis and synthesis operators.

The advantage of the method developed in this paper as compared to the one developed in [START_REF] Bastiaans | On the discrete Gabor transform and the discrete Zak transform[END_REF], is that it works with FFTs of shorter length, and does not require multiplication by complex exponentials caused by the quasi-periodicity of the Zak transform. The two methods have the same asymptotic complexity, O (N M log M ), where M is the number of channels and N is the number of time steps. A more accurate flop count is presented later in the paper.

We shall study the DGT applied to multiple signals at once. This is for instance a common subroutine in computing a multidimensional DGT. The DGT defined by (1) works on a multi-signal f ∈ C L×W , where W ∈ N is the number of signals.

II. DEFINITIONS

We shall denote the set of integers between zero and some number L by L = 0, . . . , L -1.

(

) 2
The Discrete Fourier Transform (DFT) of a signal f ∈ C L is defined by

(F L f ) (k) = 1 √ L L-1 l=0 f (l)e -2πikl/L . (3)
We shall use the • notation in conjunction with the DFT to denote the variable over which the transform is to be applied.

To denote all elements indexed by a variable we shall use the : notation. As an example, if C ∈ C M ×N then C :,1 is a M × 1 column vector, C 1,: is a 1 × N row vector and C :,: is the full matrix. This notation is commonly used in Matlab and FORTRAN programming and also in some prominent textbooks, [START_REF] Golub | Matrix computations[END_REF].

The convolution f * g of two functions f, g ∈ C L and the involution f * is given by

(f * g) (l) = L-1 k=0 f (k) g (l -k) , l ∈ L (4) f * (l) = f (-l), l ∈ L . (5)
It is well known how convolution can be computed efficiently using the discrete Fourier transform. We shall use a variant of this result

(f * g * ) (l) = √ LF -1 L (F L f ) (•) (F L g) (•) (l) . ( 6
)
The Poisson summation formula in the finite, discrete setting is given by

F M b-1 k=0 g(• + kM ) (m) = √ b (F L g) (mb), (7)
where

g ∈ C L , L = M b with b, M ∈ N. A family of vectors e j , j ∈ J of length L is called a frame if constants 0 < A ≤ B exist such that A f 2 ≤ J-1 j=0 | f, e j | 2 ≤ B f 2 , ∀f ∈ C L . ( 8
) Algorithm 1 Window factorization WFAC(g, a, M ) 1) for r = c k = p , l = q 2) for s = d 3) tmp (s) ← g (r + c • (k • q -l • p + s • p • q mod d • p • q)) 4)
end for 5)

P hi (r, k, l, :) ←DFT(tmp) 6) end for 7) return Phi

The constants A and B are called lower and upper frame bounds. If A = B, the frame is called tight. If J > L, the frame is redundant (oversampled). Finite-and infinite dimensional frames are described in [START_REF] Christensen | An Introduction to Frames and Riesz Bases[END_REF].

A finite, discrete Gabor system (g, a, M ) is a family of vectors g m,n ∈ C L of the following form

g m,n (l) = e 2πilm/M g (l -na) , l ∈ L (9)
for m ∈ M and n ∈ N where L = aN and M/L ∈ N. A Gabor system that is also a frame is called a Gabor frame. The analysis operator C g : C L → C M ×N associated to a Gabor system (g, a, M ) is the DGT given by given by ( 1). The Gabor synthesis operator D γ : C M ×N → C L associated to a Gabor system (γ, a, M ) is given by

f (l) = N -1 n=0 M -1 m=0 c (m, n) e 2πiml/M γ (l -an) . (10)
In ( 1), ( 9) and ( 10) it must hold that L = N a = M b for some M, N ∈ N. Additionally, we define c, d, p, q ∈ N by

c = gcd (a, M ) , d = gcd (b, N ) , (11)
p = a c = b d , q = M c = N d , (12)
where GCD denotes the greatest common divisor of two natural numbers. With these numbers, the redundancy of the transform can be written as L/ (ab) = q/p, where q/p is an irreducible fraction. It holds that L = cdpq. The Gabor frame operator S g : C L → C L of a Gabor frame (g, a, M ) is given by the composition of the analysis and synthesis operators

S g = D g C g .
The Gabor frame operator is important because it can be used to find the canonical dual window g d = S -1 g g and the canonical tight window g t = S -1/2 g g of a Gabor frame. The canonical dual window is important because D g d is a left inverse of C g . This gives an easy way to construct an inverse transform of the DGT. Similarly, then D g t is a left inverse of C g t . For more information on Gabor systems and properties of the operators C, D and S see [START_REF] Gröchenig | Foundations of Time-Frequency Analysis[END_REF], [START_REF]Gabor Analysis and Algorithms[END_REF], [START_REF]Advances in Gabor Analysis[END_REF].

III. THE ALGORITHM

We wish to make an efficient calculation of all the coefficients of the DGT. Using (1) literally to compute all coefficients c (m, n, w) would require 8M N LW flops.

To derive a faster DGT, one approach is to consider the analysis operator C g as a matrix, and derive a faster algorithm Algorithm 2 Discrete Gabor transform DGT(f, g, a, M )

1) P hi =WFAC(g, a, M ) 2) for r = c 3) for k = p , l = q , w = W 4) for s = d 5) tmp (s) ← f (r + (k • M + s • p • M -l • h a • a mod L) , w) 6) end for 7) P sitmp (k, l + w • q, •) ←DFT(tmp) 8)
end for 9)

for s = d 10)

G ← P hi (:, :, r, s) 11)

F ← P sitmp (:, :, s) 12)

Ctmp (:, :, s) ← G T • F 13) end for 14)

for u = q , l = q , w = W 15)

tmp ←IDFT(Ctmp (u, l + w • q, :)) 16)

for s = d 17) coef (r + l • c, u + s • q -l • h a mod N, w) ← tmp (s) 18)
end for [START_REF] Zeevi | Oversampling in the Gabor scheme[END_REF] end for 20) end for 21) for n = N ,w = W 22) coef (:, n, w) ←DFT(coef (:, n, w)) 23) end for 24) return coef through unitary matrix factorizations of this matrix. This is the approach taken by [START_REF] Strohmer | Numerical algorithms for discrete Gabor expansions[END_REF], [START_REF] Qiu | Discrete Gabor transforms: The Gabor-gram matrix approach[END_REF]. Unfortunately, this approach tends to introduce many permutation matrices and Kronecker product matrices. Another approach is the one taken in [START_REF] Bastiaans | On the discrete Gabor transform and the discrete Zak transform[END_REF] where the Zak transform is used. This approach has the downside that values outside the fundamental domain of the Zak transform require an additional step to compute. In this paper we have chosen to derive the algorithm by directly manipulating the sums in the definition of the DGT.

To find a more efficient algorithm than (1), the first step is to recognize that the summation and the modulation term in (1) can be expressed as a DFT:

c (m, n, w) = √ LF L f (•, w)g (• -an) (mb) . ( 13
)
We can improve on this because we do not need all the coefficients computed by the Fourier transform appearing in [START_REF] Janssen | Iterative algorithms to approximate canonical Gabor windows: Computational aspects[END_REF], only every b'th coefficient. Therefore, we can rewrite by the Poisson summation formula [START_REF]Advances in Gabor Analysis[END_REF]:

c (m, n, w) = √ M F M b-1 m=0 f (• + mM, w)g (• + mM -an) (m) = (F M K (•, n, w)) (m) , (14)
where

K (j, n, w) = √ M b-1 m=0 f (j + mM, w) g (j + mM -na) , (15)
for j ∈ M and n ∈ N . From ( 14) it can be seen that computing the DGT of a signal f can be done by computing K followed by DFTs along the first dimension of K.

To further lower the complexity of the algorithm, we wish to express the summation in (15) as a convolution.

We split j as j = r + lc with r ∈ c , l ∈ q and introduce h a , h M ∈ Z such that the following is satisfied:

c = h M M -h a a. ( 16
)
The two integers h a , h M can be found by the extended Euclid algorithm for computing the GCD of a and M . Using [START_REF] Qiu | Discrete Gabor transforms: The Gabor-gram matrix approach[END_REF] and the splitting of j we can express [START_REF] Prinz | Calculating the dual Gabor window for general sampling sets[END_REF] as

K (r + lc, n, w) = √ M b-1 m=0 f (r + lc + mM, w) × ×g (r + l (h M M -h a a) + mM -na) (17)
= √ M b-1 m=0 f (r + lc + mM, w) × ×g (r + ( m + lh M ) M -(n + lh a ) a) (18)
We substitute m + lh M by m and n + lh a by n and get

K (r + lc, n -lh a , w) = √ M b-1 m=0 f (r + lc + ( m -lh M ) M, w) × ×g (r + mM -na) (19)
= √ M b-1 m=0 f (r + mM + l (c -h M M ) , w) × ×g (r + mM -na) (20)
We split m = k + sp with k ∈ p and s ∈ d and n = u + sq with u ∈ q and s ∈ d and use that M = cq, a = cp and c -h M M = -h a a:

K (r + lc, u + sq -lh a , w) = √ M p-1 k=0 d-1 s=0 f (r + kM + spM -lh a a, w) × ×g (r + kM -ua + (s -s) pM ) (21)
After having expressed the variables j, m, n using the variables r, s, s, k, l, u we have now indexed f using s and g using (s -s). This means that we can view the summation over s as a convolution, which can be efficiently computed using a discrete Fourier transform. Define

Ψ f r,s (k, l + wq) = F d f (r + kM + •pM -lh a a, w) , (22) Φ g r,s (k, u) = √ M F d g (r + kM + •pM -ua) , (23)
Using [START_REF]Gabor Analysis and Algorithms[END_REF] we can now write (21) as 8M N L Eq. ( 14)

K (r + lc, u + sq -lh a , w) = √ d p-1 k=0 F -1 d Ψ f r,• (k, l + wq) Φ g r,• (k, u) (s) (24) = √ dF -1 d p-1 k=0 Ψ f r,• (k, l + wq) Φ g r,• (k, u) (s) (25)
8L Lg a + 4N M log 2 (M ) [1] L " 8q + 1 + q p " + 4L " 1 + q p " log 2 N + 4M N log 2 (M ) Alg. 2 L (8q) + 4L " 1 + q p " log 2 d + 4M N log 2 (M )
Flop counts for 4 different way of computing the DGT: By the linear algebra definition [START_REF] Bastiaans | On the discrete Gabor transform and the discrete Zak transform[END_REF], by the method based on Poisson summation [START_REF] Johnson | A Modified Split-Radix FFT With Fewer Arithmetic Operations[END_REF], by the method of Bastiaans and Geilen from [START_REF] Bastiaans | On the discrete Gabor transform and the discrete Zak transform[END_REF] and by Algorithm 2. The term Lg denotes the length of the window used so Lg/a is the overlapping factor of the window. Note for comparison that log 2 N = log 2 d + log 2 q

If we consider Ψ f r,s and Φ g r,s as matrices for each r and s, the sum over k in the last line can be written as matrix products. Algorithm 2 follows from this.

IV. RUNNING TIME

When computing the flop count of the algorithm, we will assume that a complex FFT of length M can be computed using 4M log 2 M flops. A nice review of flop counts for FFT algorithms is presented in [START_REF] Johnson | A Modified Split-Radix FFT With Fewer Arithmetic Operations[END_REF]. Table I shows the flop count for Algorithm 2 and compares it with the definition of the DGT (1), with the algorithm for short windows using Poisson summation [START_REF] Johnson | A Modified Split-Radix FFT With Fewer Arithmetic Operations[END_REF] and with the algorithm published in [START_REF] Bastiaans | On the discrete Gabor transform and the discrete Zak transform[END_REF]. The algorithm by Prinz presented in [START_REF] Prinz | Calculating the dual Gabor window for general sampling sets[END_REF] has the same computational complexity as the Poisson summation algorithm. For simplicity we assume that both the window and signal are complex valued. In the common case when both f and g are real-valued, all the algorithms will see a 2 to 4 times speedup.

The flop count for definition (1) is that of a complex matrix multiplication. All the other algorithms share the 4M N log 2 M term coming from the application of an FFT to each 'block' of coefficients and only differ in how the application of the window is performed. The Poisson summation algorithm is very fast for a small overlapping factor L g /a, but turns into an O L 2 algorithm for a full length window. In this case the other algorithms have an advantage. The term L 8q + 1 + q p in the [START_REF] Bastiaans | On the discrete Gabor transform and the discrete Zak transform[END_REF] algorithm comes from calculation of the needed Zak-transforms, and the 4L 1 + q p log 2 N term comes from the transform to and from the Zak-domain. Compared to ( 22) and (23) this transformation uses longer FFTs. Algorithm 2 does away with the multiplication with complex exponentials in the [START_REF] Bastiaans | On the discrete Gabor transform and the discrete Zak transform[END_REF] algorithm, and so the first term reduces to L (8q).

Both the Poisson summation based algorithm and Algorithm 2 can do a DGT with L ≈ 2000000 in less than 1 second on a standard PC at the time of writing. We have not created an efficient implementation of the algorithm from [START_REF] Bastiaans | On the discrete Gabor transform and the discrete Zak transform[END_REF] in C so therefore we cannot reliably time it.

V. EXTENSIONS

The algorithm just developed can also be used to calculate the synthesis operator D γ . This is done by applying Algorithm

Algorithm 3 Canonical Gabor dual window GABDUAL(g, a, M ) 1) P hi =WFAC(g, a, M ) 2) for r = c , s = d 3) G ← P hi (:, :, r, s) 4)

P hi d (:, :, r, s) ← G • G T -1 • G 5) end for 6) g d =IWFAC P hi d , a, M 7) return g d 2 in the reverse order and inverting each line. The only lines that are not trivially invertible are lines 10-12, which becomes 10) Γ ← P hi d (:, :, r, s) 11) C ← Ctmp (:, :, s) 12) P sitmp (:, :, s) ← Γ • C where the matrices P hi d (:, :, r, s) should be left inverses of the matrices P hi (:, :, r, s) for each r and s.

The matrices P hi d (:, :, r, s) can be computed by Algorithm 1 applied to a dual Gabor window γ of the Gabor frame (g, a, M ). It also holds that all dual Gabor windows γ of a Gabor frame (g, a, M ) must satisfy that P hi d (:, :, r, s) are left inverses of the matrices P hi (:, :, r, s). This criterion was reported in [START_REF] Janssen | On rationally oversampled Weyl-Heisenberg frames[END_REF], [START_REF] Janssen | The duality condition for Weyl-Heisenberg frames[END_REF].

A special left-inverse in the Moore-Penrose pseudo-inverse. Taking the pseudo-inverses of P hi (:, :, r, s) yields the factorization associated with the canonical dual window of (g, a, M ), [START_REF] Christensen | Frames and pseudo-inverses[END_REF]. This is Algorithm 3. Taking the polar decomposition of each matrix in Φ g r,s yields a factorization of the canonical tight window (g, a, M ). For more information on these methods, as well as iterative methods for computing the canonical dual/tight windows, see [START_REF] Janssen | Iterative algorithms to approximate canonical Gabor windows: Computational aspects[END_REF].

VI. SPECIAL CASES

We shall consider two special cases of the algorithm: The first case is integer oversampling. When the redundancy is an integer then p = 1. Because of this we see that c = a and d = b. This gives [START_REF] Qiu | Discrete Gabor transforms: The Gabor-gram matrix approach[END_REF] the appearance

indicating that h M = 0 and h a = -1 solves the equation for all a and q. The algorithm simplifies accordingly, and reduces to the well known Zak-transform algorithm for this case, [START_REF] Janssen | The Zak transform: a signal transform for sampled time-continuous signals[END_REF].

The second case is the short time Fourier transform. In this case a = b = 1, M = N = L, c = d = 1, p = 1, q = L and as in the previous special case h M = 0 and h a = -1. In this case the algorithm reduces to the very simple and well known algorithm for computing the STFT.

VII. IMPLEMENTATION

The reason for defining the algorithm on multi-signals, is that the multiple signals can be handled at once in the matrix product in line 12 of Algorithm 2. This is a matrix product of two matrices size q × p and p × qW , so the second matrix grows when multiple signals are involved. Doing it this way reuses the Φ g r,s matrices as much as possible, and this is an advantage on standard, general purpose computers with a deep memory hierarchy, see [START_REF] Dongarra | A set of level 3 basic linear algebra subprograms[END_REF], [START_REF] Whaley | Automated empirical optimization of software and the ATLAS project[END_REF]. The benefit of expressing Algorithm 2 in terms of loops (as opposed to using the Zak transform or matrix factorizations) is that they are easy to reorder. The presented Algorithm 2 is just one among many possible algorithms depending on in which order the r, s, k and l loops are executed. For a given platform, it is difficult a priory to estimate which ordering of the loops will turn out to be the fastest. The ordering of the loops presented in Algorithm 2 is the variant that uses the least amount of extra memory.

Implementations of the algorithms described in this paper can be found in the Linear Time Frequency Toolbox (LTFAT) available from http://ltfat.sourceforge.net. The implementations are done in both the Matlab/Octave scripting language and in C. A range of different variants of Algorithm 2 has been implemented and tested, and the one found to be the fastest on a small range of computers is included in the toolbox.


Florent Jaillet
email: florent@kfs.oeaw.ac.at

Peter Balazs
email: peter.balazs@oeaw.ac.at

Monika Dörfler

Nonstationary Gabor Frames

To overcome the limitation induced by the fixed timefrequency resolution over the whole time-frequency plane of Gabor frames, we propose a simple extension of the Gabor theory leading to the construction of frames with timefrequency resolution evolving over time or frequency. We describe the construction of such frames and give the explicit formulation of the canonical dual frame for some conditions. We illustrate the interest of the method on a simple example.

Introduction

Gabor analysis [START_REF] Feichtinger | Gabor Analysis and Algorithms -Theory and Applications[END_REF] is widely used for applications in signal processing. For some of these applications, which include processing of signals using Gabor frame multipliers [START_REF] Feichtinger | A first survey of Gabor multipliers[END_REF][START_REF] Balazs | Basic definition and properties of Bessel multipliers[END_REF], the rigid construction of the Gabor atoms results in important limitations on the signal analysis and processing ability of the associated schemes. The Gabor transform uses time-frequency atoms built by translation over time and frequency of a unique prototype function, leading to a signal decomposition having a fixed time-frequency resolution over the whole time-frequency plane. This can be very restricting when dealing with signals with characteristics changing over the time-frequency plane. For example, this led some people to prefer the use of alternative decompositions with time-frequency resolution evolving with frequency in some applications, to better fit the feature of interest of the signal. Examples of such decompositions are the wavelet transform [START_REF] Daubechies | Ten Lectures On Wavelets[END_REF] or the decompositions using filter banks based on perceptive frequency scales for processing of audio signals, as for example gammatone filters [START_REF] Hartmann | Signals, Sounds, and Sensation[END_REF]. A case for which the limitation induced by the constant time-frequency resolution of the Gabor transform can be seen is shown on the didactic example of Figure 1. On this figure, two spectrograms of the same glockenspiel signal are represented. These spectrograms are obtained by plotting the square absolute value of the Gabor coefficients using a color scale with a level coding in dB. Both spectrograms are obtained from the Gabor coefficients using the same type of window, but using two different window lengths. We see that the signal contains two very contrasting types of components:

• at the beginning of the notes, the signal presents sharp attacks which are spread in frequency, but very localized in time,

• during the resonance of the notes, the signal contains quasi-sinusoidal components which are spread in time, but very localized in frequency. We see that the use of the narrow window is well suited for the analysis and processing of the attacks, leading to a very sparse decomposition for these components, but gives an unsatisfying representation of the resonance, as the different sinusoidal components are not resolved. On the other hand, the wide window gives a good representation of the resonance part, but a blurred representation of the attacks.

For this example, it appears that if we want to build an optimised scheme for processing of both attacks and the resonances at the same time, it would be suitable to be able to adapt the time-frequency resolution locally for the different types of components. The purpose of this paper is to describe one way to achieve this goal. For this, we show that, while staying in the context of frame theory [START_REF] Casazza | The art of frame theory[END_REF][START_REF] Christensen | An Introduction To Frames And Riesz Bases[END_REF], the standard Gabor theory can be easily extended to provide some freedom of evolution of the time-frequency resolution of the decomposition in either time or frequency. Furthermore, this extension is well suited for applications as it can easily be implemented using fast algorithm based on fast Fourier transform [START_REF] Walker | Fast Fourier Transforms[END_REF]. We first describe the construction of the frames in Section 2., and then illustrate in Section 3. the potential of the approach on the preceding example of Figure 1.

Construction of the frames 2.1 Resolution evolving over time

As opposed to standard Gabor analysis, we replace time translation for the construction of atoms by the use of different windows for the different sampling positions in time. For each time position we still build atoms by regular frequency modulation. So using a set of functions {g n } n∈Z of L 2 (R), for m ∈ Z and n ∈ Z, we define atoms of the form:

g m,n (t) = g n (t)e i2πmbnt .
In practice we will choose each window g n centered around a time a n , and it will typically be constructed by translating a well localized window centered around 0 by a n , as in the standard Gabor scheme, but with the possibility to vary the window g n for each position a n . Thus the sampling of the time-frequency plane is done on a grid which is irregular over time, but regular over frequency. Figure 2 shows an example of such a sampling grid. It can be noted that some results exist in Gabor theory for semiregular sampling grids, as for example in [START_REF] Casazza | Gabor frames over irregular lattices[END_REF]. Our study here uses a more general setting, as the sampling grid is in general not separable, and more importantly, the window can evolve over time.

In this case, the coefficients of the decomposition are given by:

c m,n = f, g m,n ,
and the frame operator is given by:

Sf = m n f, g m,n g m,n .
The frame operator can be described by its kernel K given the following relation, which holds at least in a weak sense:

Sf (s) = K(t, s)f (t)dt.
Here the kernel K simplifies according to the following relations:

K(t, s) = m n g n (t)g n (s)e i2πmbn(s-t) = n g n (t)g n (s) m e i2πmbn(s-t) = n 1 b n g n (t)g n (s) k δ s -t - k b n thus, Sf (s) = k n 1 b n g n s - k b n g n (s) f s - k b n
In general, the inversion of S is not obvious. However we can identify a special case, which is analog to the "painless" case in standard Gabor analysis [START_REF] Gröchenig | Foundations of Time-Frequency Analysis[END_REF], for which the expression of S simplifies. More precisely, we suppose from now on that for every n ∈ Z, the function g n has a limited time support

supp g n = [c n , d n ] such that d n -c n < 1 bn .
Due to this support condition, the terms of the summation over k in the preceding equation are 0 for k = 0 and the frame operator S becomes a multiplication operator:

Sf (s) = n 1 b n |g n (s)| 2 f (s).
In this case the invertibility of the frame operator is easy to check and the system of functions g m,n forms a frame for L 2 (R) if and only if ∀t ∈ R, n

1 bn |g n (t)| 2 ≃ 1.
When this condition is fulfilled, the canonical dual frame elements are given by:

gm,n (t) = g n (t) k 1 b k |g k (t)| 2 e i2πmbnt
, and the associated canonical tight frame elements can be calculated by:

ġm,n (t) = g n (t) k 1 b k |g k (t)| 2 e i2πmbnt .

Resolution evolving over frequency

An analog construction is possible with a sampling of the time-frequency plane irregular over frequency, but regular over time. An example of the sampling grid in such a case is given on Figure 3.

In this case, we introduce a family of functions {h m } m∈Z of L 2 (R), and for m ∈ Z and n ∈ Z, we define atoms of the form: h m,n (t) = h m (tna m ). In practice we will choose each function h m as a well localized pass-band function having a Fourier transform centered around some frequency b n .

In this case the frame operator is given by:

Tf = m n f, h m,n h m,n ,
and the problem is completely analog to the preceding up to a Fourier transform, as we have:

Tf = m n f , h m,n h m,n ,
and h m,n = h m (ν)e -i2πnamν . So the preceding computation can be done, working on the Fourier transforms of the involved functions instead of directly on the functions. Now the "painless" case appears when we suppose that for every m ∈ Z, the function h n has a limited frequency support supp h n = [e n , f n ] such that f ne n < 1 an . Then the following expression holds:

Tf (ν) = m 1 a m | h m (ν)| 2 f (ν),
and the system of functions h m,n forms a frame of L 2 (R) if and only if ∀ν ∈ R, n

1 am | h m (ν)| 2 ≃ 1.
The associated canonical dual and tight frame can be computed as preceding, with the addition of an inverse Fourier transform.

Implementation

For the practical implementation, we have developed the equivalent theory in a finite discrete setting, that is to say working with complex vectors as signals. This theory won't be described here due to lack of space, but the construction is very similar to the one described in 2.1 and 2.2 and leads to a frame matrix which simplifies to a diagonal matrix in the "painless" case, suitable for applications. The implementation is then very similar to the implementation of the standard Gabor case and can exploit fast Fourier transform algorithms for efficiency. The only differences compared to standard Gabor implementation are due to the fact that the storage of coefficients requires more advanced storage structures due to the irregularity of the time-frequency sampling grid, and that the computation of the dual window must be performed for every time position resulting in a slight increase in computational cost.

Example

The possibility to build a decomposition with timefrequency resolution evolving over time can be exploited to solve the problem described in example of Section 1. illustrated by Figure 1. For the corresponding glockenspiel signal, as we have seen before, the use of narrow window is suitable for the attacks of the notes, while a wide window should be used for the resonances. Figure 4 shows a representation built with our approach using a narrow window of 6 ms for the attacks and a wide window of 93 ms for the resonance. The frame used for this figure is a tight frame. It should be noticed that the evolution of the window size between the two target window lengths is smoothed in order to ensure that the atoms used for the decomposition maintain a "nice" shape, in the sense of having a good time-frequency concentration. This ensures the easy interpretability of the decomposition, especially for processing using frame multipliers. This figure gives an idea of the type of decompositions that can be constructed with our approach and should be compared to the decomposition obtained using standard Gabor analysis on Figure 1. With our approach, it becomes possible to have a simultaneous good representation of both types of components of this signal while keeping the same processing ability than with standard Gabor. We see that our approach allows to build decompositions with better time-frequency localization of the signal energy. This can be helpful for many processing tasks, in particular to reduce artifacts in component extraction or denoising.

Conclusion

Our approach enables the construction of frames with flexible evolution of time-frequency resolution over time or frequency. The resulting frames are well suited for applications as they can be implemented using fast algorithms, at a computational cost close to standard Gabor frames. Exploiting evolution of resolution over time, the proposed approach can be of particular interest for applica- tions where the frequency characteristics of the signal are known to evolve significantly with time. Order analysis [START_REF] Shao | Order tracking by discrete Gabor expansion[END_REF], in which the signal analyzed is produced by a rotating machine having evolving rotating speed, is an example of such application. Exploiting evolution of resolution over frequency, the presented approach could be valuable for applications requiring the use of a tailored non uniform filter bank. In particular, it can be used to build filter banks following some perceptive frequency scale. One difficulty when using our approach is to adapt the time-frequency resolution to the evolution of the signal characteristics. If prior knowledge is available, this can be done by hand, as for the example of Figure 4. But to go further, our approach could be extended to construct an adaptive decomposition of the signal by automatically adapting the resolution to the signal. To achieve this, we plan to investigate the possibility to couple our approach with the use of sparsity criterion as proposed in [START_REF] Jaillet | Time-frequency jigsaw puzzle: adaptive multiwindow and multilayered gabor representations[END_REF]. The general idea would then be to consider time segments of the signal, and for each time segment compare the sparsity criterion obtained for Gabor transforms computed with different possible windows. We would then use in our decomposition the window corresponding to the best criterion for each time segment, leading to a decomposition optimizing the sparsity of the decomposition over time.

Frequency

Figure 1 :

 1 Figure 1: Two spectrograms of the same glockenspiel signal obtained using two different window lengths. On the top plot, a narrow window of 6 ms is used, on the bottom plot, a wide window of 93 ms is used.

TimeFrequencyFigure 2 :

 2 Figure 2: Example of sampling grid of the time-frequency plane when building a decomposition with time-frequency resolution evolving over time.

Figure 3 :

 3 Figure 3: Example of sampling grid of the time-frequency plane when building a decomposition with time-frequency resolution evolving over frequency.

FrequencyFigure 4 :

 4 Figure 4: Spectrogram of the same glockenspiel signal as in Figure 1 using a nonstationary Gabor decomposition.

Acknowledgment

This work was supported by the WWTF project MULAC ("Frame Multipliers: Theory and Application in Acoustics", MA07-025).