Keywords: Reuse, Embedded Finite-Word-Length (FWL) Controllers, Intellectual Property (IP), Linear Time Invariant (LTI) Systems, Low-Power and Speed Optimization, Proportional-Integral-Derivative (PID)

ASIC or FPGA implementation of a finite wordlength PID controller requires a double expertise: in control system and hardware design. In this paper, we only focus on the hardware side of the problem. We show how to design configurable fixed-point PIDs to satisfy applications requiring minimal power consumption, or high control-rate, or both together. As multiply operation is the engine of PID, we experienced three algorithms: Booth, modified Booth, and a new recursive multi-bit multiplication algorithm. This later enables the construction of finely grained PID structures with bit-level and unit-time precision. Such a feature permits to tailor the PID to the desired performance and power budget. All PIDs are implemented at RTL level as technologyindependent reusable IP-cores. They are reconfigurable according to two compile-time constants: set-point word-length and latency. To make PID design easily reproducible, all necessary implementation details are provided and discussed.

I. BACKGROUND AND MOTIVATION

HE PID is by far the most commonly used feedback controller due to its simple structure and robust performance [START_REF] Åström | PID Controllers: Theory, Design, and Tuning[END_REF]. An important feature of this controller is that it does not require a precise analytical model of the system that is being controlled, which makes it very attractive for a large class of dynamic systems. While PID is well adapted for linear-time-invariant (LTI) systems [START_REF] Xue | Linear Feedback Control[END_REF], it stands powerless for non-LTI ones. Nevertheless some solutions exist, such as partitioning the non-LTI control algorithm into a linear portion and a non-linear portion [START_REF] Xiaoyin | A New Motion Control Hardware Architecture with FPGA based IC-Design for Robotic Manipulators[END_REF][4] [START_REF] Qu | Design and Implementation of Intelligent PID Controller based on FPGA[END_REF]. The linear portion represents the major control loop and is computed using an integrated PID, while the non-linear portion that acts as dynamic compensation to the linear one is performed in software using a general-purposemicroprocessor or a DSP.

In embedded control applications, such as in small-scale mobile robot, the control-loop-cycle is very tight and the power budget is very limited. A low sample rate leads to poor and degraded control-performance. And high power consumption shortens the battery lifetime. To cope with these two severe and antagonistic constraints, the need for both a high-speed and low-power PID structure is of utmost importance.

Today, design-reuse [START_REF] Keating | Reuse Methodology Manual for System on a Chip Designs[END_REF] is a well established design standard that allows grasping with rapid technology changes and increasing design complexity. It consists in the use of predesigned technology-independent, generic and reconfigurable IP-cores [START_REF]Reports of the International Technology Roadmap for Semiconductors (ITRS)[END_REF], most generally implemented at register-transfer-level (RTL).

However, at RTL abstraction level, no significant optimization results can be achieved if not undertaken at architectural and especially at algorithmic level. To achieve such a goal, a deep insight into PID arithmetic is necessary. At this stage, a choice of a numeric representation format is a crucial issue. Compared to floating-point, fixed-point format is the best candidate for optimized designs as it is much simpler to implement, faster, power-efficient and requires far much less hardware resources. However, the limited dynamic range can be source of control instability. This problem, referred to as finite-word-length (FWL) effect is an active research area that aims to shorten the floating-to-fixed point conversion time while preserving control performances [8][9]. The digital implementation of PID controllers went through several stages of evolution, initially dominated by the use of commercial-of-the-shelf (COTS) components and DSP. But over the past few years, FPGAs have brought a key advantage to digital control: the inherent parallelism of FPGA architecture allows many independent control loops to run at different deterministic rates without relying on shared resources that might slow down their responsiveness as in the case of COTS and DSP [START_REF] Gretlein | DSPs, Microprocessors and FPGAs in Control[END_REF] [START_REF] Manmasson | FPGA in Industrial Control Applications[END_REF].

A survey of recent PID related works can be classified into three categories. The biggest one includes works that are straightforward FPGA implementations targeting specific applications: DC-DC converter [START_REF] Chander | FPGA-based PID Controller for DC-DC Converter[END_REF], temperature control [START_REF] Yang | The IP Core Design of PID Controller based on SOPC[END_REF], motor multi-axis control [START_REF] Lazaro | Simulink/Modelsim Simulable VHDL PID Core for Industrial SoPC Multiaxis Controllers[END_REF], liquid level control [START_REF] Fons | Custom-Made Design of a Digital PID Control System[END_REF], and Xilinx versus Altera FPGA implementation for result comparison [START_REF] Sreenivasappa | Design and Implementation of FPGA based Low Power Digital PID Controllers[END_REF]. The second category proposes methodologies that analyze the FWL effect on PID controller in order to reduce the number of hardware resources [START_REF] Lima | A Methodology to Design FPGA-based PID Controllers[END_REF] [START_REF] Urriza | Word Length Selection Method based on Mixed Simulation for Digital PID Controllers Implemented in FPGA[END_REF]. And finally the third category, paradoxically the smallest one despite the large popularity of PID, comprises architecture-optimization works. In [START_REF] Zhao | FPGA Implementation of Closed-Loop Control Systems for Small-Scale Robot[END_REF] low-power serial and parallel multiple-channel PID architectures are proposed for small mobile robots. In this work, the optimization was carried out at macro-level considering several PIDs, rather than at micro-level (optimization of the PID itself). Nevertheless, the whole architecture will deliver much more interesting results if combined with an optimized PID. The second work [START_REF] Samet | A Digital PID Controller for Real-Time and Multi-Loop Control: a Comparative Study[END_REF] proposes serial, parallel, and mixed PID architectures incorporating different number (1-3) of multiplication cores. High power consumption, even with the serial architecture, and complex control-part are the two major shortcomings of this proposal. Finally, in [START_REF] Fong | Design and Implementation of Modular FPGA-Based PID Controllers[END_REF] an attractive optimized PID structure based on distributed arithmetic (DA) is presented. Although this latter exhibits interesting results in terms of resource utilization and power consumption, it suffers from three serious drawbacks: high latency (n+1 clock-cycles for n bit set-point word-length), FPGA technology-dependent as it's essentially based upon FPGA look-up-tables (LUTs), The objective of this paper is to design optimized FWL-PID structures that overcome all above-mentioned shortcomings, and which are especially dedicated to embedded control applications. The PID cores are described at RTL level. They are highly reconfigurable and technology-independent, offering the possibility to be mapped both on FPGA and ASIC.

To reach such a goal, a special focus was put on the optimization of the inner arithmetic of PID. For that, we considered two discrete forms of PID algorithm: the commercial form [START_REF] Wittenmark | Computer control: An overview[END_REF], called also the standard or ISA form, and the incremental form. These two forms went through three successive types of FPGA implementations, using: Booth multiplication algorithm (BMA) [START_REF] Booth | A Signed Binary Multiplication Te:chnique[END_REF], modified Booth multiplication algorithm (MBMA) [START_REF] Macsorley | High-Speed Arithmetic in Binary Computers[END_REF], and a new developed version called recursive multibit recoding multiplication algorithm (RMRMA) [START_REF] Oudjida | A New Recursive Multibit Recoding Algorithm for High-Speed and Low-Power Multiplier[END_REF]. Results show gradual improvements with clear superiority over those provided in [START_REF] Fong | Design and Implementation of Modular FPGA-Based PID Controllers[END_REF]. PID control-rate and energy-consumption savings are respectively as follows: 32% and 25% with BMA, 177% and 23% with MBMA, 431% and 20% with RMRMA.

Our previous paper [START_REF] Oudjida | High-Speed and Low-Power PID Structures for Embedded Applications[END_REF] introduced a limited designspace of PID. In this paper, we extended the design-space to accommodate different application cases and provided all necessary implementation details to make the design easily reproducible.

The paper is organized as follows. In this section we outlined the main requirement specifications for embedded PID controller. Section II introduces the two mostly-used discrete versions of PID algorithm. Section III, IV and V deal with BMA, MBMA and RMRMA implementations, respectively. A discussion around the obtained results is given in section VI. Section VII describes the verification method, while Section VIII shows how the FWL-effect is tackled. And finally some concluding remarks in Section XI.

II. THE TWO MOSTLY-USED DISCRETE VERSIONS OF PID

A typical closed-loop system using a PID controller is shown in Fig. 1, where u c (k), y(k), and u(k) are the discrete signal quantities at the k th sampling instant of the reference set-point, the process-feedback measured output, and the PID controller output, respectively.

In digital control, commercial and incremental forms are the two mostly-used discrete PID versions [START_REF] Åström | PID Controllers: Theory, Design, and Tuning[END_REF] [START_REF] Wittenmark | Computer control: An overview[END_REF]. They are denoted by recurrent equations ( 1) and ( 2), respectively, and their corresponding coefficients are grouped in Table I. Equations ( 1) and ( 2) are fully detailed in the Appendix.

( ) ( ) ( ) ( )

k D k I k P k u + + = (1)
Where ( ) ( ) ( )

k y B k u A k P c ⋅ + ⋅ = ;
( ) ( ) ( )

1 1 - ⋅ + - = k e C k I k I ; ( ) ( ) ( ) k f L k D H k D ⋅ + - ⋅ = 1 . With ( ) ( ) () 1 1 1 - - - = - k y k u k e c and ( ) ( ) ( ) 1 - - = k y k y k f And ( ) ( ) ( ) ( ) ( ) 2 1 1 - ⋅ + - ⋅ + ⋅ + - = k e C k e B k e A k u k u (2)
Where ( ) ( ) ( )

k y k u k e c - = ; ( ) ( ) () 1 1 1 - - - = - k y k u k e c ; ( ) ( ) () 2 2 2 - - - = - k y k u k e c .
To satisfy different application cases, two IP versions are developed for each equation: with constant coefficients and with varying coefficients (Fig. 2). This latter requires a host side interface (HSI) to handle the runtime change of the coefficients.

The commercial version allows the three standard PID functioning modes (P, PI, PID) according to Mode input value. At the end of u(k) computation, the Done output signal toggles during one clock cycle, and the PID enters into sleep mode (whole internal activity stopped except for clocking and HSI) for maximum energy conservation.

Coefficients

Commercial PID Incremental PID A straightforward parallel implementation of PID requires an amount of 7 adders/substractors and 5 multiplication cores for equation [START_REF] Åström | PID Controllers: Theory, Design, and Tuning[END_REF], and 4 adders/substractors and 3 multiplication cores for equation [START_REF] Xue | Linear Feedback Control[END_REF]. In digital hardware, the total gate count scales linearly with word length for an adder core, while it scales quadratically for a multiplier core. Thus, any effort for a low-power optimization of PID must be focused on the implementation of the multiply-and-accumulate (MAC) function (X.Y) [START_REF] Seo | A New VLSI Architecture of Parallel Multiplirer-Accumulator Based on Radix-2 Modified Booth Algorithm[END_REF]. In this work, the optimization effort is rather concentrated on the double MAC function (X.Y+T.Z) called DMAC, considered as the main building block of our PID structures. Equations ( 1) and ( 2) are partitioned accordingly.

A b K p ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ + + s d i s p T T T T K 1 B p K - ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ + - s d p T T K 2 1 C i s p T T K - s d p T T K H s d d T N T T + _ L s d d p T N T N T K + - _ K p is
For FWL-PID, two's complement fixed-point representation is used, which is habitually expressed in Q notation as Q ni.nf The values are coded in n i bits before the point (integer word length including 1 sign bit), and n f bits after the point (fractional word length). The total word length is n=n i + n f . Booth multiplication algorithm [START_REF] Booth | A Signed Binary Multiplication Te:chnique[END_REF] belongs to the class of recoding algorithm, i.e. algorithms that recode one of the two operands to cope with signed two's complement multiplication. Let Y be the multiplier:

j n j j n n y y Y 2 2 2 0 1 1 ∑ - = - - + - = (3)
Equation ( 3) can also be expressed as follows:

( ) ∑ ∑ - = - = - = - = 1 0 1 0 1 2 2 n j j j j n j j j Q y y Y (4)
Where

0 1 = - y and } { 1 , 0 , 1 - ∈ j Q
Consequently, the multiplier Y is divided into n slices, each of 2 bits. Each pair of two contiguous slices has one bit in common. Thus, the DMAC becomes:

( )( ) ∑ ∑ - = - = + = + 1 0 1 0 2 . 2 . . . n j j j n j j j T P X Q Z T Y X (5) [ ] j n j j j T P X Q 2 1 0 ∑ - = + = . . (6)
According to (5), Booth algorithm consists in recoding the multiplier Y into a set of ternary numbers

} { 1 , 0 , 1 -
in order to generate n simple partial products which are summed subsequently. Table II summarizes the 4 possibilities that may occur. The -X can be easily formed by adding 1 to the complement of X. A direct translation of DMAC equation (5) into architecture (Fig. 3) requires one extra adder and two registers in comparison with the optimized version (Fig. 4) based on (6), called ODMAC. Additionally, one clock cycle latency is also needed in Fig. 3. The control-part responsible of producing the successive couples (y j-1 , y j ) is insignificant: just one multiplexer driven by a counter.

Based upon ODMAC as the main building block, PID architectures are constructed for both incremental (Fig. 5) and commercial (Fig. 6) forms, and their implementation results (Table III) are respectively compared to those of [START_REF] Fong | Design and Implementation of Modular FPGA-Based PID Controllers[END_REF]. Comparison was made into identical conditions using the same FPGA device (Spartan XC2S50E-7FT256), although relatively old, as well as the same synthesis-tool version (Xilinx ISE 9.1i). In [START_REF] Fong | Design and Implementation of Modular FPGA-Based PID Controllers[END_REF], only a 16-bit word-length commercial version with constant coefficients (without HSI) is implemented. PID1 and PID3 exhibits interesting results: 44%, 25%, and 32% savings and 62%, 35%, and 38% savings in terms of gate count, power, and speed, respectively. PID3 exhibits higher savings but at the expense of control-quality. Latency is rather the same [START_REF] Lima | A Methodology to Design FPGA-based PID Controllers[END_REF], which is n+1 clock cycles for all designs (PIDX).

Optimizing latency without sacrificing the three other issues is the main objective of the next two sections. 3) can also be rewritten as follows [START_REF] Macsorley | High-Speed Arithmetic in Binary Computers[END_REF]:

y j-1 y j Reg X T z j-1 z j Y Z X.Y+T.Z + Cin + Cin Fig. 4. Optimized DMAC implementation T "0" T Mux "0" X "0" X Mux "0" << j (P j .T) (Q j .X) j = 0 , n-1 ODMAC n n n n 2n+1 y j-1 y j X T z j-1 z j Y Z + Cin + Cin Fig. 3. Straightforward DMAC implementation Reg X.Y+T.Z + Reg Reg T "0" T Mux "0" X "0" X Mux "0" << j << j (P j .T) (Q j .X) j = 0 , n-1 DMAC j = 0 , n-1 MAC j = 0 , n-1 MAC n n n n 2n+1 TABLE II BOOTH ALGORITHM Y j Y j-1 Operation 0 0 + 0 0 1 + X 1 0 -X 1 1 -0 IV. MBMA BASED PID Equation (
( ) ∑ ∑ - = - = + - = - + = 1 ) 2 / ( 0 2 2 1 ) 2 / ( 0 1 2 2 1 2 2 2 2 n j j j j n j j j j Q y y y Y (7)
Where

0 1 = - y and } { 2 , 1 , 0 , 1 , 2 - - ∈ j Q
In this case, the multiplier Y is divided into n/2 slices, each of 3 bits, with one bit overlapping between adjacent slices. The proof of equation ( 7) is given in [START_REF] Rubinfield | A Proof of the Modified Booth Algorithm for Multiplication[END_REF]. Thus, the DMAC equation becomes:

[ ] j n j j j T P X Q Z T Y X 2 1 ) 2 / ( 0 2 . . . . ∑ - = + = + (8)
Likewise, n/2 simple partial products are generated (Table IV). Since ODMAC is a reconfigurable RTL block, it is parameterized to suit equation [START_REF] Hilaire | A Unifying Framework for Finite Word Length Realizations[END_REF]. The new adapted ODMAC architecture is depicted in Fig. 7. The only difference is that Mux(8:1) are used instead of Mux(4:1), and (<<2.j) hardwired shifter instead (<<1.j). Compared to BMA based PID (Table V), MBMA based one (PID1) shows much more interesting results, since latency is divided by 2 while maintaining stable power consumption and speed. Control rate is drastically improved as its equal to maximum clock frequency divided by latency. As the discrete commercial form (equation 1) can accommodate the three functioning modes, implementation of PID2 produced the following power consumption values at 47 MHz: 268 mW, 313 mW, and 366 mW for P, PI, and PID functioning modes, respectively.

With regard of these improvements, one is encouraged to pursue farther [START_REF] Macsorley | High-Speed Arithmetic in Binary Computers[END_REF] in reducing latency by considering larger slices, such as:

( ) ∑ ∑ - = - = + + - = - + + = 1 ) 3 / ( 0 3 3 1 ) 3 / ( 0 2 3 2 1 3 3 1 3 2 2 2 . 2 n j j j j n j j j j j Q y y y y Y ( 9
)
Where

0 1 = - y and } { 4 , ... , 0 , ... , 4 - ∈ j Q
But in this case, some hard partial products are required such as 3X and -3X which are not easy to generate. How to circumvent this obstacle is the purpose of the next section.

V. RMRMA BASED PID

Multiplication is a fundamental operation in digital design. Its speed and power requirements are two critical factors limiting the whole system performances (PID in our case). Since the publication of Booth's algorithm in 1951, a huge number of improvement attempts were proposed, especially after the publication of a generalized version of MBA algorithm accompanied with its proof [START_REF] Sam | A Generalized Multibit Recoding of Two's Complement Binary Numbers and its Proof with Application in Multiplier Implementation[END_REF]. Most of the proposals aimed to reduce the number of partial products either by employing digital optimization techniques [START_REF] Lamberti | Reducing the Computation Time in (Short Bit-Width) Two's Complement Multiplier[END_REF][31] [START_REF] Kang | A Simple High-Speed Multiplier Design[END_REF] or by using larger slices (higher radices) [START_REF] Crookes | Using Signed Digit Arithmetic for Low-Power Multiplication[END_REF]. However, experience showed [START_REF] Seidel | Secondary Radix Recodings for Higher Radix Multipliers[END_REF] that beyond 4-bit slices (radix 8), the complexity to generate hard partial products can not be managed in a realistic way. In [START_REF] Seidel | Secondary Radix Recodings for Higher Radix Multipliers[END_REF], three metrics are provided for comparing the tradoffs when employing higher radix Booth recodings: partial product compression factor (gain), the number of hard multiples that must be precomputed (computation complexity), and partial product generation fanin (routing complexity).

To circumvent the problem of hard partial products in higher radices, the idea proposed in [START_REF] North | β-Bit Serial/Parallel Multipliers[END_REF] is to apply a recursive Booth recoding on the r-bit slice. While the idea is interesting, it relies upon a complicated mathematical formulation, leading to a complex control circuitry and especially to an exaggerated latency (2n/r).

According to the multibit recoding algorithm presented in [START_REF] Sam | A Generalized Multibit Recoding of Two's Complement Binary Numbers and its Proof with Application in Multiplier Implementation[END_REF], a n-bit two's complement operand Y can be written as:

( ∑ - = + + - ⋅ ⋅ ⋅ + + + + = 1 ) / ( 0 2 2 1 1 0 1 2 . 2 2
0 0 0 + 0 0 0 1 + X 0 1 0 + X 0 1 1 + 2X 1 0 0 -2X 1 0 1 -X 1 1 0 -X 1 1 1 -0
) ( + - + + - - + - + - + ) 2 2 ( 2 3 4 5 2 . 2 r r rj r rj r rj y y y ) ( rj r r rj r rj r rj y y y 2 2 . 2 ) 1 2 ( 2 1 2 3 ⎥ ⎦ ⎤ - + - - + - + - + (11) ( ) ( ) ( ) rj r n j r i i i rj i rj i rj y y y 2 2 2 1 0 1 2 0 2 2 1 2 2 1 ∑∑ - = - = + + + + - ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ - + = // . ( 12
) ( ) ( ) rj r n j r i i ji Q 2 2 1 0 1 2 0 2 ∑∑ - = - = ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ = // (13) With } { 2 , 1 , 0 , 1 , 2 - - ∈ ji Q
There is no need to prove equation [START_REF] Manmasson | FPGA in Industrial Control Applications[END_REF] since it is a combination of equations ( 10) and ( 7) which were already proven in [START_REF] Sam | A Generalized Multibit Recoding of Two's Complement Binary Numbers and its Proof with Application in Multiplier Implementation[END_REF] and [START_REF] Rubinfield | A Proof of the Modified Booth Algorithm for Multiplication[END_REF], respectively. The partitioning of operand Y according to equation ( 13) is illustrated by Fig. 8.

To avoid dealing with special cases, n and r must be chosen as even numbers, with r as a divider of n. Thus, the DMAC equation becomes:

( ) rj r n j i r i ji ji T P X Q Z T Y X 2 2 . . . . 1 ) / ( 0 2 1 ) 2 / ( 0 ∑∑ - = - = ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ + = + (14)
Depending on r value ranging from 2 to n, PIDs with various levels of parallelism and latencies (n/r+1) can be automatically generated with slight control complexity. The special cases of r=n and r=2 correspond to fully-parallel and fully-sequential PID, respectively. In between (r=4,n/2), partially-parallel PIDs are obtained. The outstanding advantage of this algorithm (equation 13) is that hard partial products are generated using simple ones (2X, X) only. For a simplified hardware and lower power consumption, the step-by-step sign-propagate technique is employed [START_REF] Henlin | A 16 bit x 16 bit Pipelined Multiplier Marcrocell[END_REF].

Obviously, equation [START_REF] Yang | The IP Core Design of PID Controller based on SOPC[END_REF] does not reduce the number of partial products, but allows a modulable space-time partitioning of the multibit recoding algorithm (equation 10), where n/r sets comprising each r/2 partial products can be generated and summed either simultaneously or iteratively. Whilst the parallel implementation of equation ( 13) allows an important reduction of the critical path (using a carry-save adder CSA), it requires too much space. Therefore, only the serial implementation is retained. In this case, latency drops from (n/2+1) to (n/r+1), whereas the overhead on the total critical path, which goes through log 2 (r/2) adder levels and which is equal to D in the case of MBMA, is slightly increased D+log 2 (r/2). Note that we are using a logarithmic summation tree and not a linear one (CSA like).

An illustrative serial example with r=4 is described as follows:

( ) j n j j j j j j y y y y y Y 4 1 ) 4 / ( 0 3 4 3 2 4 2 1 4 4 1 4 2 2 2 2 ∑ - = + + + - - + + + = (15) ( ) ( )
j n ji i i j i j i j y y y 4 1 4 0 1 0 2 2 1 4 2 4 2 1 4 2 2 2 ∑∑ - == + + + + - ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ - + = / . (16)
[ ] ( )

j n j j j Q Q 4 1 4 0 2 1 0 2 2 ∑ - = + = / (17) ( ) ( ) [ ] j n j j j j j T P X Q T P X Q Z T Y X 4 1 ) 4 / ( 0 2 1 1 0 0 2 2 . . ∑ - = + + + = + ( 18
)
The mapping of equation ( 18) into a serial architecture is shown in Fig. 9. Such a case (r=4) would have required the computation of hard partial products (7X, 5X, 3X) if the simple form of equation ( 15) was used. Notice that MBMA is a special case of RMRMA for r=2. For r=1, equation [START_REF] Gretlein | DSPs, Microprocessors and FPGAs in Control[END_REF] corresponds to BMA (equation 4).

Table VI comprises the implementation results of PIDs with n=16 and r=4,8,16. For instance, PID1 with r=4 not only achieves high improvement in latency (71%), but also maintains positive savings in power (14%) and speed (13%). These important achievements are partially due to logic-trimming performed by the synthesis tool on the constant coefficients. Such an operation is impossible in the case of PID [START_REF] Fong | Design and Implementation of Modular FPGA-Based PID Controllers[END_REF] since the coefficients are stored into LUTs.

X T Y Z z 2j-1 z 2j+1 z 2j j = 0 , (n/2)-1 ODMAC Reg X.Y+T.Z + Cin Mux T "0" T 2T 2T T T "0" Mux X "0" X 2X 2X X X "0" + Cin << 2j (Q j .
Y 16+1 bits Q j 8+1 bits Q ji 2+1 bits Q 01 Q 02 Q 03 Q 10 Q 11 Q 12 Q 13 : y -1 = 0
At this stage, a key question arises: among this panoply of PIDs, which one fits the best one's application case? The answer to this question is given in the next section.

VI. DISCUSSION

In embedded control, satisfactory control-rate (without performance degradation) at minimum power consumption is the main requirement. To select the most adequate PID for a given application, it's necessary to investigate how speed, power and hardware resources scales versus r factor for a fixed word length n. Referring to equation ( 14) and aided by Fig. 9, the ODMAC architecture scales as a binary tree with one stage of r mux(8:1) followed by Log 2 (r)+1 stages of adders with a total of r adders too. Thus, the total delay cumulated by the critical path which goes through Log 2 (r)+2 stages increases with O(Log(r)) complexity, whilst latency (n/r+1) decreases linearly O(r), which makes the maximum control-rate increases as r increases. This is confirmed by implementation results shown in Table VII and VIII corresponding to PID1 and PID2, respectively. The sole exception to this general rule is PIDX_n/2 which always yields to the highest control-rate compared to PIDX_n despite the numerous tests with various n values. This is justified since they exhibit very close latencies (3 and 2, respectively) and one stage difference in the critical path (n-1 and n, respectively), but an important multiplexer fanin difference (n/4 and n/2, respectively).

In terms of resource occupation, the total complexity grows linearly O(r) as r multiplexers and r adders are required by ODMAC which is the most resource consuming block of PID architecture. This is also confirmed by the implementation results shown in Table VI. Note that each adder of each level of MAC and ODMAC as well as the two ones at the output of the PID (Fig. 5 and6) are successively extended by one bit so that the total bit size of the control output u(k) becomes 2n+log 2 (r)+2. It's necessary to do so to prevent the apparition of a possible overflow in the data-path which can cause signal clipping and instabilities in the closed loop response [START_REF] Kelly | Design and Implementation of Digital Controllers for Smart Structures Using Field Programmable Gate Arrays[END_REF].

As for power consumption, intuitively, one would expect to see PID1_16 of Table VII as being the most rapid and the most power consumer too, for the reason that it exhibits the smallest latency and the biggest total gate count! While it is almost true for the latter (13 MHz, before the first), it is quite the opposite for the former (244 mW, the smallest one). The explanation is that power consumption (

clk sw dd F C V P 2 5 . 0 =
) depends linearly on the frequency (F clk ), which is in this case 26 MHz (the smallest one) and also on the switched capacitance (C sw ) which describes the average capacitance charged during each clock period (1/F clk ). In fact, C sw depends on a number of parameter (circuit structure, logic function, input pattern dependence…) and not only on the total gate count (more precisely, not only on the total physical capacitance of the circuit). Furthermore, a study [START_REF] Shang | Dynamic Power Consumption in Virtex-II FPGA Family[END_REF] that analyzed the dynamic power consumption in Xilinx's FPGA revealed the following share: 60% by routing, 16% by logic, and 14% by clocking. The reason is that routing is intensively segmented, using pass logic and buffers.

When both high control-rate close to 13MHz and low power are required, PID1_16 (244 mW at 13MHz) stands as the best candidate compared to PID1_8 (323 mW at 13MHz). However, it's noteworthy to mention that this comparison stands valid only for the special case of 16-bit word-length PID, for a given set of coefficients, mapped on XC2S150E-7FT256 FPGA circuit and using Xilinx's XST synthesis tool, version 9.2. Results could significantly change under other conditions, especially when considering the logic trimming process which is essentially dependant on Finally, when only low power is targeted, PIDX_1 is the best candidate. We dealt here with extreme situations only, but for a given couple (cr, pc) of control-rate and power consumption, several candidates are possible. Yet, the best PID is the one which requires the smallest gate count.

So far, speed and power have been considered in isolation to area which becomes critical, and sometimes prohibitive, for large word-length n due to the fact that PID is basically built of a set of multipliers (three or five) that scale quadratically with word length. The bigger is the area, the higher is the cost. Consequently, another advantage of RMRMA algorithm is to cope also with the cost issue as an additional constraint to speed and power.

We deliberately chose Spartan2e FPGA to compare our results with those provided in [START_REF] Fong | Design and Implementation of Modular FPGA-Based PID Controllers[END_REF]. A mapping on a recent FPGA circuit (Virtex6) using XST 12.1 version of extreme PID2 delivered state-of-the-art results grouped in Table IX. Note that control-rate scaled with an average factor of 2, while power dissipation scaled with an average factor of 45. This is not surprising, since Spartan2e and Virtex6 were fabricated with two differently scaled technology processes: 150 nm and 40 nm, respectively. Therefore, the physical capacitances of the circuit in Virtex6 are relatively too much smaller. Additionally, the supply-voltages (V dd ) used for internal core (Vccint) and for output blocks (Vcco) are respectively 1.8V and 3.3V for Spartan2e, 1V and 2.5V for Virtex6. Furthermore, the efficient advances made in CAD tools (from Xilinx ISE 9.1 to 12.1 versions) as well as in FPGA architecture, such as advanced segmented-routing, much contributed to lower the power consumption [START_REF]Virtex6 FPGA: Satisfying the Insatiable Demand for Higher Bandwidth[END_REF]. Power consumption evaluation studies [START_REF] Shang | Dynamic Power Consumption in Virtex-II FPGA Family[END_REF][39] based on simulation and measurements, targeting Virtex2 and Virtex6 families revealed the following results: 5.9µW per CLB per MHz, and 1.09 mW per 100 MHz at 38% toggle rate, respectively. These studies roughly confirm our power results as proximate values are obtained.

Timing and power evaluations were performed in the following conditions. Delays were calculated for two types of paths: Clock-To-Setup and all paths together (Pad-To-Setup, Clock-To-Pad and Pad-To-Pad.) The Clock-To-Setup gives more precise information on the delays than other remaining paths, which depend in fact on I/O Block (IOB) configuration (low/high fanout, CMOS, TTL, LVDS…). Thus, all delays (frequencies) presented so far are clock-to-setup delays with the highest speed grade of the FPGA circuit. As for power, we chose the highest Vcco voltage value (3.3 for Spartan2e and 2.5 for Virex6) with a maximum toggle activity of 50%, which means that Flip-Flops (FFs) toggle one time during each clock cycle. The reason is that only simple-edge triggered FFs are used for synthesis (no double-edge FFs).

VII. VERIFICATION METHOD

The PID design verification process went through several steps. First equations ( 12) and ( 14) were tested with a random C-program. Then, a severe cycle-accurate functional verification procedure using Modelsim simulator was applied to MAC and ODMAC as they are the main building blocks of PID architecture. They were challenged against a set of special test cases (visual simulation), and then submitted to a random test for a very large number of vectors. Once tested successfully, the RTL PID module written in Verilog-2001 (IEEE 1364) was integrated into Modelsim/Simulink environment for a co-simulation. At this stage, a ZOH discrete time invariant model of a third order continuous process (G(s)=1/(s+1) 3 ) was chosen from the test set used by Åström and Hägglund [START_REF] Åström | PID Controllers: Theory, Design, and Tuning[END_REF] as examples of representative plants for the dynamics of typical industrial processes. To derive the PID parameters, a theoretical PID taken from Matlab component-library was tuned using floating-point numerical representation (IEEE 754 double format). The sampling period T s was chosen based on the magnitude of the pole time constants. For this case T s =10 ms. The following parameters were obtained: K p = 0.5913 ; T i = 0.0523 ; T d = 0.0225 for N=10 and b=1. Calculations give the following floating-point values for the coefficients of commercial PID: A=0.5913; B=-0.5913; C=0.1130; D=0.1836; E=-1.0860 To co-simulate the RTL PID, a conversion of the coefficients to 16-bit (Q 4.12 ) fixed-point representation was necessary. Variations were obtained: A=0.5911; B=-0.5911; C=0.1130; D=0.1836; E=-1.0860 Note that to represent the original parameters with fullprecision, 44 bits are needed for the fractional part. Varied simulations were performed to verify the correctness of the PID RTL code. First, to explore the precision effect on control quality, the control output of PIDs with various fractional-part sizes (Q 4.4 , Q 4.12 , Q 4.20 ) were compared to that of the Matlab floating-point PID component (Fig. 10). Simulation shows different rise-times for different precisions. The higher is the precision; the closer is the control output from the ideal model. The second simulation tests the behavior of the PID after having reached the steady state (Fig. 11). For that, two perturbations are successively exerted on control output and on the plant measure. Each time the system recovers as expected. And finally, the third simulation investigates the PID capabilities to track setpoints of arbitrary amplitudes and durations (Fig. 12). After a successful functional verification, the RTL code of PID was synthesized, placed, and routed on Xilinx's FPGA (Virtex-2). The three preceding co-simulations but with timing backannotation were performed again as a last necessary software verification step before hardware integration of the PID into an FPGA evaluation board (MEMEC V2MB1000).

Finally, as an ultimate validation step, a physical test of our PIDs is performed. We built up a classical temperature control setup (Fig. 13 and14), which consists in a tube comprising a halogen lamp (12 V, 21 W), a temperature sensor (LM35), and a DC Fan (12 V, 1.68 W). Temperature regulation inside the tube is achieved by controlling either the intensity of the lamp, or the rotation speed of the fan. This is carried out by the use of two PWMs, whose dutycycle ratios represent the PID controller output (u(k)). These two PWMs do not act directly on the fan or on the lamp but rather on transistors (IRF540) that control the power consumed by the lamp and fan.

The sensing of the actual temperature of the tube is assured by LM35 component which delivers a voltage value that grows linearly with temperature (1.5 volts corresponds to 150 °C). As the maximum voltage allowed by FPGA evaluation board (V2MB1000) is 3.3 volts, the calculation of the real temperature (T) is done as follows: T = [(val_opb_ADC * 3.3)/1023] * 100. This allows a temperature control with a minimum step of 0.32 °C.

The V2MB1000 board is connected through RS232 port to a PC running a .net application which allows a real-time display of the temperature as well as an instantaneous tuning of the set-point.

VIII. THE FINITE WORD-LENGTH (FWL) EFFECT

Fixed-point arithmetic is employed as an approximation of real numbers (floating-point), with a fixed bit-length of the word used to represent data (Finite Word-Length). This limitation leads to performance degradation (FWL effect) mainly due to quantization of coefficients (parametric errors) and roundoff errors subsequently cumulated during the computation process (numeric noise). In fact, the FWL effect is more-or-less exaggerated depending on the control algorithm used (I/O relationship, levels of parallelism, etc) as well as on the way the computations are performed (number of bits, different/unique fixed point position, round/truncation, etc). Compared to the reference floatingpoint implementation, the FWL effect can be assessed using some indicators such as transfer function sensitivity, or pole sensitivity [START_REF] Gevers | Parametrizations in Control, Estimation and Filtering Probems[END_REF][41] [START_REF] Lopez | Sum-of-products Evaluation Schemes with Fixed-Point arithmetic, and their application to IIR filter implementation[END_REF].

In fact, the objective is twofold: we need to provide an optimal ASIC/FPGA implementation of FWL PID without degrading control performances. To achieve such a goal, a double expertise is required in hardware design and control system. But usually, hardware designers do not master 11.

12.

control system design, and control system experts do not have the required skills to implement and evaluate the controllers using ASIC/FPGAs [17][43]. This is why we propose, as hardware designers, a highly reconfigurable (n, r) and technology-independent FWL PID that can systematically respond to control-engineer demands after having modelled, simulated, and evaluated the performances provided by different bit-width fixed-point representations using Matlab/Simulink environment, and finally opted for an appropriate word-length (n) of the setpoint. As for latency value (r), it depends on the application domain and intended objectives. Precise guidelines on how to choose r value were given in section VI. Now that (n, r) couple is known, the FWL problem is tackled from hardware side by simply adjusting in the RTL code the two compile-time constants: setpoint bit-size (n) and latency (r). The synthesis of such a PID generates an optimal structure that not only meets the performances specified by control-engineers, but also consumes minimum power and hardware resources. This would not have been possible without the use of the new highly serialisable multi-bit multiplication algorithm (equation 13). The incorporation of equation ( 13) [START_REF] Oudjida | A New Recursive Multibit Recoding Algorithm for High-Speed and Low-Power Multiplier[END_REF] into equations ( 1) and ( 2) as an efficient PID engine, allows the generation of PID architectures classified as regular iterative architectures (RIA) [START_REF] Rao | Regular Iterative Algorithms and their Implementation on Processor Arrays[END_REF], known for their high conformity with the principles of regularity and locality. In addition to equation ( 13), we propose in [START_REF] Oudjida | A New Recursive Multibit Recoding Algorithm for High-Speed and Low-Power Multiplier[END_REF] several new highly serialisable multiplication algorithms, offering different features in power, space and delay, depending on the operand size (n). Reader is encouraged to explore these algorithms [START_REF] Oudjida | A New Recursive Multibit Recoding Algorithm for High-Speed and Low-Power Multiplier[END_REF] to select the appropriate one that leads to best performances of its controller with regard to the size (n) of the setpoint.

Regularity and locality are two important features, highly sought in hardware design as they lead to an important gain in space and delay. Regularity is a general space feature, where the repetitiveness of just one or few elementary building-blocks (mux, adders and shifters of ODMAC, Fig. 9) and their interconnection scheme (predefined netlist) suffice to draw the whole architecture (MAC/ODMAC and then PID). In the other hand, locality is both space and time feature, in the sense where each building-block can only interact with its nearest surrounding neighbours, and any transaction from one building-block to the next is completed in one and only one unit time delay (clock period). Because of these two important features, our PIDs can be finely grained at bit level in space (setpoint bit-size n, latency r) and unit delay in time (latency r).

Experimental results depicted in Fig. 15 illustrate the FWL effects on temperature regulation. Reducing the fractional-part size of the set-point beyond a certain limit (4 bits) yields to a continuous fluctuation of the temperature inside the tube (Fig. 15.d). The best compromise is a 6-bit fractional-part (Fig. 15.c) which ensures a correct regulation while consuming less power and hardware resources. As temperature regulation system has a very slow dynamic, speed is not a concern. Therefore, the most appropriate PID in this case is PIDX_1 as it is the least power consumer. Adversely, for very fast dynamic systems, such as MEMS [START_REF] Hoover | Towards Understanding Architectural Tradeoffs in Mems Closed-Loop Feedback Control[END_REF] or microrobotics applications [START_REF] Casanova | Integartion of the Control Electronics for a mm 3sized Autonomous Microrobot into a Single Chip[END_REF], PIDX_n/2 is the most adequate option as it leads to the highest control rate.

IX. SUMMARY AND CONCLUSION

Despite the large popularity of PID controller, little attention has been paid to its optimization, either for ASIC or for FPGA integration. To break down this paradoxical situation, a series of high-speed and low-power PIDs, especially dedicated to embedded applications was proposed. They are based on two discrete forms of PID algorithm: the incremental form and the commercial form, both with constant and time-varying coefficients. The work focused more particularly on the commercial form with varying coefficients as it is the most used in industry due to the higher control-quality provided. Two types of optimizations were carried out: architectural and algorithmic optimizations. The former is a macro-level optimization, based on an efficient partitioning of PID discrete-equations, considering the double MAC (DMAC=XY+ZT) as the main building block of PID architecture. An optimized version of DMAC was developed (ODMAC) for less hardware resource occupation. As for the micro-level optimization (inner optimization of ODMAC), three multiplication algorithms were experienced: BMA, MBMA, and a new general and recursive version of MBMA called RMRMA. In addition, some low-power design techniques were incorporated, such as: sleep mode, and step-by-step sign-propagation technique.

The implementation results of PID based upon these three algorithms yielded to gradual improvements with a clear superiority over results presented in [START_REF] Fong | Design and Implementation of Modular FPGA-Based PID Controllers[END_REF]. For instance, concerning PID1_2 and PID1_4, savings of 177%, 23%, and 36%, and savings of 284%, 14%, and 26% are obtained in control-rate, power consumption, and total gate count, respectively. Additionally, analytical scaling-complexity evaluations with respect to the couple (n,r), confirmed also by software simulations, revealed useful information which is summarized as follows:

• PIDX_n/2 is the fastest PID that yields to the highest control-rate (30 MHz for PID2_8 mapped on Virtex6, with (n,r)=(16,8) ); • PIDX_1 is the most power efficient PID when speed is not a concern; • PIDX_n and PIDX_n/2 are the most efficient PIDs when both high control-rate and low-power dissipation are required. Further extension to the present work is to apply the same (or appropriate) partitioning in conjunction with RMRMA algorithm to the set of recurrent equations of an arbitrary number of multi-loop PID controllers taken as a whole.

Finally, the new recursive multiplication algorithm (RMRMA), well adapted to large word-lengths, and which was behind the drastic optimization of PID, can be efficiently applied to a variety of advanced control algorithms such as to linear-quadratic-gaussian (LQG) or sliding-mode controllers, etc.

Incremental form

The standard version of PID controller is described in a differential equation as:

() () () () ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ + ⋅ + = ∫ ⋅ t d i p dt t de T d e T t e K t u 0 1 τ τ
, where e is the system error ( () () ()

t y t u t e c - =
), u c is the command signal (setpoint), y is the process variable (measured variable). K p is the proportional gain, T i the integration time constant, and T d the derivative time constant of the controller.

Using Laplace transform, () t u is expressed in s-domain by: ()

() ( ) () ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ⋅ ⋅ + ⋅ + = s E T s T s s E s E K s U d i p .
For a small sample interval T s , the continuous time variable ( ) t u can be discretized using the following approximations: . k denotes the k th sampling instant (k.T s ). Thus, () t u can be rewritten as:

2 - ⋅ + k e C
This latter equation is called the incremental form of the controller. A drawback with the incremental algorithm is that it cannot be used for P or PD controllers.

Commercial form

For better performances of PID, two corrections are performed: limitation of the derivative gain and setpoint weighting. A pure derivative action will induce a very large amplification of measurement noise. The gain of the derivative must thus be limited. This can be done by approximating the transfer function s.T d as follows: Consequently, the discretized form of () t u d is: () ( ) ( ) ( ) ()

s d p s d T k Y k Y T K T k D k D N T k D 1 1 - - - - - ⋅ - = .
After simplification, we obtain: () ( ) () ( ) ( )

1 1 - - ⋅ + ⋅ ⋅ - - ⋅ + = k Y k Y T N T T N K k D T N T T k D

( 1 )Fig. 1 .

 11 Fig.1. Typical closed-loop control system using a PID

Fig. 2 .

 2 Fig. 2. Various PID IP-cores. (a) commercial PID with constant coefficients; (b) commercial PID with time varying coefficients; (c) incremental PID with constant coefficients; (b) incremental PID with time varying coefficients;

Q

  In this general case, the multiplier Y is divided into n/r slices, each of r+1 bits. Each pair of two contiguous slices

17 *:

 17 Fig. 6. Commercial PID architecture

Fig. 7 .

 7 Fig. 7. Optimized DMAC architecture for r=2

Fig. 9 .

 9 Fig. 9. Optimized DMAC architecture for r=4

Fig. 13 .Fig. 9 .Fig. 10 .Fig. 11 .

 1391011 Fig. 13. Synoptic scheme of the setup

Fig. 15 .

 15 Effect of the setpoint fractional length on temperature regulation (a) Floating point PID; (b) Our PID with Q ni.nf = Q 8.8 ; (c) Our PID with Q ni.nf = Q 8.6 ; (d) Our PID with Q ni.nf = Q 8.4

  After simplifications, we get the following recurrent equation:

1 , 1 U 1 . 1 .

 1111 where N is typically in the range of 3 to 20. In addition, to avoid sudden overshoots due to high variations of the setpoint, only a fraction b of u c acts on the proportional part (b.u cy). Hence, the improved PID algorithm becomes: (s) expression is discretized such that the proportional, integral and derivative terms are separately obtained, as follows: derivative term () k D, we use the differential equation representing the transfer function of By performing cross products, we get:Applying the inverse Laplace Transform to this latter equation, we obtain:

  -1 y 0 y 1 y 2 y 3 y 4 y 5 y 6 y 7 y 8 y 9 y 10 y 11 y 12 y 13 y 14 y 15

	n	n	n	n
	y 2j+1 y 2j y 2j-1
	X)		(P j .T)
		2n+1
	Q 00
			Q 1
		Y
							TABLE VI
					IMPLEMENTATION RESULT COMPARISON OF RMRMA-BASED PID
	Fig. 8. Partitioning of a 16-bit Y operand with r=8	PID Core PID [21]	Total Gate Count 16728	Power* (mW) 223	Max. Clock Freq. (MHz) 47	Latency 17
					PID1_4	12443 (+26%)	191 (+14%)	53 (+13%)	5 (+71%)
					PID1_8	15688 (+06%)	194 (+13%)	44 (-06%)	3 (+82%)
					PID1_16	23545 (-41%)	217 (+03%)	26 (-45%)	2 (+88%)
					PID2_4	22962 (-37%)	256 (-15%)	43 (-08%)	5 (+71%)
					PID2_8	26073 (-56%)	204 (+08%)	37 (-21%)	3 (+82%)
					PID2_16	40327 (-141%)	488 (-119%)	23 (-51%)	2 (+88%)

*: Dynamic power consumption at 23MHz; PIDY_X: X = r (+AB%): saving; (-AB%): overhead Q 0 y

  arrangement of the coefficients. For a minimum influence of the trimming operation on the synthesized results, appropriate coefficients were used such as all Qj terms are represented except the null one to avoid generating null partial products that greatly simplify the circuit logic. In fact, constant coefficients PIDs (PID1) are somehow unpredictable with regard to r. They are coefficient dependant. Adversely, PID2 is not involved with the trimming process since coefficients are time varying. Implementation results comprised in TableVIIIshow that PID2_8 is the best at all aspects for the same reasons cited above. In sum, when high control-rate is the ultimate objective, PIDX_n/2 is the best candidate whatever n value. But in the case where both high speed and low power are required, timing and power evaluations are necessary to decide which PID to select: either PIDX_n/2 or PIDX_n.

	the bit-
			TABLE VII
	MAXIMUM POWER-CONSUMPTION AND CONTROL-LOOP-CYCLE OF PID1
	PID Core	Power* (mW)	Max. Clock Freq. (MHz)	Latency	Max. Control Loop Cycle (MHz)
	PID [21]	456	47	17	2.76
	PID1_1	342 (+25%)	62	17	3.65 (+32%)
	PID1_2	350 (+23%)	62	9	7.66 (+177%)
	PID1_4	431 (+05%)	53	5	10.60 (+284%)
	PID1_8	365 (+20%)	44	3	14.67 (+431%)
	PID1_16	244 (+46%)	26	2	13.00 (+371%)
	*: Dynamic power consumption at maximum clock frequency;
	PID1_X: X=r; Max. control loop cycle = Max. clock frequency / Latency


Naiyu Zhang
email: naiyu.zhang@utbm.fr

Jean-Charles Creput
email: jean-charles.creput@utbm.fr

Hongjian Wang
email: hongjian.wang@utbm.fr

Cyril Meurie
email: cyril.meurie@utbm.fr

Yassine Ruichek
email: yassine.ruichek@utbm.fr

Partial Demosaicing for Stereo Matching of CFA Images on GPU and CPU

Keywords: GPU, Demaosaicing Stereovision, CFA Image, CUDA

This paper presents a GPU implementation of a partial demosaicing scheme that is specially designed for stereo matching of CFA image. This method consists of three main techniques keys: the adapted matching cost for CFA image, the estimated Second color component based on Hamilton's estimate method and a robust cost aggregation window. Experiments are carried out to explore the performance for this method on GPU both at matching quality and matching efficiency, with comparison with version on CPU. The experiments on different size image pairs from Middlebury dataset show that this method can be substantially accelerated on GPU when the image size is large and has still space for improvements in performance.

I. INTRODUCTION

In most implementation of stereo matching, such as in intelligent cars and integrated robot systems, color image pairs can be acquired by two types of cameras: the one equipped with three sensors associated with beam splitters and color filters providing the so-called full color images, of which each pixel is characterized in Red, Green and Blue levels, and the one equipped only with a single-sensor.

In the second case, the single-sensor can not provide a full color image directly but actually deliver a color filter array (CFA) image, of which every pixel is characterized by a single color component that can be one of the three color components: Red, Green and Blue. So, the missing color components has to be estimated at each pixel. This process of estimating the missing color components is usually referred to as CFA demosaicing and produces a demosaicing color image where every pixel is represented by an estimated color point [START_REF] Battiato | Recent patents on color demosaicing[END_REF]. This estimation step brings some artifacts in color density values, so, it is interesting to find one stereo matching directly based on CFA image.

As the demosaicing methods intend to produce demosaiced color images, they attempt to reduce the presence of color artifacts, such as the false colors ot zipper effects, by filtering the images [START_REF] Yang | Quality evaluation of color demosaicing according to image resolution[END_REF]. So, some useful color information for stereo matching may be lost in the color demosaiced images. As a result, the demosaiced color image pairs stereo matching quality usually suffer either from color artifacts or from the alteration of color texture caused by demosaicing schemes.

The method that is used in this paper is an alternative solution to match pixels by analyzing directly the CFA images without reconstructing the full color image by demosaicing processing, but only estimating a second color component. The method is a local dense stereo matching method that was first proposed by Halawana [START_REF] Halawana | Partial demosaicing of CFA images for stereo matching[END_REF]. Here, we study its potential for Graphics Processing Unit (GPU) implementation.

The following contents are organized in three parts. The Section II presents the partial demosaicing matching method including the estimation of the second color component and the adapted matching cost. The Section III describes the experiment platform and CUDA implantation. The Section IV details the experimental results. At the end is a conclusion.

II. PARTIAL DEMOSAICING MATCHING METHOD A. Second Color Component

The partial demosaicing matching method starts with the mosaiced CFA images. Here, the CFA images are those images obtained according to the Bayer color filter, each two-by-two submosaic contains 2 green, 1 blue and 1 red filter, each covering one pixel sensor and the mosaiced CFA image is the one whose pixel contains only one color component according to the Bayer color filter.

Different from the classic methods, which estimate all the missing color components for every pixel in the CFA images, the partial demosaicing method estimates only one color component, the Second Color Component (SCC), for every pixel. Here, the SCC is defined as the color component that is available in the same line. This means that SCC is the green color for all the red and blue pixels while for the green pixels the SCC is the red color component for even lines and the blue color component for odd lines. Summarized as [START_REF] Battiato | Recent patents on color demosaicing[END_REF].

SCC(x, y) =

   Ĝ(x, y) for red and blue pixels R(x, y) for green pixels in even lines B(x, y) for green pixels in odd lines (1)

B. Hamilton's Estimate Method

This method is an edge-adapted demosaicing method presented by Hamilton and Adams [START_REF] Hamilton | Adaptive color plan interpolation in single sensor color electronic camera[END_REF]. To select the interpolation direction, this method takes into account both gradient and Laplacian second-order values by using the green levels available at nearby pixels and red (or blue) samples located two apart.

We take the case GRG, as illustrated in Fig. 1, as example, to estimate the missing green level at the red pixels, this method uses the following algorithm: a) Approximate the horizontal ∆ x and vertical ∆ y gradients thanks to absolute differences as [START_REF] Yang | Quality evaluation of color demosaicing according to image resolution[END_REF].

b) Interpolate the green level as (3).

Fig. 1: Estimation of SCC in case of GRG.

∆ x = |G -1,0 -G 1,0 | + |2R -R -2,0 -R 2,0 | ∆ y = |G 0,-1 -G 0,1 | + |2R -R 0,-2 -R 0,2 | (2) Ĝ =          G-1,0+G1,0 2 + 2R-R-2,0-R2,0 4 if ∆ x < ∆ y G0,-1+G0,1 2
+ 2R-R0,-2-R0,2 4 if ∆ x > ∆ y G-1,0+G1,0+G0,-1+G0,1 4 + 4R-R-2,0-R2,0-R0,-2-R0,2 8 if ∆ x = ∆ y (3)
Since this method well combines two color component data in partial derivate approximations by exploiting spectral correlation in the green plane estimation, the precision is well guaranteed.

Each pixel with coordinates (x, y) in the partially demosaiced color images is characterized by a two-dimensional partial color denoted ÎP A . As shown in (4), this partial color point is composed of the available color component and the estimated second color component.

ÎP A (x, y) =        (R(x, y), Ĝ(x, y)) T if x is odd and y is even ( R(x, y), G(x, y)) T if x is even and y is even ( Ĝ(x, y), B(x, y)) T if x is even and y is odd (G(x, y), B(x, y)) T if x is odd and y is odd (4)

C. Adapted Matching Cost

The method is a local dense stereo matching method also called window-based approach. It respects the very assumption that the color information of neighbors of a left pixel is close to those of the same neighbors of its homologous right pixel in the right image. So, the matching costs are defined between the window around the left pixel and the window around the candidate right pixels in the corresponding line (epipolar line) in the right image. The window is shifted over all possible pixels so that a matching cost between the left pixel and each candidate in the right image is obtained. By the Winner-Takes-All method, the final disparity estimation is realized by selecting the window with the lowest matching cost.

The matching cost, SSD (Sum of Squared Differences cost), is adapted as [START_REF] Brad Hiebert-Treuer | Stereo Datasets[END_REF].

SSD ω (x l , y, s) = ω i=-ω ω j=-ω ÎlP A (x l + i, y + j) -ÎrP A (x l + i -s, y + j) 2 (5)
Where the • is the Euclidean norm. While s is the spatial shift along the horizontal epipolar line and ω the half-width of the (2ω + 1) × (2ω + 1) aggregation window.

Since these pixels of horizontal lines with the same parity in the left and right partially demosaiced color images are characterized by the same two color components, we can reasonably assume that the partial color points of two homologous pixels are similar. Because the partial costs compare the partial color points of left and right pixels located on the same horizontal lines, they reach an extremum when the shift is equal to the disparity.

III. EXPERIMENT

A. Experiment Platform and CUDA Implementation

We implement on GPU the partial demosaiced matching method of H. Halawana [START_REF] Halawana | Partial demosaicing of CFA images for stereo matching[END_REF]. We use the well known Middlebury databases [START_REF] Brad Hiebert-Treuer | Stereo Datasets[END_REF] to evaluate the method on both accuracy and computation time, as the image size grows. It is worth noting that the method was only implemented in a sequential computer (CPU) with no systematic evaluation of the tradeoff between quality and computation time. The ten datasets used in our experiments are entitled 'Aloe', 'Bowling1', 'Cloth1', 'Flowerpots', 'Lampshade1', 'Midd1', 'Monopoly', 'Plastic', 'Rocks1' and 'Wood1'. All the datasets of 2 views are used here (full-size (width: 1240...1396, height: 1110), half-size (width: 620...698, height: 555), and third-size (width: 413...465, height: 370) ) . In these datasets, the color stereo images are acquired by high resolution cameras equipped with one-single-sensor [START_REF] Scharstein | A taxonomy and evaluation of dense two-frame stereo correspondence algorithms[END_REF]. That's to say, the full color images are in fact color images that have been demosaiced by a specific chip integrated in the camera. They could contain artifacts caused by the demosaiced step. What's more, the work of applying a demosaicing step on CFA images, which have been generated by sampling color components from these previously demosaiced color images involves applying two successive demosaicing steps on the CFA images acquired by the camera. However, since Middlebury is the most frequently utilized database, an evaluation on this database will allow future comparative evaluation with new methods.

Our experiments are carried out both on CPU and GPU. For the CPU, we use a Intel(R) Core(TM)2 Duo CPU E8400 of 3.00GHz, with two cores having a cache of 6144KB. The GPU used in the experiments is a GPU GeForce GTX 570 of NVIDIA. It has 15 multiprocessors of 32 cores, GPU clock speed 1.54GHz, Memory clock rate 2000MHz, Memory Bus width 320 bits while the total amount of global memory is 1280 Mb (constant memory 65536 Kb, shared memory per block 49152 Kb). The system, on which the experiment is evaluated is Ubuntu 11.04 32 bits.

The programming interface we used for parallel computation on GPU is the Compute Unified Device Architecture (CUDA). The parallel computation work is realized by a kernel function, which is executed concurrently by multiple threads on data elements. All these threads are organized into a two level concepts: grid and block. A kernel has one grid, which contains multiple blocks. Every block is formed of multiple threads. The dimension of grid and of blocks can be one-dimension, two-dimension or three-dimension. The performance of GPU with CUDA is closely related to thread organization and memory accesses, which should attract much attention according to various computation works and GPU platform. Based on our experiment platform, given an image of size W×H as example and briefly lay out the CUDA settings and the parameter values mentioned in our algorithm.

In our experiments, we use the two-dimension blocks of 16 × 16 size. Every thread takes care of one pixel and the size of grid is obtained by W +16-1 16 × H+16-1 16 for a given W×H image. In the SCC estimation step, a grid of W×H threads is created and each thread takes care of one pixel in the estimation, and the shared memory is employed to for fast memory access. For the matching cost computation, a grid of W×H threads is employed to compute the matching cost for every pixel at a set of given disparity and then pick out the best homologous candidate pixel by the Winner-Takes-All method, here, the main data is stored on Global memory space while the shared memory space is used to support the computation of matching cost and the comparison in WTA processing.

B. Experiment Process

The procedures on CPU and on GPU are almost the same. As the datasets are all full color images, at the very beginning of the experiments, a simulation step for every pair of images is realized to obtain their CFA images needed by keeping only one of the three color components at every pixel. This work is done by GPU and by CPU separately . The whole evaluation is performed according to the spatial arrangement of the Bayers CFA. Then, for every CFA image, we do partial demosaicing step. Here, the partial demosaicing estimates only the missed Second Color Component (SCC) for every pixel by Hamilton's estimate method in the left image and in the right image. So, the left demosaiced color image and the right one are produced. The estimation method is the edge-adaptive demosaicing method proposed by Hamilton as presented in the precedent section.

The original full color image is shown in Fig. 2(a) and the left demosaiced color image using Hamiltons method is shown in Fig. 2(c). They look somewhat similar. However, zooming on the square areas outlined in these images as presented in Fig. 2(b) and Fig. 2(d), shows that textured areas are locally different.

Then the local stereo matching algorithm and the Winner-Takes-All method are taken into actions. In this step, as shown in Fig. 3 for every pixel in the left image, the method finds out its homologous pixel from a group of candidate pixels in the right image and compute their disparity s by computing the matching cost SSD, which is modified to adapt to the partial demosaiced color images. A matching is considered as valid when the absolute difference between the estimated disparity and the given benchmark d ω l (x l , y) is lower or equal to δ, which is the disparity error tolerance. In the experiments we set this coefficient to 1.

At the end, a disparity map, as illustrated in Fig. 2(e), is estimated for each pair of image.

IV. EXPERIMENTAL RESULTS AND DISCUSS

The experiment are executed on all the ten chosen datasets. Here, we take the 'Aloe' image group as an example. As it is shown in Fig. 4 that with the increase of the half-window, the computation time increases significantly. When the halfwindow is given as ω, for every pixel in the left image and for their every possible candidate pixel in the right image, we should compute a group of (2ω +1)×(2ω +1) pixels to obtain the matching cost. So, when the half-window ω increases, the performance. Though on CPU and on GPU, all the singleprecision floating-point computations follow the same accuracy standards. The accuracy-loss is more important on GPU owning to the accuracy problem of floating-point, that is why the results by CPU and that by GPU may have some deviations, especially when some cumulations in the computation exists. In our experiments, we use the floating point in the computation of matching cost with SSD and in the aggregation based on fixed support window, the deviation occurs at these computations, which can explain the tolerances between the RCMP by CPU and that by GPU shown in Fig.

4(b).

In addition, this method is also compared with the classic fixed windows matching methods (treating full color image and gray-level image separately) we used in former experiments, the results are shown in Table I. This method performs worse on all these four pairs both in matching quality and in computation time. That is because this method requires too many refers to the logical operations in the programming and too many branches in data treatments, which is the real weakness of GPU architecture. These branches and logical operations lead to great load on to the GPU system when loading data from the memory space and wastes GPU's CUDA cores, which are powerful in arithmetical operation, by making them do logical works.

V. CONCLUSION

This paper presents a GPU implementation of a partial demosaicing scheme specially designed for stereo matching of CFA image. By analyzing the CFA image directly, this method can handle the stereo matching works in case of single-CCD cameras usage. This method has three main techniques characters: the adapted matching cost for CFA image, the estimated Second color component based on Hamilton's estimate method, and the robust cost aggregation window. Experiments carried on show the performance of this method on GPU. The results show that this method can benefit from GPU's parallel architecture. The experiments also show that this method performs faster on GPU as the image size grows. Future improvements should come from a better setting of the memory arrangement in order to allow more coalescing accesses.

Fig. 2 :

 2 Fig. 2: 'Aloe' left image.

Fig. 3 :

 3 Fig. 3: matched pixel and its candidate ones in the right image.

Fig. 4 :

 4 Fig. 4: Compute-time and rate of correctly matched pixels (RCMP) obtained with the adapted SSD computed on the full sizeAloe stereo image pair for δ set to 1.

Fig. 5 :

 5 Fig. 5: GPU outperforms CPU in terms of computation time. The maximum, the minimum and the average acceleration factor obtained in the executions on 'Aloe' image pair based on half-window from 1 to 10 are marked. The acceleration increases along with the image dimensions.

TABLE I :

 I COMPARATIVE

EVALUATION ON THE FOUR MIDDLEBURY DATA SETS OF THE THREE DIFFERENT SIZES BY DIFFERENT GPU PROGRAMMING VERSION. 'COMPUTATION TIME (CT) ' AND 'PERCENTAGE OF BAD PIXELS (PBP)' .

Copyright (c) IARIA, 2013. ISBN: 978-1-61208-310-0 INFOCOMP 2013 : The Third International Conference on Advanced Communications and Computation

Copyright (c) IARIA, 2013. ISBN: 978-1-61208-310-0

Meanwhile, with the increase of the half-window, and so, the computation intensity, CPU's computation time increases more importantly than GPU's computation time. Moreover, as illustrated in Fig. 5, as the image dimensions augment, the acceleration factor obtained by using GPU has expanded from 1.75 to 20.75. This means that the GPU offers powerful computation capacity in intense computation thanks to the parallel organization of the blocks and the threads on GPU.

The performance of GPU is significantly influenced by the dimension of the image. The nearer to the complete load of the employed streaming multiprocessor on GPU, the higher performance we get.

The Ratio of Correctly Matched Pixel (RCMP) is the percentage of the well matched pixels in all the pixels of the image to be matched. As it is shown in Fig. 4, for 'Aloe ' image pair, the RCMP reaches the smooth peak when the half-window is between 4 and 8 (for those textureless image pairs, the half-window should be bigger to have their peak of RCMP). In fact, whatever the image type, the matching performance increases with aggregation window half-width. Small windows do not contain enough information to allow a correct matching. At the opposite, too large aggregation windows may cover image regions containing pixels with different disparities, which explains the decrease of matching