Marie-Lou Barnaud
email: marie-lou.barnaud@gipsa-lab.grenoble-inp.fr

Julien Diard

Pierre Bessière

Jean-Luc Schwartz

Assessing Idiosyncrasies in a Bayesian Model of Speech Communication

Keywords: speech development, motor learning, Bayesian modeling, idiosyncrasies

come   L'archive ouverte pluridisciplinaire

Introduction

Although speech acquisition is fast and efficient, the mechanisms underlying speech development are quite complex. If we only consider phonetic learning occurring during the first year of life, it can be decomposed in three steps [START_REF] Kuhl | Early language acquisition: Cracking the speech code[END_REF][START_REF] Kuhl | Phonetic learning as a pathway to language: New data and native language magnet theory expanded (NLM-e)[END_REF]. First, from birth, children learn to associate sounds with the phonemes of their native language. Then, from around seven months, a babbling phase occurs during which children learn to associate acoustic signals with motor gestures. Finally, around two months later, children begin to associate motor gestures with the phonemes of their native language.

These three learning steps, respectively called sensory, sensory-motor and motor learning in the following of this paper, are language specific. Indeed, the exposure to one particular language results in tuning the sensory and motor phonetic representations to this language (in the perceptual domain, this is called perceptual narrowing [START_REF] Werker | Influences on infant speech processing: Toward a new synthesis[END_REF]). As a consequence, children speaking different languages have different phonetic representations.

Conversely, we may expect children speaking the same language to have similar phonetic repertoires. However, there is also intra-language variability, called idiosyncrasies. Typically, in speech production, when two agents produce the same phoneme, acoustic results may vary extensively [START_REF] Kuhl | Early language acquisition: Cracking the speech code[END_REF][START_REF] Ménard | Invariance and variability in the production of the height feature in French vowels[END_REF].

In this paper, we focus on the development of idiosyncrasies in speech production and aim at better understanding what component of the learning process could be at their origin. Since idiosyncrasies in production concern the relationship between motor gestures and phonemes, we assume that they appear during the motor learning phase. We compare two computational models of this phase of speech development, both based on an imitation algorithm during which a computational learning agent tries to reproduce speech utterances of a master agent. In the first model, named "repetition model", the agent tries to reproduce sounds uttered by the master. In the second one, named "communication model", the agent tries to replicate phonemes produced by the master.

Our two motor learning algorithms are embedded inside a Bayesian model of speech communication called COSMO (for "Communicating Objects using Sensory-Motor Operations"), that we have been developing in the past years. COSMO is in our view an efficient framework to study and simulate various aspects of speech communication, including the emergence of sounds systems in human languages [START_REF] Moulin-Frier | Rôle des relations perception-action dans la communication parlée et l'émergence des systèmes phonologiques: étude, modélisation computationnelle et simulations[END_REF][START_REF] Moulin-Frier | Communicating about Objects using Sensory-Motor Operations"): A Bayesian modeling framework for studying speech communication and the emergence of phonological systems[END_REF] or online speech perception [START_REF] Moulin-Frier | Adverse conditions improve distinguishability of auditory, motor, and perceptuo-motor theories of speech perception: An exploratory Bayesian modelling study[END_REF][START_REF] Laurent | A computational model of perceptuo-motor processing in speech perception: Learning to imitate and categorize synthetic CV syllables[END_REF].

This paper is organized as follows: Section 2 presents the COSMO model and describes the two motor learning models. Section 3 compares results of experimental simulations with the two learning models, which are then discussed in Section 4.

COSMO, a Bayesian model of speech communication 2.1. Model description

Within a speech communication process between two agents, a speaker produces motor gestures, that result in acoustic signals perceived by a listener; this enables an exchange of linguistic information between the two agents. From this conceptual description of the communication process, the COSMO model relies on the assumption that communicative agents internalize in their brain all the involved motor, sensory and linguistic representations. In COSMO, these representations are modeled by probabilistic variables: M for motor gestures, S for sensory (acoustic) signals, OS and OL for the linguistic "objects" (in a general sense) of communication, OS relating to the object for the speaker and OL to the object for the listener, and C for the evaluation of communication success.

Based on the Bayesian Programming methodology [START_REF] Lebeltel | Bayesian robot programming[END_REF][START_REF] Bessière | Bayesian Programming[END_REF], the joint probability distribution P (C OS S M OL) is decomposed as a product of five distributions: a prior on objects P (OS), a motor system P (M | OS), a sensory-motor system P (S | M ), an auditory recognition system P (OL | S) and a communication validation system P (C | OS OL). These five distributions are the knowledge of our communicating agent.

In this study, we implement a "vowel version" of the COSMO model. It involves the use of an articulatory model of the vocal tract, VLAM [START_REF] Maeda | Compensatory articulation during speech: Evidence from the analysis and synthesis of vocal-tract shapes using an articulatory model[END_REF][START_REF] Boë | Modélisation de la croissance du conduit vocal[END_REF][START_REF] Ménard | Auditory normalization of French vowels synthesized by an articulatory model simulating growth from birth to adulthood[END_REF] in which orofacial articulators (jaw, larynx, tongue, lips) are controlled by 7 parameters.

In our model, linguistic units OS and OL correspond to the seven vowels /i u e o E O a/, which are the seven preferred vowels in human languages [START_REF] Schwartz | The dispersion-focalization theory of vowel systems[END_REF]. The motor variable M only retains three parameters of VLAM sufficient for these vowels, that are lip height LH , tongue body TB, and tongue dorsum TD, respectively monitoring vowel rounding, vowel height and vowel anterior/posterior configurations. The sensory variable S consists of formants F 1 and F 2 expressed in Barks [START_REF] Schroeder | Objective measure of certain speech signal degradations based on masking properties of human auditory perception[END_REF]. We discretize F 1 and F 2 respectively into 59 and 73 values, while M contains 15*15*15 values. C is a boolean value, expressing that OL and OS are identical or different.

We define the probability distributions of the model. P (OS) is a uniform distribution: all vowels are equiprobable. P (S | OL), P (S | M ) and P (M | OS) are conditional Gaussian distributions. To express the lack of knowledge before learning, these distributions are initially set with means in the middle of their space and large variance, approximating uniform distributions. Learning consists in providing values for objects, sensory and motor variables (e.g. o, s and m) in ways that will be explained later. From these values, parameters of the Gaussian distributions P (S | OL), P (S | M ) and P (M | OS) are updated in a straightforward manner respectively using observed data s, o , s, m and m, o . Finally, P (C | OS OL) is a "Bayesian switch" [START_REF] Gilet | Bayesian actionperception computational model: Interaction of production and recognition of cursive letters[END_REF]: when C = 1, OS and OL are constrained to the same value.

We previously showed how respectively setting OS or OL or both OS and OL as the pivot of communication enabled to switch from a motor to an auditory to a sensory-motor theory of speech communication [START_REF] Laurent | A computational model of perceptuo-motor processing in speech perception: Learning to imitate and categorize synthetic CV syllables[END_REF][START_REF] Barnaud | COSMO, a Bayesian computational model of speech communication: Assessing the role of sensory vs. motor knowledge in speech perception[END_REF]. In this paper, we keep the most general framework, that is a sensory-motor theory of speech production, so that OL and OS are always constrained (by C = 1) to be equal. Hence, to simplify notations, in the following, we note both OL and OS by a single O. This particularly concerns processes in the motor phase (see Eq. ( 3)) and processes used for the evaluation (see equations in Section 2.3).

Learning phases

Starting from scratch, we consider the three speech development stages previously introduced: a sensory learning phase associating sounds with phonemes, a sensory-motor learning phase associating motor gestures with sounds, and a motor learning phase associating motor gestures with phonemes. In agreement with other works [START_REF] Kuhl | Early language acquisition: Cracking the speech code[END_REF][START_REF] Kuhl | Phonetic learning as a pathway to language: New data and native language magnet theory expanded (NLM-e)[END_REF], we consider that these steps are consecutive and performed in interaction with a master agent.

Master agent

The master agent we use in this study disposes of a set of target motor commands for each vowel. These target sets have been defined so as to produce typical formant values for the seven considered vowels, based on data for French vowels [START_REF] Meunier | Phonétique acoustique[END_REF]. For each vowel, the master agent draws values for M according to a Gaussian distribution around the motor target, with a given variance in the articulatory space. Motor commands are then translated into acoustic values thanks to VLAM. This provides a (non Gaussian) distribution P (S|Omast) from which the master draws samples provided to the learning agent.

Sensory and sensory-motor phases

During sensory learning, the agent learns its probability distribution P (S | OL). This learning phase is straightforward: the master produces a linguistic object o resulting in an acoustic signal s, and we assume that the learning agent is able to access both s from its auditory system and o from a given parallel communication stream, e.g. deixis [START_REF] Moulin-Frier | Communicating about Objects using Sensory-Motor Operations"): A Bayesian modeling framework for studying speech communication and the emergence of phonological systems[END_REF]. The learning agent then directly updates its distribution

P ([S = s] | [OL = o]) thanks to the o, s couple.
During sensory-motor learning, the agent learns its probability distribution P (S | M ). This phase is a little more complex: as the master agent cannot directly inform the learning agent about the motor gestures it produces, the learning agent needs to infer them. We suppose that inference is based on an imitation process. As in the sensory phase, the master produces a linguistic object o resulting in an acoustic signal s. Then, the learning agent tries to imitate the master by inferring a motor gesture m thanks to the distribution P (M | [S = s]). The selection of a given motor command m results in the production of a sound s (computed thanks to VLAM). Of course, s has no reason to be equal to the target sound s provided by the master. However, the agent exploits this s , m pair to update its sensory-motor system P (

[S = s ] | [M = m]).

Motor phase

Once the sensory-motor learning phase is completed, the motor learning phase begins. During this phase, the learning agent updates its distribution P (M | OS). Although it uses an imitation process similar to the sensory-motor phase, the inference process is different.

We consider two versions of this inference process. In the first version, called "repetition model", the agent attempts to reproduce the exact sound produced by the master for a given object. For this aim, inference is based on the distribution

P (M | [O = o] [S = s]
), which means: select a motor gesture likely to be associated to the phoneme o and to result in the sound s. In the second version, called "communication model", the agent tries to select a motor gesture likely to ensure communication and hence to realize a vowel o similar to the one produced by the master. For this aim, inference is based on the distribution

P (M | [O = o] [C = 1]).
More formally, both distributions

P (M | [O = o] [S = s]) and P (M | [O = o] [C = 1]
) are computed in the COSMO model using Bayesian inference, which yields:

P (M | [O = o] [S = s]) ∝ (1) P (M | [OS = o])P ([S = s] | M ) , P (M | [O = o] [C = 1]) ∝ (

Summary of the complete learning sequence

For each motor learning model, we performed 12 simulations in which each learning phase lasted 300,000 steps. Due to random sampling, simulations differed in couples s, o given by the master and in motor gestures m (and resulting s ) selected at each learning step. This enabled to test whether different simulations would result in different final stages at the end of the whole learning process, which could possibly provide idiosyncrasies.

Model evaluation

At the end of the whole learning process, models are evaluated in two ways, assessing both communication performance and possible motor and sensory idiosyncrasies. To assess communication performance, the learning agent tries to communicate an object Oag to the master agent, by producing motor commands resulting in sounds from which the master infers Omast. We compute the confusion matrix P (Omast | Oag):

P (Omast | Oag) = S (P (Omast | S)P (S | Oag)) , ( 3
)
where P (Omast | S) is the perceptual categorization system of the master, while P (S | Oag), the sensory result of the productions of the learning agent, is computed by:

P (S | Oag) ∝ M (P (S | M )P (M | Oag)) . (4)
Here, P (S | M ) is the real motor-to-acoustic transformation provided by VLAM, and P (M | Oag) is the production process of the learning agent.

Results

Communication performance

We computed the confusion matrix P (Omast | Oag) (Eq. ( 3)), at the end of the learning process for each of the 12 simulations for each motor learning model. A global communication performance index was provided by the mean proportion of correct answers for all phonemes, that is the average value of the diagonal of the confusion matrix. The average over the 12 simulations provides 99.1 % of correct recognition in the "communication model" and 98.4 % in the "repetition model". Those two values are quite close and both indicate high performance, illustrating that both motor learning models are able to correctly learn the phoneme repertoires of their master.

To further analyze our results, let us first display the distribution P (S | Omast) of the master. Figure 1 (left) provides the classical distribution of reference acoustic data [START_REF] Meunier | Phonétique acoustique[END_REF], where each vowel covers a unique portion of the acoustic space, though with some small overlap at their boundaries.

We also display the distributions P (S | Oag) for a typical simulation of one learning agent (see Eq. ( 4)) at the end of learning. The middle and right plots of Figure 1 respectively show an instance of P (S | Oag) in the "communication" and the "repetition" models. We notice that in both cases, vowels are well defined and distinguishable. However, we notice that while the "repetition" model on the right reproduces the master

Idiosyncrasies

On Figure 2 (left), we display both P (Omast | S), i.e. the categorization regions of the master, and the means of P (S | Omast), i.e. the sensory prototypes of phonemes (see Section 2.2.1)). As expected, prototypes of the master are well centered in each categorization region. This describes the way the sensory space is structured by the distribution of vowels in the master space, acting as a reference for the learning agent.

From this basis, the other plots of Figure 2 show how the 12 simulations of the "communication model" (middle) and the "repetition model" (right) compare to the stimuli provided by the master at the end of the learning stage. These displays were obtained by computing the means of P (S | [O = o]) (see Eq. ( 4)) for each vowel o in each of the 12 simulations, for the "communication model" and "repetition model". The 12 corresponding means are shown as colored dots, keeping the master categorization regions as reference.

We observe that idiosyncrasies appear only in the "communication model". Indeed, only in this case do the 12 mean values of P (S | Oag) vary between simulations. Importantly, despite these idiosyncrasies, the means of each vowel are still in their respective categorization regions, supporting the idea that idiosyncrasies do not alter perceptual categorization, and thus do not alter communication efficiency -as indeed shown by the measured communication performance in the previous section.

In the "repetition model", in contrast, there are no idiosyncrasies: vowel means are not variable from one simulation to the other, and are concentrated around the means of stimuli provided by the master distributions P (S | Omast). Indeed, it can be mathematically shown that, in this learning algorithm, P (S | Oag) progressively converges towards P (S | Omast).

Notice that, even if there are no sensory idiosyncrasies, the "many-to-one" relation from motor to sensory spaces may generate motor idiosyncrasies, since a given sensory percept can result from various different motor gestures. As a matter of fact, we display in Figure 3 distributions P (M | OS) in the motor space for two simulations of the "repetition model". Motor distributions are clearly different. Detailed analyses of simulation results confirm that such motor idiosyncrasies appear in both the "communication" and "repetition" models, even though sensory idiosyncrasies appear only in the first case.

Discussion

In this paper, we compared two versions of the motor learning stage in speech development, to investigate idiosyncratic learning in speech production: a "communication model" and a "repetition model". For this aim, we implemented a sequence of learning steps proposed by specialists of speech development [START_REF] Kuhl | Early language acquisition: Cracking the speech code[END_REF][START_REF] Kuhl | Phonetic learning as a pathway to language: New data and native language magnet theory expanded (NLM-e)[END_REF] into the COSMO model. Our first experimental result is that, in the scope of the phonetic material considered in this paper and involving a small set of oral vowels, COSMO is able to correctly produce learned phonemes whatever the version used.

The second and main result of this study is that idiosyncrasies are only obtained in the "communication model" of motor learning. Since idiosyncratic behaviors are a commonly observed phenomenon, we infer that speech development likely involves some motor learning process guided by a communicative goal, during which children would try to replicate perceived phonemes rather than perceived sounds. Such learning process based on a communicative goal could actually take a wide variety of forms, including communication scenarios based on inverse imitation games (see, e.g. [START_REF] Messum | Creating the cognitive form of phonological units: The speech sound correspondence problem in infancy could be solved by mirrored vocal interactions rather than by imitation[END_REF]).

The sequence of learning stages within speech development that we considered in the present study could be embedded within a more general scenario based on hierarchical learning, with a first stage guided by sensory representations (our sensory and sensory-motor phases), followed by a second, higher-level stage guided by phonetic representations (our motor phase).

Our model has several limitations. Just to mention one, we only considered learning interaction with a single master, which is unrealistic for child speech development. Simulations with several masters are likely to provide idiosyncrasies also in the "repetition model". However, such idiosyncrasies would be centered on the average of the different masters' productions, and iteration of this process over generations would likely gradually reduce the spread of idiosyncrasies. It is not sure whether that would reflect realistic idiosyncrasies.

Whatever the obvious limitations of this initial study, we believe that the proposed strategy -based on the comparison of different computational architectures within a single computational framework -is promising, in order to assess the role of specific components of the general speech communication model we are aiming at here. The specific component that is tested here, that is, the existence of a learning process based on efficient communication, appears to be a potentially important ingredient in the future developments of COSMO.

  2) P (M | [OS = o]) S (P (S | M )P ([OL = o] | S)) , where P (M | OS), P (S | M ), and P (OL | S) are probability distributions of the learning agent. In both versions, the inferred motor gesture m is used to update parameters of the motor system P ([M = m] | [OS = o]).

Figure 1 :

 1 Figure 1: Vowel distributions P (S | O). Plots use the classical view of the acoustic space, with F1 on the y-axis, F2 on the x-axis, both reversed. Axes values are in Barks. High probabilities are in red, low probabilities in blue. Each region with a color scale from green-yellow to red represents a vowel. Left: Distributions P (S | Omast) of the master; Middle: Learned distributions P (S | Oag) in the "communication model"; Right: Learned distributions P (S | Oag) in the "repetition model".

Figure 2 :

 2 Figure 2: Acoustic space is represented as in Figure 1, with values in Barks. Ellipses in the three plots correspond to the categorization regions of the master (distribution P (Omast | S)). Points respectively correspond to the means of: Left: the master distributions P (S | Omast); Middle: the distributions P (S | Oag) in the "communication model" for the 12 simulated agents; Right: the distributions P (S | Oag) in the "repetition model" for the 12 simulated agents.

Figure 3 :

 3 Figure 3: Comparison of the motor distribution P (M | OS) for two simulations of learning agents in the "repetition model": tongue body (TB) on the x-axis, tongue dorsum (TD) on the y-axis and lip height (LH ) on the z-axis. Axes values are based on VLAM values. Points in the same color correspond to the same vowel.

Acknowledgements

Research supported by a grant from the European Research Council (FP7/2007-2013 Grant Agreement no.

339152, "Speech Unit(e)s").


Zhouye Chen

Adrian Basarab

Denis Kouamé

A simulation study on the choice of regularization parameter in ℓ 2 -norm ultrasound image restoration*

Ultrasound image deconvolution has been widely investigated in the literature. Among the existing approaches, the most common are based on ℓ2-norm regularization (or Tikhonov optimization) or the well-known Wiener filtering. However, the success of the Wiener filter in practical situations largely depends on the choice of the regularization hyperparameter. An appropriate choice is necessary to guarantee the balance between data fidelity and smoothness of the deconvolution result. In this paper, we revisit different approaches for automatically choosing this regularization parameter and compare them in the context of ultrasound image deconvolution via Wiener filtering. Two synthetic ultrasound images are used in order to compare the performances of the addressed methods.

I. INTRODUCTION

Ultrasound (US) medical imaging has the advantages of being noninvasive, harmless, cost-effective and portable over many other imaging modalities such as X-ray computed tomography or Magnetic Resonance Imaging [START_REF] Szabo | Diagnostic ultrasound imaging: inside out[END_REF]. However, the limited bandwidth of the imaging transducer, the characteristics of the ultrasound propagation such as diffraction and the imaging system tend to degrade the quality of the US images. To deal with this problem, significant efforts have been made in the last few decades.

Under some weak assumptions (first order Born approximation and weak scattering), the recorded US radiofrequency image can be modeled as the result of a 2D convolution between the tissue reflectivity function and the point-spread function (PSF) [START_REF] Jensen | Deconvolution of in vivo ultrasound b-mode images[END_REF]. As a consequence, deconvolution methods have been intensively considered to enhance the spatial resolution of US images. Among the existing approaches, Wiener filtering (and variants) was one of the most explored tracks in US imaging [START_REF] Michailovich | Blind deconvolution of medical ultrasound images: A parametric inverse filtering approach[END_REF], [START_REF] Jirik | Two-dimensional blind bayesian deconvolution of medical ultrasound images[END_REF]. However, it is well-known that the results of the Wiener filter or its variants largely depends on the choice of the regularization parameter (RP) that provides a compromise between data fidelity and smoothness of the deconvolution result. In US imaging, the choice of the RP is either done manually (empirically or by trial-and-error) [START_REF] Michailovich | Blind deconvolution of medical ultrasound images: A parametric inverse filtering approach[END_REF], or is simply related to the signal-to-noise ratio (the RP is considered the inverse of the SNR) which is further estimated from the data [START_REF] Jirik | Two-dimensional blind bayesian deconvolution of medical ultrasound images[END_REF]. In the general literature related to Wiener filter deconvolution, several approaches have been proposed to find an optimal value of RP in an automatic manner. In this paper, we propose to evaluate and compare their performances in the framework of ultrasound imaging. To do so, two simulated US images are used.

The remainder of the paper is organized as follows. The considered model and the Wiener filtering approach are presented in Section II. Section III reviews the existing approaches of optimal RP choice. The simulation setup, the deconvolution results and the comparison in US imaging are provided in Section IV.

II. PROBLEM FORMULATION

Under the assumption of weak scattering and using the first order Born approximation, the interaction between the tissues and the propagating ultrasound waves is classically modeled by the following 2D convolution model [START_REF] Jensen | Deconvolution of in vivo ultrasound b-mode images[END_REF]:

y = Hx + n (1)
where y ∈ R M N ×1 is the vertical concatenation of M acquired radio-frequency (RF) signals of length N (also known as "lexicographical notation" of the RF image), x ∈ R M N ×1 is the tissue reflectivity function using the same lexicographical notation, H ∈ R M N ×M N is a block circulant with circulant block (BCCB) matrix related to the 2D PSF of the system and n ∈ R M N ×1 is a zero-mean additive white Gaussian noise with variance σ 2 . The purpose of non-blind deconvolution methods is to recover the tissue reflectivity function x from the recorded US RF image y considering the PSF known (or provided a previous estimate of the PSF). A common way to solve the non-blind deconvolution problem is the well-known Wiener/Tikhonov approach, that provides the following analytical estimation of x:

x = (H t H + λQ t Q) -1 H t y = F (λ)y (2)
where x is the estimate of x, and Q ∈ R M N ×M N is a BCCB matrix representing the regularization operator, classically considered the identity matrix or the 2D Laplacian operator. λ is a hyperparameter providing compromise between data fidelity and smoothness of the estimate (addressed as regularization parameter (RP) hereafter). The value of λ that guarantees an optimal compromise between data fidelity and smoothness depends on the variance of the noise σ 2 and on the properties of H, Q and y [START_REF] Golub | Generalized cross-validation as a method for choosing a good ridge parameter[END_REF].

III. METHODS FOR REGULARIZATION PARAMETER OPTIMAL CHOICE

In the general Wiener filtering literature, several approaches have been proposed to automatically fix the value of the RP λ [START_REF] Galatsanos | Methods for choosing the regularization parameter and estimating the noise variance in image restoration and their relation[END_REF]. We review hereafter the main approaches, that we further evaluate in the context of US images. Following [START_REF] Galatsanos | Methods for choosing the regularization parameter and estimating the noise variance in image restoration and their relation[END_REF], we classify them into two categories. The first includes methods that need the knowledge of the variance of the noise (or an estimate of it) to estimate the RP, while the second one does not have this constraint.

A. Methods with knowledge of the variance of the noise 1) Constrained least squares (CLS):

The principle of this method, originally proposed in [START_REF] Hunt | The application of constrained least squares estimation to image restoration by digital computer[END_REF], is based on the residual between the data y and H x for a given λ, expressed as:

φ(λ) = y -H x 2 2 = y -HF (λ)y 2 2 (3)
The optimal λ in the sense of CLS is obtained in the case where the residual in ( 3) is equal to the variance of the noise, or more precisely to M N σ 2 . From this it results an equation with the unknown λ, which was shown in [START_REF] Hunt | The application of constrained least squares estimation to image restoration by digital computer[END_REF] to have an unique solution. This optimal value of λ, that we will denote by λ CLS in this paper, is obtained by solving the above equation using a numerical iterative procedure.

2) Degree of freedom (EDF): Similar to CLS approach, the equivalent degree of freedom method [START_REF] Wahba | Bayesian" confidence intervals" for the cross-validated smoothing spline[END_REF] estimated the RP λ from an equation relating the residual and the variance of the noise. Based on the linear relation between y and x, the residual in ( 3) is further expressed as:

φ(λ) = σ 2 trace[I M N -HF (λ)] (4)
where

I M N ∈ R M N ×M N is the identity matrix.
The optimal RP given by EDF approach will be denoted hereafter by λ EDF .

3) Mean square error (MSE): The main idea behind this method is to find the value of λ that minimizes the mean square error between x and its estimate, denoted by ǫ(λ):

E( ǫ(λ) 2 2 ) = E[ x-x 2 2 ] = x 2 2 +E[ x 2 2 ]-2E[x t x]
(5) In order to eliminate the cross-term 2E[x t x], which is impossible to compute in practice because of the nonknowledge of x, it is replaced in [START_REF] Galatsanos | Methods for choosing the regularization parameter and estimating the noise variance in image restoration and their relation[END_REF] by a term depeding on the data y and on the noise variance σ 2 , based on a circulant assumption in DFT domain. In this case, it is shown in [START_REF] Galatsanos | Methods for choosing the regularization parameter and estimating the noise variance in image restoration and their relation[END_REF] that minimizing [START_REF] Golub | Generalized cross-validation as a method for choosing a good ridge parameter[END_REF] with espect to λ is equivalent to soving the following equation:

Q -1 (I M N -HF (λ)) 3 2 y 2 2 = σ 2 trace[Q -2 (I M N -HF (λ)) 2 ] (6)
We denote by λ M SE the value of λ verifying [START_REF] Galatsanos | Methods for choosing the regularization parameter and estimating the noise variance in image restoration and their relation[END_REF] and obtained in practice with numerical optimization techniques.

4) Predicitve mean square error (PMSE): Similar to the previous approach, this method is also based on the minimization of the MSE between x and its estimate. More precisely, a weighted ℓ 2 error norm E[ Hǫ(λ) 2 2 ], also called the predictive mean square error, is minimized [START_REF] Golub | Generalized cross-validation as a method for choosing a good ridge parameter[END_REF], [START_REF] Hall | Common structure of techniques for choosing smoothing parameters in regression problems[END_REF], [START_REF] Ramani | Regularization parameter selection for nonlinear iterative image restoration and mri reconstruction using gcv and sure-based methods[END_REF]. Following the same assumptions as the previous approach, this minimization turns out in solving the equation hereafter:

E[ Hǫ(λ) 2 2 ] = (I M N -HF (λ))y 2 2 +2σ 2 trace[HF (λ)] (7)
The solution of [START_REF] Hunt | The application of constrained least squares estimation to image restoration by digital computer[END_REF], obtained by numerical optimization, is called λ P M SE in what follows.

5) Generalized Stein's unbiased risk estimate (GSURE):

The main idea behind this method is to find the RP that minimizes the following error measure between x and its estimate [START_REF] Eldar | Generalized sure for exponential families: Applications to regularization[END_REF], [START_REF] Giryes | The projected gsure for automatic parameter tuning in iterative shrinkage methods[END_REF]:

e(λ) = P (x -x) 2 (8)
where P = H t (HH t ) † H is the projection operator and (•) † represents the pseudo-inverse. Keeping in mind that the Wiener filter, given a value of λ, provides an analytical relation between x and y, we denote by f a function from R M N to R M N with x = f (y). As shown in [START_REF] Ramani | Regularization parameter selection for nonlinear iterative image restoration and mri reconstruction using gcv and sure-based methods[END_REF], [START_REF] Eldar | Generalized sure for exponential families: Applications to regularization[END_REF] and [START_REF] Giryes | The projected gsure for automatic parameter tuning in iterative shrinkage methods[END_REF], an unbiased estimator for the projected MSE in [START_REF] Wahba | Bayesian" confidence intervals" for the cross-validated smoothing spline[END_REF], depending on the variance of the noise and separating the true x from its estimation f (y), is given by: e(λ) = P x 2 + P f (y) 2 2 -2f t (y)H t (HH t ) † y +2σ 2 trace(P

(H t H + λQ t Q) -1 ) (9)
Minimizing ( 9) with respect to λ provides the optimal value of the RP, in the sense of GSURE approach, denoted by λ GSU RE .

B. Methods without knowledge of the variance of the noise

The methods presented in the previous section are all based on the knowledge of the variance of the noise. Depending on the application, this may be an issue and a bad estimation of σ 2 may cause severe errors in the optimal choice of the RP.

In the following, we briefly describe two existing approaches that do not use σ 2 to provide an automatical choice of λ.

1) Generalized cross-validation (GCV): Generalized cross-validation (GCV) is one of the most popular methods of choosing optimal RP and does not require the knowledge of the noise variance σ 2 . Introduced in [START_REF] Golub | Generalized cross-validation as a method for choosing a good ridge parameter[END_REF], its main idea is based on the "leave-one-out" principle. In the case of linear algorithms, the GCV method leads to a simple function to minimize depending on λ:

GCV (λ) = (I M N -HF (λ))y 2 2 [trace(I M N -HF (λ))] 2 (10)
The value obtained by numerical minimization of ( 10) is denoted by λ GCV in this paper.

2) Marginal likelihood (ML): This approach is based on the minimizing of the marginal likelihood (ML) function given hereafter [START_REF] Galatsanos | Methods for choosing the regularization parameter and estimating the noise variance in image restoration and their relation[END_REF]:

M L(λ) = y t (I M N -HF (λ))y (det[I M N -HF (λ)]) 1/M N (11)
Hereafter, we denote the λ minimizing ( 11) by λ M L .

IV. RESULTS AND DISCUSSION

A. Experimental Setting

To evaluate the performance of these different methods, we used two synthetic ultrasound images. The two images were simulated by 2D convolution between realistic PSFs and the tissue reflectivity functions shown as the first images in Fig. 1 and Fig. 2. Independent and identically distributed (IID) zero-mean Gaussian noise was added to the data, yielding different SNRs. This leads to the original B-mode images, as shown in Fig. 1 and Fig. 2. 1) Cyst phantom image: For the first simulated image, the PSF was generated with Field II [START_REF] Jensen | A model for the propagation and scattering of ultrasound in tissue[END_REF] and corresponds to a 3.5 MHz linear probe, sampled in the axial direction at 20 MHz. The tissue refelctivity function was obtained by generating scatterers at uniform random positions with random amplitudes following gaussian distribution with zero mean and variance depending on their spatial position. The medium consists in five hyperechoic circular cysts, five hypoechoic circular cysts and five point reflectors.

2) Simulated cardiac image: The PSF was also generated with Field II [START_REF] Jensen | A model for the propagation and scattering of ultrasound in tissue[END_REF] and corresponds to a sectorial probe with the central frequency equal to 4 MHz and a axial sampling frequency of 40 MHz. The scatterer positions were uniformly random distributed. In order to obtain an ultrarealistc simulation, the amplitudes of the scatterers were related to the amplitude of an in vivo cardiac image, as suggested in [START_REF] Alessandrini | Simulation of realistic echocardiographic sequences for ground-truth validation of motion estimation[END_REF].

B. Deconvolution Results

For both US simulations, the deconvolution results were obtained with the Wiener filtering approach, considering the true PSF known. The RP was chosen with one of the approaches in Section III or using the classical choice of λ equal to the inverse of the SNR. For the methods that need the knowledge of the SNR, the true value was employed. Two classical cases were considered using two different regularization operators Q: the identity matrix and the 2D Laplacian operator.

TableI shows the deconvolution results for each approch of RP optimal choice. The deconvolved images are compared to the true reflectivity function using the structural similarity measure (SSIM) [START_REF] Wang | Image quality assessment: from error visibility to structural similarity[END_REF]. For each simulation and for a given SNR, we highlight in bold fonts the best result. Note that the figures we give in Table I are obtained with a regularization operator equal to the identity matrix. The results follow the same trend for a Laplacian operator.

The different results may be appreciated from a qualitative viewpoint in Fig. 1 and2, highlighting the original refelctivity function, the data and the deconvolved images. All the processing was done in the RF domain. However, for visualisation reasons, we plot the corresponding B-mode images, obtained after log-compression of the corresponding enveloppe images.

C. Discussion

We may first remark that the different strategies of RP optimal choice provide very different values of λ and thus influence the quality of the deconvolution results. We should note that this conclusion is less true for higher SNRs (larger than 40 dB).

Second, as predicted in [START_REF] Golub | Generalized cross-validation as a method for choosing a good ridge parameter[END_REF], the results provided by PMSE and GCV approaches are very similar. Similarly, the values of λ M SE and λ GSU RE are also very close. This is due to the fact that both MSE and GSURE methods are based on the minimization of the MSE. Despite the different ways of minimizing the MSE, they still provide similar results. Third, we may remark that the approaches that are not based on the knowledge of the SNR (GCV and ML) are still providing interesting results in US image deconvolution compared to the other approaches. The quantitative and qualitative results that we obtained point out that the ML approach is more adapted to US imaging than GCV. For all the cases, ML is among the best three methods (or very close to the best three), with CLF and the inverse of the SNR, both of which require the knowledge of the variance of the noise.

V. CONCLUSIONS

In this paper, we have compared eight methods of choosing optimal RP in the context of US deconvolution via Wiener filtering. The comparison was performed on two synthetic US images. An interesting conclusion is that the approaches not using the knowledge of the variance of the noise are competitive against the others and are thus very attractive in practical situations. In a future work, these results should be confirmed on experimental data, combined with approaches of PSF and noise variance estimation.

Fig. 1 :

 1 Fig. 1: SNR=30dB. From left to right, the images are Cyst phantom tissue reflectivity function, Cyst phantom B-mode US image and its deconvolution results (B-mode visualisation) for Q = I.

Fig. 2 :

 2 Fig. 2: SNR=30dB. From left to right, the images are Cardiac tissue reflectivity function, Cardiac B-mode US image and its deconvolution results (B-mode visualisation) for Q = Q 2DL .

TABLE I :

 I Deconvolution Results of Synthetic Ultrasound Images

	SNR		10dB			20dB			30dB
	Image	Cyst Phantom	Cardiac Image	Cyst Phantom	Cardiac Image	Cyst Phantom	Cardiac Image
		Lambda	SSIM	Lambda	SSIM	Lambda	SSIM	Lambda	SSIM	Lambda	SSIM	Lambda	SSIM
	1/SNR	1,00E-01	82,05	1,00E-01	63,50	1,00E-02	86,88	1,00E-02	67,29	1,00E-03	91,61	1,00E-03	76,15
	CLS	2,12E-01	82,95	8,32E-01	65,32	1,60E-02	86,53	4,13E-02	67,93	1,81E-03	91,76	2,66E-03	75,24
	EDF	2,73E-02	76,67	3,48E-02	61,05	2,12E-03	83,74	1,91E-03	60,24	3,18E-04	89,49	2,22E-04	67,49
	MSE	2,37E-02	75,74	2,49E-02	40,46	2,25E-03	84,04	1,79E-03	59,69	4,09E-04	90,55	2,89E-04	70,56
	PMSE	3,41E-02	77,98	5,84E-02	63,73	2,02E-03	83,46	2,12E-03	61,02	2,06E-04	86,59	1,79E-04	64,27
	GSURE	2,37E-02	75,74	2,49E-02	40,46	2,25E-03	84,04	1,79E-03	59,69	4,09E-04	90,55	2,89E-04	70,56
	GCV	3,46E-02	78,06	6,05E-02	63,80	2,02E-03	83,46	2,13E-03	61,08	2,06E-04	86,59	1,79E-04	64,27
	ML	4,57E-02	79,46	2,79E-01	63,59	3,90E-03	86,04	1,51E-02	67,78	2,80E-04	88,65	7,00E-04	70,24


Designing and Assessing Interactive Systems Using Task Models

Keywords: User Interaction Design, Task description and modelling ACM Classification, H.5.m. Information interfaces and presentation (e.g., HCI)

Introduction

Task analysis is meant to identify user goals and tasks when using an interactive system. In the case of users performing real life work, task analysis can be a cumbersome process gathering a huge amount of unorganized information. Task Models (TM) provide a mean for the analyst to store information gathered in an abstract way that can be further detailed and analyzed if needed. A task model allows HCI researchers and practionners to record in a systematic, complete and unambiguous way the set of user goals and the way those user goals can be performed on an interactive system. Reasoning about the Task Models produced supports the assessment of effectiveness of an interactive system (which is one of the most difficult dimension of usability to assess). Task models have also proven being of great help for structuring user documentation, designing and assessing a training program, assessing the complexity of the users' work. If used for analysis, they can also provide support for identifying types, location and likelihood of human errors. When used for design they

Abstract

This two-part course takes a practical approach to introduce the principles, methods and tools in task modelling. Part 1: A non-technical introduction demonstrates that task models support successful design of interactive systems. Part 2: A more technical interactive hands-on exercise of how to "do it right", such as: How to go from task analysis to task models? How to assess (through analysis and simulation) that a task model is correct? How to identify complexity of user tasks … also provide precious support for identification of good candidates for task migration towards automation.

Contribution and benefit

This course intends to provide newcomers with background in task modeling. It provides an overview on how the recent advances in task description techniques can be exploited to design and assess interactive systems. As task models can be large, it is important to provide the analyst with computer-based tools for editing task models and for analyzing them. To this end, this course provides attendees with the HAMSTERS task modeling tool that can be directly applied in practice.

Objectives

On completion of this tutorial, attendees will:

-Know the benefits of using task modeling techniques to design, structure and assess UIs, -Be able to describe users' activities in a systematic and structured way, -Have experience in analyzing an interactive systems focusing on the tasks users have to perform with it, -Know how to use the HAMSTERS tool suite for editing, analyzing and simulating task models.

Instructors' background

The instructors have applied task modeling techniques to several industrial projects such as the design of collaborative environments to manage collision risks between satellite and space objects, to improve interactive satellite Ground Segment Systems, to assess impact of dependability on usability for interactive cockpits and in Air Traffic Management (http://www.irit.fr/recherches/ICS-site).

  Philippe Palanque is Professor in Computer Science at University of Toulouse 3. The main driver of Philippe's research has been to address in an even way Usability, Safety and Dependability in order to build trustable safety critical interactive systems. As for conferences he is a PC member of conferences in these domains such as DSN 2014 (44th conference on Dependable Systems and Networks), and was co-chair of CHI 2014. Célia Martinie (PhD) is Assistant Professor in Computer Science at University of Toulouse 3. She has been working on task modeling techniques for the design and development of interactive systems since 2009 and led the development of HAMSTERS tool.


Hugues Cassé
email: casse@irit.fr

Haluk Ozaktas
email: ozaktas@irit.fr

Christine Rochange
email: rochange@irit.fr

A Framework to Quantify the Overestimations of Static WCET Analysis *

Keywords: 1998 ACM Subject Classification C.3 Computer Systems Organization -Special-Purpose and Application-Based Systems -Real-Time and Embedded Systems static WCET analysis, uncertainty, overestimation, cache analysis Digital Object Identifier 10.4230/OASIcs.WCET.2015

To reduce complexity while computing an upper bound on the worst-case execution time, static WCET analysis performs over-approximations. This feeds the general feeling that static WCET estimations can be far above the real WCET. This feeling is strengthened when these estimations are compared to measured execution times: generally, it is very unlikely to capture the worstcase from observations, then the difference between the highest watermark and the proven WCET upper bound might be considerable. In this paper, we introduce a framework to quantify the possible overestimation on WCET upper bounds obtained by static analysis. The objective is to derive a lower bound on the WCET to complement the upper bound.

Introduction

Tasks in hard real-time systems are subject to strict deadlines and any failure to meet a deadline might have severe consequences. Many approaches to hard real-time task scheduling have been proposed in the literature [START_REF] Davis | A survey of hard real-time scheduling for multiprocessor systems[END_REF]. Based on the knowledge of the worst-case execution time (WCET) of individual tasks, they compute a task schedule that is valid in the sense that all deadlines are met.

The correctness of a task schedule relies on the confidence in the estimated WCETs of tasks. An under-estimated WCET might lead to an unsafe schedule and possibly to a violation of deadlines at runtime. Therefore, considering reliable techniques to derive those WCET estimates is crucial. The flexibility of measurement-based techniques is appealing but they cannot bring strict guarantees that the worst case has been observed or can be extrapolated. Static analysis techniques instead provide reliable estimates as soon as any detail on the hardware architecture of the target processor is known and correctly modelled 1 .

Static analysis approaches determine an upper bound on the WCET based on an overapproximation of the processor state at any point in the program. For example, usual strategies to analyse the dynamic behaviour of a cache memory are based on abstract interpretation techniques that build abstract states of the cache at any program point.

By nature, WCET upper bounds are pessimistic and static WCET analysis techniques are often blamed for generating excessive overestimation. This feeling is exacerbated when the largest observed (measured) execution time is compared to the statically-determined WCET. The difference between the two is often exhibited as an indicator of overestimation although without any evidence (even any conviction) that the largest observed execution time was measured on the longest possible path. In the same time, it can generally not be shown that the estimated WCET value could be really reached. The goal of this paper is to show how additional numbers could be delivered together with the trustworthy WCET estimation in order to give better insight into how imprecise it might be.

Three factors can contribute to overestimating the WCET: (a) the analysis identifies a longest path that is not feasible in practice due to restrictions on input data values which are unknown (unspecified) to the analysis and/or to the incapacity of the analysis to take such restrictions correctly into account; (b) the analysis assumes a worst-case initial hardware state which cannot be observed in practice; (c) the way the analysis is performed generates overestimation (generally to reduce complexity). In this paper, we focus on factor (c) and leave the study of the two others for future work.

Our objective is to isolate the part of the estimated WCET that does not result of any over-approximation. We refer to this value as a lower bound on the real WCET : it represents an execution time that is feasible (assuming that both the sets of possible paths and initial hardware states are not over-approximated). The difference between the lower and upper bounds on the WCET gives an insight into the accuracy of the estimated WCET.

The paper is organised as follows: Section 2 introduces a framework to determine lower bounds along with upper bounds on the WCET of a task. A possible use of this framework for the cache analysis is presented in Section 3. Experimental results are reported in Section 4 and Section 5 concludes the paper.

2

A Framework to Quantify Pessimism of Static WCET Estimations

Static WCET Analysis Techniques

The usual flow for static WCET analysis is a sequence of several steps. The input is the program to analyse, more precisely its binary code required for low-level analyses. Its source code can be useful for path analyses.

In the first step (Flow/Path Analysis), possible execution paths are identified and encoded as a CFG (Control Flow Graph) and a set of flow facts (loop bounds, infeasible paths, etc.).

The next step (Global Low-level Analysis) determines the behaviour of history-based hardware schemes, such as (instruction and data) caches or branch predictors. These schemes make use of limited-capacity storage which may engender conflicts between distant parts of the program. In addition, they have been designed to improve the average execution time of the program by exploiting the execution history and hence may exhibit highly dynamic behaviours. These effects are particularly interesting to handle with static analysis approaches because time-expensive behaviours are relatively rare (but must be considered for safety reasons) and are particularly difficult to capture through measurements.

The third step (Local Low-level Analysis) takes the global and local behaviours of the hardware components into account to compute the possible execution times of each basic block in the CFG.

Finally, an ILP (Integer Linear Programming) system is built [START_REF] Li | Efficient microarchitecture modelling and path analysis for real-time software[END_REF] for WCET Computation (IPET method). Its objective is to maximize the task execution time expressed as the sum of the basic blocks execution times weighted by their respective execution counts. These execution counts are bounded by a collection of linear constraints that express: (a) restrictions on the possible execution paths in the code, and (b) restrictions on the possible occurrences of some hardware-level behaviours.

Sources of Overestimation

The overestimation of static WCET analyses is due to three kinds of reasons: imprecise flow facts: in the absence of full information on possible execution paths, the WCET analysis might consider infeasible paths and those might happen to exhibit longer execution times than valid paths. The goal of flow analysis [START_REF] Suhendra | Efficient detection and exploitation of infeasible paths for software timing analysis[END_REF] is to eliminate such infeasible paths but it might fail to identify all of them for several reasons: (a) restrictions on input data values are under-specified by the user; (b) some information given by the user cannot be translated into flow facts (this depends on the annotation format considered by the WCET analysis tool, complex scenarii might be difficult to describe); (c) some flow facts derived from user annotations or automatically extracted from the source or binary code cannot be exploited during the analysis (e.g. they cannot be turned into linear constraints when the WCET is computed using the IPET method). In this paper, we assume that the applications under analysis do not suffer such issues and we consider that every relevant information of possible execution flows is known and taken into account when computing WCETs. We leave the analysis of imprecision due to incomplete information on input data for future work. imprecise information on the initial hardware state: the state of the processor (pipeline) and of the memory hierarchy (cache memories) has an impact on the (worst-case) execution time of a task. For example, if part of the task code already resides in the instruction cache, the overall latency of instructions fetches is shorter than if the cache is empty.

In the absence of such information, the WCET analysis systematically assumes the worst-case situation and this might lead to overestimating execution times [START_REF] Nemer | Inter-Task WCET computation for A-way Instruction Caches[END_REF]. We ignore this issue in this paper and we assume that the tasks under analysis can effectively start with the worst-case hardware state. static analysis techniques do not unroll any possible execution path. Instead, they derive at each point in the code some properties that hold for any possible execution. This keeps the complexity of determining the longest path tractable. A consequence is that these properties might be pessimistic for a particular execution. Here, the pessimism is generated by the analysis technique and not by incomplete information on possible execution scenarii. The contribution of this paper is to quantify this overestimation which is inherent to static WCET analysis.

Proposed Framework

According to the discussion above, the execution time of an application code running on a given platform that ensures full isolation (i.e. the execution time of a task cannot be impacted by any other task or system operation) is a function of the initial hardware state and of the followed execution path:

ET : H × P → N
where H is the set of possible initial platforms states and P is the set of possible execution paths in the code. The WCET of the program is defined as:

WCET = max (h,p)∈H×P ET (h, p) .
Static WCET analysis aims at determining an upper bound on the execution time. Instead of running the code with different (h, p) ∈ {H × P} pairs, as measurement-based approaches would do, it performs abstractions to derive, at each program point, some properties that hold for any execution path. These abstractions may lead to some undecided local behaviours which are then over-approximated to compute the final WCET estimation. For example, it can happen that the latency of an access to a cache is unknown because the cache analysis can not determine whether it will be a hit or a miss (because this depends on the execution history). One of these options must be chosen to compute the WCET. We use the term scenario to refer to a combination of selected options for all the unpredictable behaviors in the program and C denotes the set of all the possible scenarii.

The estimated (abstract) WCET can then be defined as follows:

WCET # : H × P × C → N .
Each set D, where D stands for H, P or C, can be considered in several flavours: D ⊤ is the set of theoretically possible values. For example, P ⊤ is the set of possible paths expressed by the program control flow graph.

D + contains all the values that may (but are not guaranteed to) be feasible.

In other words, D + is a subset of D ⊤ that excludes the values that can be proven infeasible. For example, P + includes all the execution paths that fulfil the flow constraints (both those computed by the flow analysis step and those specified by the user). Similarly, C + is the set of all the scenarii that may occur and is considered by static WCET analysis to compute a trustworthy upper bound on the WCET.

D r is the set of all the values that are really feasible.

Usually, this set is unknown or too large to be exhaustively explored.

D -is a set of values that can be guaranteed to be feasible.

To illustrate this, let us consider an access to the cache that cannot be classified as a hit or a miss by the static cache analysis. If the processor is free of timing anomalies, the optimistic option is to consider the access as a hit, while the pessimistic assumption is a miss. Any scenario in C -specifies a hit, and any scenario in C + specifies a miss. We have:

D -⊆ D r ⊆ D + ⊆ D ⊤ .
Usual static WCET analysis techniques produce an upper bound on the real WCET: WCET + : H + × P + × C + → N. Provided H + and P + are trustworthy, and assuming that the low-level analyses are correct (then C + is trustworthy too), we must have WCET r ≤ WCET + where WCET r is the real (but unknown) WCET of the application.

On the other hand, given that all the values in H -, P -and C -are feasible, it must be that WCET -≤ WCET r with:

WCET -: H -× P -× C -→ N.
To summarize, WCET -≤ WCET r ≤ WCET + . The possible overestimation on the real WCET is then quantified by U = (WCET + -WCET -)/WCET -. We refer to WCET -as the lower bound on the WCET2 in contrast to WCET + , the upper bound on the WCET. Only WCET + is a reliable estimation of the program's WCET.

Note that computing the WCET from sets of different flavours, e.g. H -× P -× C + , does not provide any useful estimate: it cannot be ordered with respect to the real WCET.

In this paper, we focus on the computation of C -and C + for various low-level analyses.

Upper-Bounding the Possible Overestimation of Low-Level Analyses

Background on Cache Analysis

Category-Based Approaches

Global low-level analyses account for the impact of history-based hardware schemes. Such schemes include instruction or data caches, or dynamic branch predictors. They offer either a fast access (a hit in the usual cache terminology), or a slow access (a miss), and are designed to maximize the number of hits in the average case.

A common approach to support such schemes is to assign a category to each access to the device. It has been successfully applied to LRU3 instruction [START_REF] Ferdinand | Applying compiler techniques to cache behavior prediction[END_REF] and data [START_REF] Ferdinand | On predicting data cache behavior for realtime systems[END_REF] caches but also to branch prediction history tables [START_REF] Colin | Worst case execution time analysis for a processor with branch prediction[END_REF] and flash memory prefetchers [START_REF] Holsti | WCET Tool Challenge[END_REF].

Instruction Cache Analysis

In [START_REF] Ferdinand | Applying compiler techniques to cache behavior prediction[END_REF] an analysis for LRU instruction caches introduces three categories: Always-Hit (AH ) -on each access, the instruction block is present in the cache, resulting in a fast access; Always-Miss (AM ) -the instruction block is never in the cache, always causing a slow access; and Not Classified (NC ) -the behaviour cannot be predicted at analysis time. These categories have later been extended in [START_REF] Ferdinand | A fast and efficient cache persistence analysis[END_REF][START_REF] Ballabriga | Improving the first-miss computation in setassociative instruction caches[END_REF] with the Persistence category (PERS) that is assigned when the first occurrence of an access belonging to a loop body cannot be classified but the following accesses (in the next iterations) are hits.

The computation of categories is based on Abstract Interpretation (AI) techniques. A cache state is derived at each point of the program, as a function that assigns an age to each cache block: S : B → A. The age represents the position of the block in a cache set according to the LRU replacement policy: it is an integer value between 0 and A, where A is the associativity of the cache -an age equal to A means that the block is not in the cache. The Update function expresses the behaviour of the cache: U : S × B → S (the accessed block gets age 0, the age of blocks younger than it is incremented by 1, the age of other blocks remains unchanged). To avoid an explosion of the number of possible cache states at a program point, abstract interpretation merges states using a Join function J # : S # × S # → S # that operates on an abstract representation S # of the concrete cache states.

A category is assigned to an access to cache block b ∈ B at any program point p as follows: (i) if any state s ∈ S contains b whatever the path leading to p, the access always hits (AH ); (ii) if no path leading to p produces an s ∈ S that contains b, then the access always misses (AM ). Two abstract interpretation analyses are needed to draw such conclusions: (i) the MUST analysis computes the worst-case age of block b (if the worst-case age of b is less than A at point p, this means that b is always in the cache at point p) with

J # (s 1 , s 2 ) = s such that ∀b ∈ B, s[b] = max(s 1 [b], s 2 [b]); and (ii) the MAY analysis computes the best-case age of b (if the best-case age of b is A at point p, this means that b is never in the cache at point p) with J # (s 1 , s 2 ) = s such that ∀b ∈ B, s[b] = min(s 1 [b], s 2 [b]).
If none of the two previous predicate holds, the access is considered as NC 4 . This category is the main source of overestimation in the WCET estimation.

Data Cache Analysis

The analysis of data caches is very similar to that for instruction caches. The additional part consists in performing a data flow analysis to discover the addresses accessed by load and store instructions and to determine which cache blocks are used [START_REF] Ferdinand | On predicting data cache behavior for realtime systems[END_REF].

While an instruction fetch addresses a single memory block (known at analysis time), an access to data performed within a loop may reference several memory blocks across iterations. The data flow analysis may find that the target of a memory access has a single value, or a set or an interval of values. A special value, ⊤, represents the cases where the analysis fails to determine the possible address values. Such a failure could reflect input-data dependent addresses or might be due to an imprecise address analysis that makes over-approximations.

The imprecision in the address analysis is also a source of imprecision in the computation of the abstract cache states used to assign the categories. When a single address value has been found for a load/store instruction, the cache state is updated in the same way as for an instruction cache. If the instruction may reference several different addresses, the analysis accounts for a possible impact on several cache sets (on all cache sets if the predicted address is ⊤. For the MUST analysis, all blocks are aged by one (whatever the number of possible address values, only one is accessed at a time); in the MAY analysis, the age of all accessed blocks is unchanged while other blocks are aged. In the end, an imprecision in the address analysis directly induces over/underestimation (MUST/ MAY) of ages in the cache states and hence on the categories: memory access with imprecise address sets are more likely to be categorized as NC .

Exploitation of Cache Categories to Compute the WCET

A straightforward way to account for the cost of instruction and data cache misses in the estimated WCET of a program is to add it to the expression of the execution time so that the ILP solver maximises the number of misses together with the execution time. However, this solution is reliable only if the execution time of a sequence of instructions experiencing a cache miss is less than the adding the the cost of that miss to the execution time with a hit. This is not the case when the program runs on an architecture that enables timing anomalies [START_REF] Lundqvist | Timing anomalies in dynamically scheduled microprocessors[END_REF], i.e. for which the local worst case for the execution time of the sequence might be a hit.

A more reliable approach is to compute the local WCET of each basic block in the program CFG as many times as possible combinations of hits and misses for NC accesses exist. Then the execution count of each of these different WCET values is bounded by a combination of category-related execution counts in the ILP formulation. This approach has been implemented in OTAWA [START_REF] Ballabriga | OTAWA: An open toolbox for adaptive wcet analysis[END_REF], the static WCET analysis tool used to carry out the experiments reported in this paper. In OTAWA, a local WCET is computed for each sequence of two basic blocks (then related to an edge in the CFG) based on contextual execution graphs [START_REF] Rochange | A context-parameterized model for static analysis of execution times[END_REF]. It considers an optimistic or pessimistic context (availability of hardware resources) for the first and second block in the sequence, respectively, thus ensuring that an upper bound is found.

Contribution to the Lower Bound on the WCET (WCET -)

Approximation in category-based instruction cache analysis clearly comes from the join operator J # : even if the input cache states are precise, the result may be imprecise. In data cache analysis, another source of approximation is the address analysis that sometimes fails to derive precise address values, which leads to additional inaccuracy in the computation of A simplistic approach to estimate WCET -would consist in considering each NC access to the cache as a hit (similarly, each PERS access as a hit in the first loop iteration). However, to account for possible timing anomalies, both options must still be considered (as explained in Section 3.1). More precisely, an opposite assumption is taken for each basic block in a sequence. For example, all NC accesses in the first block are accounted for as misses, while those in the second block are considered as hits. This way, the local WCET of the second block is minimized (assuming no timing anomaly can occur). At the end, when all possible scenarii have been explored, the lowest local WCET is kept for each basic block.

4

Experimental Study and Discussions

Methodology

All the experiments reported here have been done using our OTAWA toolset [START_REF] Ballabriga | OTAWA: An open toolbox for adaptive wcet analysis[END_REF] in which we implemented support for the analysis of the possible overestimation.

We considered a microprocessor architecture with a 7-stage in-order pipeline similar to that of the MPC5554 from Freescale. We assumed a single level of separated instruction and data caches with the following parameters: 4KB of size, 4-way associative with LRU replacement, with 16B cache lines and a write-allocate policy. A memory latency (for cache misses) of 20 cycles was considered. We considered an infinite-capacity write buffer, then only loads are considered in the experiments.

We have analysed several benchmarks from the MediaBench suite [START_REF] Lee | Mediabench: a tool for evaluating and synthesizing multimedia and communicatons systems[END_REF]: cjpeg_jpeg6b_wrbmp is a JPEG image bitmap writing code; gsm, gsm_decode and gsm_encode are respectively a GSM 06.10 provisional standard codec, decoder and encoder; h264dec_ldecode_block and h264dec_ldecode_macroblock are H.264 block and macroblock decoding functions. The dynamic numbers of accesses to the instruction cache and loads of these benchmarks (i.e. counts on the worst-case path) are given in Table 1.

We assume that each execution path that fulfils the flow facts specified to the WCET analysis tool and converted into constraints in the ILP formulation is feasible (i.e. no overestimation comes from the flow analysis). Similarly, we assume that any initial hardware state is feasible.

Results

Table 2 shows the possible overestimation on (WCET + ). For cjpeg_jpeg6b_wrbmp), the maximum overestimation on WCET + is small (36.25%), meaning that the analysis is precise in this case. The other benchmarks exhibit much larger potential overestimation, ranging from 96.33% to 216%. As we will see below, this is mainly due to the weakness of the data cache analysis we considered, which fails to compute the targets of load instructions.

In the following, we try to isolate the respective contributions of the instruction and data cache analyses to the possible overestimation of the WCET upper bound. For that purpose, we compute U i$ (resp. U d$ ) considering pessimistic results for the other analysis.

Table 3 shows the contribution to the overestimation by the instruction cache analysis (U i$ ). These low numbers (less than 5%) show how instruction caches can be taken into account with much accuracy in WCET analysis. There are two reasons for this. First, the instruction cache analysis generates relatively precise results: columns 2 and 3 show how many accesses are classified as PERS (Persistent, generating a hit for all loop iterations but the first one) or NC (Not Classified, then considered as a hit to be optimistic and as a miss to be pessimistic). Except for one benchmark, the number of instruction fetches that cannot be classified is less than 20%. Second, the impact of a miss on an instruction fetch can sometimes be partly hidden by stalls generated by instructions inter-dependencies, and this is captured by our pipeline analysis.

The major part of the possible overestimation on the static WCET comes from the data cache analysis, as shown in Table 4. Contrarily to the instruction cache, many accesses (from 46% to more than 80%) are categorized as NC . This is due to the basic address analysis implemented in OTAWA, that only handles global and stack data. Addresses of accesses to array elements are not determined by this analysis. As reported in column 4, the precise target address could not be determined by the data analysis for a considerable amount of load instructions (from 22% up to 64%). An unknown target address impacts the load instruction itself, since it cannot be determined whether it hits or misses in the cache, but it also impacts subsequent loads by producing a destructing effect on the abstract cache state. Clearly, the address analysis can be identified as the weakest link here.

Possible Interpretations of the Lower Bound on the WCET

Several uses of the overestimation metric (U ) -difference between the upper and lower bounds on the WCET -can be envisioned: An obvious application is the qualification of a particular estimated WCET. If U is small, then the estimated WCET can be considered as precise. This means that WCET + is not far above the real WCET (even if measured execution times look much lower). Then the estimated WCET will not drive to oversize the hardware required to meet deadlines, nor break irrelevantly the real-time constraints verification. Instead, if U is large, this may mean either that the WCET computation approach does not fit the particular profile of application under analysis, or that the way the application is coded drives too many non predictable behaviours.

Local contributions to overestimation could also be used to implement an adaptive WCET analysis. For example, a first fast analysis could be performed with a very simple (but imprecise) method. Then the program parts exhibiting the largest overestimation could be processed again with a more precise but also more time-consuming method. Abstract interpretation is well adapted to support such an approach: to enhance accuracy, the analysis could choose which states would be merged and which would be kept separate by the Join operator. Finally, the overestimation metric could also be considered as a basis to compare different analyses with respect to the precision they can achieve.

Conclusion

Static WCET analysis has the reputation of overestimating WCET, which is not completely false since it is designed to compute upper bounds on the real WCET. However the size of the gap between the real and estimated WCETs is unknown (since the real WCET cannot usually be determined). Comparison of static WCET estimates to measured execution times often shows large discrepancies but nobody can tell whether the upper bound is too pessimistic or if the measured execution time was observed under particularly favourable conditions.

Our objective in this paper was to propose a framework to evaluate the possible overestimation of static WCETs. We defined a lower bound on the real WCET (WCET -). This value can been seen as a complement to the traditional upper bound, WCET + . It represents an execution time that can be reached, provided all information on infeasible paths and initial hardware states were known and taken into account at analysis time. The real WCET is guaranteed to be between the lower and upper bounds. We envision that these information could be useful not only to increase confidence in the static WCET analysis but also to guide code optimisations to enhance precision. It could also be used to compare analysis techniques. Note that any observed execution time that would be greater than WCET -could be used to tighten the interval in which the real WCET falls, and thus to improve precision.

As future work, we plan to apply the proposed approach to other hardware components in single-and multi-core architectures. We will also investigate how it can be used to evaluate the precision of flow analysis techniques.

Table 1

 1 Benchmarks. Imprecise cache states lead to imprecise categories (NC , and PERS to a lesser extent) then to possibly overestimating WCET + .

	benchmark # accesses to I$ # loads
	cjpeg_jpeg6b_wrbmp	37,328	6,187
	gsm	4,2319,41	895,605
	gsm_decode	1,739,854	211,854
	gsm_encode	2,378,890	681,917
	h264dec_ldecode_block	11,942	1,802
	h264dec_ldecode_macroblock	51,373	51,373
	cache states.

Table 2

 2 Possible overestimation on WCET + .

	benchmark U = (WCET +	-WCET -)/WCET -
	cjpeg_jpeg6b_wrbmp		36.25%
	gsm	161.95%
	gsm_decode		96.33%
	gsm_encode	216.00%
	h264dec_ldecode_block	126.11%
	h264dec_ldecode_macroblock	144.41%

Table 3

 3 Possible overestimation in instruction cache analysis.

	benchmark PERS	NC	U i$
	cjpeg_jpeg6b_wrbmp 25.0%	5.1%	0.54%
	gsm 34.9% 19.6% 2.36%
	gsm_decode 38.4% 29.3% 1.83%
	gsm_encode 31.5% 18.7% 2.29%
	h264dec_ldecode_block 52.6% 18.8% 4.84%
	h264dec_ldecode_macroblock 57.7%	1.8%	0.32%

Table 4

 4 Possible overestimation in data cache analysis.

	benchmark PERS	NC	unknown @	U d$
	cjpeg_jpeg6b_wrbmp	9.2%	46.1%	33.8%	35.29%
	gsm	0.0%	55.8%	38.4%	66.25%
	gsm_decode	0.0%	69.4%	52.6%	143.66%
	gsm_encode	0.0%	51.7%	34.1%	190.01%
	h264dec_ldecode_block	0.0%	80.5%	64.4%	82.51%
	h264dec_ldecode_macroblock	0.0%	54.4%	22.1%	141.90%

Note that this lower bound on the WCET (WCET -) is something different from the best-case execution time (BCET).

LRU stands for Least Recently Used and refers to the cache block replacement policy. We consider this policy in the remainder of the paper.

As mentioned above, another category for blocks that remains in the cache across the iterations of a loop but cannot be classified for their first access can be assigned. This analysis is not described here for the sake of simplicity. The reader may refer to[START_REF] Ferdinand | A fast and efficient cache persistence analysis[END_REF][START_REF] Ballabriga | Improving the first-miss computation in setassociative instruction caches[END_REF] for further details.

* This work was partially supported by the French research foundation (ANR) through the W-SEPT project (ANR-12-INSE-0001). See http://wsept.inria.fr. strong assumption that is difficult to achieve but this issue is beyond the scope of this paper.


Zhouye Chen

Ningning Zhao

Adrian Basarab

Denis Kouamé

ULTRASOUND COMPRESSIVE DECONVOLUTION WITH ℓ P -NORM PRIOR

Keywords: ultrasound imaging, compressive sampling, deconvolution, generalized Gaussian distribution, alternating direction method of multipliers

It has been recently shown that compressive sampling is an interesting perspective for fast ultrasound imaging. This paper addresses the problem of compressive deconvolution for ultrasound imaging systems using an assumption of generalized Gaussian distributed tissue reflectivity function. The benefit of compressive deconvolution is the joint volume reduction of the acquired data and the image resolution improvement. The main contribution of this work is to apply the framework of compressive deconvolution on ultrasound imaging and to propose a novel ℓ p -norm (1 ≤ p ≤ 2) algorithm based on Alternating Direction Method of Multipliers. The performance of the proposed algorithm is tested on simulated data and compared with those obtained by a more intuitive sequential compressive deconvolution method.

INTRODUCTION

Ultrasound (US) medical imaging has the advantages of being noninvasive, harmless, cost-effective and portable over many other imaging modalities such as X-ray Computed Tomography or Magnetic Resonance Imaging [START_REF] Szabo | Diagnostic ultrasound imaging: inside out[END_REF]. US images have a textural appearance caused by the presence of a speckle noise resulting from the interferences between the US waves backscattered by the randomly-located scatterers of the tissue being imaged. However, the limited bandwidth of the imaging transducer, the characteristics of the US propagation such as the diffraction and the imaging system tend to degrade the resolution of the US images. Thus, one of the research tracks extensively explored in the litterature is the deconvolution of US images [START_REF] Taxt | Two-dimensional noise-robust blind deconvolution of ultrasound images[END_REF][START_REF] Michailovich | Blind deconvolution of medical ultrasound images: A parametric inverse filtering approach[END_REF][START_REF] Morin | Semi-blind deconvolution for resolution enhancement in ultrasound imaging[END_REF][START_REF] Zhao | Restoration of ultrasound images using a hierarchical bayesian model with a generalized gaussian prior[END_REF][START_REF] Alessandrini | A restoration framework for ultrasonic tissue characterization[END_REF]. Based on the first order Born approximation, these methods assume that the US radiofrequency (RF) images follow a 2D convolution model between the point spread function (PSF) and the tissue reflectivity function to be recovered [START_REF] Jensen | Deconvolution of in vivo ultrasound b-mode images[END_REF].

More recently, a few research teams (e.g. [START_REF] Achim | Compressive sensing for ultrasound rf echoes using a-stable distributions[END_REF][START_REF] Quinsac | Frequency domain compressive sampling for ultrasound imaging[END_REF][START_REF] Chernyakova | Fourier-domain beamforming: the path to compressed ultrasound imaging[END_REF][START_REF] Liebgott | Compressive sensing in medical ultrasound[END_REF][START_REF] Liebgott | Pre-beamformed rf signal reconstruction in medical ultrasound using compressive sensing[END_REF][START_REF] Schiffner | Pulse-echo ultrasound imaging combining compressed sensing and the fast multipole method[END_REF]) evaluated the application of compressive sampling (CS) in the context of 2D and 3D US imaging. CS [START_REF] Donoho | Compressed sensing[END_REF][START_REF] Candès | Robust uncertainty principles: Exact signal reconstruction from highly incomplete fre-quency information[END_REF] routines, an image from few linear measurements (below the limit imposed by the Shannon-Nyquist theorem). For this, two conditions must be fulfilled: the image must have a sparse representation in a known basis or frame and the measurement and sparifying basis must be incoherent [START_REF] Candès | An introduction to compressive sampling[END_REF]. The benefits of CS in US imaging reported by the existing approaches are the reduction of the acquired data (very useful for instance in 3D or in Doppler imaging) or/and of the acquisition time.

The direct models of both deconvolution and CS are linear models leading to ill-posed inverse problems. Thus, combining these two frameworks, leading to the so-called compressive deconvolution problem [START_REF] Amizic | Compressive blind image deconvolution[END_REF][START_REF] Spinoulas | Simultaneous bayesian compressive sensing and blind deconvolution[END_REF], is a very interesting and recent research track. The benefit of compressive deconvolution is the joint volume reduction of the acquired data and the image quality improvement.

In this paper, we propose a novel scheme of compressive deconvolution and its application to US imaging. The resulting acquisition model is as follows:

y = ΦHx + n (1)
where Φ ∈ R M ×N corresponds to the CS acquisition matrix composed for example by M random Gaussian vectors with (M << N ), H ∈ R N ×N is a block circulant with circulant block (BCCB) matrix related to the 2D PSF of the system, x ∈ R N and n ∈ R M represent the lexicographically ordered tissue reflectivity function and the zero-mean additive white Gaussian noise, respectively. y ∈ R M corresponds to the M compressed measurements of one US RF image.

The aim of this work is to estimate the tissue reflectivity function x from the compressed and blurred measurements y. The solution proposed herein is based on the alternating direction method of multipliers (ADMM) [START_REF] Boyd | Distributed optimization and statistical learning via the alternating direction method of multipliers[END_REF][START_REF] Zhao | Alternating direction method of multipliers for nonlinear image restoration problems[END_REF] and uses two constraints. The first one is related to CS and imposes via an l 1 -norm the sparsity of the RF images in the 2D Fourier domain [START_REF] Quinsac | Frequency domain compressive sampling for ultrasound imaging[END_REF]. The second one is related to an l p -norm [START_REF] Antoniadis | Wavelet thresholding for some classes of non-Gaussian noise[END_REF]. Based on the assumption of Generalized Gaussian Distributed (GGD) x, we employ the minimization of an l p -norm of the tissue reflectivity function [START_REF] Zhao | Restoration of ultrasound images using a hierarchical bayesian model with a generalized gaussian prior[END_REF][START_REF] Alessandrini | A restoration framework for ultrasonic tissue characterization[END_REF]. The simulation results show that our method outperforms a more intuitive sequential scheme estimating firstly the blurred image and processing the deconvolution in a second step.

The remainder of the paper is organized as follows. In section 2 we formulate the compressive deconvolution problem. Section 3 details our proposed compressive deconvo-lution algorithm based on ADMM. Simulation results are shown in section 4 before drawing the conclusions in section 5.

PROBLEM FORMULATION

An intuitive idea to invert the direct model in ( 1) is to proceed in two sequential steps. The aim of the first step is to recover the blurred US RF image r = Hx from the compressed measurements y by solving the following optimization problem:

min a∈R N a 1 + 1 2µ y -ΦΨa 2 2 ( 2
)
where µ is a parameter, a is the sparse representation of the blurred image r in the transformed domain Ψ. That is, a = Ψ -1 Hx, where Ψ ∈ R N ×N represents an orthonormal basis. Different basis have been shown to provide good results in the application of CS in US imaging, such as the wavelets, the wave atoms or the 2D Fourier basis [START_REF] Liebgott | Compressive sensing in medical ultrasound[END_REF] .

Once the blurred RF image recovered by minimizing the convex problem in (2), one can use an US dedicated deconvolution method to restore the tissue reflectivity function x.

While the sequential approach represents the most intuitive way to solve the compressive deconvolution problem, dividing a single problem into two separate subproblems will inevitably generate larger estimation errors. Therefore, we propose herein a method to solve the CS and deconvolution problem simultaneously.

Following recent results in the US deconvolution literature [START_REF] Zhao | Restoration of ultrasound images using a hierarchical bayesian model with a generalized gaussian prior[END_REF][START_REF] Alessandrini | A restoration framework for ultrasonic tissue characterization[END_REF], we assume that the tissue reflectivity function x follows a GGD. Here, we are interested in shape parameters ranging from 1 to 2. Taking into account this statistical prior information for x, we reformulate our compressed deconvolution problem in:

min x∈R N Ψ -1 Hx 1 +α x p p + 1 2µ y -ΦHx 2 2 (3)
where α, µ are parameters, p is related to the shape parameter and 1 ≤ p ≤ 2. In the next section, we propose an ADMM-based scheme to solve the inverse problem above.

PROPOSED ULTRASOUND COMPRESSIVE DECONVOLUTION ALGORITHM

Basics of Alternating Direction Method of Multipliers

Alternating Direction Method of Multipliers (ADMM) has been extensively studied in the areas of convex programming and variational inequalities [START_REF] Boyd | Distributed optimization and statistical learning via the alternating direction method of multipliers[END_REF]. The general optimization problem is as follows:

min u,v f (u) + g(v) s.t. Bu + Cv = b, u ∈ U, v ∈ V (4)
where U ⊆ R s and V ⊆ R t are given convex sets, f : U → R and g : V → R are closed convex functions; B ∈ R r×s and C ∈ R r×t are given matrices; b ∈ R r is a given vector.

By attaching the Lagrangian multiplier λ ∈ R r to the linear constraint, the Augmented Lagrangian (AL) function of (4) is

L(u, v, λ) = f (u) + g(v) -λ T (Bu + Cv -b) + β 2 Bu + Cv -b 2 2 (5)
where β > 0 is the penalty parameter for the linear constraints to be satisfied. The standard ADMM framework follows the iterative process:

       u k+1 ∈ argmin u∈U L(u, v k , λ k ) v k+1 ∈ argmin v∈V L(u k+1 , v, λ k ) λ k+1 = λ k -β(Bu k+1 + Cv k+1 -b) (6)
The advantage of ADMM is that it can split awkward intersections and objectives to easy subproblems and the iterations are comparable to those of other first-order methods. Moreover, it is relatively easy to implement.

Proposed ADMM parameterization for Ultrasound Compressive Deconvolution

In this subsection, we propose an ADMM method for solving the ultrasound compressive deconvolution problem in [START_REF] Michailovich | Blind deconvolution of medical ultrasound images: A parametric inverse filtering approach[END_REF].

Using a trivial variable change, the minimization problem in (3) can be rewritten as:

min x∈R N w 1 +α x p p + 1 2µ y -Aa 2 2 (7)
where a = Ψ -1 Hx. Let us denote z = w x , w = a and A = ΦΨ. The reformulated problem in [START_REF] Jensen | Deconvolution of in vivo ultrasound b-mode images[END_REF] can fit the general ADMM framework in (4) by choosing:

f (a) = 1 2µ y -Aa 2 2 , g(z) = w 1 +α x p p , B = I N ×N Ψ , C = -I N ×N 0 0 -H and b = 0.
The augmented Lagrangian function of ( 7) is given by

L(a, z, λ) = f (a) + g(z) -λ T (Ba + Cz) + β 2 Ba + Cz 2 2 ( 8
)
where λ ∈ R 2N stands for λ =

λ 1 λ 2 , λ i ∈ R N (i = 1, 2).
According to the standard ADMM iterative scheme, the minimization with respect to a and z will be performed alternatively, followed by the update of λ. More precisely, the numerical scheme for solving (8) is described in Algorithm 1.

The algorithm stops when the convergence criterion x kx k-1 / x k-1 < 5e -3 is satisfied.

Algorithm 1 The Alternating Direction Method of Multipliers for Solving (3) Input: a 0 , λ 0 , α, µ, β Output: x 1: while not converged do 2:

Step 1. Given a k-1 , update w k by using the shrinkage formula: 3:

w k = Shrink 1 β (a k-1 - λ k-1 1 β ) 4:
Step 2. Given a k-1 , update x k by using the proximal operator: 5:

x k = prox αγ • p p /β {x k-1 -γh ′ (x k-1 )} 6:
Step 3. Given z k = w k

x k and λ k-1 , update a k by solving the following linear system 7:

a k = ( 1 µ A T A + βB T B) -1 ( 1 µ A T y + B T λ k-1 - βB T Cz k ) 8:
Step 4. Update λ k by 9:

λ k = λ k-1 -β(Ba k + Cz k ) 10: end while
We solve the subproblems in Step 1 and Step 2 by using proximal operator as proposed in [START_REF] Pesquet | A parallel inertial proximal optimization method[END_REF][START_REF] Pustelnik | Parallel proximal algorithm for image restoration using hybrid regularization[END_REF][START_REF] Pustelnik | Relaxing tight frame condition in parallel proximal methods for signal restoration[END_REF]. The wsubproblem in Step 1 is equivalant to:

min w∈R N w 1 + β 2 a k-1 -w - λ k-1 1 β 2 2 (9)
The solution of ( 9) is given explicitly in the literature (see e.g. [START_REF] Combettes | Proximal splitting methods in signal processing[END_REF])

w k (i) = max |a k-1 (i) - λ k-1 1 (i) β | - 1 β , 0 • sign(a k-1 (i) - λ k-1 1 (i) β ) (10)
where i ∈ {1, 2, ..., N }, w i denotes the intensity value of pixel i. The x-subproblem in Step 2 is as follow:

min x∈R N α x p p -λ k-1 2 T (Ψa k-1 -Hx)+ β 2 Ψa k-1 -Hx 2 2 (11) Denoting by h(x) = 1 2 Ψa k-1 -Hx - λ k-1 2 β 2 2 , we approximate h(x) by h ′ (x k-1 )(x -x k-1 ) + 1 2γ x -x k-1 2 2 (12)
where γ > 0 is a proximal parameter and h ′ (x k-1 ) is the gradient of h(x) when x = x k-1 , which is equal to

h ′ (x k-1 ) = H T (Hx k-1 -Ψa k-1 - λ k-1 2 β ) (13)
We should note that the proximal operator associated to • p p has been given in [START_REF] Combettes | Proximal splitting methods in signal processing[END_REF]. Herein, we use Newton's method to obtain its numerical solution.

In Step 3, the update of a is written in the form of solving an N × N linear system or inverting an N × N matrix. However, since the sparse basis Ψ is orthogonal, it can be reduced to solving a smaller M × M linear system or inverting an M × M matrix by Sherman-Morrison-Woodbury formula [START_REF] Deng | Group sparse optimization by alternating direction method[END_REF]:

(β 1 I N +β 2 A T A) -1 = 1 β 1 I N - β 2 β 1 A T (β 1 I M +β 2 AA T ) -1 A (14)
where I N ∈ R N ×N and I M ∈ R M ×M are identity matrices. In this paper, we considered the compressive sampling matrix Φ as Structurally Random Matrix (SRM) [START_REF] Do | Fast and efficient compressive sensing using structurally random matrices[END_REF]. Therefore, A was formed by randomly taking a subset of rows from orthonormal transform matrices, that is, AA T = I M . As a consequence, there is no need to solve a linear system and the main computational cost becomes two matrix-vector multiplications per iteration.

SIMULATION RESULTS

In this section we evaluate the performance of the proposed compressive deconvolution algorithm for synthetic ultrasound data against a sequential resolution processing separately the CS reconstruction and the deconvolution. For the sequential approach, YALL1 [START_REF] Yang | Alternating direction algorithms for l1problems in compressive sensing[END_REF] was used to process the CS reconstruction from (2), while the deconvolution was done using an ℓ pnorm minimization. To solve this minimization problem, we also used one of the parameterizations of ADMM to split the ℓ p -norm term and the quadratic term, then solved them separately and iteratively.

The ultrasound data in Figure 1 and Figure 2 were simulated by 2D convolution between known PSFs and reflectivity images. For the first simulated image, the PSF was generated with Field II [START_REF] Jensen | A model for the propagation and scattering of ultrasound in tissue[END_REF] corresponding to a 3.5 MHz linear probe, sampled in the axial direction at 20 MHz. For the simulated kidney image, the PSF was also generated with Field II [START_REF] Jensen | A model for the propagation and scattering of ultrasound in tissue[END_REF] Fig. 1: The 1 st row represents the tissue reflectivity function, the ultrasound RF image and the corresponding B-mode image (obtained after envelope detection and log-compression from the RF image) from left to right; The 2 nd row depicts the compressive deconvolution results obtained with the proposed method for SNR=40dB and for CS ratios of 80%,60% and 40% respectively; The 3 rd row shows the corresponding B-mode images of the results in the 2 nd row.

Lateral distance [mm] Axial distance [mm]

Fig. 2: Simulated kidney image and its compressive deconvolution results for SNR=40dB. From left to right, the tissue reflectivity function, the B-mode image and the compressive deconvolution results for CS ratios of 80%, 60% and 40% respectively.

corresponding to a 4 MHz central frequency and a axial sampling frequency of 40 MHz. The amplitudes of the scatters in these two reflectivity images are random variables distributed according to GGD with the shape parameter set to 1.5. The number of scatterers was considered sufficiently large to ensure fully developed speckle.

The variables a 0 and λ 0 were initialized to zeros vec-tor and image. The rest of the hyperparameters, α, µ, β and γ, were set to their best possible values by crossvalidation. For the results in Figure 1, the parameters were set to {α, µ, β, γ} = {1e-6, 1e-3, 1e-3, 1e-1} and for those in Figure 2, {α, µ, β, γ} = {1e -5, 5e -3, 1e -4, 2e -4}. Figure 1 and Figure 2 show a series of compressive deconvolution results for the two tested data sets, for various CS ratios and for a SNR of 40dB. Quantitative results are shown in Table 1 and Table 2 and allow the comparison between the proposed method and the sequential approach. Three metrics were employed: the Blurred-PSNR (B-PSNR) (defined and used in [START_REF] Amizic | Compressive blind image deconvolution[END_REF]), the Blurred-Structural Similarity (B-SSIM, defined in the same manner as B-PSNR) and the Structural Similarity (SSIM) . Ten experiments were conducted for each CS ratio. We may remark that in all cases our approach provides better reults than the sequential scheme. We should remark that in both cases the same priors were used for the CS reconstruction and for the deconvolution.

CONCLUSIONS

In this paper, we firstly presented a compressive deconvolution framework for ultrasound imaging systems. Then we developed a compressive deconvolution algorithm based on Alternating Direction Method of Multipliers (ADMM). Based on the assumption of generalized Gaussian distributed tissue refelectivity function, our method uses an ℓ p -norm minimization for p between 1 and 2. The results obtained on simulated ultrasound images are promising and clearly show that the proposed simutaneously compressive deconvolution scheme performed better than the typical sequential approach. Future work will include the study of medical CS reconstruction algorithms such as [START_REF] Florescu | A majorize-minimize memory gradient method for complex-valued inverse problems[END_REF], the automatic estimation of p for the ℓ p -norm, the comparison with other existing methods such as [START_REF] Amizic | Compressive blind image deconvolution[END_REF] and experiments on clinical data.

  is a mathematical framework allowing to recover, via non linear optimization This work is partially supported by CSC (Chinese Scholarship Council).

Table 1 :

 1 Quatitative results for US image in Figure1PSNR 52.33 47,28 40,65 28,97 B-SSIM 99.41 98.27 93.24 48.72 SSIM 72.52 65.06 43.74 22.40 sequential B-PSNR 49,25 40,98 31,90 22,34 B-SSIM 98.91 93.37 64.20 23.84 SSIM 70.93 55.72 36.82 20.75

		CS ratio	80%	60%	40%	20%
	proposed	B-

Table 2 :

 2 Quatitative results for US image in Figure 2 PSNR 73.75 66.98 52.38 52.12 B-SSIM 100.0 99.99 99.64 79.11 SSIM 60.25 59.06 52.80 26.23 sequential B-PSNR 36.46 29.60 25.90 23.87 B-SSIM 88.36 64.31 44.88 32.37 SSIM 39.89 40.89 28.02 18.32

		CS ratio	80%	60%	40%	20%
	proposed	B-