Damien Hardy
email: damien.hardy@irisa.fr

Benjamin Lesage
email: benjamin.lesage@irisa.fr

Isabelle Puaut
email: isabelle.puaut@irisa.fr

Scalable Fixed-Point Free Instruction Cache Analysis

Keywords: Hard real time systems, worst-case execution time estimation, instruction caches

Estimating worst-case execution times (WCETs) for architectures with caches requires the worst-case number of cache misses to be upper bounded. Most existing static cache analysis methods use fixed-point computation and do not scale well with large code sizes. To address this scalability issue, we propose in this paper a new fast and scalable instruction cache analysis technique. In contrast to existing work, neither fixedpoint computation nor heavyweight interprocedural analysis are required. Thus, code sizes too long to analyze with existing techniques are then analyzable with lower analysis time and memory consumption, and with only a slight degradation of the analysis precision. Experimental results show a reduction of the analysis execution time of a factor 5 in average (with a peak near 30 for the largest and most complex code) with only a degradation of the analysis precision of 0.5% in average. The proposed technique is intended to be used in situations where the need for fast analysis outweighs the need for very tight results: early software development phases, timing analysis of large pieces of software, or iterative WCET-oriented compiler optimizations.

I. INTRODUCTION

Caches aim at filling the ever broadening gap between fast micro-processors and slow memories. By taking advantage of the temporal and spatial locality inherent to most applications, caches improve average-case memory access latencies. This positive impact on access latencies, and thereby systems' throughputs, led to the wide adoption of architectures comprehending caches even in the context of embedded real-time systems.

A requirement of hard real-time systems is to prove that tasks meet their deadlines in all execution scenarios, including the worst-case. To attain this level of confidence, the worst-case execution time (WCET) has to be estimated for any task in the system. WCET estimates must offer both safety and tightness properties. Safety is the guarantee that the estimated WCET is greater than or equal to any possible execution time. Tightness means that the estimated WCET is as close as possible to the actual WCET to correctly dimension the resources required by the system.

WCET estimation methods have to account for all possible flows in a program and determine the longest execution path (so-called high-level analysis). They also have to account for the hardware the program is running on, through a low-level analysis. A number of static analysis methods have been designed in the last two decades at both levels [START_REF] Wilhelm | The worst-case execution time problem-overview of the methods and survey of tools[END_REF].

The presence of caches in real-time systems makes the low-level analysis harder due to the dynamic behavior of caches. Safely estimating WCETs on architectures with caches requires knowledge of all possible cache contents at every program point, which in turn requires knowledge of the cache replacement policy.

Much research has been undertaken, during the last two decades, with the objective of predicting WCETs in architectures equipped with caches. In most static cache analysis methods, fixed-point analysis is applied on the control flow graph of tasks to determine all possible cache contents in every execution point. The analysis result is then a classification of the outcome of every memory reference in the worst-case execution scenario (e.g. always-hit, alwaysmiss, first-miss, etc.). Although shown to be tight, fixed-point based methods can be very time and memory consuming. To overcome these scalability issues, we propose in this paper a fixed-point free instruction cache analysis. This analysis uses standard information extracted from the control flow graph and commonly used in compilers [START_REF] Muchnick | Advanced compiler design and implementation[END_REF]: dominator information, and loops entries and back-edges. Thanks to this information, the proposed analysis identifies from the control flow graph the most common reference patterns resulting in cache hits (among other reference patterns, sequential execution of instructions stored to the same memory block, execution within loops of code small enough to stay resident in the cache). Our method is modular in the sense that each reference pattern can be applied independently of the others depending on the analysis requirements in terms of tightness and analysis cost.

A fast and simple basic analysis is first presented. Two extensions are then proposed to tighten WCET estimates, at the price of a reasonable extra computation time.

The goal of our method is not to provide results as tight as existing analyses. Rather, our goal is to allow the user to trade analysis tightness for analysis efficiency. Experimental results show that WCET estimates are close to an existing state-of-the-art analysis [START_REF] Theiling | Fast and precise WCET prediction by separated cache and path analyses[END_REF], whereas on the largest piece of code analyzed, the analysis run-time is only one-thirtieth of the time required by the state-of-the-art approach.

The proposed technique is thus intended to be used in situations where the need for fast analysis outweighs the need for very precise results. A typical situation is early software development phases, during which the software code is modified frequently and tightness of WCET estimates is not that crucial. Another case is the analysis of large and complex software, for which the tightest analyses are too time and memory consuming. Finally, some WCET-oriented optimizing compilers, such as WCC [START_REF] Falk | A compiler framework for the reduction of worst-case execution times[END_REF], may need to iteratively estimate WCETs in order to apply optimizations along the longest execution path. Hence, they require to estimate WCETs as fast as possible even if these estimates are slightly less tight.

The remainder of the paper is organized as follows. Section II first compares our work with related work. Assumptions and notations are given in Section III. Our cache analysis method is presented in Section IV. Experimental results are given in Section V. Finally, Section VI concludes the paper and gives directions for future work.

II. RELATED WORK

Many WCET estimation methods have been designed in the last decades (see [START_REF] Wilhelm | The worst-case execution time problem-overview of the methods and survey of tools[END_REF] for a survey). Among them, the mostly used WCET computation technique, called IPET (Implicit Path Enumeration Technique) [START_REF] Malik | Performance analysis of embedded software using implicit path enumeration[END_REF], estimates the WCET through the resolution of an Integer Linear Programming (ILP) problem constraining the execution frequencies of the program's basic blocks.

Static WCET estimation techniques need a hardware analysis phase to determine the worst-case timing behavior of the microarchitectural components: pipelines and out-oforder execution [START_REF] Engblom | Pipeline timing analysis using a trace-driven simulator[END_REF], [START_REF] Li | Modeling out-oforder processors for WCET estimation[END_REF], branch predictors [START_REF] Colin | Worst case execution time analysis for a processor with branch prediction[END_REF] and caches. Regarding caches, two main classes of approaches have been designed: integrated approaches analyse the cache behavior simultaneously with WCET computation, whereas decoupled approaches perform a separate cache analysis whose results are subsequently used by the WCET computation step.

Among integrated approaches, the methods presented in [START_REF] Li | Efficient microarchitecture modeling and path analysis for real-time software[END_REF], [START_REF]Cache modeling for real-time software: beyond direct mapped instruction caches[END_REF] generate an ILP system that both computes the WCET and captures cache effects, through variables in the ILP computation that model for each memory reference if the reference hits the cache or not. Such a fine-grain cache modeling, as reported in [START_REF]Cache modeling for real-time software: beyond direct mapped instruction caches[END_REF], might result in very long computation times for set-associative caches (e.g. more than 15 hours to analyze the jpeg benchmark (5.5KB) for a 2-way set-associative instruction cache). In [START_REF] Lundqvist | Timing anomalies in dynamically scheduled microprocessors[END_REF], a cycleaccurate processor simulator is extended to enable symbolic execution with unknown input data values; it uses alternative instruction semantics to handle unknown operands. While shown to provide tight WCET estimates, a limitation of the approach is its computational complexity, mainly due to the use of a processor simulator and a (limited) exploration of feasible paths.

For the sake of separation of concerns and analysis efficiency, decoupled approaches analyse the cache behav-ior separately. Static cache simulation [START_REF] Mueller | Static cache simulation and its applications[END_REF], [START_REF]Timing analysis for instruction caches[END_REF], based on dataflow analysis, and the methods described in [START_REF] Ferdinand | Cache behavior prediction for real-time systems[END_REF], [START_REF] Theiling | Fast and precise WCET prediction by separated cache and path analyses[END_REF], [START_REF] Ferdinand | Reliable and precise WCET determination for real-life processor[END_REF], based on abstract interpretation, both use fixedpoint computation to determine the cache contents at all program points, considering all execution paths altogether. These two classes of methods provide a classification of the outcome of every memory reference in the worst-case execution scenario (e.g. always-hit, always-miss, first-miss, etc.). These methods were originally designed for code only, and for direct-mapped or set-associative caches with a Least Recently Used (LRU) replacement policy. They have been later extended to other replacement policies [START_REF] Heckmann | The influence of processor architecture on the design and the results of WCET tools[END_REF], data and unified caches [START_REF] White | Timing analysis for data and wrap-around fill caches[END_REF], [START_REF] Chattopadhyay | Unified cache modeling for WCET analysis and layout optimizations[END_REF], and caches hierarchies [START_REF] Hardy | WCET analysis of instruction cache hierarchies[END_REF]. Such techniques were shown to be very tight, especially when considering instruction caches only. However, the need for fixed-point computation makes these techniques prohibitively expensive for large code sizes. In contrast, our method, while providing slightly less tight results, is shown to have a much lower time complexity.

A similar effort, to reduce the cost of cache analyses, was initiated in [START_REF] Patil | Compositional static instruction cache simulation[END_REF]. Separate analyses of each function are performed using static cache simulation. Afterward, the results are composed in a final step. Instead, we completely abstract ourselves from fixed-point based cache analyses.

III. ASSUMPTIONS AND NOTATIONS

The static cache analysis presented in this paper focusses on instruction caches. An architecture with one level of cache (L1) is assumed. The cache is set-associative, with an associativity degree of k, and implements the LRU (Least Recently Used) replacement policy. The address space of every program is logically partitioned into a set of memory blocks. Upon a miss in the cache, a memory block is loaded as a whole in the cache line of equal size. Each instruction is located in a single memory block.

The outcome of the proposed static cache analysis is very similar to the one of existing static cache analysis methods. For each reference within a memory block, a Cache Hit/Miss Classification (CHMC) defines its worst-case behavior:

• always-hit (AH ): the reference will always result in a cache hit, • first-miss (FM ℓ ): the reference could not be classified as hit in the first time it occurs inside loop ℓ, but will result in cache hits in later iterations, • not-classified (NC ): in all other cases. The proposed analysis is contextual: it categorizes each reference in each of its call contexts. The call context of a function f is defined as a path in the program call graph (assumed to be acyclic), and is represented as a list of call basic blocks traversed from the program entry function up to f . The program entry point has a single, empty call context. CTX (i) denotes the set of call contexts attached to instruction i.

The analyses focus on natural loops [START_REF] Muchnick | Advanced compiler design and implementation[END_REF], loops possessing a single entry block, their entry basic block, and possibly several back-edges leading back to the aforementioned entry block.

IV. FAST AND SCALABLE INSTRUCTION CACHE

ANALYSIS

The proposed analysis, illustrated in Figure 1, is divided and introduced in two parts. The first one, described in paragraph IV-A and called basic instruction cache analysis, focuses on the easy-to-capture reference patterns: references appearing inside each basic block in isolation; references in loops, by determining if the memory blocks accessed inside the loop can be ensured to stay in the cache during its execution.

Then, two kinds of additional (and optional) extensions are proposed, in paragraph IV-B, to tighten the estimates of the basic analysis, depending on the user-selected tradeoff between analysis precision and execution time. The first extension focuses on each memory block accessed by different basic blocks. The second one targets memory blocks reused between successive calls to the same function. The two extensions can be applied separately or jointly, and in the latter case, with no imposed order of application. While our analysis produces a cache classification for every call context, the analysis is applied on each function separately, to avoid heavyweight interprocedural analysis.

A. Basic instruction cache analysis

In this section, we define a basic instruction cache analysis divided in two different parts. The first one, called intra basic block analysis, captures a part of the spatial locality, whereas the second one, called loop analysis, captures a part of the temporal locality due to program loops.

Intra basic block analysis

Caches are very efficient at capturing spatial locality in programs by loading information at the memory block granularity. A memory block contains several consecutive instructions. Upon a cache miss, the whole memory block containing the missing instruction is loaded into the cache, and consecutive accesses to the same memory block result in cache hits, as illustrated in Figure 2. In the figure, when instruction i 1 's fetch results in a cache miss, the whole memory block will be loaded from the memory, and thus consecutive instructions i 2 and i 3 will hit. According to this behavior, a very simple intra basic block analysis is defined as follows. Inside a basic block, instructions are executed sequentially in their address order, resulting in consecutive accesses to the same memory block. Thus, every access to a memory block in a basic block can be classified as AH, except the one with the lower address. Formally:

∀i ∈ Instr , ∀ctx ∈ CTX (i), CHMC (i, ctx) = AH if ¬ first access (i, mb(i), bb(i)) ( 1
)
where Instr is the set of all program's instructions, bb(i) is the basic block which contains i, mb(i) is the memory block containing i and first access is a boolean function returning true if i is the first access to mb(i) inside bb(i).

This analysis is applied to each basic block of the program in all call contexts. Only first accesses to memory blocks inside basic blocks remain not-classified (N C) after this analysis step. Such accesses are called i mb f irst hereafter, and the set of all such accesses is denoted IMBFirst . The other analyses defined in this paper focus only on accesses in the IMBFirst set.

Loop analysis

During the execution of a loop, a memory block, denoted mb, accessed inside the loop may be kept in the cache after it has been first loaded. More precisely, assuming a LRU replacement policy, block mb is guaranteed to stay in the cache if the number of other distinct memory blocks possibly accessed within the loop and mapped to the same cache set as mb, is strictly lower than the cache associativity degree. In that case, the access to this memory block can result in a miss or a hit the first time it is referenced, and in hits afterwards (classification FM).

After an identification of loops using standard compiler techniques [START_REF] Muchnick | Advanced compiler design and implementation[END_REF], the loop analysis determines for each loop ℓ of the program, all memory blocks possibly accessed inside ℓ, denoted MB loop (ℓ). MB loop (ℓ) is computed in a straightforward manner by determining all memory blocks referenced by the basic blocks executed within the scope of ℓ, denoted BB loop (ℓ). This includes all basic blocks of functions called, directly or indirectly, by basic blocks in the loop, as illustrated in Figure 3. In the figure, when analyzing loop ℓ in function main, MB loop (ℓ) contains the memory blocks of all basic blocks shown in the figure, due to the direct call to foo and the indirect call to bar.

Then, for a given memory block mb mapped to cache set s, if the number of memory blocks of MB loop (ℓ), mapped to s and different from mb, is strictly lower than the cache's associativity degree, i mb f irst accesses to mb are classified as FM. Formally:

∀i ∈ IMBFirst, ∀ctx ∈ CTX (i), ∀ℓ ∈ Loop, CHMC (i, ctx) = FM ℓ if bb(i) ∈ BB loop (ℓ) ∧ | {c ∈ MB loop (ℓ) : c = mb(i) ∧ set(c) = set(mb(i))} |< k (2)
where Loop is the set of all the program loops, BB loop (ℓ) is the set of all basic blocks referenced when executing ℓ, set is a function returning the cache set where a memory block is mapped and k is the cache associativity degree.

B. Extensions of the instruction cache analysis

The basic instruction cache analysis presented in IV-A is local (to a basic block for the intra basic block analysis, to a loop for the loop analysis). Consequently, it does not detect that a memory block accessed by a basic block might have been loaded by its predecessor(s) in the control flow graph. Similarly, the basic analysis does not detect that a memory block referenced in a function might have been loaded in the cache by previous calls to the same function. The analysis steps defined below detect these two situations, and are presented respectively in paragraphs B.1 and B.2.

Let us first introduce two conditions used in these analyses: the Not evicted condition and the Dominance condition.

• Not evicted condition.

This condition is introduced to capture the fact that a memory block mb cannot be evicted between two of its accesses, once it has been set to the most recently used position by the first one. Thanks to the properties of the LRU replacement policy, mb cannot be evicted if less than k different memory blocks, in conflict with mb, are accessed between the two accesses, with k the cache associativity. If BB denotes the set of memory blocks possibly referenced between the two accesses to mb, a boolean function not evicted (mb, BB ), that returns true if mb is ensured not to be evicted between the two accesses, and false otherwise, is defined as follows:

not evicted (mb, BB ) = | {c ∈ Blocks(BB) : c = mb ∧ set(c) = set(mb)} |< k (3)
In the formula, set(mb) is a function which returns the cache set where mb is mapped in the cache, and Blocks(BB ) is a function returning the sets of memory blocks accessed by each basic block of BB . The not evicted condition is used in the analysis extensions for different parameters BB , whose computation are detailed in due course.

• Dominance condition.

This condition is introduced to ensure that given two basic blocks, any path reaching the second basic block traverses the first one before. This allows us to ensure that a memory block accessed by the dominated basic block is set to the most recently used position into the cache set by the dominating basic block.

The verification of the dominance condition relies on the computation of dominators [START_REF] Muchnick | Advanced compiler design and implementation[END_REF] within each function.

Inter procedural dominator information is not required. Each analysis defines its own dominance condition detailed in due course.

1) Inter basic blocks analysis (IB)

Our inter basic blocks analysis takes into account three different cases where a given basic block bb accesses a memory block mb previously loaded by other basic block(s): (i) the direct predecessors of bb have loaded mb ; (ii) an indirect predecessor which dominates bb has loaded mb ; (iii) the direct predecessors of bb in a loop have loaded mb in a loop.

It can be noticed that the set of basic blocks which are accessing the same memory block mb is easily determined by scanning all memory blocks accessed by each basic block. Then, only basic blocks accessing the same memory blocks have to be matched against the following patterns.

• (i) Direct predecessors analysis.

Intuitively, when two basic blocks p and bb access the same memory block mb and p is the only direct predecessor of bb in a function, the i mb f irst access to mb in bb will result in a cache hit. Indeed, we are sure that any path reaching bb traverses p before, and sets mb to the most recently used position in the cache set. The analysis focuses on the most widespread case: the memory block in common is accessed by the last instruction (i last ) of p and the first instruction (i f irst ) of bb.

Based on this observation, and generalizing it to basic blocks with multiple predecessors, the analysis classifies as AH, i mb f irst accesses to memory blocks ensured to be in the cache because they have been loaded in the cache by all their direct predecessors in the function. A memory block is guaranteed to be in the cache at bb if two conditions are met: (i) all its direct predecessors have accessed this memory block and (ii) any possible execution path between each of the predecessors and bb cannot evict it. This latter condition is required in case some of the predecessors call a function. Formally:

∀i ∈ IMBFirst, ∀ctx ∈ CTX (i), CHMC (i, ctx) = AH if ∀p ∈ predecessors(bb(i)), mb(i last (p)) = mb(i) ∧ i f irst (bb(i)) = i ∧ caller(p) ⇒ not evicted (mb(i), BB callee (callee(p)) (4)
where i last (p) is the last instruction of p, i f irst (b) is the first instruction of b, caller(p) a boolean function which returns true if p calls a function, callee(p) returns the function called by p and BB callee (f ) the set of all basic blocks possibly reached by function f . BB callee (f ) contains all the basic blocks of function f and transitively called functions, as illustrated in Figure 4. In the figure, BB callee (bar) contains all basic blocks of function bar, called by p 2 in function main.

• (ii) Indirect predecessor analysis

In some cases, an indirect predecessor p of a basic block bb, with p dominating bb, can access the same memory block mb as bb. As exemplified in Figure 5, this case may occur when the loop's entry basic block is put at the end of the loop in the binary code. In such a case, bb is inside a loop ℓ, p is ℓ, and they are both accessing mb. In this loop layout, the memory block in common is accessed by the last instruction (i last ) of p and the first instruction (i f irst ) of bb. The analysis focuses on this common case. Then, the i mb f irst access to mb by bb can be classified as AH if mb cannot be evicted by accesses inside the loop (determined by the loop analysis presented in Section IV-A).

To verify that mb is not evicted between p and bb, the set of basic blocks possibly executed between p and bb, and noted BB 1 p,bb , has to be determined. As we know that p is a dominator of bb, this is easily done by gathering the basic blocks (except those coming from back-edges of loop ℓ) met during a reverse graph traversal starting from bb (not included) up to p (not included) and then, by adding all basic blocks of functions called directly or indirectly by p or by basic blocks reached during the graph traversal. Formally:

∀i ∈ IMBFirst, ∀ctx ∈ CTX (i), ∀ℓ ∈ Loop if ∃p ∈ dom(bb(i)), bb(i) ∈ BB loop (ℓ) ∧ p ∈ BB loop (ℓ) ∧ mb(i) = mb(i last (p)) ∧ i f irst (bb(i)) = i ∧ not evicted (mb(i), BB 1 p,bb(i) ) then if CHMC (i, ctx) = FM ℓ then CHM C(i, ctx) = AH (5)

• (iii) Direct predecessor(s) of loop entry

The last case of our inter basic blocks analysis also results from the mapping, in the binary code, of the loop's entry basic block after the loop's body. Let us assume, as illustrated in Figure 6, that the loop's entry basic block bb and its direct predecessor p, inside the loop, are accessing the same memory block mb. In such a case, when bb is executed, we are sure that p was traversed before bb inside the loop. Thus, i mb f irst accesses to mb in bb result in hits after the first iteration of the loop (FM classification). In this situation, p is linked to bb by the back-edge of the loop. Again, the memory block in common is accessed by the last instruction (i last ) of p and the first instruction (i f irst ) of bb. Based on this observation, and generalizing it to loops with multiple back-edges, the dominance condition is verified if each in-loop predecessor of the loop's entry basic block is accessing mb. Regarding the not evicted condition, it is ensured in a straightforward manner if none of bb's predecessors calls a function, because in that case, no access to other memory blocks can occur between the i last instruction of each predecessor and the i f irst access of bb. Otherwise, the not evicted condition is ensured by considering the basic blocks of functions called directly or indirectly by each caller predecessor. Formally:

∀i ∈ IMBFirst , ∀ctx ∈ CTX (i), if ∃ℓ ∈ Loop, bb(i) = entry(ℓ) ∧ i f irst (bb(i)) = i ∧ ∀p ∈ predecessors(bb(i)) ∩ BB loop (ℓ), mb(i) = mb(i last (p)) ∧ caller(p) ⇒ not evicted (mb(i), BB callee (callee(p)) then if CHMC (i, ctx) = AH CHMC (i, ctx) = FM ℓ (6)

2) Inter call analysis (IC)

Between two executions of a given function, some memory blocks loaded by the first execution may be kept in the cache when the second execution occurs. In such a situation, a memory block mb is kept in the cache between these runs, and can be classified as AH in the second execution, if the following conditions are met:

• (p1): whenever the second instance occurs, the first instance must have taken place before; • (p2): mb must be accessed by the function, each time this function is called; • (p3): mb cannot be evicted between its access in the first instance and its access in the second one. We can notice that properties p1 and p2 form the dominance condition while p3 is the not evicted condition. The concordance of two instances of a function to properties p1, p2 and p3 is detailed below.

Each instance of a function is represented by its call context, ctx1 and ctx2 respectively for the first and the second instance of the studied function. The call contexts are denoted [c 1 , c 2 , ..., c n ] for ctx1 and [c ′ 1 , c ′ 2 , ..., c ′ m ] for ctx2, and contain the list of call basic blocks in the acyclic call graph, traversed from the program entry function up to the studied instance.

As ctx1 and ctx2 are different call contexts, their paths share at least a common prefix, possibly empty. Hence, ctx1 (resp. ctx2) can be expressed as

[c 1 , c 2 , ..., c d-1 , c d , ..., c n ] (resp. [c 1 , c 2 , ..., c d-1 , c ′ d , ..., c ′ m ]) such that c d = c ′ d ; (c d , c ′ d )
c d ∈ dom(c ′ d ) (7) ∧ ∀i > d ∧ i ≤ n, ∀bb ∈ exit(function(ci)), ci ∈ dom(bb) (8)
where exit(f ) returns all the exit basic blocs of function f and function(c i ) returns the function containing basic bloc c i .

Intuitively, c d and c ′ d belong to the same function since they have been reached by going through the same call basic blocks path. c d is ensured to be traversed before c ′ d according to the dominator relation (equation 7). Then according to equation 8, each call basic block of ctx1 located after c d is ensured to be traversed. Each of them dominates its function's exit basic blocks.

Property p2: Given a memory block mb, mb is ensured to be accessed each time function f is called if:

∃i ∈ Instr , function(bb(i)) = f ∧ mb(i) = mb ∧ ∀bb ∈ exit (function(bb(i))), bb(i) ∈ dom(bb) (9)
Intuitively, we check that mb is accessed by at least one basic block which dominates function f exit basic blocks. This implies the traversal of this basic block each time f is called.

Property p3: Given a memory block mb, mb cannot be evicted between its access in ctx1 and its access in ctx2 of a same function f if:

not evicted (mb, BB ctx1,ctx2) (10)
where BB ctx1,ctx2 is the set of basic blocks possibly reached between the call contexts ctx1 and ctx2 of f . The computation of BB ctx1,ctx2 , is achieved in three steps. The first step gathers all the basic blocks contained in f . Then, based on the divergence point (c d , c ′ d ) of ctx1 and ctx2, and property p1, basic blocks met in the following graph traversals are gathered: The last step adds to BB ctx1,ctx2 , all basic blocks of functions called directly or indirectly by a basic block determined by the two previous steps.

Function f may have more than two call contexts. In such a situation, we only have to consider in the analysis the immediate dominating context in accordance with the dominance condition. An immediate dominating context is derived from its set of dominators similarly to the derivation of a basic block's immediate dominator from its set of dominators. Finally, for each function f , and each of its call contexts ctx, all i mb f irst accesses to memory blocks which verify the not evicted condition and the dominance condition between the immediate dominating context of ctx and ctx are classified AH.

V. EXPERIMENTAL RESULTS

In this section, we compare the tightness and cost (in terms of computation time and memory consumption) of the proposed method with those of a state-of-the-art instruction cache analysis [START_REF] Theiling | Fast and precise WCET prediction by separated cache and path analyses[END_REF] used as a baseline. We first describe the experimental conditions ( § V-A) and then we give and analyze results of our experiments ( § V-B).

A. Experimental setup

Analyzed codes: The experiments were conducted on seven benchmarks from the Mälardalen WCET benchmark suite 1 , and on the Debie software, developed by the European Space Agency. Debie monitors space debris and micrometeoroids in situ by detecting impacts using mechanical and electrical sensors. We have analyzed two tasks of Debie, as well as the test suite of Debie, which is to the best of our knowledge the largest real-time application publicly available. Table I summaries the applications' characteristics. The analyzed programs include both small loop intensive programs (jfdctint, minver, crc, qurt, fft, adpcm) and control programs, containing less loops and more complex flows of control (statemate, Debie tasks and Debie test suite).

Cache analysis and WCET estimation:

The experiments were conducted on MIPS R2000/R3000 binary code compiled with gcc 4.1 with no optimization. The default linker memory layout is used: functions are represented in memory sequentially, and no alignment directive is used. The WCETs of tasks are computed by the Heptane timing analyzer [START_REF] Colin | A modular and retargetable framework for tree-based WCET analysis[END_REF],

1 http://www.mrtc.mdh.se/projects/wcet/benchmarks.html more precisely its Implicit Path Enumeration Technique (IPET), solving linear equations generated from the program control flow graph, as originally proposed in [START_REF] Malik | Performance analysis of embedded software using implicit path enumeration[END_REF].

Only the outermost loops are analyzed in our experiments to underline the impact on the tightness of referenced patterns not captured by our loop analysis. Moreover, to separate the effect of caches from those of other parts of the processor micro-architecture, the WCET estimation only accounts for the contribution of caches to the WCET. The effects of other architectural features are not considered. In particular, we do not take into account timing anomalies caused by interactions between caches and pipelines, as defined in [START_REF] Lundqvist | Timing anomalies in dynamically scheduled microprocessors[END_REF]. The cache classification NC is thus assumed to have the same worst-case behavior as a miss during the WCET computation.

The state-of-the-art baseline instruction cache analysis used in our experiments is an implementation of the cache analysis of [START_REF] Theiling | Fast and precise WCET prediction by separated cache and path analyses[END_REF], based on abstract interpretation. For each of its composing fixed-point analyses, respectively Must and Persistence, abstract cache states (ACSs) are computed upon each basic block entry and exit. ACSs for the Must analysis model memory blocks known to be present in cache, whereas for the Persistence analysis they model memory blocks persistent in loops. Until a fixpoint on ACSs' contents is reached, a basic block's entry ACS is computed from its predecessors' exit ACSs; similarly, an exit ACS is computed from said basic block's input and its sequence of memory accesses. The cache analysis starts with an empty cache state, which is safe in the assumed context (absence of timing anomalies, use of the LRU replacement policy). The safety issue in the original persistence analysis of [START_REF] Theiling | Fast and precise WCET prediction by separated cache and path analyses[END_REF], recently identified and corrected in [START_REF] Huynh | Scope-aware data cache analysis for WCET estimation[END_REF] is solved using the method given in [START_REF] Huynh | Scope-aware data cache analysis for WCET estimation[END_REF].

The three most interesting configurations of our cache analysis framework are compared to each others and to the state-of-the-art fixed-point based analysis:

• the basic analysis only, denoted BA;

• the basic analysis with the inter basic blocks analysis, denoted BA+IB; • the analysis with the inter basic blocks analysis and inter call analysis, denoted BA+IB+IC. Cache configurations: The experiments were conducted using two different cache sizes. For the Mälardalen WCET benchmarks and for the two tasks of Debie, a 1KB cache is used. A 8KB cache is used for the test suite of Debie. In both configurations, a cache associativity of 4 and a 32B cache line size are used.

Metrics: To evaluate the precision of our approach, we compare the predicted hit ratio obtained using the different configurations of our analysis framework with the hit ratio predicted by the state-of-the-art analysis. The predicted hit ratio is the percentage of instructions classified as hits, computed on the worst-case execution path identified by the IPET WCET computation. Regarding the evaluation of the cost of our analysis, the computation time evaluation has been realized on an Intel Core 2 Duo 3 GHz with 8 GB of RAM and the memory consumption is measured with Valgrind.

B. Experimental results

Analysis tightness: Table II compares the tightness of our analysis with the one of a state-of-the-art analysis, for the three studied configurations of our analysis framework (BA, BA+IB, BA+IB+IC).

As shown in Table II, the hit ratios predicted by the stateof-the-art analysis and those predicted by our analysis with all extensions enabled are always close to each other. The difference in percentage between the baseline state-of-theart method and the BA+IB+IC analysis is computed by the following formula: hit ratio baseline -hit ratioBA+IB+IC hit ratio baseline * 100.

In average for all analyzed codes, the difference is 0.53%.

For most analyzed small codes, the two predicted hit ratios are equal (jfdctint, minver, crc, qurt, acquisition task, and tc execution task). In some cases, our analysis provides slightly less tight results than the state-of-the-art method (fft, adpcm, statemate, with a difference in predicted ratio of respectively 0.31%, 0.05% and 0.55%). A careful examination of the structure of these codes allowed us to identify two different scenarios implicitly captured by the state-ofthe-art analysis only: (i) memory blocks shared between different functions due to function frontiers not aligned on memory block boundaries, which are not captured by our analysis; (ii) memory blocks kept between successive calls to the same function in different loop iterations, which is not investigated by our analysis. Those two cases constitute interesting, tightness-improving, future extensions of our method. The maximum difference between the state-of-theart method and BA+IB+IC is 4.40% for the largest analyzed code (Debie test suite); due to the code complexity, further experiments would be needed to pinpoint the exact sources of pessimism of our analysis.

As one would expect, for all analyzed codes, the hit ratio of the fixed-point free analysis with all extensions (BA+IB+IC) is always strictly higher than the hit ratio obtained with the basic analysis (BA). Concerning the individual impact of the analysis extensions on tightness, we observe that the IB extension benefits to all of analyzed codes. The IC extension benefits to some analyzed codes only (qurt, fft, acquisition task, tc execution task, Debie test suite); jfdctint, minver, crc, adpcm and statemate are unaffected by the IC extension, because their code does not reuse memory blocks between function calls. A good knowledge of the analyzed application, or an initial analysis run during an iterative process, may guide the activation of the IC extension during WCET analyses.

Analysis cost (time and memory): Table III compares the cost of the proposed cache analysis, with all extensions enabled, with the cost of the baseline analysis. The left part of the table gives the analysis run-time, while the right part gives the analysis memory consumption. The costs are evaluated for the cache analysis only, and exclude costs for loading the control flow graph in memory and computing the WCET from the cache classification.

Regarding the analysis execution time, the proposed analysis is always faster than the baseline analysis. In average for all analyzed codes, the proposed analysis is more than 5 times faster than the state-of-the-art. The most interesting feature is that the larger speedup is obtained with the largest and most complex analyzed code (Debie test suite), which is the objective of our analysis. For that code, the proposed analysis is close to 30 times faster than the baseline.

Similarly, the memory consumption of the proposed analysis is much lower than the one of the baseline analysis (3.9 times lower in average for all the benchmarks). And interestingly again, the code with the best improvement in terms of memory consumption is the largest and most complex code we have analyzed (Debie test suite), which required 12 times less memory to be analyzed using our method than using the baseline analysis.

Table IV gives a detailed view of the analysis cost for the Debie test suite, for the different studied analysis configurations (BA, BA+IB, BA+IB+IC, state-of-the-art analysis). The basic analysis BA detects a part of the application locality captured by the cache. Nevertheless, the improvements in terms of precision attained by the use of the IB extension clearly outweighs the negligible extra computation time required to perform this analysis. This is an argument in favor of the inclusion of the IB extension in the default analysis framework. Higher precision can still be reached by applying the IC extension at the price of extra computation time, larger than the execution time of IB, but still reasonable.

In summary, depending on the context of WCET estimation (early software development phases or iterative WCEToriented compiler optimizations), and the user knowledge about the application structure, the user can selected the best trade-off between analysis precision and analysis execution time.

VI. CONCLUSION

In this paper, we have introduced a novel fixed-point free instruction cache analysis, to abstract cache analysis from costly fixed-point computation. A basic analysis and two optional extensions have been proposed based on classical graph information commonly used in compilers. This modularity allows for trade-offs between analysis tightness and efficiency, depending on the analysis context such as early software development phases, timing analysis of large pieces of software, or iterative WCET-oriented compiler optimizations.

The experimental results show that the proposed analysis yields results close to the ones obtained using a state-ofthe-art analysis (average loss of precision of 0.5%). Further-more, the analysis with all the extensions selected shows a significant reduction in terms of computation time and memory consumption; in average, the analysis is 5 times faster and consumes 4 times less memory than one state-ofthe-art analysis.

In future research, we plan to improve the tightness of our method by defining additional extensions, as an example to capture reuse between calls to functions performed inside loops. We also plan to extend our analysis to support other replacement policies. We conjecture that this can be achieved by using the minimum life-span bound defined in [START_REF] Reineke | Timing predictability of cache replacement policies[END_REF]. Another direction for future work is the extension of this approach to analyze data caches and cache hierarchies, for uni-core and multi-cores architectures.

Figure 1 .

 1 Figure 1. Analysis framework

Figure 2 .

 2 Figure 2. Intra basic block analysis

Figure 3 .

 3 Figure 3. Loop analysis.

Figure 4 .

 4 Figure 4. Direct predecessors analysis

Figure 5 .

 5 Figure 5. Indirect predecessor analysis

Figure 6 .

 6 Figure 6. Direct predecessor(s) of loop entry

  is the divergence point of ctx1 and ctx2. Property p1: Given two different call contexts, ctx1 and ctx2, represented by their paths [c 1 , c 2 , ..., c d-1 , c d , ..., c n ] and [c 1 , c 2 , ..., c d-1 , c ′ d , ..., c ′ m ], and (c d , c ′ d ) their divergence point. ctx1 is ensured to occur before ctx2 if:

• ( 1 )

 1 from each call basic block reached after c d (not included) in ctx1, down to their respective function exit basic blocks; • (2) from c ′ d up to c d (not included); • (3) from each call basic block reached after c ′ d in ctx2, up to their respective function entry basic block. As illustrated in Figure 7, (1) considers all the basic blocks from the end of the first instance of f to the successors of c d , (2) considers all the basic blocks from the successors of c d to c ′ d , and (3) considers all the basic blocks from c ′ d to the second instance of f .

Figure 7 .

 7 Figure 7. Inter call analysis for the function bar

Table II PREDICTED

 II HIT RATIO OF THE FIXED-POINT FREE ANALYSIS COMPARED TO THE STATE-OF-THE-ART ANALYSIS. CACHE SIZE OF 1KB (8KB FOR DEBIE TEST SUITE).

	Analyzed code	BA	BA+IB	BA+IB+IC	State-of-the-art
						analysis
	jfdctint	0.9073	0.9098	0.9098	0.9098
	minver	0.9583	0.9658	0.9658	0.9658
	crc	0.9991	0.9993	0.9993	0.9993
	qurt	0.9671	0.9755	0.9770	0.9770
	fft	0.8360	0.8739	0.8971	0.8999
	adpcm	0.9701	0.9788	0.9788	0.9793
	statemate	0.7612	0.8294	0.8294	0.8340
	acquisition task	0.9940	0.9958	0.9975	0.9975
	tc execution task	0.9938	0.9958	0.9959	0.9959
	Debie test suite	0.7120	0.8184	0.8758	0.9162
	Analyzed code	Execution time (s)		Memory consumption (MB)
	BA+IB+IC	State-of-the-art	BA+IB+IC State-of-the-art
			analysis			analysis
	jfdctint	0.003	0.004		0.15	0.30
	minver	0.003	0.009		0.24	0.59
	crc	0.003	0.006		0.15	0.44
	qurt	0.003	0.009		0.28	0.81
	fft	0.009	0.022		0.54	1.65
	adpcm	0.016	0.046		0.94	4.91
	statemate	0.018	0.069		1.51	6.60
	acquisition task	0.015	0.022		0.97	2.37
	tc execution task	0.036	0.060		4.08	5.47
	Debie test suite	8.07	239.17		250.70	3088.07

Table III EXECUTION

 III TIME AND MEMORY CONSUMPTION OF THE FIXED-POINT FREE ANALYSIS COMPARED TO THE STATE-OF-THE-ART ANALYSIS. CACHE SIZE OF 1KB (8KB FOR THE DEBIE TEST SUITE).Table IV COMPUTATION TIME AND MEMORY CONSUMPTION FOR THE DEBIE TEST SUITE COMPARED TO THE STATE-OF-THE-ART ANALYSIS. CACHE SIZE OF8KB.

	Metrics	BA	BA+IB	BA+IB+IC	State-of-the-art
					analysis
	Analysis	2.31	2.85	8.07	239.17
	time (s)		(2.31 + 0.54)	(2.31 + 0.54 + 5.22)
	Memory	152.66	171.38	250.70	3088.07
	consumption (MB)

ACKNOWLEDGMENT

The authors would like to thank André Seznec and the anonymous reviewers for their fruitful comments on earlier versions of this paper.


Soukayna Ghandour-Haidar

Laurent Ros
email: laurent.ros@gipsa-lab.grenoble-inp.fr

Jean-Marc Brossier
email: jean-marc.brossier@gipsa-lab.grenoble-inp.fr

On the Use of First-order Autoregressive Modeling for Rayleigh Flat Fading Channel Estimation with Kalman Filter

Keywords: Channel estimation, Autoregressive model, Kalman Filter, Jakes's spectrum, Rayleigh channel, Flat fading, Bayesian Cramér-Rao Bounds BCRB

. This letter gives analytic justification by providing approximated closed-form expressions of the estimation variance for the CM and MAV criteria, and of the optimal AR(1) parameter.

Introduction

The Rayleigh fading channel model with Jakes's Doppler spectrum is the most accepted random model to represent temporal variations (fading) of the equivalent baseband channel complex gain (CG) in wireless communication. However, this model is not always directly convenient for highly important tasks such as channel estimation or equalization. Usually, the autoregressive model AR(p) with some Gaussian assumptions (most often used with order p = 1) is used to approach this channel and facilitate its manipulation. In [START_REF] Baddour | Autoregressive modeling for fading channel simulation[END_REF], the authors demonstrated that the AR model can be considered for the computer simulation of correlated fading channels insisting that low orders are appropriate for narrowband Doppler fading processes. Moreover, [START_REF] Wang | On verifying the first-order Markovian Assumption for a Rayleigh Fading Channel Model[END_REF] showed that a first-order model is enough to capture most of the channel tap dynamics. This approximation has been widely used to track the true Jakes's spectrum channel by a Kalman Filter (KF) in various wireless communication systems ([4, 5, 6, 7, 8, 9, 10, 11], ...). This letter deals with the choice of the AR(1) coefficient. In many papers, e.g., the aforementioned works, the AR(1) coefficient calculation for a given normalized Doppler frequency ( f d T ) is based on the same criterion, called Correlation Matching (CM) in this letter. The CM criterion 1 consists in imposing that the autocorrelation coefficients R α [n] of the approximated AR(1) process α perfectly match the sampled autocorrelation function (which is a Bessel function) of the true CG α for lags n ∈ {-1, 0, 1} [START_REF] Baddour | Autoregressive modeling for fading channel simulation[END_REF] (see also [START_REF] Komninakis | Multi-Input Multi-Output fading channel tracking and equalization using Kalman estimation[END_REF], Fig. 1). But in some situations the results were disappointing (in terms of the distance from the Bayesian Cramér-Rao lower Bound (BCRB)), as pointed out recently in [START_REF] Abeida | An EM algorithm for path delay and complex gain estimation of a slowly varying fading channel for CPM signals[END_REF] (Fig. 3) and in [START_REF] Ros | Paths complex gain tracking algorithms for OFDM receiver in slowly-varying channels[END_REF] (Fig. 2) for the special case of a slow fading scenario where the channel is theoretically easier to estimate. Note also that some authors propose adding a positive ε to the zeroth autocorrelation lag, which modifies slightly the CM constraint to better approximate the original process in some sense [START_REF] Baddour | Autoregressive modeling for fading channel simulation[END_REF]. In the recent work [START_REF] Barbieri | On the ARMA Approximation for Frequency Channels Described by the Clarke Model with Applications to Kalmanbased Receivers[END_REF], the search for the AR(1) parameter is based on another method, consisting in the minimization of the steady-state estimation error variance, called the Minimum Asymptotic Variance (MAV) criterion in this letter. This method seems very effective compared to the CM method for a Clarke (i.e., Rayleigh-Jakes) Model channel, but the study in [START_REF] Barbieri | On the ARMA Approximation for Frequency Channels Described by the Clarke Model with Applications to Kalmanbased Receivers[END_REF] does not provide analytic results nor closed-form expressions about the choice of the AR(1) parameter (nor the associated Mean Square Error (MSE)) for a given channel state ( f d T , SNR). In this letter, we first provide a general theoretical frequency-domain analysis of the estimation error (in terms of static and dynamic contributions) that explains analytically the previous disappointing results with the CM-criterion based method. Moreover, the analytic approach is used to calculate the optimal AR(1) coefficient under the MAV criterion, without the CM constraint. For a given Doppler and SNR scenario, we give approximate closedform expressions of the optimal coefficient and of the corresponding MSE.

Mathematical Model

We consider the estimation of a flat Rayleigh fading channel. The discrete-time observation is2 :

y k = α k + n k ( 1
)
where k is the (symbol) time index, n k is a zero-mean additive white circular complex Gaussian noise with variance σ 2 n , and α k is a zero-mean circular complex Gaussian CG with variance σ 2 α = 1. The normalized Doppler frequency of this channel is f d T , where T is the symbol period. A Jakes's Doppler spectrum is assumed for this channel:

Γ α ( f ) =    σ 2 α π f d √ 1- ( f f d ) 2 if | f | < f d 0 if | f | > f d (2)
The autocorrelation matrix R k α of a block of k consecutive CG is then defined by its [i, j]th entry:

[ R k α ] i, j = R α [i -j] = σ 2 α J 0 (2π f d T.(i -j)) 1 ≤ i, j ≤ k (3)
where J 0 is the zeroth-order Bessel function of the first kind.

The time-varying CG α k is approached by a first-order autoregressive (AR(1)) model αk :

αk = a. αk-1 + e k ( 4
)
where e k is a white circular complex Gaussian noise with variance σ 2 e = (1a 2 )σ 2 α . The observation equation ( 4) is then approximated by:

y k = αk + n k (5)
The AR(1) coefficient verifies

a = R α [1] R α [0] (6)
If the correlation matching (CM) criterion [START_REF] Baddour | Autoregressive modeling for fading channel simulation[END_REF][START_REF] Komninakis | Multi-Input Multi-Output fading channel tracking and equalization using Kalman estimation[END_REF][START_REF] Liu | Space-time coding and Kalman Filtering for time selective fading channels[END_REF][START_REF] Mcguire | Low-order Kalman filters for channel estimation[END_REF][START_REF] Al-Naffouri | An EM-based forward-backward Kalman for the estimation of time-variant channels in OFDM[END_REF] is imposed, the two first autocorrelation coefficients of the AR(1) process α coincide with those of the exact CG

α (i.e., R α [0] = R α [0], and R α [1] = R α [1]
), and using equation ( 6) and (3) the AR(1) coefficient a noted a CM becomes

a CM = J 0 (2π f d T ) (7)

MSE Analysis, Optimization, and Simulation Results

Given the model and observation equations (4 and 5), we use a KF to get an on-line unbiased estimate, αk , of the true α k . We compare the variance of the error ε k = α k -αk to the on-line BCRB (see [START_REF] Hijazi | Bayesian Cramér-Rao Bounds for complex gain parameters estimation of slowly varying Rayleigh channel in OFDM systems[END_REF] for details about the computation of this BCRB):

E { |ε k | 2 } ≥ BCRB(k) with BCRB(k) = [ J -1 ] k,k
3 the last element of the inverse of the Bayesian Information

Matrix J : J = 1 σ 2 n I k + [R k α ] -1
and I k is the k × k identity matrix. The KF is given by [START_REF] Chui | Kalman Filtering with real-time applications[END_REF] (see also p. 436 in [START_REF] Kay | Fundamentals of Statistical Signal Processing -Estimation Theory[END_REF]):

Prediction αk|k-1 = a αk-1 (i) Prediction MSE P k|k-1 = a 2 P k-1 + σ 2 e (ii) Kalman Gain K k = P k|k-1 σ 2 n +P k|k-1 (iii) Correction αk = αk|k-1 + K k (y k -αk|k-1 ) (iv) MSE P k = (1 -K k )P k|k-1 (v)
Replacing (ii) in (iii), (v) in (ii), and (i) in (iv), we obtain the simplified system:

K k = [a 2 P k-1 + σ 2 e ]/[a 2 P k-1 + σ 2 e + σ 2 n ] (8)
P k = (1 -K k )(a 2 P k-1 + σ 2 e ) (9) αk = a αk-1 + K k (y k -a αk-1 ) (10)
where K k is the Kalman gain at iteration k and P k is the estimation error variance. Since the linear system (1) and ( 4) is observable and controllable, an asymptotic regime for which K k = K k+1 = K ∞ and P k = P k+1 = P ∞ is quickly reached ( [START_REF] Chui | Kalman Filtering with real-time applications[END_REF]).

Based on equations ( 8) and ( 9), and given

a 2 = 1 -σ 2 e σ 2 α
, we deduce

P ∞ = σ 2 n (a 2 -1) -σ 2 e + √ ∆ 2a 2 K ∞ = a 2 P ∞ + σ 2 e a 2 P ∞ + σ 2 e + σ 2 n ( 11
) with ∆ = (σ 2 n +σ 2 e -a 2 σ 2 n ) 2 +4a 2 σ 2 n σ 2 e .
From now on, we assume that the noise model variance satisfies σ 2 e << σ 2 n . Note that this assumption means that K ∞ << 1, which seems reasonable for low normalized Doppler (i.e., f d T ≤ 10 -2 ) and the usual range of SNR (i.e., 0.01 ≤ σ 2 n σ 2 α ≤ 1). We then have the following approximations (see Appendix A): √ ∆ ≈ 2σ e σ n and a 2 P ∞ ≈ σ e σ n [START_REF] Hijazi | Bayesian Cramér-Rao Bounds for complex gain parameters estimation of slowly varying Rayleigh channel in OFDM systems[END_REF] leading to the following approximate expression of the steady-state Kalman gain:

K ∞ ≈ σ e σ n ( 13
)
On the other hand, we see that in the asymptotic regime (or tracking mode), equation [START_REF] Abeida | An EM algorithm for path delay and complex gain estimation of a slowly varying fading channel for CPM signals[END_REF] of the KF reduces (using ( 1)) to a time-invariant linear filter, expressed in the Z-Transform domain by:

α(z) = az -1 α(z) + K ∞ (α(z) + N(z) -az -1 α(z))
Then α(z) = T (z)(α(z) + N(z)) and the error is

ε(z) = α(z) -α(z) = (1 -T (z))α(z) -T (z)N(z) (14)
where T (z) is the transfer function of the steady-state KF:

T (z) = K ∞ 1 -a(1 -K ∞ )z -1 (15)
Using ( 14), the asymptotic MSE is divided into two parts:

MSE = MSE1 + MSE2 ( 16
)
• MSE1 is the dynamic error variance, due to the CG variations α k filtered by the high pass filter 1 -T (z):

MSE1 def = ∫ + 1 2T -1 2T |1 -T (e j2π f T )| 2 Γ α ( f )d f (17)
Under the assumption 1 2

σ 2 e σ 2 α << 2π f d T < 2π f c T << 1,
which is reasonable for a low-Doppler scenario, we obtain (see Appendix B):

MSE1 ≈ ∫ + f d -f d 2π f T 2π f c T 2 Γ α ( f )d f = σ 2 α (2π f d T ) 2 2(2π f c T ) 2 (18)
where f c is the cut-off frequency of the low-pass filter T(z), defined such that 2π

f c T = 1-a(1-K ∞ ) a(1-K ∞ ) . Still assuming that σ 2 e << σ 2 n , it results that 2π f c T ≈ σ n -a(σ n -σ e )
a(σ n -σ e ) ≈ σ e σ n , which means that the (normalized) cut-off pulsation is approximately equal to the steady-state Kalman gain. The approximate closed-form expression of MSE1 is then

MSE1 ≈ σ 2 α (2π f d T ) 2 σ 2 n 2σ 2 e = ( 2π
f d T ) 2 σ 2 n 2(1 -a 2 ) ( 19
)
• MSE2 is the static error variance, due to the additive noise n k filtered by the low pass filter -T (z):

MSE2 def = σ 2 n T ∫ + 1 2T -1 2T |T (e j2π f T )| 2 d f = σ 2 n K 2 ∞ |1 -a 2 (1 -K ∞ ) 2 | (20)
Inserting ( 13) in (20), using again

a 2 = 1 -σ 2 e σ 2 α
and the previous assumption σ 2 e << σ 2 n , it follows that

MSE2 ≈ σ 2 n . σ e 2σ n = √ (1 -a 2 )σ n σ α 2 (21)
These approximate closed-form expressions ( 19) and (21) will be useful in the perspective of an optimization (MAV criterion). But first we apply them to the case of the CM criterion in order to explain the performance of CM-based methods.

By imposing the CM criterion [START_REF] Chen | Kalman-filter channel estimator for OFDM systems in time and frequency-selective fading environment[END_REF], and using

[J 0 (ε)] 2 ≈ 1 -ε 2 2 for ε << 1, we have (1 -a 2 ) ≈ 2(π f d T ) 2 and then MSE1 CM ≈ σ 2 n (22) MSE2 CM ≈ π f d T σ n σ α √ 2 (23)
Note that these approximate closed-form expressions assumed that f d T < f c T (for a valid computation of MSE1 as ( 18) and ( 19)) and then are valid only for σ 2 n σ 2 α < 1 2 , since with the CM criterion we get σ e ≈ √ 2π f d T , and then

f c T ≈ K ∞ ≈ σ e σ n ≈ √ 2π f d T
σ n . However, for the MAV criterion, we can find the coefficient of the AR(1) model by minimizing the MSE in ( 16) with respect to a, assuming ( 19) and ( 21):

a * = 1 -4 3 √ (π f d T ) 4 σ 2 n σ 2 α (24)
and the theoretical corresponding MSE is

MSE * = 3 2 ( π f d T σ 2 n σ α ) 2 3 (25)
In Fig. B.1, we compare the MSE obtained with CM and MAV criteria. We also plot the BCRB (for k = 2000) as reference.

• With the CM criterion, the MSE is approximately constant with respect to the Doppler frequency and at low Doppler, the error seems to be far from the BCRB. This agrees with equations ( 22) and ( 23 • With the MAV criterion, MSE * seems to be approximately the same as the MSE computed by Monte Carlo simulation for all the usual range of SNR between 0 and 20 dB, so we validate the closed-form expression (25).

• The most important observation is that with the optimization (MAV criterion), we are closer to the bounds than with the CM criterion, especially for low frequencies ( f d T = 10 -4 ) and low SNRs (this corroborates the results of [START_REF] Barbieri | On the ARMA Approximation for Frequency Channels Described by the Clarke Model with Applications to Kalmanbased Receivers[END_REF]).

In Fig. B.2, we complete our comparison with the AR(2) KF (under the CM criterion) of [START_REF] Komninakis | Multi-Input Multi-Output fading channel tracking and equalization using Kalman estimation[END_REF], and the AR(1) KF of [START_REF] Enescu | On the estimation of state matrix and noise statistics in state space models[END_REF] where the AR(1) coefficient a is estimated on-line from observations under the assumption that the channel is a true AR(1) process (instead of a Jakes process). The performance of the first algorithm is close to that of the AR(1) KF under the CM criterion (that corroborates the results in [START_REF] Barbieri | On the ARMA Approximation for Frequency Channels Described by the Clarke Model with Applications to Kalmanbased Receivers[END_REF] and [START_REF] Komninakis | Multi-Input Multi-Output fading channel tracking and equalization using Kalman estimation[END_REF]) while the second one is poorer, and then the AR(1) KF optimized by the MAV criterion performs better than both of them. Although the knowledge of SNR = σ 2 α σ 2 n is required to design the KF in [START_REF] Al-Naffouri | An EM-based forward-backward Kalman for the estimation of time-variant channels in OFDM[END_REF], its value is also needed for MAV-based optimization while it is not for the CM criterion. For the sake of fairness, we depict the sensitivity to imperfect SNR knowledge: Fig. B.2 plots the MSE versus the true SNR for wrong values of a (obtained from (24) using arbitrary values of the SNR, fixed to SNR ′ = 10, 15, or 20 dB), it shows that the sensitivity to the SNR knowledge is quite weak and that MAV-based optimization even with an SNR mismatch outperforms other CM-based methods.

Conclusion

This paper addresses the problem of estimating a Rayleigh channel using a first order AR model. An analytic study clearly shows that the most widely used choice for the AR(1) pole estimation (the CM criterion) is not accurate for low SNR and low Doppler f d T . Therefore, switching to an estimation error variance criterion as already proposed in [START_REF] Barbieri | On the ARMA Approximation for Frequency Channels Described by the Clarke Model with Applications to Kalmanbased Receivers[END_REF], we carry out the optimization of the AR(1) model and the calculation of its performance. We provide an approximate expression of the MSE for the CM and MAV criteria first, and of the AR(1) (MAV) parameter for a given SNR and Doppler scenario. It is demonstrated that the MSE of the AR(1) KF (MAV) is proportional to the (2/3) power of the product ( f d T × σ 2 n ), where σ 2 n is the observation noise variance.

Appendix A. Approximations of ∆ in [START_REF] Hijazi | Bayesian Cramér-Rao Bounds for complex gain parameters estimation of slowly varying Rayleigh channel in OFDM systems[END_REF] With the assumption σ 2 e << σ 2 n , it is straightforward to see that

∆ = (σ 2 n + σ 2 e -a 2 σ 2 n ) 2 + 4a 2 σ 2 n σ 2 e = σ 2 n [ σ 2 n σ 4 e σ 4 α + 4σ 2 e -4 σ 4 e σ 2 α ] ≈ 4σ 2 n σ 2 e
Note that for an AR(1) model with CM criterion, we have where

σ 2 e σ 2 α = 1-a 2 CM = 1-[J 0 (2π f d T )] 2 = 1 -[1 -(π f d T ) 2 + (π f d T ) 4 2 -....] 2 ≈ 2(π f d T ) 2 ≤ 2π
A = K ∞ a(1-K ∞ ) , 2π f c T = 1-a(1-K ∞ ) a(1-K ∞ )
, and 2π f V LF T = 2π f c T -A = 1-a a . Note that for proper tracking, the cut-off frequency of the steady-state filter, f c , should be chosen higher than the Doppler frequency f d , but much lower than the sampling rate 1/T . The two conditions are quite compatible in the low Doppler scenario, and are summarized by

f d T < f c T << 1.
Note also that 2π f V LF T = 1 a -1, and given

1 a = [ 1 -σ 2 e σ 2 α ] -1 2 ≈ 1 + 1 2 σ 2 e σ 2 α
, we have

2π f V LF T ≈ 1 2 σ 2 e σ 2 α
. And according to Appendix A, we have for the AR(1) model

with CM criterion, σ 2 e σ 2 α ≈ 2(π f d T ) 2 << 4π f d T .
From now on, our complete assumptions for the computation of the integral (17) will be:

f V LF T << f d T < f c T << 1 2.
The function 1 -T (e j2π f T ) can be (asymptotically) divided into three parts (as can be checked by the plot of 1 -T (e j2π f T ) 2) KF (CM), and AR(1) KF of [START_REF] Enescu | On the estimation of state matrix and noise statistics in state space models[END_REF], and the effect of a mismatch on the knowledge of the SNR for the AR(1) KF (MAV) with an AR(1) coefficient a computed by ( 24) with (false) fixed values of SNR, SNR ′ = 10, 15 and 20 dB (solid lines), and (true) SNR (circles).

2 in Fig. B.3): (a) for 0 < f T < f V LF T , a constant part approximately equal to 2π f V LF T 2π f c T (b) for f V LF T < f T < f c T , a linear part approximately equal to 2π f T 2π f c T (c) for f T > f c T , a constant part equal to 1, not considered in the integral (17), since Γ α ( f ) = 0 for | f | > f d , with f d T < f c T . Then MSE1 ≈ 2 + f V LF ∫ 0 [ 2π f V LF T 2π f c T ] 2 Γ α ( f )d f + 2 + f d ∫ + f V LF [ 2π f T 2π f c T ] 2 Γ α ( f )d f (B.2) = 2I 1 [0; + f V LF ] + 2I 2 [+ f V LF ; f d ] (B.3) but I 2 [0; + f V LF ] < I 1 [0; + f V LF ] since f < f V LF in this computation, and I 2 [+ f V LF ; + f d ] < I 2 [0; + f d ], then 2I 2 [0; + f d ] < MSE1 < 2I 1 [0; + f V LF ] + 2I 2 [0; + f d ] but given the form of Γ α , it results that I 1 [0; + f V LF ] < f V LF f d I 2 [0; + f d ] << I 2 [0; + f d ] since f V LF << f d . In conclusion, MSE1 ≈ 2I 2 [0; + f d ] = 2 + f d ∫ 0 [ 2π f T 2π f c T ] 2 Γ α ( f )d f = (2π f d T ) 2 2(2π f c T ) 2 σ 2 α .

2 n σ 2 α∈

 22 ) since at low f d T ≤ 10 -2 and medium SNR (i.e., σ [0.01; 0.5]), MSE2 CM is small and the mean square error MSE CM ≈ MSE1 CM ≈ σ 2 n .

2 α≤ 1 )

 21 2 10 -4 (due to the low Doppler assumption, i.e., f d T ≤ 10 -2 ). And due to the additional usual range SNR assumption (i.e., 0 dB ≤ SNR ≤ 20 dB and thus 0.01 ≤ σ 2 n σ , we verify the assumption σ 2 e << σ 2 n .

Figure B. 1 :Figure B. 2 :

 12 Figure B.1: MSE comparison of the AR(1) KF estimator based on the literature CM criterion[START_REF] Chen | Kalman-filter channel estimator for OFDM systems in time and frequency-selective fading environment[END_REF] or on the MAV criterion (24), for different SNR, as functions of f d T .

Figure B. 3 :

 3 Figure B.3: The high-pass filter 1 -T (z) for various f d T , SNR = 20 dB, and σ 2 e σ 2 α = f d T 1000

And for a more general AR(p) model (with more degrees of freedom than the AR(1) model), the coefficients of the AR(p) process are calculated[START_REF] Baddour | Autoregressive modeling for fading channel simulation[END_REF] by still imposing the correlation matching constraint for lags n ∈ {-p, ...., -1, 0, 1, ...., p}, but followed by the resolution of the Yule-Walker equations in order to minimize the prediction error (model noise).

[START_REF] Baddour | Autoregressive modeling for fading channel simulation[END_REF] In general, the observation is rather y k = s k α k + n k where s k denotes the transmitted symbols. Our simplified model (1) assumes the symbols are known (or decided). This model relies also on the assumption that the CG variation can be neglected during one symbol period. It should be noted that even if it is simple, this observation model is the same as the one obtained through multipath frequency selective channel in an OFDM system, for each subcarrier.

The Cramér-Rao inequality is verified at every iteration k but the BCRB converges to a constant value after a certain number of iterations ([START_REF] Hijazi | Bayesian Cramér-Rao Bounds for complex gain parameters estimation of slowly varying Rayleigh channel in OFDM systems[END_REF],[START_REF] Ghandour | Estimation de canal à évanouissements lents pour les communications radio-mobiles[END_REF]), (i.e., BCRB ≈ BCRB(k) for a sufficiently large value of k).

Appendix B. Approximation of T (z) and Approximation of MSE1 in (18)

The approximations made to obtain (18) are:

1. The steady-state KF, T(z), can be approximated in the frequency domain using z -1 = e -j2π f T ≈ 1 -j2π f T for low frequencies compared to the sampling-rate, i.e., for f << 1/T , by: