Amira Zaylaa

Sébastien Ménigot

Jean-Marc Girault

Jamal Charara

Empirical Optimization Of Frequency Parameters In Chirp Inversion Imaging

Contrast ultrasound imaging has become a significant imaging technique for the last ten years. Physiological and pathological information are accessible thanks to microbubble injection and nonlinearity detection. To improve resolution and penetration, frequency modulated signals called chirps combined with compression filter have been utilized by radar technology. Recent studies applied the use of chirps in harmonic imaging. And since the setup parameters of the chirp are important to enhance the microbubble nonlinearities and thus the contrast; yet their full potential is not tapped. The present study aimed to empirically optimize the harmonic energy backscattered by microbubbles through a polynomial frequency modulation law of order three. Two imaging techniques were considered: harmonic chirp imaging and chirp inversion imaging. In simulation, the frequency parameters of chirps were sought for microbubbles of different sizes (1 µm and 2.5 µm). It was beneficial to increase the order of the polynomial chirp to improve the contrast. Chirp Inversion Imaging was preferable over harmonic chirp imaging using a polynomial chirp for contrast enhancement. Moreover, optimal frequency parameters were selected manually. Optimization of polynomial chirps was empirically capable of enhancing the contrast in images however, this study is paramountly conceived by an automatic technique.

Introduction

Contrast Ultrasound Imaging (CUSI) is a significant imaging technique requiring the injection of a gas-containing contrast agent; it is possible because microbubbles resonate when impinged upon by ultrasound and emit ultrasound at harmonic frequencies. Ultrasound harmonic images are produced by anatomic structures that contain these agents. Tissues that do not contain gaseous contrast agents presumably produce negligibal harmonic signals. The contrast between the medium infused with microbubbles and the surrounding tissues is given by the contrast-to-tissue ratio (CTR) [START_REF] Misaridis | Use of modulated excitation signals in medical ultrasound. part I:basic concepts and expected benefits[END_REF][START_REF] Hoskins | Diagnostic Ultrasound: Physics and Equipment (Cambridge Medicine)[END_REF].

Ultrasonic Contrast agents (UCA) are used in ultrasound imaging to increase the sensitivity of the imaging technique by altering the image contrast between different structures. A dramatic increase in scattering which can be achieved with this microbubble has made it useful both in difficult-to-image patients and in areas for which the ultrasound signal is reduced. Injection of the contrast microbubbles and observation of their real-time passage in an ultrasound image can significantly improve the sensitivity and diagnostic power of ultrasonic imaging. Basically, Harmonic Imaging can be achieved by various aquisition techniques such as second harmonic aquisition and pulse inversion.

Most applications of harmonic imaging implicates that the spectrum of the backscatter contains not only the fundamental but also the harmonics. Therefore, by second harmonic imaging a bandpass filter selects the second harmonic, however, the bandpass filter cannot avoid totally the spectral overlap between the fundamental and the harmonics. Conversely, by pulse inversion imaging the f o lobe is suppressed totally and we are left with the harmonics and/or subharmonics [START_REF] Shung | Diagnostic Ultrasound: Imaging and Blood Flow Measurements[END_REF]. Regardless of the imaging technique, whether second harmonic or pulse inversion imaging, one of the major problems that ultrasonic images suffer from is bad contrast, low resolution and small penetration. To improve both res-olution and penetration, frequency modulated signals known as "Chirps" often utilized in Radar technology were recommended in CUSI [START_REF] Misaridis | Use of modulated excitation signals in medical ultrasound. part I:basic concepts and expected benefits[END_REF][START_REF] Chiao | Coded excitation for diagnostic ultrasound: a system developer's perspective[END_REF]. Although the use of these chirps is capable of providing clinicians with images having good resolution and penetration in the future, the question remains what about the CTR?

Chirp imaging constitutes of sending frequency modulated signal, chirp excitation, into the medium under exploration, receiving the backscattered signal and filtering it in order to produce the image (Figure 1). The chirp could be linear or nonlinear; it depends on the frequency modulation profile or law. Traditional harmonic chirp imaging (HCI) has been utilized to reconstruct images at the second harmonic after transmitting chirps to the medium in order to improve the resolution and contrast.

Recent studies have utilized pulse inversion imaging to improve the CTR by using conventional pulses [START_REF] Pasovic | Second harmonic inversion for ultrasound contrast harmonic imaging[END_REF], however, conventional pulses have led a tradeoff between resolution and penetration. Others have used chirp inversion imaging (CII) based on a linear chirp to improve the signal-to-noise ratio (SNR) and solve the tradeoff between resolution and penetration [START_REF] Misaridis | Use of modulated excitation signals in medical ultrasound. part I:basic concepts and expected benefits[END_REF][START_REF] Song | Coded tissue harmonic imaging with nonlinear chirp signals[END_REF]. Exponential chirp was utilized to investigate nonlinear systems in audio signaling [START_REF] Novak | Nonlinear System Identification Using Exponential Swept-Sine Signal[END_REF] and hasn't yet been investigated in CUSI. Nonlinear chirp of order two has been optimized by HCI [START_REF] Pasovic | Second harmonic inversion for ultrasound contrast harmonic imaging[END_REF] merely. For that, optimization of the chirp's frequency law became profoundly essential generally due to the relative safety of UCA and particularly due to the robustness of chirp imaging in enhancing penetration, resolution and the CTR of medical images. The optimization of such chirps in diagnostic and theraputic ultrasound imaging is still flourishing and not until now did anyone investigated their full potential of chirps.

In this paper we aimed to empirically optimize the frequency modulation parameters of a third order polynomial chirp using Chirp Inversion Imaging (CII) in contrast media, i.e. instead of pulses in [START_REF] Simpson | Pulse inversion doppler: a new method for detecting nonlinear echoes from microbubble contrast agents[END_REF] we used chirps, with a small overview on harmonic chirp imaging (HCI) in order to enhance the contrast of contrast ultrasound images. We proposed an algorithm for CII and implemented HCI, in addition, we employed the compression filter in the two imaging techniques. Furthermore, since the energy or the power of the signal is related to the CTR, we computed the energy of the signal backscattered by 1 µm and 2.5 µm micro-bubbles during their nonlinear response and focused our study on improving these energy outcomes i.e. improving indirectly the CTR. Under this condition we anticipated that we not only get an image with good resolution and penetration but also with enhanced contrast that is capable of providing amended clinical diagnosis in the future. We tended to optimize a polynomial chirp of order higher than currently optimized chirps. The novelty of our approach was based on utilizing a polynomial chirp of order three and employing chirp inversion imaging in CUSI to obtain primary ultrasonic medical images of high contrast and consequently discriminate between regions that were not easily differentiated.

Methods And Materials

We simulated chirp imaging with Matlab (Mathworks, Natick, MA, USA). First, frequency modulated (FM) signals having the form of Eq. ( 1) were transmitted into the medium infused with contrast bubbles by means of a transducer using HCI and/or CII.

C(t) = G(t). cos((2π( f i (t))t), (1)
where f i is the instantaneous frequency, t the time and G(t) the gaussian envelope defined by:

G(t) = A exp (-( (α) 2 2(σ) 2 )) . ( 2
)
where A is the normalized amplitude of the chirp, α the shift of the envelope from the center and σ the width of the gaussian's envelope.

The instantaneous frequency we used in Eq. ( 3) also known as FM law was governed by a polynomial of order three:

f i (t) = f o t + β 1 t 2 + β 2 t 3 (3)
where f o is the starting frequency of the chirp in (MHz) and (β 1 , β 2 ) were the modulation indices in (THz/s) and (THz/s 2 ) respectively. In particular, when (β 2 =0) the chirp is said to be a Quadratic FM, otherwise it is a Polynomial FM [START_REF] Barlow | A theoretical investigation of chirp insonification of ultrasound contrast agents[END_REF][START_REF] Sboros | Absolute measurement of ultrasonic backscatter from single microbubbles[END_REF][START_REF] Misaridis | Potential of coded excitation in medical ultrasound imaging[END_REF]. Second, chirps made their way towards the medium. This resulted in backscattered signals which were further detected by the same transducer and fed into a compression filter to be finally interpreted and transformed into an energy image. Third, for different microbubble sizes and chirp types several simulations were carried out to optimize the frequency parameters of chirps. We implemented the compression filter as proposed in [START_REF] Borsboom | Harmonic chirp imaging method for ultrasound contrast agent[END_REF]. As compared to a regular compression filter, this harmonic compression filter has twice the instantaneous frequency at every time point. The actual chirp C(t) was given by Eq. ( 1), the filters' impulse response was defined by :

h(t) = [G(t)] 4 . cos((2π(2 f i (t))t) (4)
with f i (t) given by Eq. ( 3) and the output of the filter was given by y(t):

y(t) = C(t) * h(t) (5)
which described the compressed echo with good axial resolution and a central frequency situated at the second harmonic of the transmitted chirp. Also we looked up for the highest energy values in order to enhance the contrast since the energy is directly proportional to the CTR and is more accessible

Harmonic Chirp Imaging

It is a classical technique that permits the selection of the backscattered signal whose central frequency lies within the bandwidth of the transducer and is equal to the second harmonic value. This is equivalent to bandpass filtering. The backscattered signal is bandpass filtered such that the second harmonic is selected and the image is ready to be reconstructed at 2 f o [START_REF] Hoskins | Diagnostic Ultrasound: Physics and Equipment (Cambridge Medicine)[END_REF].

In HCI algorithm we:

• Specified the bandwidth of the transducer

• Fired a quadratic chirp C(t) given by Eq. ( 1)

• Simulated the backscatter e(t) = Echo{C(t)}, by solving the Modified Rayleigh-Plesset differential equation

• Fed e(t) into the compression filter and computed the energy image of the compression filter's output y(t) for different frequency values f o as a function of β 1

• Selected manually the optimum f o and its corresponding β 1

We further extended the algorithm to send a polynomial chirp (β 2 0), receive the backscattered signal, feed it into the compression filter and compute the energy image of the backscattered signal, but here, for different β 2 values as a function of β 1 . We finally selected the best β 2 and its corresponding β 1 manually from the image.

Chirp Inversion Imaging

It is the process of transmitting two pulses 180 0 out of phase sequentially into the medium infused with microbubbles. The received echoes are summed up such that echoes at the fundamental frequency will cancel out yet those at the harmonic frequencies will show up [START_REF] Hoskins | Diagnostic Ultrasound: Physics and Equipment (Cambridge Medicine)[END_REF].

In CII algorithm we proposed that:

• Two quadratic chirps C 1 (t) and C 2 (t) having the form of Eq. ( 1) with C 2 satisfying equation ( 6) are sent in such a way that a chirp was transmitted into the medium of micro-bubbles and, after a suitable delay t 0 , a second pulse was transmitted that is an inverted copy of the first pulse.

C 2 (t) = -C 1 (t -t 0 ) (6)
where C 1 (t) and C 2 (t) represents the acoustic pressure of the first and second pulses as a function of time, t.

• These signals were backscattered separately by microbubbles. The echo of each firing was then simulated:

e 1 (t) = Echo{C 1 (t)} e 2 (t) = Echo{C 2 (t)} (7)
• Then the following pulse inversion detection signal was manipulated:

s(t) = e 1 (t -t 0 ) + e 2 (t) (8)
• s(t) was fed into the compression filter and the energy of its output y(t) was computed for various central frequency values as a function of β 1 • The optimum f o is selected manually Moreover, the previous quadratic chirp was extended to a polynomial one having (β 2 0) and a fixed f o value, as previously described, the resultant backscattered signal was computed and led to a specific energy image and optimum β 2 and β 1 values after filtering.

Results and Discussion

In this part we will show the simulation results of the emitted quadratic and polynomial FM, particularly, the energy of the backscattered signal during the nonlinear response of micro-bubbles by traditional HCI and CII. We will show also the effect of the chirp's frequency parameters as well as discuss the results of the energy backscattered from 1 µm and 2.5 µm bubbles in HCI and CII .

Simulated Emitted Chirp

The chirp excitation given by Eq. ( 1) was depicted in Figure 2, the upper plot manifested the temporal representation of the third order polynomial chirp insonating the medium of micro-bubbles, the image exhibited its spectrogram with the energy indicated by the colorbar in an ascending order and the left plot represented merely its spectrum. The spectrogram topologically represented a small region having significant energy values. The highest energy value occured at t=1 µs when the chirp attained its maximum amplitude.

The Influence of Chirp's Frequency Parameters

In concern of the effect of the frequency parameters we simulated Figure 3 and Figure 4. Figure 3 represents the energy of the signal backscattered from 2.5 µm bubble by HCI and Figure 4 represeants the energy of the resultant signal backscattered by the same bubble by CII. For 2.5 µm bubble when insonifying the medium by a polynomial chirp of order 3 with (β 1 =0) it was vivid that the backscattered energy relied on both β 2 and f 0 in HCI and CII. Table 1: The variation of the maximum energy by CII and HCI of a= 2.5µm bubble as a function of the frequency parameters β 1 and β 2 in the frequency range f=[0.5 MHz, [START_REF] Song | Coded tissue harmonic imaging with nonlinear chirp signals[END_REF] MHz].

Frequency Parameters

β 1 =0 ∀ β 1 β 1 =0 ∀ β 1 =0 β 2 =0 β 2 =0 ∀ β 2 ∀ β 2
Max(E)-CII (dB) -31.6 -31.5 -30 -29.2

Max(E)-HCI (dB) -37.0 -36.4 -37.0 -36.3

In the range tested, optimal values of β 2 and f 0 existed in the regions where the best backscattered energy occured (dark red color). These optimal values were different in HCI and CII. The highest backscattered energy was attained in CII. The same remark was valid when transmitting a quadratic chirp, and studying the variation of ( f 0 ) as a function of β 1 in HCI and CII techniques. Moreover, we noticed that lines β 1 =0 and β 2 =0 in (Figures 3,4) were not the axis of symmetry this implied that the frequency parameters were asymmetric.

Many believed that nonlinear chirp signals were designed to minimize the peak range sidelobe level by optimizing its frequency sweep rate and time duration [START_REF] Song | Coded tissue harmonic imaging with nonlinear chirp signals[END_REF]. The major finding of our study manifested the strong impact of the FM law which could be quadratic or polynomial on the energy of the backscattered signal. Our data (see Table 1) exhibited that the higher the degree of the polynomial in FM law, the higher the energy was carried by the backscattered signal. In other words, as the FM law was more curved the energy received was dramatically higher.

From Table 1 we inferred that by increasing the degree of the FM law from: 1 to 2 an improvement of 0.1 dB for CII and 0.6 dB for HCI was observed, from: 1-3 an improvement of 1.6 dB for CII yet 0.0 dB for HCI was obtained and from: 2-3 (β 1 0 , β 2 0) 2.3 dB for CII and 0.1 dB for HCI was attained. The best results were those governing a third order polynomial governed by equation (2) using CII. However, in HCI the improvement of the Energy-to-complexity ratio was low as compared to that obtained for CII.

The Influence of Microbubble's Size

We impinged two bubbles of radius 2.5 µm and 1 µm respectively by a third oder polynomial chirp and summarized the results in Table 2. Table 2 represented the optimum frequency parameters f 0 * , (β 1 ) * and (β 2 ) * in the two imaging modalities utilized and for 2.5 µm and 1 µm bubbles. In summary, the maximum energy values, when the signal was backscattered by 1µm bubble, were much imperative than those when the signal was backscattered by 2.5 µm bubble. In addition, it was advantageous to use CII over CHI, since the backscattered energy obtained with CII is always superior of 7 dB as compared to CHI technique.

The optimal values depicted in Table 2 led to the best backscattered energy that were not only dependent on the degree of the FM law but also on the imaging technique. Consequently, searching for the best imaging technique and the best FM law are crucial. Note however that the optimal value of the frequency parameters varied with both the transmits' frequency ( f 0 ) and the size of the microbubble. Nevertheless, it was not empirically demonstrated in this study, it was obvious that the optimal value of each parameter depended not only on the transmits' frequency ( f 0 ) and the size of the microbubble but also on other parameters such as the bandwidth of the transducer and the initial pressure level.

Conclusion

Our study manifested that accounting for both the chirp type and the size of the micro-bubble it was necessary to investigate the frequency parameters that govern the behaviour of time-frequency transmit signal to optimize the energy reflected by microbubbles and indirectly the CTR. The third order Polynomial chirp insonation using CII have provided the optimal and most promising results as compared to CHI during the nonlinear response of microbubbles. Using the robust CII algorithm and a third order polynomial several optimum frequency parameters that tended to enhance the CTR were selected manually. Indeed the improvement of the CTR opened the door for high quality images and consequntly amanded diagnosis. Future step could be the validation of our findings experimentally for improving clinical diagnosis. And since the selection of the optimum parameters in our empirical study was done manually, we ought to consider an advanced automatic technique.

Figure 1 :

 1 Figure 1: Block Diagram of Contrast Ultrasound Imaging using chirp excitation.

Figure 2 :

 2 Figure 2: The transmitted polynomial Chirp.

Figure 3 :

 3 Figure 3: The energy of the compressed polynomial FM backscatter in HCI for 2.5 µm bubble.

Figure 4 :

 4 Figure 4: The energy of the compressed and backscattered polynomial FM in CII for 2.5 µm bubble.

Table 2 :

 2 The variation of the maximum energy and the modulation indicies as a function of the radius (2.5µm,1µm) of bubbles and the imaging method (HCI,CII) in the frequency range f=[0.5 MHz,6 MHz].

	Radius (µm)	2.5			1
	ImagingMethod	CII	HCI	CII	HCI
	Max(E) (dB)	-29.2	-36.3 -24.9	-28.1
	f 0	*	2.75	3.3	3.1	3.1
	(β 1 ) *	1.2	34.3	0.8	4.8
	(β 2 ) *	-20410 -0.8	2.4	-25010

Acknowledgements

This research has been supported by the National Council for Scientific Research in Lebanon and the Lebanese University.


Dumitru Potop-Butucaru

Isabelle Puaut

Analyse de pire temps de réponse intégrée d'applications multi-coeurs non-préemptives

Keywords: WCET, WCRT, analyse d'ordonnançabilité, système sur puce Integrated Worst-Case Response Time Evaluation of Multicore Non-Preemptive Applications WCET, WCRT, schedulability analysis, network-on-chip, system-on-chip

Les techniques d'analyse de pire temps d'exécution (WCET) ont atteint une très bonne précision dans l'analyse de programmes séquentiels s'exécutant sur architectures monoprocesseurs. Dans cet article, nous étendons une technique récente d'analyse WCET et l'outil associé pour permettre de calculer des estimations précises du temps de réponse (WCRT) d'applications parallèles non-préemptives s'exécutant sur des plates-formes multi-coeurs. La technique proposée est intégrée dans le sens où elle calcule en même temps les estimations WCET des fragments de code séquentiel et le WCRT global. L'utilisation de cette methode produit des estimations plus précises que les approches découplées plus classiques où le calcul du WCRT est réalisé à partir de valeurs WCET précalculées séparément pour chacun des fragments de code séquentiel. Sur 2 exemples d'applications de contrôle embarqué notre technique d'analyse améliore les estimations WCRT de 21% en moyenne.

Introduction

Multi-core systems are becoming prevalent in both general purpose systems and embedded systems. Their adoption is driven by scalable performance (speed, power, etc.) arguments, but this scalability comes at the price of increased software complexity. Indeed, multi-core systems run parallel software involving potentially complex synchronizations between the sequential programs executed on the various cores. To ease the temporal validation of real-time multi-core systems, we address in this paper the issue of evaluating the worst-case response time (WCRT) of parallel applications. We focus on the case of non-preemptive, statically scheduled code.

Worst-case response time (WCRT) estimation consists in providing safe upper bounds for the real-time duration between the trigger and the termination of some finite computation. The worst-case bounds are derived through static analysis techniques applied on a model of both hardware and software of the system. The fundamental trade-off in WCRT estimation is between model precision, computational tractability of the analysis, and generality of the technique in terms of the class of systems where it applies. Two main classes of techniques, usually applied sequentially, contribute to this trade-off :

-WCET (Worst-case execution time) estimation, which works on sequential programs, -WCRT estimation, that use WCET values as inputs.

WCET analysis emphasizes the importance of hardware micro-architecture. Indeed, in its double quest for execution speed and simplicity of the programming model, modern hardware architectures include user-transparent performance enhancing features such as pipelining, cacheing and arbitrated access to shared resources. The presence of these architectural elements significantly impacts execution time and complicates WCET estimation. Limiting generality to sequential code running on mono-processors (the classical form of WCET analysis) and selecting moderately complex hardware allows the preservation of computational tractability while the hardware micro-architecture is precisely modeled. This allows the computation of very tight execution time bounds.

WCRT analysis emphasizes the system-level complexity, by taking into account aspects such as the distribution of computations, task communication/synchronization, and interaction with the environment. The objective here is usually to provide execution time bounds on execution flows involving several tasks, possibly running on multiple processors, and their communications and synchronizations.

To limit computational complexity of WCRT analysis, hardware and software are usually represented in much less detail than in WCET estimation techniques. Typical objects at this level are sequential tasks with a simple read-compute-write execution model and characterized by large-grain functional and non-functional properties, such as : inputs and outputs, worst-case duration, period, execution conditions.

Ensuring the safety of the bounds computed through WCRT analysis necessarily relies on the safety of this large-grain characterization. One aspect here is the correctness of the quantitative characterization. For instance, durations must be computed through a WCET analysis of the task code (which takes into account micro-architectural details).

WCRT techniques are available either as stand-alone tools [START_REF] Kim | An ilp-based worst-case performance analysis technique for distributed real-time embedded systems[END_REF], or included in offline allocation and scheduling tools [START_REF] Grandpierre | From algorithm and architecture specification to automatic generation of distributed real-time executives[END_REF].

Contribution Isolating micro-architecture and system-level timing analysis has advantages, such as a reduction in the complexity of each phase, and a clear separation of concerns between system-level aspects, driven by WCRT estimation, and task-level aspects, driven by WCET estimation.

However, the major drawback of this isolation is a loss in precision which can be significant.

To consider only one aspect, micro-architectural analysis information (e.g. cache information) is lost at the frontier between tasks. This implies that to be safe, the analysis of each task must start from an undetermined micro-architecture state, which may result in overestimated WCETs, and consequently on pessimistic WCRTs. The contribution of this paper is an integrated WCRT analysis approach that handles jointly the task-level and system-level information. This is achieved by extending a state-of-the-art WCET estimation technique and tool to manage system-level information. Such a holistic approach is made possible by the use of deterministic and composable software and hardware architectures (multi-cores without cache sharing, off-line scheduling) as detailed later in this paper. We demonstrate the interest of the approach using an adaptive differential pulse-code modulation (adpcm) encoder where the integrated approach provides significantly tighter response time estimations than the more classical decoupled approaches.

Paper organization

The rest of this paper is organized as follows. Related work is briefly surveyed in Section 2. Section 3 presents the task model and defines more formally what is meant by worst-case response time. Section 4 details and motivates the class of multi-core architectures considered in this study. Section 5 defines our WCRT estimation method. Experimental results are given in Section 6. We conclude and discuss future work in Section 7.

Related work

Much research effort has been spent in the past in estimating Worst-Case Execution Times (WCETs) and Worst-Case Response Times (WCRTs).

Research on WCET estimation has mainly targeted software running on single-core architectures (see [START_REF] Wilhelm | The worst-case execution-time problem overview of methods and survey of tools[END_REF] for a survey of methods and tools). A lot of effort has been put on so-called low-level analysis, allowing architectures with caches and in-order pipelines to be analyzed precisely. The research presented in this paper is not a new WCET estimation technique, but rather takes benefit of state-of-the-art low level analysis to produce tight WCRT estimates of parallel applications.

Many WCRT estimation methods estimate end-to-end response times of distributed applications communicating using message passing [START_REF] Tindell | Holistic schedulability analysis for distributed hard real-time systems[END_REF][START_REF] García | Schedulability analysis of distributed hard real-time systems with multiple-event synchronization[END_REF] on multiprocessor systems [START_REF] Bertogna | Response-time analysis for globally scheduled symmetric multiprocessor platforms[END_REF][START_REF] Kuo | Efficient wcrt analysis of synchronous programs using reachability[END_REF]. To our best knowledge, all these methods can be qualified as decoupled , in the sense that they use as input WCET estimations of code snippets (computed before the WCRT analysis). By comparison, we argue that an integrated analysis that computes WCETs and WCRTs altogether allows to produce tighter WCRTs than a decoupled approach, because it allows effects between code portions, like for instance cache effects, to be captured accurately.

The research we found to be closer to our approach is described in [START_REF] Rochange | Wcet analysis of a parallel 3d multigrid solver executed on the merasa multicore[END_REF][START_REF] Nemer | Inter-task wcet computation for away instruction caches[END_REF][START_REF] Kim | An ilp-based worst-case performance analysis technique for distributed real-time embedded systems[END_REF].

Paper [START_REF] Rochange | Wcet analysis of a parallel 3d multigrid solver executed on the merasa multicore[END_REF] is devoted to WCET estimation of a parallel application running on a predictable multi-core architecture. Similarly to our work, emphasis is put on predictability of the hardware and software architectures (Merasa predictable multicore architecture at the hardware level, statically-defined synchronizations and communications at the software level). However, in contrast to [START_REF] Rochange | Wcet analysis of a parallel 3d multigrid solver executed on the merasa multicore[END_REF] that provides formulas to combine WCETs of code snippets to obtain the WCET of the parallel application, in our work the application is analyzed as a whole. As a consequence, we are able to exploit knowledge of the hardware state between code snippets and thus can provide tighter estimations, especially for fine-grain parallelism.

In [START_REF] Nemer | Inter-task wcet computation for away instruction caches[END_REF], a method to determine residual cache states after a task execution in mono-core platforms, is provided. The method allows to obtain tighter WCETs in case of repetitive task executions. Using our method, we obtain the same benefits, but without needing a specific analysis.

Inria

This benefit comes as a side product of our method because WCRT computation is integrated into a WCET estimation tool.

Paper [START_REF] Kim | An ilp-based worst-case performance analysis technique for distributed real-time embedded systems[END_REF] proposes an ILP formulation for WCRT computation of parallel applications running on multi-core systems. The method computes the application WCRT given a task mapping, architecture and scheduling policy, with contentions when accessing shared resources. Unlike [START_REF] Kim | An ilp-based worst-case performance analysis technique for distributed real-time embedded systems[END_REF], we currently rule out resource contentions, that is left for future work. However, contrary to [START_REF] Kim | An ilp-based worst-case performance analysis technique for distributed real-time embedded systems[END_REF], our analysis of the hardware is expected to be tighter because of our integrated approach, through the integration of WCRT computation into a WCET analysis tool.

[ i ]= tqmf [ i -2 ] ; } tqmf [ 1 ] = x i n 1 ; tqmf [ 0 ] = x i n 2 ; // R

Task model

The simplest embedded control systems running on mono-processor architectures follow a socalled simple control loop paradigm. In such systems, the software is simply a loop whose body is the sequence of calls to the various input sampling, processing, and actuation functions/tasks. Scheduling (i.e. order of the sequence of calls) is non-preemptive and fixed in an off-line fashion.

In this paper, we consider the multi-processor equivalent of simple control loops. More precisely, we consider systems where each processor/core executes a simple control loop. We shall denote with τ i the program that forms the body of the loop executed by processor/core CPU i , 1 ≤ i ≤ n. Tasks τ i , 1 ≤ i ≤ n can communicate with each other through a set of messagepassing channels C = {c 1 , . . . , c m } which are bounded FIFO buffers that do not loose, duplicate, or corrupt messages. Communication is done using send and receive primitives that can be used at any known position in the programs τ i . The two primitives are blocking (send on full channel, receive on empty channel), which means that channels can be used for synchronization.

For the scope of this paper, we make the following assumptions concerning the channels :

-Each channel connects exactly two processors (one sender and one receiver).

-Task τ i contains a single send or receive primitive call for each channel connected to CPU i .

-Each channel allows the storage of only one message. No assumption is made on how the message passing channels are implemented on the execution platform. For instance, the execution platform we used for experimentation implements the message passing primitives over a distributed shared memory physical architecture.

Finally, for the definition of our analysis technique, we make the following assumption concerning scheduling : We assume that the loops of the various processors/cores advance in lockstep, meaning that iteration m of τ i starts at the same date on all processors/cores. This can be ensured using channel-based synchronization. Note that this hypothesis means that system execution features system-wide iterations involving all CPUs and with a defined start date.

Illustrating example

An illustrating toy application, that will be used all along the paper, is depicted in Figure 1. The application is a portion of a bi-processor parallel version of the adpcm (adaptive pulse code modulation algorithm) from the Mälardalen WCET benchmark suite [START_REF] Gustafsson | The Mälardalen WCET benchmarkspast, present and future[END_REF]. We emphasized with arrows the two send/receive pairs associated with channel1 and channel2 respectively.

WCRT estimation problem formulation

Assuming the previously-defined task model, the worst-case response time analysis problem we solve in this paper is the following.

Definition 1 (WCRT problem formulation) Compute the worst-case duration of an iteration of the system, from start to the start of the next iteration.

To allow timing analysis, we shall assume that each program τ i , considered in isolation, satisfies the classical requirements allowing WCET analysis. In particular, all loops it contains (except the main loop of each task, that is by definition infinite) have statically bounded numbers of iterations. We also assume that inter-task communications is free of deadlocks by construction.

Execution platform

In the context of multi-processor systems, we believe that obtaining precise timing information is only possible if the hardware meets some basic assumptions ensuring that timing interference during concurrent access to shared resources can be tightly bounded. By shared resource we mean here all components of the hardware system that may introduce functional and temporal interferences between components. The shared resources we consider in our work are the memory subsystem, the buses and networks (including I/O), DMA controllers, and the synchronization subsystem. Multi-processor systems with shared caches, although amenable to WCET estimation [START_REF] Li | Timing analysis of concurrent programs running on shared cache multi-cores[END_REF][START_REF] Hardy | Using bypass to tighten wcet estimates for multicore processors with shared instruction caches[END_REF], may yield pessimistic WCET estimates, because the state of these caches becomes difficult to approximate in the presence of concurrent requests. Our choice is to consider architectures where each processor has its own cache subsystem, independent from the ones of other processors. This general hypothesis covers caches with partitioning mechanisms, but in the experimental section of this paper we only consider cores with private L1 instruction and data caches, and no other caches.

For the scope of this paper, we make the following assumptions concerning the caches :

-Each core has separate and private instruction and data caches. We consider such architectures because separate caches are analyzable more precisely than unified caches by WCET estimation techniques. -All caches have a Least Recently Used (LRU) replacement policy. LRU is assumed because it was shown to be the most predictable cache replacement policy [START_REF] Reineke | Timing predictability of cache replacement policies[END_REF]. Aside from shared caches, another significant source of WCET estimation imprecision is the presence of shared memory banks and shared communication busses. Our choice here is to consider architectures where the duration of all memory accesses and data transmissions can be precisely determined. As pictured in Fig. 2, the timing precision can be ensured :

-Fully by hardware mechanisms, for instance through the use of time division (TDM) memory controllers. -Through a mix of software and hardware mechanisms. In these cases, software and/or hardware synchronization mechanisms (semaphores, locks) are used to guarantee the absence of contentions due to access to RAM banks or communication lines. In this paper we consider architectures of the second type, as this case covers classical distributed bus-based architectures, shared memory architectures featuring multiple RAM banks, but also mixes of the two, such as the Network-on-Chip (NoC) based architectures proposed by various vendors [START_REF] Harrand | Network on chip with quality of service[END_REF][START_REF] Benini | Programming heterogeneous many-core platforms in nanometer technology : the p2012 experience[END_REF]. In particular, our NoC-based experimentation platform falls in this last case, and the aforementioned synchronization mechanisms are also used to control access to DMA controllers.

WCRT computation

Our approach to WCRT computation consists in extending a state-of-the-art WCET estimation method to compute WCRTs of parallel applications. The state-of-the-art WCET estimation technique we use as base for this extension is presented in Section 5.1 ; Section 5.2 then presents the extension itself.

x1 = x f 1 + x n 1 x f 1 ≤ 1 ... //

Existing state-of-the-art WCET estimation technique

Static WCET estimation techniques are organized in three phases performing different analyses [START_REF] Wilhelm | The worst-case execution-time problem overview of methods and survey of tools[END_REF] : Control-flow analysis, Hardware-level analysis and WCET calculation.

Control-flow analysis This phase extracts information about possible execution paths from the program source or binary. The output of this phase is a data structure representing a superset of possible flows, represented by either a control flow graph (CFG) or a syntax tree. For the scope of this paper, we will represent the program control flow using a CFG, extracted from the program binary, and annotated with additional flow information such as maximum number of loop iterations. The CFG of the loop body of function core1 (from our sample application) is given in Fig. 3.a. In this figure, basic blocks of code are represented with circles, and control flow between basic blocks with arrows. Loops, annotated with their maximum number of iterations, are represented as rectangles.

Hardware-level analysis This step, also called low-level analysis, estimates the worst-case execution times of portions of sequential code, typically basic blocks. The difficulty during this task is to take into account micro-architectural components of (all) modern processors, such as caches, pipelines, branch predictors. In the presence of such components, the execution time of a statement is dependent on the context it is called in.

The hardware-level analysis is often made of several sub-analyses, one for each considered hardware component (instruction and data caches, pipelines, branch predictors, ...). The overall Inria typical outcome of hardware-level analysis is a maximum execution time per basic block in two different contexts to cope with cache effects : the first execution of the basic block, denoted t f and its subsequent executions, denoted t n ; (negative) execution times may also be associated to edges to account for pipeline effects between basic blocks.

WCET calculation

The purpose of this final phase is to determine an estimate for the WCET, based on the flow and timing information derived in the previous phases. The most widespread calculation method, that will be adopted in this paper, is called implicit-path enumeration (IPET). In IPET, program flow and basic-block execution time bounds are combined into sets of arithmetic constraints. Each entity (basic block or program flow edge) in the task is assigned two values :

-A time coefficient, denoted t entity , which expresses the upper bound of the contribution of that entity to the total execution time every time it is executed. -A count variable (x entity ), corresponding to the number of times the entity is executed. An upper bound of the WCET of the program is determined by maximizing the sum of products of the execution counts and times ( i∈entities x i * t i ), where the execution count variables are subject to constraints reflecting the structure of the task and possible flows. The result of an IPET calculation is an upper timing bound and a worst-case count for each execution count variable.

Fig. 3.b illustrates the constraints and formulas generated by an IPET-based bound calculation method on the loop body of core1, assuming it is the program entry point. The start constraint states that the code is executed once. The structural constraints reflect the possible program flows, meaning that each basic block must be entered the same number of times as it is exited. The loop bounds constrain the number of executions of basic blocks inside loops. Cache induced constraints express that basic blocks have different execution times, one for their first execution, another for the next ones. In the WCET expression to be maximized, there are two execution durations per basic block to model cache effects. For instance, in the formula, t f 2 = 120, whereas t n 2 = 14.

WCRT computation

Starting from the WCET estimation method sketched above, we build our WCRT estimation technique for parallel applications. We were able to do so by only modifying the control flow analysis phase to make it :

-Add new edges between basic blocks to model synchronization/communications between tasks. -Incorporate in the CFG the low-level information we need concerning inter-processor communications. The 3 phases of the modified analyzer run as follows :

1. A new step dubbed Modelling of communications, is attached at the end of Control flow analysis. Once the classical control flow analysis runs (separately) on the functions/tasks that are put in parallel, the new step adds new edges between the control flow graphs of these functions, the result being a single CFG. One new edge is added for each inter-task communication.

2. Hardware-level analysis. Application of the hardware-level analysis of the WCET estimation method runs unmodified on each function/task taken in isolation. During this first step, each task is analyzed as if it was not communicating with the other tasks executed on the other cores. 3. WCRT computation. Even though the CFG corresponds to a parallel application, and has slightly different topological properties, the WCET calculation can be run unmodified. This is due to the fact that the analysis works by finding the critical path in a directed acyclic graph (under certain constraints). That the graph represents purely sequential behaviors, or parallel ones (including the new edges that model communications) is not important.

The method is illustrated on our toy application of Fig. 1. Fig. 4 depicts the control flow graphs of the loop bodies of each of the two functions core1 and core2. In the figure, the shaded areas labeled st1-st6 correspond to sub-tasks, which are by definition parts of the two functions that are separated by inter-task communication primitives. For instance, node N4 of the core1_body CFG is the basic block containing the send call on channel1.

After the classical control flow analysis phase runs unmodified on the various sequential functions of the parallel program, the modelling of communications adds new edges in the application CFG to model communications. This is depicted by bold red arrows in Fig. 5. These new edges correspond to :

-Message passing between tasks (edges N 4 → N 11 and N 19 → N 9).

-Parallel launching of tasks on the different cores (edges to and from nodes N 21 and N 22 in the application entry point main). Every new edge is associated a duration to model its execution time (message transmission time for communications). During the hardware analysis phase, our WCRT analysis method applies instruction and data cache analysis and pipeline analysis on the two CFGs core1_body and core2_body considered in isolation. This allows to benefit from the tightness of hardware-level analysis on each subtask. For instance, in CFG core1_body, it allows to detect that array tqmf is still in the data cache after calling primitive send . This would not have been possible if a decoupled approach was used (WCET estimation of sub-tasks followed by WCRT estimation). If a decoupled method was used, conservative assumptions would have to be taken such that the analysis of sub-tasks is safe (assuming the worst-case hardware state, i.e. empty cache at WCET analysis start). Using an integrated approach, the first step of the analysis is able to capture hardware effects between sub-tasks (instruction caches, data caches, pipeline) naturally.

Finally, the WCET computation step is applied unmodified. Thanks to the introduction of the new edges, new constraints are automatically added in the WCET calculation equations (see new constraints below for the example), and communication delays are automatically taken into account. Given that our claims mainly concern the precision of the timing analysis, we considered an evaluation platform allowing us to perform very precise (cycle-accurate) estimations and measurements of execution time, in both single-processor and multi-processor cases.

We achieved this by using the SocLib library [START_REF]SoClib : an open platform for virtual prototyping of multi-processors system on chip[END_REF] for virtual prototyping of multi-processor systems-on-chips (MPSoC). The hardware components we use are of cycle-accurate, bit accurate type, written in SystemC [START_REF] Grotker | System Design with SystemC[END_REF]. Synthesizable VHDL models exist for these components, allowing FPGA prototyping (although in this paper we only used the SystemC compiled simulators).

The precise architecture model we worked on using SocLib is a scaled-down version of that of [START_REF] Djemal | Programmable routers for efficient mapping of applications onto noc-based mpsocs[END_REF]. While the original platform scales up to 4096 cores, we have only used single-, double-, and quad-core configurations for our tests. Our MPSoC is organized as a set of 1, 2, or 4 singleprocessor tiles, each one with the architecture of Fig. 6. As mentioned before, to facilitate cache analysis, each core has separate L1 instruction and data caches, and separate instruction and data memory. In the experimental setting of this paper all caches implement a Least Recently Used (LRU) cache replacement policy, and feature 32 sets, 4 ways, and 8 4-byte words per cache line. All CPU cores are of the same type, using the MIPS32 instruction set.

Processing elements are interconnected in a 2-D mesh using the DSPIN Network-on-Chip provided by SocLib. The quad-core configuration, for instance, is obtained by connecting 4 tiles in a 2x2 mesh, as pictured in Fig. 7 (the arrows represent communication lines of the Networkon-Chip).

The resulting hardware platform follows a distributed shared memory model, where all memory banks are assigned addresses in a global address space, and can be addressed by all CPU cores. Of course, memory access is of NUMA (Non-Uniform Memory Access) type, meaning that access time depends on the distance between the CPU and the addressed RAM bank. In building our tiles, we followed classical techniques for reducing access contentions. Thus, the local interconnect of each tile is a full crossbar, and each tile provides multiple RAM banks to allow non-interferent accesses to program text and data by the CPU cores (program and data) and DMA units. A set of hardware locks allows low-overhead synchronization between the CPU cores.

Inria

WCET estimation tool

The Heptane static WCET estimation tool [START_REF] Colin | A modular and retargetable framework for tree-based wcet analysis[END_REF] was used for the experiments. Heptane implements the its Implicit Path Enumeration Technique (IPET) WCET calculation technique. Regarding hardware-level analysis, Heptane supports instruction caches, data caches and pipeline analysis. Heptane's pipeline analysis module was customized to analyze the pipeline structure of each core.

Modifying Heptane to calculate WCRTs mainly consisted in adding a new analysis pass to Heptane, interposed between hardware-level analysis and WCET calculation, for a total volume of 200 lines of C++ code. Inter-task communications are detected by Heptane through annotations in the source code of the analyzed application, specifying at each communication point the recipient of the message and the communication latency.

Studied applications

The proposed WCRT estimation method was experimented on two signal processing applications. The first one is the complete parallel version of the adaptive differential pulse-code modulation (adpcm) from the Mälardalen WCET benchmark suite [START_REF] Gustafsson | The Mälardalen WCET benchmarkspast, present and future[END_REF]. The global dataflow of the code executed at each iteration of the modulation applications is depicted in Fig. 8, where boxes represent subtasks and arrows communications/synchronizations between subtasks.

When parallelized for a 2-core architecture, the mapping of subtasks to cores is depicted in Fig. 9. Only the arrows crossing CPU boundaries are coded as communications ; the sequencing of QMF, low-band encoder and multiplexer is implemented simply by calling successively the When parallelized for a 4-core architecture, every subtask is assigned to a different core, and software pipelining is used to authorize more parallelism from the application. The communication latency of every inter-core communication was determined by an analysis of the hardware platform as a formula dependent on the volume of data to be transferred.

The second application is a simple load balancing example, where two processors are needed to improve the throughput of a simple image filter. In this bi-processor application, processor 0 successively receives image lines in a buffer. The buffer content must be stored elsewhere to allow a new line to arrive, and this new line will be sent to processor 1, cyclically.

Application code was compiled using a standard GNU Mips compilation toolchain.

For the scope of this performance evaluation, application code was parallelized manually. Automatic code parallelization software like [START_REF] Bondhugula | A practical automatic polyhedral parallelizer and locality optimizer[END_REF] or offline real-time scheduling tools like [START_REF] Grandpierre | From algorithm and architecture specification to automatic generation of distributed real-time executives[END_REF] that generate efficient parallel code could have been used instead.

Experimental results

Demonstrating the tightness of our WCRT estimation method requires to : (i) demonstrate that the accuracy of the hardware model used in the base timing analysis tool with respect to the analyzed architecture ; (ii) compare WCRTs obtained using our integrated approach against those obtained using a decoupled estimation method ; (iii) compare WCRTs obtained with our tool with actual execution durations of the parallel code.

Accuracy of hardware model To show the accuracy of the hardware model used in the analysis, we have compared the number of cycles returned by the SocLib simulation software with the one returned by the Heptane timing analysis tool. Experiments were conducted on randomly generated sequential code, starting with known contents of the instruction of data caches. After a careful, extensive and always fastidious comparison of the analyzer and simulator cycle counts, both tools returned exactly the same number of cycles for all considered code.

Comparison with baseline WCRT estimation method

To evaluate the tightness of WCRT estimates, we have compared them with a baseline approach that estimates WCETs and WCRT separately. The baseline method :

-computes WCET of all sub-tasks (sequences of code between communications) separately ; to be safe, the worst-case hardware state is assumed by the static analyzer at the start of every subtask ; -WCRTs are computed in an ad hoc manner according to the synchronization pattern of each task ; this turned out to be very easy for the considered applications, that have simple and regular communications, never more complex that the ones illustrated in Figure 1.

Inria Table 1 -Worst-case response times of the considered applications using our approach (Integrated) and the baseline approach (Isolated). The improvement over the baseline approach is defined as The computed response time bounds are given in Table 1. The estimations produced by our integrated approach are always tighter than using the isolated method (21% in average on the three applications). The gain varies depending on the dependencies between subtasks of tasks running on the same core. When the amount of dependencies is high, like in application adpcm, that features intensive code and data reuse between subtasks, the gain is significant. When the amount of dependencies is smaller like in filter (no reuse of data, modest reuse of code between subtasks), the gain is much smaller.

Comparison with observed response times The pessimism of our WCRT evaluation method is evaluated by comparing estimated WCRTs with observed ones, estimated using the SocLib simulation software. Regarding simulation results, due to time constraints, we made no attempt to identify the worst-case input data and execute the code with typical input data, not necessarily representative of the worst-case situation. The estimated pessimism is thus an upper bound of the method pessimism. Results are reported in Table 2.

The numbers show that even without executing the code using its worst-case input data, the results are encouraging : estimated and measured response times are close to each other (of 7.7% in average). Further experiments need to be conducted to identify the actual overestimation and not only an upper bound of the overestimation.

Conclusion

We have presented in this paper a method to compute the WCRTs of parallel applications running on multicore platforms. WCRTs computation is integrated into a WCET estimation method such that hardware effects across portions of the application are dealt with naturally. We have demonstrated that our approach produces WCRTs that are tighter than using decoupled computations of WCETs and WCRTS (of 21% in average on three signal processing applications). Preliminary experiments show that the WCRT overapproximation is below 7.7% in average. The modification of a state-of-the-art WCET estimation tool turned out to be easy to achieve. We believe that our method can be integrated easily in other WCET estimation tools using the implicit path enumeration techniques to the extent that the analysis framework is sufficiently modular (hardware-level analysis and WCET computation are clearly separated, and new analysis passes can be inserted in between).

In this paper, assumptions have been made regarding the software structure in order to demonstrate the validity of our approach on simple but yet realistic setting. In our future work, our first objective will be to relax as much as possible these assumptions to broaden the scope of application of the approach. Another area for future research will be to use obtained WCRTs to refine task placement and scheduling.

Figure 1 -

 1 Figure 1 -Toy illustrating example : Code snippets of a parallel version of the adpcm benchmark from the Mälardalen WCET benchmark suite[START_REF] Gustafsson | The Mälardalen WCET benchmarkspast, present and future[END_REF]

Inria

Figure 2 -

 2 Figure 2 -Examples of execution platforms allowing simple contention avoidance : (a) Shared multi-bank RAM (, (b) Shared single-bank RAM with TDMA arbiter, (c) Distributed, NUMA, or network-on-chip architecture (computing cluster = a or b).

  CFG of the loop body of core1 // Start constraint x1 = 1 // Structural constraints x1 = x1,2 x2 = x1,2 + x2,3 = x2,3 + x2,4 x3 = x2,3 = x3,2 x4 = x2,4 = x4,5 ... // Loop bound constraints x3 ≤ 12 x7 ≤ 22 // Cache-induced constraints

Figure 3 -

 3 Figure 3 -State-of-the-art WCET calculation on the loop body of core1 in our illustrating example.

Figure 4 -

 4 Figure 4 -Original control flow graph for the toy example. Shaded areas represent task portions between communications (subtasks)

Figure 5 -

 5 Figure 5 -Control flow graph of our small example after adding communication arcs between sequential tasks.

Figure 6 -

 6 Figure 6 -A processing element

Figure 7 -Figure 8 -

 78 Figure 7 -The 4-core 2x2 mesh configuration

Figure 9 -

 9 Figure 9 -Allocation and scheduling of adpcm on two processors

Table 2 -

 2 Worst-case response times of the considered applications using our approach (Integrated) compared with observed response times (Measured). The pessimism is defined as

	Integrated-Measured Measured	* 100
	Name	Integrated Measured	Pessimism
			(cycles)	(cycles, typical input) (%)
	adpcm -2 cores 73563	64944	13.3%
	adpcm -4 cores 44568	41468	7.5%
	filter -2 cores	110825	108296	2.3%

RR n°8234

Potop & Puaut

RR n°8234RESEARCH CENTRE PARIS -ROCQUENCOURTDomaine de Voluceau, -Rocquencourt B.P. 105 -78153 Le Chesnay Cedex Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr ISSN 0249-6399