F Lima Kastensmidt

L Sterpone
email: luca.sterpone@polito.it

L Carro
email: carro@eletro.ufrgs.brfglima@inf.ufrgs.br

M Sonza Reorda
email: matteo.sonzareorda@polito.it

On the Optimal Design of Triple Modular Redundancy Logic for SRAM-based FPGAs

Triple Modular Redundancy (TMR) is a suitable fault tolerant technique for SRAM-based FPGA. However, one of the main challenges in achieving 100% robustness in designs protected by TMR running on programmable platforms is to prevent upsets in the routing from provoking undesirable connections between signals from distinct redundant logic parts, which can generate an error in the output. This paper investigates the optimal design of the TMR logic (e.g., by cleverly inserting voters) to ensure robustness. Four different versions of a TMR digital filter were analyzed by fault injection. Faults were randomly inserted straight into the bitstream of the FPGA. The experimental results presented in this paper demonstrate that the number and placement of voters in the TMR design can directly affect the fault tolerance, ranging from 4.03% to 0.98% the number of upsets in the routing able to cause an error in the TMR circuit.

Introduction

Triple Modular Redundancy (TMR) is a well-known fault tolerant technique for avoiding errors in integrated circuits. The TMR scheme uses three identical logic blocks performing the same task in tandem with corresponding outputs being compared through majority voters. TMR is especially suitable for protecting designs synthesized in SRAM-based Field Programmable Gate Arrays (FPGAs) [START_REF] Carmichael | Triple Module Redundancy Design Techniques for Virtex® Series FPGA: Application Notes 197[END_REF], due to the peculiar effect of an upset in the user's combinational logic design.

In a SRAM-based FPGA, both the user's combinational and sequential logic are implemented in programmable complex logic blocks (CLBs), which are customizable by writing in the SRAM cells of the program memory. When a charged particle strikes one of the sensitive nodes of a cell in such a memory, as a drain in an off state transistor, it generates a transient current pulse that can turn the gate of the opposite transistor on. The effect can produce an inversion in the stored value (in other words, a bit flip in the memory cell). This is called Single Event Upset (SEU), or simply upset. The value of an SRAM cell affects either the CLB logic function (by implementing it in a Lookup Table or LUT), or the connections among CLBs (by customizing the routing). Changing the value of such a cell has a transient effect followed by a permanent effect on the controlled logic or routing [START_REF] Lima Kastensmidt | Designing and Testing Fault-Tolerant Techniques for SRAM-based FPGAs[END_REF]. The upsets in the routing represent the main concern, as 90% of the SRAM cells inside the FPGA control the routing. The main effects of an upset in the routing are open wires and shortcuts between distinct wires. When an upset in the routing affects signals from two distinct redundant logic parts, for example, provoking an undesirable connection between redundant logic part 0 and redundant logic part 1, the TMR approach may not vote the correct output. The probability of this occurrence depends on the logic placement and the number of majority voters in the design. Previous results from bitstream fault injection [START_REF] Lima | A fault injection analysis of Virtex FPGA TMR design methodology[END_REF][START_REF] Fuller | Radiation Test Results of the Virtex FPGA and ZBT SRAM for Space Based Reconfigurable Computing[END_REF] and radiation ground testing [START_REF] Rebaudengo | Simulation-based Analysis of SEU effects of SRAM-based FPGAs[END_REF] showed that there are few upsets in the routing that can generate an error in the output. This work focuses on how to avoid these types of upsets.

The majority voters perform a very important task in the TMR approach, because they are able to block the effects of an upset through the logic at the final output. In this way, the voters can be placed in the end of the combinational and sequential logic blocks, creating barriers for the upset effects. The problem is to determine the optimal partition of the TMR logic that must be voted inside the circuit in order to reduce the probability that upsets in the routing affect two distinct redundant parts that are voted by the same voters. A small size block partition requires a large number of majority voters that may be too costly in terms of area and performance. On the other hand, placing only voters at the last output increases the probability of an upset in the routing affecting two distinct redundant logic parts overcoming the TMR.

The goal of this paper is to investigate the optimal partition of the TMR logic that must be voted in such a way that the minimal number of voters is used and the maximal upset tolerance is achieved. The considered case study circuit is a digital filter. This design has a large amount of combinational and sequential logic corresponding to dedicated multipliers and adders and registers present in each tap of the filter. To evaluate the TMR robustness a previously developed upset injection tool [START_REF] Bernardi | On the evaluation of SEU sensitiveness in SRAM-based FPGAs[END_REF] was used to randomly inject faults in the FPGA bitstream. Results show that there is an optimal solution in terms of number and placement of voters and robustness, which in the case of the filter can reduce from 4.03% to 0.98% (four times smaller!) the number of upsets in the routing able to cause an error in the TMR.

Related Work

The correct implementation of a TMR circuitry into a SRAM-based FPGA depends on the type of data structure to be protected. As stated in [START_REF] Carmichael | Triple Module Redundancy Design Techniques for Virtex® Series FPGA: Application Notes 197[END_REF], the logic may be grouped into four different structure types: Throughput Logic, State-machine Logic, I/O Logic, and Special Features (Select block RAM, DLLs, etc.). The Throughput Logic is a logic module of any size or functionality, synchronous or asynchronous, where all of the logic paths flow from the inputs to the outputs of the module without ever forming a logic loop. In this case, it is only necessary to triplicate the logic, creating three redundant logic parts (tr0, tr1 and tr2), as presented in fig. 1.

If an upset occurs in one of the redundant combination logic parts (LUTs or routing), its effect will remain until the load of the next bitstream. The constant reconfiguration of the device avoids the accumulation of upsets in the programmable matrix. This continuous loading of the bitstream is called scrubbing, and it does interrupt the application. It is important to notice that in Throughput Logic structures composed by registers, the only way to correct an upset in a register is by loading a new data in the input of the register, or by implementing this refreshing structure with voters (see fig. 2). In the case of the registers, it is not possible to load the configuration bitstream without interrupting the application because the correct state of the register cannot be saved and loaded by the bitstream.

State-machine logic is any structure where a registered output, at any register stage within the module, is fed back into any prior stage within the module, forming a registered logic loop. This structure is used in accumulators, counters, or any custom state-machine or state-sequencer where the given state of the internal registers is dependent on its own previous state. In this case, it is necessary to triplicate the logic and have majority voters in the outputs. The register cannot be locked in a wrong value, and for this reason there is a voter for each redundant logic part in the feedback path, making the system able to recover by itself. The same structure presented in fig. 2 can be used. One majority voter can be implemented by one LUT. Because the LUT can be upset (permanent effect), the voters are also triplicated. In this way, if one voter is upset, there are still two voters working properly.

The primary purpose of using a TMR design methodology is to remove all single points of failure from the design. This begins with the FPGA inputs. If a single input was connected to all three redundant logic paths within the FPGA, then a failure at that input would cause these errors to propagate through all the redundancies, and thus the error would not be mitigated. Therefore, each redundant logic part of the design that uses FPGA inputs should have its own set of inputs, as shown in fig. 1. Thus, if one of the inputs suffers a failure, it will only affect one of the redundant logic parts. The outputs are the key to the overall TMR strategy. Since the full triple module redundancy generates every logic path triplicated, the TMR output majority voters, inside the output logic block, allow converging the output again to one signal outside the FPGA, as presented in fig. 1.

There are mainly four types of upsets that may occur in designs synthesized into the FPGA matrix, whose effects are summarized in table 1. They can be classified by the upset location: upsets in the logic (LUT), upsets in the customization routing bits inside the CLB (customization logic in general), upsets in the routing connecting CLBs and pins (routing), and upsets in the CLB flip-flops (flip-flops). The percentages of each type of cells in the whole set of customizable elements in the programmable matrix of our filter are as follows: the LUTs represent 7.4%, the flip-flops represent 0.46%, the customization bits in the CLB represent 6.36% and the general routing represents 82.9%. Based on these values, one can conclude that the probability of an upset affecting the registers is very low, compared to the probability of this same upset affecting the routing, which is our main concern as the other type of upsets are protected by TMR. The probability that an upset in the routing overcomes the TMR is related to the routing density and logic placement. In fig. 1, there are two examples of upsets in the routing. Upset "a" connects two signals from the same redundant part, which does not generate an error in the TMR output, because the upset effect will be voted by the outermost voters. However, upset "b" may provoke an error in the TMR output, because it connects two signals from distinct redundant logic parts affecting two out of three redundant parts of the TMR.

Dedicated floorplanning for each redundant part of the TMR can reduce the probability of upsets in the routing affecting two or more logic modules, but it may not be sufficient, since placement can be too complex in some cases. Remember that each time it is necessary to include voters, there are connections between the redundant parts (see fig. 2), which make impossible to place the redundant logic parts very far away from each other with no connections at all.

If the redundant logic parts tr0, tr1 and tr2 (represented in fig. 1 after the TMR register with voters and refresh) were partitioned in smaller logic blocks with voters, a connection between signals from distinct redundant parts could be voted by different voters. This logic partition by voters is represented in fig. 3. Notice that now the upset "b" can not provoke an error in the TMR output, which increases the robustness of the TMR in the presence of routing upsets without being of concern to floorplanning. The problem is to evaluate the best size of the logic to achieve the best robustness. If the logic is partitioned in very small blocks, the number of voters will increase dramatically, causing an overly costly TMR implementation. The objective is finding the best partition in terms of area cost, performance and robustness.

The Case Study Design: Digital Fir Filter

The designed test case to evaluate the robustness of the TMR according to the logic partition is an 11 taps 9-bit digital low-pass filter. The original coefficients calculated by Matlab [START_REF] Mathworks | Matlab and Simulink Documentation[END_REF] were multiplied by the constant 512. The final multiplier coefficients are: 1, -1, -9, 6, 73 and 120. In this filter design, there are eleven dedicated 9-bit multipliers, ten 18-bit adders and ten 9-bit registers.

An upset can affect the registers (in this case causing a transient effect), or can affect the logic (multipliers, adders, voters), causing a permanent effect. Five different versions of the filter were implemented. The first is the standard one with no protection at all (filter). The second version was protected by TMR using the maximum logic partition (TMR_p1). In this case, each combinational logic component, such as an adder or a multiplier, was triplicated and majority voters were inserted in the output as presented in fig. 4(a). The third implementation is a TMR filter, where each partition logic block has one adder, one multiplier and voters placed in the output (TMR_p2), as shown in fig. 4(b). The fourth TMR filter version has only majority voters in the outermost output signals (TMR_p3), as presented in fig. 4(c). All the TMR implementations so far have TMR registers with voters and refresh. The last implementation differs from the fourth version because the registers in this case are not voted (TMR_p3_nv). They are just triplicated and the final majority voters are responsible for voting the upsets in the combinational logic (adders and multipliers) and in the sequential logic (registers), as well as the routing upsets.

The comparison between the main characteristics of the four implemented versions of the FIR digital filter are reported in Table 2, where slice reports the number of FPGA slices that the circuit uses, bitstream is composed of the routing bits, LUTs bits and CLB flip-flops. The routing bits is the number of configuration memory bits for interconnections and multiplexers used for routing the signals of the circuits, LUTs bits is the number of configuration memory bits used to program the LUTs within the CLBs, and CLB flip-flops is the number of configuration memory bits used to program the flip-flops within the CLBs.

The amount of routing bits in the unprotected design corresponds roughly to 80% of the total customizable bits. In the TMR version, the design with the maximum partition has 77% of the total customizable bits corresponding to routing bits, whereas the minimum partition presents 81% of routing bits. This result is very interesting, because it shows that increasing the logic partition with voters actually decreases the percentage of routing bits compared to the logic. This can happen because of the amount to logic that is inserted by the voters.

Table 2. Comparison between TMR partitioned designs in XC2S200E-PQ208

Fault Injection Experiments

The robustness of the TMR technique implemented in a high-level description language synthesized in the Spartan ® FPGA was evaluated by injecting faults in the configuration bits of the matrix (LUTs and configuration routing cells). We used a fault-injection system we developed [START_REF] Bernardi | On the evaluation of SEU sensitiveness in SRAM-based FPGAs[END_REF] which is composed of three modules: 1. Fault List Manager: it generates the list of faulty bitstreams to be injected within the circuit under analysis. 2. Fault Injection Manager: it manages the fault injection process, by selecting one fault from the fault list, performing its injection in the circuit under analysis, and then observing and analyzing the obtained results to provide the fault classification. 3. Golden device: it is a copy of the DUT without the usage of TMR. The fault injection system consists of the Fault List Manager implemented as a software process that runs on a PC, a Fault Injection Manager running in part on a PC as a software process, and in part on the same FPGA where the circuit under analysis is placed, and finally a Golden device placed on the same FPGA, and linked with the Fault List Manager and with the DUT. Every clock cycle their output signals are compared with those coming from the DUT.

In order to improve the fault-injection process, we designed the Fault List Manager so that it is able to identify the configuration memory bits that are actually programmed to implement the DUT and generate the bitflips only for them. To implement this solution we first generated the file storing the configuration memory file for the whole design; secondly the Fault List Manager identified those bits related to the DUT. Please note that faults were selected in such a way that common-mode errors were not possible. This process is possible thanks to a data base of the programmed resources (LUTs and configuration routing cells) we developed by decoding the Xilinx bitstream. The list of faults for the DUT is finally computed and stored. Each element of the fault list is a faulty bitstream for the FPGA, where one bit at a time is modified emulating the effect of a SEU.

In the used fault-injection system the DUT and the Fault Injection Manager implemented on the FPGA share the same FPGA device. However, to guarantee that the modules are placed on the FPGA device in such a way that any fault injected in the DUT does not interfere with the Fault Injection Manager, we constrained the place and route tools to organize the FPGA-resource floorplanning allocation. The DUT and the Fault Injection Manager are placed separately on the floorplanning by at least 10 slice columns to avoid the link between the two units generated by an undesirable interconnection due to an upset.

The FPGA device used in our experiment is a Spartan ® XC2S200EPQ208 from Xilinx [START_REF]XILINX, INC. Spartan®™ 2.5 V Field Programmable Gate Arrays: Datasheet DS003[END_REF], whose configuration memory is composed of 1,442,016 bits organized in 2,501 frames of 576 bits each. The configuration memory controls an array of 28 x 42 slices. The robustness of the TMR technique implemented in a high-level description language was observed while running five fault-injection campaigns, one for each circuit.

The fault injection process took about 13 seconds for each fault injected. The test pattern generation and the output analyzer placed in the same FPGA need a negligible time for applying all the input stimuli and classifying the faults (5.7 ms on the average). Most of the time for evaluating each injected fault was due to the time for downloading the faulty bitstream within the SRAMbased FPGA. The results we gained performing the faultinjection campaigns are reported in table 3, where Injected Faults reports the number of SEUs injected, and Wrong Answer reports the percentage of SEUs provoking a fault within the TMR DUT.

Roughly 10% of the whole configuration memory bits were injected. This set was randomly chosen from the fault list. It would take approximately 3 months for testing the whole configuration memory bits of the five circuits. The results confirm our original hypothesis, since they show that placing voters just at the final output is not sufficient to avoid errors. There is an optimum size of triplicated logic that must be voted inside the design, creating voter barriers along the logic that help block the effect of upsets in the routing to propagate to the output. The results experimentally demonstrate that the logic partition done in the TMR_p2 version reduce the percentage of uncovered faults from 4.03% to 0.98%, corresponding to a four times reduction. This solution is also efficient in term of estimated performance, as shown in Table 2.

Table 3. Fault injection campaign results

Wrong 4 reports a classification of the effects of the injected upsets that caused an error in the TMR. The effect analysis of each faulty bitstream was done using the classification tool developed in [START_REF] Bellato | Evaluating the effects of SEUs affecting the configuration memory of an SRAMbased FPGA[END_REF]. The results confirm that the routing resources are the most sensitive to upsets. No upsets in the LUTs could provoke an error in the TMR design, as it was expected. Upsets in the CLB multiplexers (mux) and in the initialization that are related to the routing resources represent less than 8% of the total effects. The main problem is the general routing that connects CLBs and I/O pads, corresponding to more than 90% of routing resources. The general routing effects that can provoke an error in the TMR output are classified as Open, Bridge, Input-Antenna, Conflict and others. The Open effect corresponds to a Programmable Interconnection Point (PIP) in open state (open connection). The Bridge effect means that a new PIP is permitted. Both Bridge and Open effects may influence the behavior of the circuit. The Input-Antenna effect also represents a new permitted PIP. In this case, this new PIP connects an unused input node to a used output node; this may influence the behavior of the circuit, since the CLB or output pad are driven to an unknown logic value. The Conflict effect is when a new PIP links both used input and output nodes, creating a conflict and a propagation of an unknown value along the TMR. These effects are a real challenge for those designers who are involved in devising solutions for hardening their FPGA-based circuits, and these effects must be avoided or at least minimized.

It is important to notice that one upset can provoke more that one effect in the programmable matrix. The same can be said about different upsets in the bitstream that can present the same effect. This explains why the number of table 3 and 4 are not exactly the same.

Conclusions

The results presented in this paper suggest that there is a trade off between the logic partition of the throughput logic (and consequently between the number of voters) and the number of routing upsets that could provoke an error in the TMR. In contrary to what was expected, a large number of voters does not always mean larger protection against upsets. There is an optimal logic partition for each circuit that can reduce the propagation of the upset effect in the routing. For the case study circuit, the best partition is the medium partition (TMR_p2). This version of the TMR design presents a small sensitivity to routing upsets (0.98%, a four times improvement over normal TMR) and small performance degradation (about 10%) compared to the standard version (not protected).

Future works will consider analyzing the effect of logic partition in power dissipation, and resolving the problem of the uncovered routing upsets by means of a combination of logic partition and dedicated floorplanning.

Figure 1 .

 1 Figure 1. Triple Modular Redundancy (TMR) Scheme in the FPGA

Figure 2 -

 2 Figure 2-TMR register with voters and scrubbing

Figure 3 .

 3 Figure 3. Triple Modular Redundancy (TMR) scheme with logic partition in the FPGA

Figure 4 .

 4 TMR digital filter schemes used to evaluate the optimal logic partition of TMR

Table 1 . Upset analysis in the Triple Modular Redundancy approach

 1

	Upset Location	Upset Effect	Consequences	Upset Correction
	LUT	Modification in the	-Error in one redundant part with no error	By scrubbing
		Combinational logic	in the TMR design output
	Routing	Connection or disconnection	-Error in the redundant part with no error	By scrubbing
		between any two or more signals	in the TMR design output
		in the design	-Error in more than one redundant part
			with error in the design output
	Customization	Connection or disconnection	-Error in the redundant part with no error	By scrubbing
	logic in	between any two signals in the	in the TMR design output
	general	same CLB	-Error in more than one redundant part
			with error in the design output
	Flip-flops	Modification in the sequential	-Error in the redundant part, no error in	By design
		logic	the TMR design output	modification

Table 4 . Effects induced by the injected upsets

 4

			Standard	TMR_p1	TMR_p2	TMR_p3	TMR_p3_nv
			[#]	[%]	[#]	[%]	[#]	[%]	[#]	[%]	[#]	[%]
	CLBs	LUT	852	16	0	0	0	0	0	0	0	16
	Logic and	MUX	123	2	16	2	1	0	15	4	367	2
	Routing	Initialization	174	3	13	1	0	0	11	3	400	3
		Open	1321	25	276	38	82	40	126	37	1672	25
	General	Bridge	427	8	62	9	41	20	42	12	403	8
	Routing	Input-Antenna	76	1	33	5	7	3	14	4	73	1
		Conflict	1342	25	26	4	13	6	6	2	185	25
		Others	1006	20	301	41	66	31	128	38	756	20
		Total	5321		727		210		342		3856

Proceedings of the Design, Automation and Test in Europe Conference and Exhibition (DATE'05) 1530-1591/05 $ 20.00 IEEE

Proceedings of the Design, Automation and Test in Europe Conference and Exhibition (DATE'05)


Jean-Francois Crouzet
email: jean-francois.crouzet@univ-montp2.fr

Jean-François Crouzet

Multi

Keywords: Fuzzy transform, fuzzy projection, signal representation, image representation, stability, robustness

published or not. The documents may come

Introduction

Fuzzy and inverse fuzzy transforms for the representation of signal and images have been initiated by Irina Perfilieva in 2004 ([? ],[? ]) and has led to many applications in their processing ( [? ],[? ],[? ],[? ],[? ],[? ],[? ]). In their continuity the author of the present study introduced in 2010 the operator of fuzzy projection for signals ([? ],[? ]). Independently, Giuseppe Patanè ([? ]) proposed in 2011 a study of the fuzzy transform via leastsquares approximation that brought out the same operator. Fuzzy projection (F-projection) proved to perform a good interpolation of the samples of a signal obtained by fuzzy transform (F-transform). In particularly it handles correctly high frequencies, even in presence of an important noise ([? ]).

In this paper, we present and discuss of F-projection in higher dimension of space. We propose to show it is a simple tool, very easy to implement numerically, especially for image analysis. (We will deal in a future work ([? ]) with numerical experiments such as compression, denoising and restauration of images, particularly a comparison with wavelet representation.)

Moreover, we point out that a key parameter permits to control the robustness of fuzzy projection, i.e. how much it is numerically well conditioned. It also indicates its stability, say its sensitivity, when choosing different partitions of the space, i.e. different generic kernels of the fuzzy partition: the effect of this choice does not modify significantly the representation process. This parameter is simply the squared L 2 -norm of the generic kernel µ, generating the fuzzy partition in dimension 1: a = µ [START_REF] Crouzet | The fuzzy projection: a simple adaptive representation of signals[END_REF] 2 . As it measures, in any dimension of space, how the generic kernel of the fuzzy partition is concentrated around the center of its supporting cell, we name it the 'localization parameter'. Thanks to the parameter a, we discuss the linear representation of the considered object. More precisely, we suggest that choosing a linear fuzzy partition (a triangular one for signals, or a pyramidal one for images) is a reasonable compromise between numerical robustness, stability of the representation, and the requirements of a fuzzy representation of a signal or a multi-dimensional object.

In section 2, we recall definitions and results necessary to present F-projection in dimension 1 of space that are been presented and proved in [? ]. In section 3, we highlight the natural analogy of these results in dimension 2 of space. In section 4, a rapid generalization is made in any dimension. Section 5 adresses a notion of stability for this representation, highlighting it is controlled by the 'localisation parameter'. We conclude with discussing the relevance of a linear fuzzy partition.

Fuzzy projection in dimension 1

Following the seminal work of Irina Perfilieva ([? ],[? ]) on direct and inverse F-transform, we collect in this section the results on F-projection that were proved in [? ] and [? ]. This will permit to introduce easily multidimensional F-projection, and the role of the 'localization parameter'.

Uniform fuzzy partition

A fuzzy partition P = {µ k , k = 0 . . . p + 1} of an interval is defined as follows.

Definition 2.1.1. (fuzzy partition P)

Let S = {ω k = k∆ , k = 0 . . . p + 1} be an uniform subdivision of [0, d] with diameter ∆ = d/(p + 1). Let µ ∈ C 0 ([-1, 1]
), the generic kernel of the partition, be an even continuous function such that µ(0) = 1, µ(1) = 0 and µ(1

-x) + µ(x) = 1 , ∀x ∈ [-1, 1]. Kernels µ k of the partition P are defined by µ k (x) = µ x-ω k ∆
, and have thus for support the cells

Ω k = [ω k-1 , ω k+1 ] (with the convention that ω -1 = ω 0 = 0 and ω p+2 = ω p+1 = d.)
Then, it is immediate to check that P is a partition of unity on [0, d]:

Property 2.1.1. (partition of unity) i) ∀x ∈ [0, d] , p+1 k=0 µ k (x) = 1 (1)
ii)

∀k = 0 . . . p + 1 , 1 ∆ d 0 µ k = 1 ∆ Ω k µ k = 1 -1 µ = 1 (2)
∀k = 0 . . . p + 1, (F f ) k := d 0 f µ k d 0 µ k = 1 ∆ (f |µ k ) (3)
Fuzzy transform thus consists in taking a sequence of means, weighted by kernels µ k , of the signal f . Note that it maps constant functions to a constant sampling with the same value. We also need to express its adjoint. Consider, for u = (u 0 , . . . , u p+1 ) and v = (v 0 , . . . , v p+1 ), the scalar

product in R p+2 : u.v = ∆ p+1 k=0 u k v k . For γ = (γ 0 , . . . , γ p+1 ), γ.F f = ∆ p+1 k=0 γ k (F f ) k = p+1 k=0 γ k (µ k |f ) = f | p+1 k=0 γ k µ k yields the expression of F * : R p+2 → L 2 ([0, d]):
Proposition 2.2.1. (adjoint fuzzy transform F * : interpolation) For γ = (γ 0 , . . . , γ p+1 ), the adjoint fuzzy transform is

F * γ = p+1 k=0 γ k µ k (4)
and can be viewed as an interpolation of the finite sequence γ 0 , . . . , γ p+1 by kernels µ k of the partition.

Remark 2.2.1. By 'interpolation' we mean in this study the passage from a discrete sampled representation of a signal (say F f in our context) to a continuous representation by mean of kernels functions µ k of the partition (say for instance F * F f for the inverse fuzzy transform introduced by I. Perfilieva). It is directly related to fuzzy projection that will be treated in the following section, where the interpolation functions are different.

Fuzzy projection

Let us first define the p+2-dimensional sub-space Π of functions generated by partition P, or, equivalently, the image of F * :

Definition 2.3.1. (projection sub-space Π) Π = p+1 k=0 γ k µ k , γ ∈ R p+2 (5)
For f ∈ L 2 ([0, d]), F-projection is defined by: Definition 2.3.2. (F-projection P )

P f = Argmin g∈Π f -g 2 (6)
Computing P f is thus equivalent to find

α = Argmin γ∈R p+2 f -F * γ 2 (7)
Its solution is given by the normal equation F F * α = F f , so that F-projection is the the classical pseudo-inverse projector ([? ]):

Proposition 2.3.1. (pseudo-inverse)
The expression of the F-projection P : L 2 ([0, d]) → Π is the following:

P f = F * (F F * ) -1 F f (8)
It is established in [? ] uniform and quadratic convergence (when p → ∞, i.e. ∆ = d p+1 → 0):

Theorem 2.3.1. (convergence) i) Let f ∈ C 0 ([0, d]). Then ∀ε > 0 , ∃∆ | f -P f ∞ ≤ ε. ii) Let f ∈ L 2 ([0, d]). Then ∀ε > 0 , ∃∆ | f -P f 2 ≤ ε.
Let us now define the announced control parameter with which we will deal till the end of the present study.

Definition and Property 2.3.1. (localization parameter a)

We define

a = µ 2 2 = 1 -1 µ 2 (9)
as the 'localization parameter' for µ. The fact that a = 1/2 + 4

1/2 0 (µ -1/2) 2 shows that it represents the degree of localization of µ in the central interval [-1/2, 1/2] of its support [-1, 1]. It varies from 1/2 (limit case of the step function with maximum overlap, µ = 1 2 χ [-1,1] ) to 1 (limit case of the basic step function with no overlap, µ = χ [-1/2,1/2] ):

1 2 < a < 1 (10)
Remark 2.3.1. There exists, of course, an infinity of generic kernels providing a same parameter a. For instance, the triangle kernel, corresponding to linear fuzzy sampling, verifies a = 2/3, and it is easy to construct other ones with the same parameter a, but presenting oscillating shapes (even if they strictly decrease on [0, 1], as demanded). Nevertheless, it is reasonable to keep in view this most natural and simple kernel.

Let us now identify F F * to its matrix M in the canonical basis of R p+2 . Then, as proved in [? ], we can state the following result: Theorem 2.3.2. (expression of F-projection and role of a) i) F F * corresponds to a (p + 2) × (p + 2) (almost bistochastic) symmetrical tridiagonal matrix

M =           a 1-a 2 1-a 2 a 1-a 2 (0) 1-a 2 a 1-a 2 . . . . . (0) 1-a 2 a 1-a 2 1-a 2 a           (11)
depending only on the localization parameter a.

ii) M is invertible and has p + 2 simple eigenvalues

λ k = a + (1 -a) cos kπ p + 1 ( 12
)
verifying 0 < 2a -1 ≤ λ k ≤ 1, k = 0 . . . p + 1.
iii) The associated normalized eigenvectors are denoted v k and defined by

(v k ) j = 2 d 1 2 cos jkπ p + 1 ( 13
)
They do not depend on the generic kernel.

iv) The set of functions

B = f k = 1 √ λ k F * v k , k = 0 . . . p + 1 (14)
forms an orthonormal basis of the projection subspace Π. Thus, for any f ∈ L 2 ([0, d]) we have

P f = p+1 k=0 (f |f k ) f k (15)
v) The functions

ν k (x) = p+1 j=0 M -1 j,k µ j (x) (16)
make P f reads equivalently as the following interpolation -natural continuous representation -of the fuzzy transform F f :

P f = p+1 k=0 (F f ) k ν k (17) Remark 2.3.2.
1) The first point of the theorem comes from a computation that the reader is invited to read in [? ]. Points ii) and ii) are then a classical property of tridiagonal matrixes (see for instance [? ]). As F F * is identified to the matrix M with eigenvalues λ k , point iv) is immediate.

2) The length d of the interval [0, d] where signals f are considered appears in the formula

(v k ) j = 2 d 1 2 cos jkπ p+1
. This is due to the normalization of those vectors with the weighted scalar product u.v = ∆ p+1 k=0 u k v k , where the diameter of the subdivision verifies ∆ = d/(p + 1).

3) Basis functions f

k = 1 √ λ k F * v k are the interpolations by F * of the finite sequence defined by (v k ) j = 2 d 1 2 cos jkπ p+1 , j = 0 . . . p + 1.
Thus, they look like a Fourier basis, but with the shape of the generic kernel µ appearing in F * (for instance a piecewise linear sinus-like function for linear interpolation with a triangular kernel.)

4) P f (x) = F * M -1 F f (x) = p+1 j=0 µ j (x) p+1 k=0 M -1 j,k (F f ) k provides formula (??) in point v).
Two important numerical requirements are that matrix M must be well conditioned and that functions {ν k } (that bring the discrete samples to a continuous interpolation representation) have well localized supports: that is what we call the 'robustness' of the fuzzy representation.

c := cond 2 (M ) = max |λ k | min |λ k | = 1 2a -1 (18) ii) Each function ν k (x) = p+1 j=0 µ j (x)M -1 j,k
, that operates the interpolation from the discrete sampling F f to its continuous approximation P f , is exponentially decreasing in the sense that |M -1 j,k | ≤ |M -1 k,k |r |k-j| for all indexes j and k. The speed of decrease around the diagonal 0 < r < 1 is given by

r = 1 -a 3a -1 (19) Remark 2.3.3. 1) The fact that that |M -1 j,k | ≤ |M -1 k,k |r |k-j|
, for all indexes j and k, means that the coefficient of the bounded and compactly supported kernels functions µ j in the expression of ν k are exponentially decreasing. Thus, function ν k itself is also exponentially decreasing. An illustration is given in Figure ?? for triangular kernels.

2) In the limit case of µ = 1 2 χ [-1,1] , we have a = 1/2, so that c = ∞ and r = 1: this is the worst numerical situation.

3) In the limit case of µ = χ [-1/2,1/2] , we have a = 1, so that c = 1 and r = 0: this is the best numerical situation, but we will discuss this case in Section ??. 4) In the case of the classical triangular kernel µ = Tri [-1,1] , we have a = 2/3, so that c = 3 and r = 1/3: this is a good numerical situation; we will also discuss this case in the last Section ??, as it is the situation we want to bring out in this study. 5) In the case of the classical sinusoidal kernel : µ(x) = 1 2 (1 + cos(πx)), we have a = 3/4, so that c = 2 and r = 1/5: this is an even better numerical situation we will also comment Section ??.

Thus, while r = r(a) controls the localization of the interpolating functions, the condition number c = c(a) controls the numerical sensibility of the representation. This shows that the localization parameter a plays then a key role in basic numerical requirements. Moreover, we can also note that a governs the rapidity of the uniform convergence in Theorem ??, as shown by its proof in [? ].

Fuzzy projection in dimension 2

This section is very analogous to the previous one. However, as it is useful for image representation, we present all the results, following exactly the same way. Most of their proofs are let to the reader, as they are nearly the same that those in dimension 1 (presented in [? ]).

) Let T = (ω k , ω l ) = (k∆, l∆) , (k, l) ∈ [0 . . . p + 1] 2 be an uniform subdi- vision of [0, d] 2 . Still, ∆ = d/(p + 1)
is called the diameter of the subdivision. Recall that µ denotes the generic kernel in the 1D Definition ??. Kernels M k,l of the partition T are then defined by

M k,l (x, y) = µ k (x)µ l (y) := µ x-ω k ∆ µ y-ω l ∆
, and are supported in the cells

Ω k,l = [ω k-1 , ω k+1 ] × [ω l-1 , ω l+1 ]. The generic kernel is M(x, y) = µ(x)µ(y), sup- ported in the unit square [-1, 1] 2 .
Again, it is easy to verify that Q is a partition of unity on [0, d] 2 , precisely:

Property 3.1.1. (2D partition of unity) i) ∀(x, y) ∈ [0, d] 2 , p+1 k=0 p+1 l=0 M k,l (x, y) = 1 (20) ii) ∀(k, l) ∈ [0 . . . p + 1] 2 , 1 ∆ 2 [0,d] 2 M k,l = 1 ∆ 2 Ω k,l M k,l = 1 -1 µ 2 = 1 (21) Remark 3.1.1.
As an illustration, the classical 1D triangular generic kernels yields a partition with square pyramids.

2D F-transform

The 2D fuzzy transform

F : L 2 ([0, d] 2 ) → R p+2 × R p+2 is defined as follows: Definition 3.2.1. (2D F-transform F ) ∀(k, l) ∈ [0 . . . p + 1] 2 , (F f ) k,l := [0,d] 2 f M k,l [0,d] 2 M k,l = 1 ∆ 2 (f |M k,l ) (22)
Considering on R p+2 × R p+2 the scalar product u.v = ∆ 2 p+1 k=0 p+1 l=0 u k,l v k,l , we have similarly the adjoint fuzzy transform

F * : R p+2 × R p+2 → L 2 ([0, d] 2 ): Proposition 3.2.1. (2D adjoint F-transform F * ) F * γ = p+1 k=0 p+1 l=0 γ k,l M k,l (23)

2D fuzzy projection

The (p + 2) 2 -dimensional sub-space Π generated by the partition Q is now:

Definition 3.3.1. (projection sub-space Π) Π = p+1 k=0 p+1 l=0 γ k,l M k,l , γ ∈ R p+2 × R p+2 (24)
This leads to the same definition, expression and convergence results for F-projection:

Definition 3.3.2. (2D F-projection P ) P f = Argmin g∈Π f -g 2 (25) Proposition 3.3.1. (pseudo-inverse)
The expression of the 2D F-projection P : L 2 ([0, d] 2 ) → Π is the following:

P f = F * (F F * ) -1 F f (26) Theorem 3.3.1. (2D convergence) i) Let f ∈ C 0 ([0, d] 2 ). Then ∀ε > 0 , ∃∆ | f -P f ∞ ≤ ε. ii) Let f ∈ L 2 ([0, d] 2 ). Then ∀ε > 0 , ∃∆ | f -P f 2 ≤ ε.
Let us now identify R p+2 × R p+2 to R (p+2) 2 by the bijective map (k, l) → k + l(p + 2). Then, the components of the fuzzy transform (F f ) k,l will be stored in a (p+2) 2 -dimensional vector. Thus, operator F F * can be identified to a matrix N in the canonical basis of R (p+2) 2 . This allows to get the following result. (Notations are the same as in 1D Theorem ??; in particular for a, λ k , v k and matrix M , given in formula (??) in the 1D case.)

N =           aM 1-a 2 M 1-a 2 M aM 1-a 2 M (0) 1-a 2 M aM 1-a 2 M . . . . . (0) 1-a 2 M aM 1-a 2 M 1-a 2 M aM           (27)
depending only on the localization parameter a.

ii) F F * is invertible and has (p + 2) 2 simple eigenvalues

Λ k,l (a) = λ k (a)λ l (a) = a + (1 -a) cos kπ p + 1 a + (1 -a) cos lπ p + 1 (28) verifying 0 < (2a -1) 2 ≤ Λ k,l ≤ 1, ∀(k, l) ∈ [0 . . . (p + 1)] 2 .
iii) The associated normalized 'eigenvectors' are V k,l = v k • v l with:

(V k,l ) j,m = (v k • v l ) j,m := (v k ) j (v l ) m = 2 d cos jkπ p + 1 cos mlπ p + 1 (29)
Again, they do not depend on the generic kernel.

iv) The set of functions

B = f k,l = 1 Λ k,l F * V k,l , (k, l) ∈ [0 . . . p + 1] 2 (30)
forms an orthonormal basis of the projection subspace Π. Thus,

P f = p+1 k=0 p+1 l=0 (f |f k,l ) f k,l (31)
v) The functions

N k,l (x, y) = p+1 i=0 p+1 j=0 (F F * ) -1 (i,j),(k,l) M i,j (x, y) (32)
make P f reads as the following interpolation of the fuzzy transform F f :

P f (x, y) = p+1 k=0 p+1 l=0 (F f ) k,l N k,l (x, y) (33)
Basic numerical requirements are the same: operator F F * (associated, for numerical purposes, to the matrix N ) must be well conditioned and interpolation functions N k,l well localized. Here is a slightly similar result, point ii) being however more technical.

C := cond 2 (F F * ) = max |Λ k,l (a)| min |Λ k,l (a)| = 1 (2a -1) 2 (34) ii) Interpolating function, N k,l (x, y) = p+1 i=0 p+1 j=0 (F F * ) -1
(i,j),(k,l) M i,j (x, y) are exponentially decreasing in the sense that, given the matrix N , there exists a constant α > 0 such that:

a) |N -1 i,j | ≤ αR i , if i > j and i is large enough, b) |N -1
i,j | ≤ αR j , if j > i and j is large enough. The speed of decrease away from the diagonal 0 < R < 1 is given by

R = 1 -a a + √ 2a -1 (35) Proof of ii).
It is proved in [? ] that a block tridiagonal block toeplitz matrix of the form

T =       A B B A B (0) . . . (0) . . B B A      
with the same dimensions as N verifies a)

|T -1 i,j | ≤ α|ρ p+1 | -i , if i > j and i is large enough, b) |T -1 i,j | ≤ α|ρ p+1 | -j
, if j > i and j is large enough, where {ρ k } are the 2(p + 1) eigenvalues, sorted in absolute value increasing order, of the 2(p + 1) × 2(p + 1) matrix

U = 0 B -B -1 B -1 A
In our case, we have (Theorem ??) A = aM and B = 1-a 2 M . Denote λ k (A) and λ k (B) the respective eigenvalues of A and B, and recall that v k are their corresponding common eigenvectors. When searching for eigenvectors of the form

V k = α k v k v k we have that U V k = 0 B -B -1 B -1 A α k v k v k = λ k (B)v k -α k λ k (B) + λ k (A) λ k (B) v k . A necessary condition for V k to be an eigenvector is then that α k = λ k (B)/ λ k (A)-α k λ k (B)
. This yields a second degree equation for α k whose 2(p + 1)solutions are, for k = 0 . . . p + 1: -α

(1) k = 1 2 λ k (A) + λ 2 k (A) -4λ 2 k (B) -α (2) k = 1 2 λ k (A) -λ 2 k (A) -4λ 2 k (B)
Then, as ρ k , the 2(p + 1) eigenvalues of U , verify by definition

U V k = U α k v k v k = ρ k α k v k v k , we have that ρ k α k = λ k (B). It
follows that, sorted in absolute value increasing order, they write, for k = 0 . . . p:

-ρ k = 2λ k (B)/ λ k (A) + λ 2 k (A) -4λ 2 k (B) -ρ p+1+k = 2λ k (B)/ λ k (A) -λ 2 k (A) -4λ 2 k (B) In the case of A = aM and B = 1-a 2 M , the computation simply comes down to ρ p+1 = 2λ 0 ( 1-a 2 M )/ λ 0 (aM ) -λ 2 0 (aM ) -4λ 2 0 ( 1-a 2 M ) = 2 1-a 2 (2a - 1))/ a(2a -1) -(a(2a -1)) 2 -4( 1-a 2 (2a -1)) 2 = 1-a a- √ 2a-1 = a+ √ 2a-1 1-a .
Then, as R = 1/ρ p+1 , we have the expected result. Remark 3.3.1. Let us examine again the four limit and classical cases of Section ??:

1) When µ = 1 2 χ [-1,1]
, a = 1/2 still yields C = ∞ and R = 1: nothing can be done numerically. This case corresponds to M(x, y) = 1 4 χ [-1,1] (x)χ [-1,1] (y), the limit generic constant kernel with maximal 2D overlap. 2) When µ = χ [-1/2,1/2] , a = 1 also still yields C = 1 and R = 0: this optimal numerical situation will be discussed in Section ??. This case corresponds to M(x, y) = χ [-1/2,1/2] (x)χ [-1/2,1/2] (y), the classical limit plateau kernel with no overlap.

3) When µ = Tri [-1,1] , we have a = 2/3, so that C = 9 and R = 1/(2 + √ 3) 0.26 (instead of r 0.33 in dimension one): this situation remains numerically very good, and the interpolating functions are more localized than in 1D. This case will be examined in Section ??; it corresponds to M(x, y) = Tri [-1,1] (x)Tri [-1,1] (y), the classical square pyramid. 4) For M(x, y) = 1 4 (1 + cos(πx))(1 + cos(πy)), corresponding to a = 3/4, C = 4 and R = 1/(3 + 2 √ 2) 0.17 (instead of r 0.2 in dimension one),

we have a 2D sinus-like kernel which is still better numerically. 5) We observe that R < r: the speed of decrease of the matrix elements around its diagonal is faster for

N in 2D than for M in 1D. (R < r ⇔ 2a -1 < √ 2a -1, which is true because 1/2 < a < 1 ⇔ 2a -1 < 1.)
Thus, the localization parameter a, defined in the 1D case, plays the same key role in the 2D case. It describes the F-transform in a simple manner, as well as it governs the numerical robustness of the projection : how much matrix N of the projection numerically well conditioned and interpolation functions N k,l (x, y) have a 'small' support.

Fuzzy projection in higher dimension

The generalization to any dimension q of space is immediate, following the link between formulas in dimension 1 and 2. With the same formalism we establish, skipping all identical expressions and in the same way avoiding all proofs, the analogous results. We chose, for simplicity, the same notations as in one dimension of space, denoting in particular x = (x 1 , . . . , x q ) and k = (k 1 , . . . , k q ).

A fuzzy partition P = {µ k , k ∈ [0 . . . p + 1] q } of the hypercube [0, d] q has the following definition and property: Definition 4.1. (qD fuzzy partition) Let S = {ω k = k∆ , k ∈ [0 . . . p + 1] q } be an uniform subdivision of [0, d] q . Kernels µ k of the partition P are then defined similarly by µ k (x) = Π i=1...q µ x i -ω k i ∆ and are supported in the q-dimensional hypercube

Ω k = Π i=1...q [ω k i -1 , ω k i +1 ].
The generic kernel µ is supported in the unit square [-1, 1] q and is the natural tensor product of the generic kernel µ defined in dimension 1.

Property 4.1. (qD partition of unity) i) ∀x ∈ [0, d] q , k∈[0...p+1] q µ k (x) = 1 (36) ii) ∀k ∈ [0 . . . p + 1] q , 1 ∆ q [0,d] q µ k = 1 ∆ q Ω k µ k = 1 -1 µ q = 1 (37)
The q-dimensional fuzzy transform corresponds to the operator F :

L 2 ([0, d] q ) → (R p+2 )
q such that:

Definition 4.2. (fuzzy transform) ∀k ∈ [0 . . . p + 1] q , (F f ) k := [0,d] q f µ k [0,d] q µ k = 1 ∆ q (f |µ k ) ( 38
)
With the scalar product u.v = ∆ q k∈[0...p+1] q u k v k on (R p+2 ) q we have the same expression for the adjoint fuzzy transform F * . Then, the (p + 2) q -dimensional sub-space Π, F-projection on it, and also convergence results are defined and established identically.

Identifying (R p+2 ) q to R (p+2) q by the bijective map k → k 1 + k 2 (p + 2) + k 3 (p+2) 2 +• • •+k q (p+2) q
, components of the fuzzy transform (F f ) k,l will be now stored in a (p + 2) q -dimensional vector. Operator F F * is then identified to a matrix M in the canonical basis of R (p+2) q , and this yields:

Theorem 4.1. (expression of qD F-projection and role of a)) i) F F * corresponds to a (p +2) q × (p +2) q almost bistochastic matrix which is symmetrical block tridiagonal with symmetrical tridiagonal blocks, themselves defined in the same way by an induction process with length q-1, and starting at M 1 = M :

M q =           aM q-1 1-a 2 M q-1 1-a 2 M q-1 aM q-1 1-a 2 M q-1 (0) 1-a 2 M q-1 aM q-1 1-a 2 M q-1 . . . . . (0) 1-a 2 M q-1 aM q-1 1-a 2 M q-1 1-a 2 M q-1 aM q-1           (39)
It still depends only on the 1D localization parameter a.

ii) F F * is invertible and has (p + 2) q simple eigenvalues λ k such that

λ k = Π q j=1 a + (1 -a) cos k j π p + 1 (40) verifying 0 < (2a -1) q ≤ λ k ≤ 1.
iii) The associated normalized "eigenvectors" are the v k such that

(v k ) l = 2 d q 2 Π q j=1 cos l j k j π p + 1 (41)
Again, they do not depend on the generic kernel.

iv) The set of functions

B = f k = 1 √ λ k F * v k (42)
is an orthonormal basis of Π, so that

P f = k (f |f k ) f k (43)
v) The functions

ν k (x) = i (F F * ) -1 (i),(k) ν i (x) (44)
make P f reads as the following interpolation of the fuzzy transform F f :

P f (x) = k (F f ) k ν k (x) (45)
We have the analogous obvious result for the condition number of M q . Nevertheless, as shown in Theorem ??, the induction process generates a matrix M q whose blocks consist in a q -1 times repeated structure of the 1D case matrix M . Even in the 3D case, we didn't succeed in computing the speed of decrease r of the interpolating functions ν k .

All these results show that the overlap parameter a = 1 -1 µ 2 , even if defined in dimension 1, plays the same central role in any dimension of space. Both for describing the F-projection as for controlling its numerical robustness. In the next section, we will see that it is also an indicator for its stability.

Stability of multi-dimensional fuzzy projection and localization parameter a

A question arising naturally is the following: when performing a representation of a signal by fuzzy projection with different generic kernel µ 1 and µ 2 can we measure the risk of deviation when choosing one or another? In other terms, can we measure the stability, say the sensitivity, of the representation, balancing from a projection P 1 to another one, P 2 ? This question is subject to interpretation because it is not really well-posed. The sole purpose of the following section is to propose tracks for this question.

Match between projection subspaces

Here is a way to characterize, in a general manner, the match between two finite dimensional subspaces of an Hilbert space.

Let us consider two finite subspaces F 1 and F 2 of an Hilbert space, with respective orthonormal basis

B 1 = φ (1) i , i = 1 . . . n 1 and B 2 = φ (2) i , i = 1 . . . n 2 .
Call G, the n 1 × n 2 Gram matrix of scalar products:

G = φ (1) i |φ (2) j (47)
and its singular values {σ k }, sorted in decreasing order. Definition 5.1.1. (match between subspaces) We call match between F 1 and F 2 the lowest singular value of the Gram matrix of their respective basis:

Match (F 1 , F 2 ) = min k {σ k } ≤ 1 (48)
It has been proved that singular values correspond to the cosines of 'principal angles' between the two subspaces ([? ]), so that the lowest one is the cosine of the highest angle. When the latter is (at worst) π/2, F 1 and F 2 will be considered as 'distants' in this sense, and it corresponds to Match (F 1 , F 2 ) = 0. On the contrary, we are interested in the situation where Match (F 1 , F 2 ) approches 1, but for somehow distant generic kernels.

Remark 5.1.1. There are other ways to model a criterion for this match, say angle or distance. None is perfect, each depends on the nature of the problem. For example finding a cut-off in the graph of {σ k } will indicate an incompatibility between the two subspaces, or a loss of robustness in the model with which they are related.

Consider now, with the same notations and context as in section ??, two fuzzy transforms F 1 and F 2 with kernels µ [START_REF] Boffi | Characterizing the inverses of block tridiagonal[END_REF] and µ [START_REF] Crouzet | The fuzzy projection: a simple adaptive representation of signals[END_REF] , and their corresponding F-projections P 1 and P 2 on projection subspaces Π 1 and Π 2 , i.e. the subspaces that are used to provide the projection of the signal: their representation by fuzzy sampling (F f ) followed by inverse inverse fuzzy transform: interpolation by the adjoint transform F * . Denote λ

(i) k = a i + (1 -a i ) cos kπ p+1 , the eigenvalues of M i = F i F *
i , and a i = µ (i) 2 2 , for i = 1, 2. Recall finally that they are associated to the same eigenvectors v k (Theorem ??).

As in section ??, we define a localization parameter a 1,2 := µ (1) |µ [START_REF] Crouzet | The fuzzy projection: a simple adaptive representation of signals[END_REF] (49) called the 'cross localization parameter' and that measures the overlap between the two kernels. Then, we have the following result:

G k,k = σ k = λ (1,2) k λ (1) k λ (2) k , k = 0, . . . , p + 1 ( 50
)
where λ

(1,2) k = a 1,2 + (1 -a 1,2
) cos kπ p+1 are the eigenvalues of the symmetric tridiagonal matrix F 1 F * 2 with diagonal term a 1,2 = µ (1) |µ [START_REF] Crouzet | The fuzzy projection: a simple adaptive representation of signals[END_REF] and sub or super diagonal term a 1,2 -1 2 (similarly to the case of M ). Moreover, σ k form a decreasing sequence, so that

Match (Π 1 , Π 2 ) = σ p+1 = 2a 1,2 -1 (2a 1 -1) (2a 2 -1) (51)
Proof. Using respective basis of Π 1 and Π 2 in Theorem ??, we compute

G k,l = f (1) k |f (2) l = 1 λ (1) k F * 1 v k | 1 λ (2) l F * 2 v l = 1 λ (1) k λ (2) l F 2 F * 1 v k .v l . As the matrix M 1,2 = F 2 F * 1 = F 1 F * 2 has for elements µ (1) k |µ (2) l
, it is also a symmetric tridiagonal matrix (as matrix M ) with the property that its diagonal terms are a 1,2 and its sub or super diagonal terms are a 1,2 -1 2 (as both µ k also define partitions of unity). Thus, its orthogonal eigenvectors are the same as M , and its eigenvalues reads λ

(1,2) k = a 1,2 + (1 -a 1,2 ) cos kπ p+1 . Then, G k,l = 1 λ (1) k λ (2) l λ (1,2) k v k .v l proves that G is diagonal and provides the first result. Now, we have to show that σ k = a 1,2 +(1-a 1,2 ) cos kπ p+1 (a1+(1-a1) cos kπ p+1 )(a2+(1-a2) cos kπ p+1 )
form a decreasing sequence. Then, as Match (F 1 , F 2 ) := min k {σ k } = σ p+1 , the second result will be immediate. An equivalent condition is that the function ϕ

(x) = a 1,2 +(1-a 1,2 )x √ (a 1 +(1-a 1 )x)(a 2 +(1-a 2 )x) is increasing on [-1, 1]. Its logarithmic derivative is ψ(x) = 1-a 1,2 a 1,2 +(1-a 1,2 )x - 1 2 1-a 1 a 1 +(1-a 1 )x + 1-a 2 a 2 +(1-a 2 )x . As a 1,2 ≤ a i for i = 1, 2, we immediately check that 1-a i a i +(1-a i )x ≤ 1-a 1,2 a 1,2 +(1-a 1,2 )
x for i = 1, 2. Thus ψ(x) ≥ 0 on [-1, 1] and ϕ is increasing, which concludes the proof.

Remark 5.1.2. (relations between localization parameters) 1) Direct computation yields the upper bounds a 1,2 ≤ a 1 and a 1,2 ≤ a 2 . 2) Cauchy-Schwartz inequality provides a 1,2 ≤ √ a 1 a 2 .

3) Note that if each a i verifies 1/2 < a i < 1, the 'cross localization parameter' a 1,2 can be lower than 1/2. For example, we have a 1,2 = 1/3 for the choice of the overlapping step and triangular generic kernels.

Again, we see that these quantities, a 1 , a 2 , and a 1,2 , play a central role in analyzing the F-projection. We have now to take advantage of formula (??) to highlight a form of stability in the representation by fuzzy projection. But choosing two parameter a 1 and a 2 is not enough for computing a 1,2 = µ (1) |µ [START_REF] Crouzet | The fuzzy projection: a simple adaptive representation of signals[END_REF] , as the only relation between them is the Cauchy-Schwartz inequality a 1,2 ≤ √ a 1 a 2 : we have to choose first µ [START_REF] Boffi | Characterizing the inverses of block tridiagonal[END_REF] and µ [START_REF] Crouzet | The fuzzy projection: a simple adaptive representation of signals[END_REF] . The natural question is how to make this choice in a wide range of generic kernels that also let the different computations and interpretations possible. To do this, These power functions thus seem to fit in a simple and efficient way classical generic kernels. Moreover, a simple computation yields

a [α] := µ [α] 2 2 = 1 - 1 1 + α + 1 2(1 + 2α) (52)
and

a [α,β] := µ [α] |µ [β] = 1 - 1 2(1 + α) - 1 2(1 + β) + 1 2(1 + α β) (53)
Thus, it becomes possible to deal with

Match (Π 1 , Π 2 ) = 2a 1,2 -1 (2a 1 -1) (2a 2 -1)
using a wide class of kernels.

In Figure ??, we show the match computed this way. We can estimate that the zone indicating a good stability between two fuzzy partitions roughly starts when a 1 , a 2 ≥ 2 3 . This corresponds to the triangular kernel. Thus, as those represented in Figure ??, generic kernels µ [α] seem to lead to a good stability for α ≥ 1 (as a [START_REF] Boffi | Characterizing the inverses of block tridiagonal[END_REF] = 2/3). In other terms, changing the partition for the fuzzy projection does not imply an important deviation in the representation when we consider kernels near to the triangular one.

Also, we point out that the match formula is the same in any dimension of space q. According to equation(??), formula (??) becomes:

Proposition 5.2.1. (match in dimension q) Match (Π 1 , Π 2 ) = σ (p+1,...,p+1) = 2a 1,2 -1 (2a 1 -1) (2a 2 -1) q (54)
In particular, 2-dimensional kernels that should lead to a certain stability should be near the classical square pyramidal kernel. This stands in any dimension of space, according a favor in the linear interpolation, when considering that a fuzzy overlap between cells Ω k is necessary for a good description or representation of the uncertainty of the considered signal or multi-dimensional object.

Conclusion

Representation of signal and images by fuzzy and inverse fuzzy transforms has led to interesting numerical applications, as mentioned in the introduction ([? ], [? ], [? ], [? ], [? ], [? ], [? ]). Fuzzy projection itself, introduced in [? ],[? ], is a direct simple tool directly related to these techniques, and, in addition to noise removal, has the advantage to take better into account high frequency of signals.

In this study, we have shown that multi-dimensional fuzzy projection is a very simple mathematical tool for representing multi-dimensional objects, and is analogous to the 1D case for signals, so that it should present the same numerical advantages. For instance, our next work ([? ]) will be to use the 2D the fuzzy projection for simultaneous denoising and enhancement of very noisy images, in concurrence with classical and well-tested methods. A project is also to reconstruct 3D objects, taking benefit of the easy implementation of the F-projection in any dimension.

In any case, F-projection is in practice easier to handle than a wavelet or related transform representation, because it only demands the choice of a generic kernel and its associated partition. Moreover, we have highlighted the fact that a key parameter, the localization of the generic kernel in one dimension of space, a = µ 2 2 , permits to illustrate and control the representation in any dimension of space. It also governs its numerical robustness (condition number of the matrix playing the central role in the projection) and its stability -say its sensitivity when modifying the sampling partition. Moreover, as studied in [? ],[? ], F-projection offers an interesting flexibility, as soon as the kernel and the diameter of the partition can be locally modified in order to better represent some regions of the signal having different behaviors (more oscillating ones for instance).

Another track for a future work is to show that, in a more precise manner, or in a more specific context of application, the linear -triangular or pyramidal -fuzzy projection is the most simple and natural one. Indeed, we have highlighted the fact that the sensitivity of the representation with respect to modifications of its generic kernel is low for the choice of the linear kernel.

For all these reasons, F-projection should be a good numerical compromise between easy implementation, numerical robustness and stable representation in the requirements of fuzzy representations.

2. 2 .f 2 1/ 2 and

 222 Fuzzy transform Denote L 2 ([0, d]) the space of finite energy signals with norm f 2 = d 0 scalar product (f |g) = d 0 f g. Defined by Irina Perfilieva ([? ],[? ]) on the subspace C 0 ([0, d]), the more general fuzzy transform we consider in this study is the operator F : L 2 ([0, d]) : → R p+2 such that: Definition 2.2.1. (fuzzy transform F : sampling)

Theorem 2 . 3 . 3 .

 233 (numerical robustness and role of a) i)

Figure 1 :

 1 Figure 1: Example of function ν k for a triangular generic kernel µ = Tri [-1,1] (d = 10)

Theorem 3 . 3 . 2 .

 332 (expression of 2D F-projection and role of a) i) F F * corresponds to a (p + 2) 2 × (p + 2) 2 (almost bistochastic) matrix which is symmetrical block tridiagonal with symmetrical tridiagonal blocks:

Theorem 3 . 3 . 3 .

 333 (2D numerical robustness and role of a) i)

Proposition 4 . 1 .

 41 (qD numerical robustness and role of a) c(a) := cond 2 (F F * ) = 1 (2a -1) q

Theorem 5 . 1 . 1 .

 511 (match between subspaces and role of localization parameters) The (p + 2) × (p + 2) Gram matrix G related to Π 1 and Π 2 is diagonal, with diagonal terms

Figure 3 :

 3 Figure 3: Match between projection subspaces according the localization parameters a 1 and a 2 , using the power generic kernels µ [α] .

let us introduce a simple reasonable model for parametrizing kernels of fuzzy partitions. Remark 5.2.1. 1) µ [0] is the overlapping step kernel: [START_REF] Boffi | Characterizing the inverses of block tridiagonal[END_REF] is the triangular kernel: µ [START_REF] Boffi | Characterizing the inverses of block tridiagonal[END_REF]

) is very well approached by µ [1.75] , as µ [1.75] -µ [cos] 2 / µ [cos] 2 = 0.0077.


K C Shashidhar

Maurice Bruynooghe

Francky Catthoor
email: catthoor@imec.be

Gerda Janssens

Functional Equivalence Checking for Verification of Algebraic Transformations on Array-Intensive Source Code

Development of energy and performance-efficient embedded software is increasingly relying on application of complex transformations on the critical parts of the source code. Designers applying such nontrivial source code transformations are often faced with the problem of ensuring functional equivalence of the original and transformed programs. Currently they have to rely on incomplete and time-consuming simulation. Formal automatic verification of the transformed program against the original is instead desirable. This calls for equivalence checking tools similar to the ones available for comparing digital circuits. We present such a tool to compare array-intensive programs related through a combination of important global transformations like expression propagations, loop and algebraic transformations. When the transformed program fails to pass the equivalence check, the tool provides specific feedback on the possible locations of errors.

Introduction

Source code transformations come into play in situations where a designer wants much better optimizations than those a compiler can provide. Such a situation is common for designers of mobile computing and communicating systems. They are required to program complex signal processing algorithms for complex platform architectures and yet meet stringent constraints on the energy consumption and performance of the final implementation. Research has shown that application of source-to-source code transformations on an original implementation of the algorithm can greatly help in meeting such constraints (cf. [START_REF] Brandolese | The impact of source code transformations on software power and energy consumption[END_REF][START_REF] Catthoor | Custom Memory Management Methodology: Exploration of Memory Organization for Embedded Multimedia System Design[END_REF][START_REF] Catthoor | Data Access and Storage Management for Embedded Programmable Processors[END_REF][START_REF] Wolf | Memory system optimization of embedded software[END_REF]).

Every stage in an implementation activity brings forth an associated verification problem. Source code transformations are no exception. The problem here is to ensure that the transformed program preserves the functionality of the original program. Designers are at present using simulation of the transformed program to gain confidence in the correctness of the transformed program before forwarding it to the synthesis stage. But simulation is both incomplete and time-consuming. Also, when the transformed program is in error, it is hard to detect the exact cause with simulation. Clearly, formal automatic verification of the transformed program against the original, with support for error diagnostics, is desirable here. A pragmatic approach to this problem is to separate the two concerns, viz., applying transformations and verifying that they preserve the functional equivalence. This implies an ex post facto solution that requires a program equivalence checking tool. Our work addresses this requirement. Since, in general, the program equivalence problem is undecidable, we target the most important transformations applied on a decidable, and yet relevant, class of programs.

Code transformations considered. We are interested in verification of source code transformations that reduce the accesses to the data memory hierarchy. Broadly, there are two kinds of such transformations viz., global loop transformations and global data-flow transformations. Global loop transformations are applied to reorder and restructure the for-loops in the complete program in order to minimize the data transfers between different layers of the hierarchy by improving the temporal and spatial locality of the accessed data. On the other hand, global data-flow transformations are applied either to remove repeated computation or to break bottlenecks caused by data-flow dependencies in the program. They comprise of expression propagations that introduce or eliminate temporary variables that hold intermediate values and global algebraic transformations that take advantage of algebraic properties of the operators in transforming the dataflow. The need for verification support for these transformations is rather high because they invariably involve error prone manipulation of the index expressions of the array variables, especially when applied manually.

We do not distinguish between the transformations as long as they are only from the above categories. The transformed program can be under a combination of the transformations. The equivalence checking is done oblivious of any information about the particular instances of the above transformations that were applied and the order of their application.

An example problem. Suppose that we are given program functions, as in

k ∈ [0 ...N -1] : C[k] = B[2*k] + B[k] + A[2*k] + A[k]. Program (d) is erroneously obtained. It is inequivalent to them ∀ even k ∈ [0 ...N -1], where C[k] = A[k] + B[k] + A[k] + B[k], but equivalent ∀ odd k ∈ [0 ...N -1].
That is, the functions are input-output equivalent under the applied transformations. We have developed a tool to check such equivalences fully automatically. Because of an erroneous transformation, (d) is not equivalent to (a), (b) and (c). It is helpful if the reason for nonequivalence of the function can be ascertained and debugged. To this end, our tool provides diagnostics when it fails to show an equivalence.

Related Work

Motivated by pragmatics, we are interested in a fully automatic, push-button style, a posteriori solution. This precludes discussion of the vast research on formal verification of the transformation tool or the library of transformations.

Undecidability of the program equivalence problem enforces that any effort start by defining a decidable class of programs that is of interest. Hence, the problem has been addressed by various researchers for different program classes with different applications in mind. Without enumerating the methods, to the best of our knowledge, none of the methods is able to show equivalence of program functions as in our example, in a scalable way. The problem we address is distinct by its central requirement to represent and maintain the relationships among elements of the array variables in the programs in closed form. Unrolling deeply nested loops with large bounds is clearly infeasible for real-life signal processing programs. To add to this, algebraic transformations will require a prohibitive search for normalization on the unrolled statements. Hence, we restrict our discussion of related work to methods that do not propose unrolling of loops.

Translation validation [START_REF] Goldberg | Into the loops: Practical issues in translation validation for optimizing compilers. 3rd Compiler Optimization Meets Compiler Verification[END_REF][START_REF] Necula | Translation validation for an optimizing compiler[END_REF] and fractal symbolic analysis (FSA) [START_REF] Mateev | Fractal symbolic analysis[END_REF], both present methods which show semantic equivalence of two versions of programs. In the case of the former, the comparison is between the source and the target code. These methods are distinct from ours in that they essentially try to heuristically infer a sequence of legal transformations that can relate the two programs. Instead, we are able to directly check for equivalence of programs that are in a suitable language class. Also, their methods do not handle algebraic transformations. The work most related to ours, because we address the same class of programs, is the algorithm recognition method presented in [START_REF] Barthou | On the equivalence of two systems of affine recurrence equations[END_REF]. Again, algebraic data-flow transformations are not handled by them. Another distinction is that, all these methods do not stress on debugging support which is very important in the context of source code transformations.

With respect to the equivalence problem of algebraic expressions, we would like to point out that methods from symbolic algebra do not suffice due to the presence of loop transformations on array variables. The quest here is for an analysis that deals with algebraic transformations involving array variables in the expressions.

We have presented a preliminary method that checks equivalence under only expression propagations and loop transformations in [START_REF] Shashidhar | Automatic functional verification of memory oriented global source code transformations[END_REF]. In this paper, based on an improved program representation, we present an extended method that is able to handle algebraic data-flow transformations, in combination with the other two transformations, in a single pass.

Program Representation

In this section we present a representation that captures the computation and the relationships between elements of array variables in the program. The representation is possible only when the programs belong to a class that we describe first.

Class of allowed programs

The class of programs that we consider is based on the recurring features in signal processing programs, and the availability of some tools (mentioned below) to convert programs that are not in the class. We assume that a program is first subject to a source-to-source code preprocessing phase in order to ensure that the original program has the properties of this class before beginning to apply the code transformations. Restriction of programs to this class is primarily to ease the analyses required to identify optimization opportunities and apply transformations. Apparently, what eases transformations also eases their verification.

The following properties distinguish the class: ① Singleassignment form: Programs have been converted to a form (called the dynamic-single assignment form) where every memory location is written only once (methods exist to automate this, for example, [START_REF] Feautrier | Array expansion[END_REF]); ② Static control-flow: Datadependent while-loops have been converted to for-loops with worst-case bounds by moving the conditions inside the loops, and the data dependent if-conditions are simple enough to be handled by if-conversions; ③ Affine indices: All expressions in the index and the loop bounds are either affine or piece-wise affine, and ④ No pointer references: Pro-grams have been converted to a form where all references to the memory are with explicit indexing (for example, using a method as in [START_REF] Van Engelen | An efficient algorithm for pointer-to-array access conversion for compiling and optimizing DSP applications[END_REF]).

Array data dependence graphs

Given that a program has the above mentioned properties, we can extract the complete data-flow in the program. This extracted data-flow can be represented in the form of an array data dependence graph (ADDG). It is a directed graph, where, nodes represent the variables and occurrences of operators/functions in the program functions, and edges represent the data dependence (in the direction opposite to the flow of data). For example, Fig. 2 gives the ADDGs of the program functions in Fig. 1. The edges outgoing from operator nodes are labeled by the position of the operand for the computation and the thick edges outgoing from the variables are labeled by the labels of assignment statements where they are defined. The distinction when compared to a standard data dependence graph is that, in an ADDG, the data dependence, denoted by a reverse directed edge, refers not just to a single value, but to a set of values. Since the program is required to be in single-assignment form, the values are guaranteed to be assigned to different elements of the array variable defined in the statement. The dependence relation between the sets of values of defined variables and the operand variables is given by the so-called dependency mappings, one for each operand variable of the statement. For example, for the statement s2 in the original function (a), the two dependency mappings viz., from buf[] to the first A[] and from buf[] to the second A[] are as defined below, where D :

= {[k] | 1 ≤ k ≤ 1024 ∧ k ∈ Z}. M buf,A 1 := {[x] → [y] | x = 2k -2 ∧ y = 2k -2 ∧ k ∈ D} M buf,A 2 := {[x] → [y] | x = 2k -2 ∧ y = k -1 ∧ k ∈ D}.
When there are multiple occurrences of a variable in the same statement we distinguish them with subscripts denoting their position as operands. In an ADDG, the root nodes and the leaf nodes correspond to the output and the input variables, respectively, of the program function. The rest of the variable nodes in the ADDG are intermediate variables. Intermediate variable reduction. For a given path, reduction of an intermediate variable is an operation that we use as a primitive in our method. This involves updating the dependency mapping from the predecessor variable node to the intermediate variable node being reduced. The new dependency mapping will then be from the predecessor variable node to the successor variable node on the path in question. This is obtained by the composition of the two mappings.

For example, let us consider reduction of the intermediate variable node tmp on the leftmost path (path 1) in the ADDG of the original function (a) from the output variable C to the input variable B, that is,

[ C s3 -→ + 1 -→ tmp s1 -→ + 1 -→ B ].
The predecessor variable node to tmp is C and the successor variable node is B. Reducing tmp on the path involves updating the dependency mapping from C to tmp (M C,tmp ) into dependency mapping from C to B (M

C 1 B ). It is computed as below,
where is the natural join operation on two relations [START_REF] Kelly | The Omega Calculator and Library[END_REF] and

D := {[k] | 0 ≤ k < 1024 ∧ k ∈ Z}. M C 1 B := M C,tmp M tmp,B 1 := {[k] → [k] | k ∈ D} {[k] → [2k] | k ∈ D} := {[k] → [2k] | k ∈ D}.
On a given path, if all the intermediate variables between the current variable and the output variable are reduced, we obtain the output-current mapping. When the current variable is an input variable, this mapping gives the relation between the elements of the output variable to the elements of the input variable for that path. Then it is called the output-input mapping for that path of computation.

In the example path we mentioned above, tmp is the only intermediate variable, hence the output-input mapping from C to B for that path is just the M C 1 B that was computed above.

Algebraic Transformations

Algebraic data-flow transformations take advantage of the properties of the operators or user-defined functions and modify the data-flow such that the semantics of the original function are preserved (modulo overflow). The algebraic transformations are not restricted to the expression in a statement, but can have a global scope. This can be seen in our simple example in Fig. 1, where the algebraic transformations applied on (a) to obtain (c) and (d) are across expressions of multiple assignment statements at an algorithmic level. The ADDGs of the functions, as shown in Fig. 2, also reflect this.

Typically, most of such global transformations just rely on the associativity and/or commutativity properties of the operators like addition and multiplication on a fixed-point datatype like integer. Hence in what follows, we restrict our discussion to only these transformations. Other algebraic properties related to identity, inverse, distributivity and evaluation of constants are less common in practice and can be handled in a way similar to what we present.

The effect of such algebraic transformations on an ADDG is shown in Fig. 3, where, operators, ⊕ is associative, ⊗ is commutative and is both commutative and associative. Associativity. As shown in Fig. 3(a), the end-nodes are regrouped with respect to the chain of ⊕-nodes (the associative chain) by the associative transformation, while maintaining their order.

Commutativity. As shown in Fig. 3(b), the effect of a commutative transformation is to permute the positions of the outgoing edges of the ⊗-node.

Combination of associativity and commutativity. As shown in Fig. 3(c), the effect of the transformation is to maintain the same end-nodes with any possible tree of -nodes between them and the root -node. The designers, based on the knowledge of the overall computation, are often able to apply algorithmic or big-step transformations. For example, as can be seen from the ADDGs in Fig. 2, the transformed functions (c) and (d) are the result of several applications of the basic algebraic transformations above. Here, the transformation is motivated by the observation that they perform N/2 integer additions less when compared to functions (a) and (b) (i.e., 3N -5N/2).

Equivalence Checking Method

The method shows equivalence of the original and the transformed program functions by checking a sufficient condition. The condition is that, in the ADDGs extracted from the two program functions, every corresponding data dependence paths have identical ① computation and ② output-input mappings. The method checks these two parts of the sufficient condition based on a traversal of the two ADDGs.

We first briefly explain the basic method that is able to handle only expression propagations and loop transformations in Section 5.1. In Section 5.2, we extend it to handle algebraic transformations along with the other two transformations.

Basic method

Central to our equivalence checking is a synchronized depth-first traversal of the two ADDGs. The traversal begins from each of the corresponding root nodes (output variables with same names) and proceeds in lock-step on the two ADDGs. Initially the output-current mapping is set to identity mapping from an output variable to itself, for all its elements.

When an intermediate variable is encountered on the path in either of the two ADDGs, it is reduced and the output-current mappings are updated for each of the outgoing paths in the sub-ADDG rooted at the node being reduced. When an operator node is reached on one of the ADDGs, the same operator node must be reached next on the other ADDG, possibly after a sequence of intermediate variable reductions. Whenever multiple outgoing branches are present, the paths on the two ADDGs are paired for further traversal. If the branching is at an operator node, the paths with the same labels on the outgoing edges on either side are paired. If the branching is at an intermediate variable, the pairing is based on the outputcurrent mapping. At any given point during the traversal, the paths already traversed on the two ADDGs are both guaranteed to have the same operator nodes appearing in the same sequence on them. This implies that, when a path ends at a leaf node, the same computation is guaranteed on the corresponding paths traversed in the two ADDGs. This satisfies the first part of the sufficient condition.

When the output-current mappings are updated at the leaf nodes on the two corresponding paths we have their outputinput mappings. The second part of the condition is satisfied when they are checked to be identical. This implies that corresponding paths supply the same operators with the same values in any execution of the function. When the traversal has exhausted all the paths, satisfying the sufficient condition for each path, the two ADDGs, hence the two programs, are both guaranteed to apply the same computation on the same input values and hence assign the same output values.

For example, consider the ADDGs of functions (a) and (b) in Fig. 2, where (b) has been obtained by applying only expression propagations and loop transformations. There are 4 paths in total from the output variable C to the input variables in (a). But in (b), assignment to C is distributed among statements t3 and t4, as a result, it has 8 paths. For both (a) and (b), if we number the paths from left to right, the traversal corresponds, path 1 in (a) to paths 1 and 5 in (b), path 2 in (a) to paths 2 and 6 in (b) and so on. This satisfies first part of the sufficient condition since corresponding paths are found without any mismatch. It can also be checked that the updated output-input mappings on all pairs of corresponding paths are identical. For instance, for path 1 in (a) and (b), we have the following identical mappings.

a M C 1 B ⇔ b M C 1 B ⇔ {[k] → [2k] | 0 ≤ k < 512 ∧ k ∈ Z}
Note that since the branching at C in (b) divides the elements of C into two groups of assignments, output-input mapping is split for all paths of both (a) and (b). The remaining pairs of output-input mappings on the corresponding paths are similarly identical. This satisfies the second part of the sufficient condition.

Extended method

Let us now consider the situation when algebraic transformations are allowed. As discussed in Section 4, they may shuffle the paths and/or redistribute some operator nodes. Clearly, the traversal as described above will not suffice anymore. For example, consider the ADDGs (a) and (c) in Fig. 2, where (c) has been obtained by additionally applying algebraic transformations. In path 1 of the two ADDGs, a mismatch occurs upon reaching the input variable (leaf node) which is B in the original ADDG (a), whereas it is A in the transformed ADDG (c). The mismatch prevents the equivalence proof of program pairs under algebraic transformations.

In order to handle algebraic transformations, the traversal has to do additional work to know which paths to pair up. This requires that upon reaching an operator which permits algebraic transformations, a specific normal form be established before continuing the traversal. Such a normalization on two operations viz., flattening and matching, invoked depending on the properties that hold for the operator. When the operator is associative, flattening is invoked and when the operator is commutative, matching is invoked. When the operator is both associative and commutative, a flattening operation is followed by a matching operation.

Flattening for an associative operator. Suppose that an associative operator (⊕) is reached on an ADDG. The flattening operation involves a lookahead traversal of the sub-ADDG of the associative chain rooted at the ⊕-node. It constructs an ordered list of nodes such that each node in the list is either a leaf node or the first node for an outgoing path which is an operator node different from ⊕. Any intermediate variables that exist on the path between the nodes in the list and the root ⊕-node are reduced. The effect of flattening is that it brings all operands of the chain to the same level as successor nodes of the root ⊕-node. Fig. 4 illustrates this. The order in which the operands are reached during the traversal is maintained by labeling the edges accordingly. Matching for a commutative operator. Suppose that a commutative operator (⊗) is reached on both the ADDGs. Since any permutation of the outgoing edges of ⊗-nodes on the two ADDGs is a valid transformation, the labels on them are of no consequence in pairing the corresponding paths. The correct pairing is then provided by the matching operation. Separately on both the ADDGs, the sub-ADDG rooted at the ⊗-node is traversed and any intermediate variables that are present are recursively reduced until all the successor nodes of the ⊗-node are either leaf nodes or operator nodes. This yields two lists of successor nodes, one for each ADDG. If the lists have unique nodes, then the pairing of the paths is one-to-one. Otherwise, if the non-unique nodes are input variables, output-input mappings are checked to pair them. If the non-unique nodes are operators, a lookahead traversal is recursively employed to reveal more successor nodes until a unique pairing is obtained. Going back to our example function pair (a) and (c), we see that the first operator starting from the output variable is an addition operator. Therefore the flattening operation is applied at the node. This results in the ADDGs shown in Fig. 5. A subsequent matching operation here has to deal with the non-unique input variables as the successor nodes. Nodes B and A each appear twice as successor nodes on each of the ADDGs. This requires that the matching be established by checking the output-input mappings to the nodes. The check reveals that the following equalities hold, where

D := {[k] | 0 ≤ k < 1024 ∧ k ∈ Z}. a M C p B ⇔ c M C z B ⇔ {[k] → [2k] | k ∈ D} a M C q B ⇔ c M C x B ⇔ {[k] → [k] | k ∈ D} a M C r A ⇔ c M C y A ⇔ {[k] → [2k] | k ∈ D} a M C s A ⇔ c M C w A ⇔ {[k] → [k] | k ∈ D}
This results in the matching: {(p, z), (q, x), (r, y), (s, w)}.

Here, a leaf node has been reached now on each of the paths and there exists an identical output-input mapping on the corresponding paths of the two ADDGs. Our sufficient condition therefore implies that they are equivalent.

To summarize, with the help of the flattening and matching operations, the synchronized traversal can be continued with a correct pairing of the branching paths at the operators which permit algebraic transformations. On each path, the traversal culminates at the matching leaf nodes, at which point, the second part of the sufficient condition is checked. That is, the output-input mapping computed on the path to that input variable during the traversal must be identical to the one on the corresponding path on the other ADDG. If this check fails, the traversal stops, reporting a failure and generating diagnostics. If it succeeds, the traversal continues until all the paths of the ADDGs are exhausted.

Before concluding, a brief remark is in order related to presence of any cycles in an ADDG. A cycle implies that the data-flow has recurrences, that is, a set of statements that read values (possibly computed from values -) written by themselves in earlier iterations. The method efficiently deals with cycles by employing the computation of transitive closure of the total dependence mapping of the cycle. This is computable only under certain conditions that usually hold in most real-life programs that we have checked on.

Verification and Debugging

Our prototype transformation verification tool implements the scheme shown in Fig. 6. The tool accepts original and transformed functions in the C language. Our sufficient condition assumes that the code is correctly scheduled. Therefore, it is required to check separately that all the reads for values follow their writes, that is, the def-use order is correct in the two programs. This can be checked by standard array data-flow analysis (cf. [START_REF] Allen | Optimizing Compilers for Modern Architectures[END_REF]). If the order is correct, the tool extracts ADDGs based on a source code analysis. The equivalence checker takes these ADDGs and applies the method we have presented. It relies on the OMEGA calculator [START_REF] Kelly | The Omega Calculator and Library[END_REF] for an efficient implementation of all the required operations on integer sets and tuple relations. If the checking succeeds, the two program functions are guaranteed to be functionally equivalent. If it fails, it generates error diagnostics.

Error diagnostics

The simplest case of an error is when there is a mismatch in the expected operator node or the leaf node. The information about the mismatch and the statement number is output, which suffices to correct the error. The more complicated errors involve the index expression. For example, we earlier mentioned that function (d) is in error. Let us see the diagnostics generated for it. It is similar to (c), but has the assignment to C interleaved among statements v3 and v4. Therefore flattening ADDG of function (d) results in the same ADDG as shown in Fig. 5(c), for each of the two branches at C. As explained earlier, the presence of non-unique leaf nodes requires that the matching, between the flattened ADDGs (a) and (d), be based on the dependency mappings. The matching succeeds for two paths, viz., {(q, x), (s, w)}, but fails for the other two paths, viz., {(p, z), (r, y)}. This is shown below.

( a M C p B ⇔ {[x] → [2x]}) ( d M C z B ⇔ {[x] → [x]}) ( a M C r A ⇔ {[x] → [2x]}) ( d M C y A ⇔ {[x] → [x]}) where, x ∈ {[k] | (∃ j | 2 j = k ∧ 0 ≤ k < 1023 ∧ k ∈ Z)}.
The above mismatch in dependency mappings implies that paths z and y are in error, which correspond to statements v3 and v1 in the program text. The diagnostic points the user to these two statements, displays the index expressions of variables C, buf 2 , A and B in the statements as possible places of error and the difference in the output-input mappings. A further heuristic on this information deduces that variable buf 2 is common to the two paths and hence its index expression is most likely to be in error. This is indeed the case in statement v3 of function (d), where, it should have been buf[2*k].

When desired, the designer can also limit the checker to focus on only certain parts of the input programs. This can be done by specifying the subsets of output and input variables, or by declaring a correspondence of intermediate variables in the two programs. This helps not only in reducing the checking time but also in generating better error diagnostics.

Experience

The method is based on the depth-first traversal of the ADDGs and it uses tabling of established equivalences to avoid reworking on any overlapping sub-ADDGs. Therefore the complexity of traversal is linear in the size of the larger of the two ADDGs. Also, the supposedly expensive operations on the integer sets and tuples can be safely assumed to be bound by a small constant as the lengths of the formulae describing them are usually small enough in practice.

We have earlier reported verification times taken by our tool that implemented the basic method to be in the order of only few seconds [START_REF] Shashidhar | Automatic functional verification of memory oriented global source code transformations[END_REF]. In experiments with the new tool implementing the extended method, on realistic examples involving algebraic transformations, we have observed no significant degradation in performance. On problem instances that we experimented on, where we used source codes whose control complexity and ADDG sizes were comparable to reallife application kernels, verification consistently took less than 100 seconds on a desktop. This shows that our verification tool can be conveniently used to increase designer productivity while applying source code transformations.

Conclusions

We have presented a verification tool as required by designers applying source code transformations on signal processing programs. The tool is fully automatic, fast and able to provide useful error diagnostics. It is based on equivalence checking of the original and the transformed programs and is able to handle important transformations like expression propagations, loop and algebraic transformations that have been widely reported in the literature. Allowing algebraic transformations has significantly increased the class of programs that can be shown equivalent by the tool and hence its applicability in practice.

Fig. 1 .Fig. 1 :

 11 Fig. 1: Program functions (a), (b) and (c), equivalent under the considered transformations, compute ∀k ∈ [0 ...N -1] : C[k] = B[2*k] + B[k] + A[2*k] + A[k]. Program (d) is erroneously obtained. It is inequivalent to them ∀ even k ∈ [0 ...N -1], where C[k] = A[k] + B[k] + A[k] + B[k], but equivalent ∀ odd k ∈ [0 ...N -1].

Fig. 2 :

 2 Fig. 2: The ADDGs of the program functions in Fig. 1.

Fig. 3 :

 3 Fig. 3: Effect of algebraic transformations on an ADDG.

Fig. 4 :

 4 Fig. 4: Illustration of the flattening operation.

Fig. 5 :

 5 Fig. 5: The ADDGs of functions (a) and (c) after flattening.

Fig. 6 :

 6 Fig. 6: The verification and debugging scheme.

Proceedings of the Design, Automation and Test in Europe Conference and Exhibition (DATE'05)

Proceedings of the Design, Automation and Test in Europe Conference and Exhibition (DATE'05) 1530-1591/05 $ 20.00 IEEE

Acknowledgment. The authors gratefully thank anonymous referees for their detailed constructive comments.