Calculation Scheme Based on a Weighted Primitive: Application to Image Processing Transforms

This paper presents a method to improve the calculation of functions which specially demand a great amount of computing resources. The method is based on the choice of a weighted primitive which enables the calculation of function values under the scope of a recursive operation. When tackling the design level, the method shows suitable for developing a processor which achieves a satisfying trade-o ﬀ between time delay, area costs, and stability. The method is particularly suitable for the mathematical transforms used in signal processing applications. A generic calculation scheme is developed for the discrete fast Fourier transform (DFT) and then applied to other integral transforms such as the discrete Hartley transform (DHT), the discrete co-sine transform (DCT), and the discrete sine transform (DST). Some comparisons with other well-known proposals are also provided.


INTRODUCTION
Mathematical notation aside, the motivation behind integral transforms is easy to understand.There are many classes of problems that are extremely difficult to solve or, at least, quite unwieldy from the algebraic standpoint in their original domains.An integral transform maps an equation from its original domain (time or space domain) into another domain (frequency domain).Manipulating and solving the equation in the target domain is, ideally, easier than manipulating and solving it in the original domain.The solution is then mapped back into the original domain.Integral transforms work because they are based upon the concept of spectral factorization over orthonormal bases.Equation (1) shows the generic formulation of a discrete integral transform where f (x), 0 ≤ x < N, and F(u), 0 ≤ u < N, are the original and the transformed sequences, respectively.Both have N = 2 n values, n ∈ N and T(x, u) is the kernel of the transform ( The inverse transform can be defined in a similar way.Table 1 shows some integral transforms ( j = √ −1 as usual).The Fourier transform (FT) is a reference tool in image filtering [1,2] and reconstruction [3].A fast Fourier transform (FFT) scheme has been used in OFDM modulation (orthogonal frequency division multiplexing) and has shown to be a valuable tool in the scope of communications [4,5].The most relevant algorithm for FFT calculation was developed in 1965 by Cooley and Tukey [6].It is based on a successive folding scheme and its main contribution is a computational complexity reduction that decreases from O(N 2 ) to O(N • log 2 N).The variants of FFT algorithms follow different ways to perform the calculations and to store the intervening results [7].These differences give rise to different improvements such as memory saving in the case of in-place algorithms, high speed for self-sorting algorithms [8] or regular architectures in the case of constant geometry algorithms [9].These improvements can be extended if combinations of the different schemes are envisaged [10].The features of the different algorithms point to different hardware trends.The in-place algorithms are generally implemented by pipelined architectures that minimize the latency between stages and the memory [11] whereas the constant geometry algorithms  have an easier control because of their regular structure based on a constant indexation through all the stages.This allows parallel data processing by a column of processors with a fixed interconnecting net [12,13].The Hartley transform is a Fourier-related transform which was introduced in 1942 by Hartley [14] and is very similar to the discrete Fourier transform (DFT), with analogous applications in signal processing and related fields.Its main distinction from the DFT is that it transforms real inputs into real outputs, with no intrinsic involvement of complex numbers.The discrete Hartley transform (DHT) analogue of the Cooley-Tukey algorithm is commonly known as the fast Hartley transform (FHT) algorithm, and was first described in 1984 by Bracewell [15][16][17].The transform can be interpreted as the multiplication of the vector (x 0 , . . ., x N−1 ) by an N ×N matrix; therefore, the discrete Hartley transform is a linear operator.The matrix is invertible and the DHT is its own inverse up to an overall scale factor.This FHT algorithm, at least when applied to power-of-two sizes N, is the subject of a patent issued in 1987 to the University of Stanford.The University of Stanford placed this patent in the public domain in 1994 [18].The DHT algorithms are typically slightly less efficient (in terms of the number of floatingpoint operations) than the corresponding FFT specialized for real inputs or outputs [19,20].The latter authors published the algorithm which achieves the lowest operation count for the DHT of power-of-two sizes by employing a split-radix algorithm, similar to that of the FFT.This scheme splits a DHT of length N into a DHT of length N/2 and two real-input DFTs (not DHTs) of length N/4.A priori, since the FHT and the real-input FFT algorithms have similar computational structures, none of them appears to have a substantial speed advantage [21].As a practical matter, highly optimized realinput FFT libraries are available from many sources whereas highly optimized DHT libraries are less common.On the other hand, the redundant computations in FFTs due to real inputs are much more difficult to eliminate for large prime N, despite the existence of O(N • log 2 N) complex-data algorithms for that cases.This is because the redundancies are hidden behind intricate permutations and/or phase rotations in those algorithms.In contrast, a standard prime-size FFT algorithm such as Rader's algorithm can be directly applied to the DHT of real data for roughly a factor of two less computation than that of the equivalent complex FFT.This DHT approach currently appears to be the only way known to obtain such factor-of-two savings for large prime-size FFTs of real data [22].A detailed analysis of the computational cost and specially of the numerical stability constants for DHT of types I-IV and the related matrix algebras is presented by Arico et al. [23].The authors prove that any of these DHTs of length N = 2 t can be factorized by means of a divide-andconquer strategy into a product of sparse, orthogonal matrices where in this context sparse means at most two nonzero entries per row and column.The sparsity joint with orthogonality of the matrix factors is the key for proving that these new algorithms have low arithmetic costs and an excellent normwise numerical stability.
DCT is often used in signal and image processing, especially for lossy data compression, because it has a strong "energy compaction" property: most of the signal information tends to be concentrated in a few low-frequency components of the DCT [24,25].For example, the DCT is used in JPEG image compression, MJPEG, MPEG [26], and DV video compression.The DCT is also widely employed in solving partial differential equations by spectral methods [27] and fast DCT algorithms are used in Chebyshev approximation of arbitrary functions by series of Chebyshev polynomials [28].Although the direct application of these formulas would require O(N 2 ) operations, it is possible to compute them with a complexity of only O(N • log 2 N) by factorizing the computation in the same way as in the fast Fourier transform (FFT).One can also compute DCTs via FFTs combined with O(N) pre-and post-processing steps.In principle, the most efficient algorithms are usually those that are directly specialized for the DCT [29,30].For example, particular DCT algorithms resemble to have a widespread use for transforms of small, fixed sizes such as the 8×8 DCT used in JPEG compression, or the small DCTs (or MDCTs) typically used in audio compression.Reduced code size may also be a reason for using a specialized DCT for embedded-device applications.However, even specialized DCT algorithms are typically closely related to FFT algorithms [22].Therefore, any improvement in algorithms for one transform will theoretically lead to immediate gains for the other transforms too [31].On the other hand, highly optimized FFT programs are widely available.Thus, in practice it is often easier to obtain high performance for generalized lengths of N with FFTbased algorithms.Performance on modern hardware is typically not simply dominated by arithmetic counts and optimization requires substantial engineering effort.
As DCT which is equivalent to a DFT of real and even functions, the discrete sine transform (DST) is a Fourierrelated transform using a purely real matrix [25].It is equivalent to the imaginary parts of a DFT of roughly twice the length, operating on real data with odd symmetry.As for DCT, four main types of DST can be presented.The boundary conditions relate the various DCT and DST types.
The applications of DST are similar to those for DCT as well as its computational complexity.The problem of reflecting boundary conditions (BCs) for blurring models that lead to fast algorithms for both deblurring and detecting the regularization parameters in the presence of noise is improved by Serra-Capizzano in a recent work [32].The key point is that Neumann BC matrices can be simultaneously diagonalized by the fast cosine transform DCT III and Serra-Capizzano introduces antireflective BCs that can be related to the algebra of the matrices that can be simultaneously diagonalized by the fast sine transform DST I.He shows that, in the generic case, this is a more natural modeling whose features are both, on one hand a reduced analytical error, since the zero (Dirichlet) BCs lead to discontinuity at the boundaries, the reflecting (Neumann) BCs lead to C • continuity at the boundaries, while his proposal leads to C 1 continuity at the boundaries, and on the other hand fast numerical algorithms in real arithmetic for deblurring and estimating regularization parameters.This paper presents a method that performs function evaluation by means of successive iterations on a recursive formula.This formula is a weighted sum of two operands and it can be considered as a primitive operation just as computational usual primitives such as addition and shift.The generic definition of the new primitive can be achieved by a two-dimensional table in which the cells store combinations of the weighting parameters.This evaluation method is suitable for a great amount of functions, particularly when the evaluation needs a lot of computing resources, and allows implementation schemes that offer a good balance between speed, area saving, and error containing.This paper is focused on the application of the method for the discrete fast Fourier transform with the purpose to extend the application to other related integral transforms, namely DHT, DCT, and DST.
The paper is structured in seven parts.Following the introduction, Section 2 defines the weighted primitive.Section 3 presents the fundamental concepts of the evaluation method based on the use of the weighted primitive, outlining its computational relevance.Some examples are presented for illustration.In Section 4, an implementation based on look-up tables is discussed and an estimation of the time delay, area occupation, and calculation error is developed.Section 5 is entirely devoted to the applications of our method for digital signal processing transforms.The calculation of the DFT is developed as a generic scheme and other transforms, namely the DHT, the DCT, and the DST are considered under the scope of the DFT.In Section 6 some comparisons with other well-known proposals considering operation counts, area, time delay, and stability estimations are presented.Finally, Section 7 summarizes results and presents the concluding remarks.

DEFINITION OF A WEIGHTED PRIMITIVE
The weighted primitive is denoted as ⊕ and its formal definition is as follows: (2) The operation ⊕ can also be defined by means of a twoinput table.Table 2 defines the operation for integer values in binary sign-magnitude representation; k stands for the number of significant bits in the representation.
In Table 2 the arguments have been represented in binary and decimal notation and the results are referred to in a generic way as combinations of the parameters α and β.The operation ⊕ is performed when the arguments (a, b) address the table and the result is picked up from the corresponding cell.The first argument (a) addresses the row whereas the second (b) addresses the column.
The same operation can be represented for greater values of k (see Table 3, for k = 2).Central cells are equivalent to those of Table 2.
The amount of cells in a table is (2 (k+1) − 1) 2 and it only depends on k.These cells are organized as concentric rings centred in 0. It can be noticed that increasing k causes a growth in the table and therefore the addition of more peripheral rings.The number of rings increases 2 k when k increases one unit.The smallest table is defined for k = 1 but the same information about the operation ⊕ is provided for any k value.When the precision of the arguments n is greater than k, these must be fragmented in k-sized fragments in order to perform the operation.So, t double accesses are necessary to complete t cycles of a single operation (if n = k • t).A single operation requires picking up from a table so many partial results as fragments are contained in the argument.The overall result is obtained by adding t partial results according to their position.
As the primitive involves the sum of two products, the arithmetic properties of the operation ⊕ have been studied with respect to those of the addition and multiplication.

Commutative
As shown, the commutative property is only verified when a = b or when α = β.
As noticed, the operation ⊕ is not associative except for a particular case given by αa( The lack of associative property obliges to fix arbitrarily an order in calculations execution.We assume that the operations are performed from left to right: (5) No neutral element can be identified for this operation.

Symmetry
Spherical symmetry can be proved by looking at the table: Proof So, a⊕b and −[a⊕b] are stored in diametrically opposite cells.
The primitive ⊕ does not fulfill the properties that allow the definition of a set structure.

A FUNCTION EVALUATION METHOD BASED ON THE USE OF A WEIGHTED PRIMITIVE
This section presents the motivation and the fundamental concepts of the evaluation method based on the use of the weighted primitive, outlining its computational relevance.

Motivation
In order to improve the calculation of functions which demand a great amount of computing resources, the approach developed in this paper aims for balancing the number of computing levels with the computing power of the corresponding primitive.That is to say, the same calculation may get the advantages steaming from the calculation at a lower computing level by other primitives than the usual ones whenever the new primitives assume intrinsically part of the complexity.This approach is considered as far as it may be a way to perform a calculation of functions with both algorithmic and architectural benefits.
Our inquiry for a primitive operation that bears more computing power than the usual primitive sum points towards the operation ⊕.This new primitive is more generic (usual sum is a particular case of weighted sum) and, as it will be shown, the recursive application of ⊕ achieves quite different features that mean much more than the formal combination of sum and multiplication.This issue has crucial consequences because function evaluation is performed with no more difficulty than applying iteratively a simple operation defined by a two-input table.

Fundamental concepts of the evaluation method
In order to carry out the evaluation of a given function Ψ we propose to approximate it through a discrete function F defined as follows: The first value of the function F is given by (F 0 ) and the next values are calculated by iterative application of the recursive equation (9).The approximation capabilities of function F can be understood as the equivalence between two sets of real values: on one hand {F i } and on the other hand {Ψ(i)} which is generated by the quantization of the function Ψ.The independent variable in function Ψ is denoted by z = x + ih, where x ∈ R is the initial value, h ∈ R is the quantization step, and i ∈ N can take successive increasing values.The mapping implies three initial conditions to be fulfilled.They are Table 4: Approximation of some usual generic functions by the recursive function F.

Usual function Ψ Mapping parameters for
(b) the successive samples of function Ψ are mapped to successive F i values irrespectively to the value of the quantization step, h; (c) the two previous assumptions allow not having to discern between i (index belonging to the independent variable of Ψ) and i (iteration number of F), that is to say: The mapping of function Ψ by the recursive function F succeeds in approximating it through the normalization defined in (a), (b), and (c).It can be noticed that the function F is not unique.Since different mappings related to different values of the quantization step h can be achieved to approximate the same function Ψ, different parameters α and β can be suited.
Table 4 shows the approximation of some usual generic functions.The first column shows different functions Ψ that have been quantized.The next four columns present the mapping parameters of the corresponding recursive functions F. All cases are shown for x = 0.
Any calculation of {F i } is performed with a computational complexity O(N) whenever {G i } is known or whenever it can be carried out with the same (or less) complexity.It can be outlined that the interest of the mapping by the function F is concerned with the fulfillment of this condition.This fact draws at least two different computing issues.The first develops new function evaluation upon the previous; that is to say, when function F has been calculated, it can play the role of G in order to generate a new function F. This spreading scheme provides a lot of increasing computing power, always with linear cost.The second scheme deals with the crossed paired calculation of functions F and G; that is to say, G is the auxiliary function involved in the calculation of F as well as F is the auxiliary function for calculation of G.In addition to the linear cost, the crossed calculation scheme provides time delay saving as both functions can be calculated simultaneously.
Figure 1: Arithmetic processor for the spreading calculation scheme.

PROCESSOR IMPLEMENTATION
As mentioned in Section 3, the two main computing issues lead to different architectural counterparts.The development of a new function evaluation upon the previous one in a spreading calculation scheme is carried out by the processor presented in Figure 1 that requires function G to be known.The second scheme deals with the crossed paired calculation of the F and G functions.The corresponding processor is shown in Figure 2.
The implementation proposed uses an LRA (acronym for look-up table (LUT), register, reduction structure, and adder).The LUT contains all partial products αA k + βB k ; A k , B k are portions of few bits of the current input data F i and G i .On every cycle, the LUT is respectively accessed by A k and B k coming from the shift registers.Then, the partial products are taken out of the cells (partial products in the LUT are the hardware counterpart of the weighted primitives presented in Tables 1 and 2).The overall partial product αF i +βG i is obtained by adding all the shifted partial products corresponding to all fragment inputs A k , B k of F i and G i , respectively.In the following iteration, both the new calculated F i+1 value and the next G i+1 value are multiplexed and shifted before accessing the LUT in order to repeat the addressing process.The processor in Figure 2 is different from Figure 1 in what concerns function G.The G values are obtained in the same way as for F but the LUT for G is different from the LUT for F.

Area costs and time delay estimation
In order to have the capability to make a comparison of computing resources, an estimation of the area cost and time delay of the proposed architectures is presented here.The model we use for the estimations is taken from the references [33,34].The unit τ a represents the area of a complex gate.The complex gate is defined as the pair (AND, XOR) that provides a meaningful unit, as these two gates implement the most basic computing device: the one bit full-adder.The unit τ t is the delay of this complex gate.This model is very useful because it provides a direct way to compare different architectures, without depending on their implementation features.As an example, the area cost and time delay for 16 bits one-bit fragmented data are estimated for both processors, as shown in Table 5.
If the fragments of the input data are greater than one bit, then the occupied area and the time delay access of the LUT vary.The relationship between area, time delay, and fragment length k for 16 bits data is shown in Table 6 for processor 2.
Table 6 outlines that the LUT area increases exponentially with k, and represents an increasing portion of the overall area as k increases.The access time for the LUT decreases as 1/k.The percentage of access time versus overall processing time decreases slowly as 1/k.The trade-off between area and time has to be defined depending on the application.
The proposed architecture has also been tested in the XS4010XL-PC84 FPGA.Time delay estimation in usual time units can also be provided assuming τ t ≈ 1 ns.

Algorithmic stability
A complete study of the error is still under consideration and numerical results are not yet available except for particular cases [35].Nevertheless, two main considerations are presented: on one hand, the recursive calculation accumulates the absolute error caused by the successive round-off which is performed as the number of iterations increases, on the other hand, if round-off is not performed, the error can become lower as the length in bits of the result increases, but the occupied area as well as the time delay increase too.In what follows, both trends are analyzed.

Round-off is performed
The drawback of the increasing absolute error can be faced by decreasing the number of iterations, that is to say the number of calculated values, with the corresponding loss of accuracy of the mapping.A trade-off between the accuracy of the approximation (related to the number of calculated values) and the increasing calculation error must be found.Parallelization provides a mean to deal with this problem by defining more computing levels.The N values of function F that are to be calculated can be assigned to different computing levels (therefore different computing processors) in a tree-structured architecture, by spreading N into a product as follows: -1st computing level: F 0 is the seed value that initializes the calculation of N 1 new values, -2nd computing level: the N 1 obtained values are the seeds that initialize the calculation of N 1 •N 2 new values (N 2 values per each N 1 ).
And so on until achieving the pth computing level: the N p−1 obtained values are the seeds that complete the calculation of If the error for one value calculation is assumed to be ε, the overall error after N values calculation is The parallelized calculation decreases the overall error without having to decrease the number of points.The minimum value for the overall error is obtained when the sum (N 1 + N 2 + • • • + N p ) is minimized, that is to say when all N i in the sum are relatively prime factors.
It can be mentioned that the time delay calculation follows a similar evolution scheme as the error.Considering T as the time delay for one value calculation, the overall time delay is The minimization of the time delay is also obtained when the N i are relatively prime factors.
For the occupied area, the precise structure of the tree in what concerns the depth (number of computing levels) and the number of branches (number of calculated values per processor) is quite relevant for the result.The distribution of the N i is crucial in the definition of some improving tendencies.The number of processors P in the tree-structure can be bounded as follows: P increases at the same rate as the number of computing levels p, but the growth can be contained if N p is the maximum value of all N i , that is to say in the last computing level p − 1, the number of calculated values per processor is the highest.It can be observed that the parallel calculation involves much more processors than sequential one processor.

Summarizing the main ideas
(i) The parallel calculation provides benefits on error bound and time delay whereas sequential calculation performs better in what concerns area saving.
(ii) A trade-off must be established between the time delay, the occupied area, and the approximation accuracy (through the definition of the computing levels).

Round-off is not performed
As explained in Section 2, we assume the first input data length is n, the data have been fragmented (n = kt), and the partial products in the cells are p bits long.If t accesses have been performed to the table and t partial products have to be added, the first result will be p + t + 1 bits long (t bits represent the increase caused by the corresponding shifts plus one bit for the last carry).The second value has to be calculated in the same way so that the p + t + 1 bits of the feedback data is k-fragmented and the process goes on.This recursive algorithm can be formalized as follows: Initial value n bits = A 0 bits 1st calculated value and so on.Table 7 presents the data length evolution and the corresponding error for n = p = 16, 32, and 64 bits data, as well as the number of calculated values that lead to the maximum data length achievement.
It can be noticed that the increase of the number of bits is bounded after a finite and rather low number of calculated values that decreases as k grows.As usual, the error decreases as the number of the data bits increases and the results are improved in any case by small fragmentation (k = 2).When round-off is not performed, time delay and area occupation increase because of the higher number of bits involved, so Tables 5 and 6 should be modified.It can be outlined that small fragmentation makes error to decrease, but time delay would increase too much.By increasing the fragment length value, time delay improves but the error and the area cost would make this issue infeasible.The trade-off between area, time delay, and error must be set regarding to the application.

GENERIC CALCULATION SCHEME FOR INTEGRAL TRANSFORMS
In this section, a generic calculation scheme for integral transforms is presented.The DFT is taken as a paradigm and some other transforms are developed as applications of the DFT calculation.

The DFT as paradigm
Equation ( 13) is the expression of the one-dimensional discrete Fourier transform.Let us have The Cooley and Tukey algorithm segregates the FT in even and odd fragments in order to perform the successive folding scheme, as shown in (14): For any u ∈ [0, M[, the Cooley and Tukey algorithm starts by setting the M initial two-point transforms.In the second step M/2 four-point transforms are carried out by combining the former transforms and so on till to reach the last step, where one M-point transform is finally obtained.
For values of u ∈ [M, N[ no more extra calculations are required as the corresponding transforms can be obtained by changing the sign, as shown by the second row in (14).
Our method enhances this process by adding a new segregation held by both real (R) and imaginary (I) parts in order to allow the crossed evaluation presented at the end of Section 3. Due to the fact that two segretations are considered (even/odd, real/imaginary) there will be, for each u, four transforms, which are R p,q even , R p,q odd , I p,q even , and I p,q odd where p, q denote the step of the process and the number of the transform in the step, respectively, p ∈ [0, n − 1], and Equations ( 15), (16), and (17) show the first, the second, and the last steps of our process, respectively, for any u ∈ [0, M[.Parameters α p (u) = cos pπu/M and β p (u) = sin pπu/M define the step p.The u argument has been omitted in ( 16) and ( 17) in order to clarify the expansion.In the first step, M two-point real and imaginary transforms are set in order to start the process.In the second step M/2 real and imaginary transforms are carried out following the calculation scheme shown in (9).At the end of the process, one real and one imaginary M-point transform are achieved and, without any more calculation, the result is deduced for u ∈ [M, N[.As observed in ( 16) and ( 17), each step involves the results of R and I obtained in the two previous steps; therefore, in each step the number of equations is halved.After the first step, a sum is added to the weighted primitive.This could have an effect on the LUT as the parameter set becomes (α, β, 1), R 1,0 even = R 0,0 even + α 1 R 0,1 odd − β 1 I 0,1 odd = R 0,0 even + R 0,1 odd ⊕ I 0,1 odd , I 1,0 even = I 0,0 even + β 1 R 0,1 odd + α 1 I 0,1 odd = I 0,0 even + R 0,1 odd ⊕ I 0,1 odd , R 1,1 odd = R 0,2 even + α 1 R 0,3 odd − β 1 I 0,3 odd = R 0,2 even + R 0,3 odd ⊕ I 0,3 odd , I 1,1 odd = I 0,2 even + β 1 R 0,3 odd + α 1 I 0,3 odd = I 0,2 even + R 0,3 odd ⊕ I 0,3 odd , The number of operations has been used as the main unit to measure the computational complexity of the proposal.The operation implemented by the weighted primitive has been denoted as weighted sum WS, and the simple sum as SS.The calculations take into account both real and imaginary parts for any u value.The initial two-point transforms are assumed to be calculated.An inductive scheme is used to carry out the complexity estimations.From these results two induced calculation formulas can be proposed referring to the count of needed weighted sums and simple sums, Proof.Starting from WS(1) = 2 and SS(1) = 0, for any n, n > 1, it may be assumed that By the application of the inductive scheme, after substituting n by n + 1 the formulas become Comparing the expressions for n and n + 1, it can be noticed that The proposed formulas (see (19)) have been validated by this proof.
Comparing with the Cooley and Tukey algorithm, where M(n) is the number of multiplications and S(n) the number of sums, we have The contribution of the weighted primitive is clear as we compare (19) and (23).The quotient M(n)/ WS(n) increases linearly versus n.The same occurs with the quotient S(n)/ SS(n) but with a steeper slope.So, the weighted primitive provides best results as n grows.

Other transforms
This calculation scheme can be applied to other transforms.As DHT and DCT/DST are DFT-related transforms, a common calculation scheme can be presented after we perform some mathematical manipulations.

Hartley transform
Let H(u) be the discrete Hartley transform of a real function f (x): H(u) is the transformed sequence that can split into two fragments: R(u) corresponds to the cosine part and I(u) to the sine part.The whole previous development for the DFT can be applied but the last stage has to perform an additional sum of the two calculated fragments, The number of simple sums increases as one last sum must be performed per each u value.Nevertheless, (19) suits because only the initial value varies, SS(1 Cosine/sine transforms Let C(u) be the discrete cosine transform of a real function f (x): C(u) is the transformed sequence that can split into two fragments as follows: So that ( 27) leads to ( 29) Then, cos[πu/2N] and − sin[πu/2N] are constant values for each u value and can lay outside the summation: Both fragments, R(u) (for the cosine part) and I(u) (for the sine part), can be carried out under the DFT calculation scheme and combined in the last stage by an additional weighted sum: A similar result could be inferred for sine transform with the following parameter values: cos(πu/2N) = α u , sin(πu/2N) = β u .
The number of weighted sums increases because of the last weighted sum that must be performed, see (31).The equation has been modified as the constant value in WS(n) varies.The reason is that the initial value WS(1) = 3,

Summarizing
The calculation based upon the DFT scheme leads to an easy approach for the calculation of the DHT and the DCT/DST, as expected.This scheme can be extended to other integral transforms with trigonometric kernel.

COMPARISON WITH OTHER PROPOSALS AND DISCUSSION
In this section, some hardware implementations for the calculation of the DFT, DHT, and DCT are presented in order to provide a comparison for the different performances in terms of area cost, time delay, and stability.

DFT
The BDA proposal presented by Chien-Chang et al. [36] carries out the DFT of variable length by controlling the architecture.The single processing element follows the Cooley and Tukey algorithm radix-4 and calculates 16/32/64 points transform.When the number of points N grows, it can split out into a product of two factors N 1 × N 2 in order to process the transform in a row-column structure.Formally, the four terms of the butterfly are set as a cyclic convolution that allows performing the calculations by means of distributed arithmetic based on blocks.The memory is partitioned in blocks that store the set of coefficients involved in the multiplications of the butterfly.A rotator is added to control the sequence of use of the blocks and avoids storing all the combinations of the same elements as in conventional distributed arithmetic.This architecture improves memory saving in exchange for increasing the time delay and the hardware because of the extra rotator in the circuit.This proposal substitutes the ROM by a RAM in order to make more flexible the change of the set of coefficients when the length of the Fourier transform varies.The processing column consists of an input buffer, a CORDIC processor that runs the complex multiplications followed by a parallel-serial register and a rotator.Four RAM memories and sixteen accumulators implement the distributed arithmetic.At last, four buffers are where N 1 is the length of the transform, M = 4 in the design, and W L is the data length.When the transform is longer as 64 points, N 1 is substituted by the N 1 × N 2 .Table 8 shows the results obtained by the synopsis implementation of the circuit that has been described in Verilog HDL.
In order to compare the performance of our architecture and that of the BDA, an estimation of the occupied area and time delay is provided.The devices for both implementations are listed in Table 9 and evaluated in terms of τ t and τ a in Table 10.For the crossed evaluation scheme, the architecture is double because of the two segregations (even/odd and real/imaginary); 64 cells LUTs are assumed as the parameter set is (α, β, 1).Data is 16 bits long for any proposal.In Table 10, neither the rotator nor the CORDIC processor has been considered in the BDA implementation because the reference does not facilitate any detail upon their structure.The estimations of the time delay are based on the author's indications and presented in terms of τ a and τ t units.
It can be observed that the BDA architecture is worse than the crossed one in what concerns the occupied area because the BDA hardware needs to be increased stepwise when the number of points of the transforms increases.The time delay is lower for the crossed architecture than for the BDA for the values of N that have been considered and will remain lower for any N, because it achieves a linear growing in both implementations.
Table 11 summarizes the hardware cost as well as the time delay of proposals for the Fourier transform calculation presented by different authors [13,[37][38][39][40].The four proposals in the beginning of the list have based their design on systolic matrices, the following one on adders and the others on distributed arithmetic (the DA is a generic distributed arithmetic approach).At the end of the list appears our proposal.Average computation time is indicated as It appears that our proposal is the best in what concerns the hardware resources but time delay has a linear growth with respect to N (number of points of the transform) and with the data precision.It can be remembered that parallel architecture may present a better performance for this case.

DHT
As mentioned in Section 1, the DHT algorithms are typically less efficient (in terms of the number of floating-point for the DHT and the DFT for power-of-two sizes, as achieved by the split-radix Cooley-Tukey FHT/FFT algorithm in both cases.Notice that, depending on DFT and DHT implementation details, some of the multiplications can be traded for additions or vice versa.The third column of the table estimates the operation counts (weighted sums + simple sums) to be performed by our proposal, following (19).As expected, our proposal behaves better in what concerns the operation counts than both the DHT algorithm and the corresponding DFT algorithm specialized for real inputs or outputs.With respect to the particular hardware implementations, as the DFT has already been compared above with our proposal, the concluding remarks related to the DHT have to be deduced.
A detailed analysis of the computational cost and especially of the numerical stability constants for DHT is presented by Arico et al. in [23].The authors base their research on the close connection existing between fast DHT algorithms and factorizations of the corresponding orthogonal Hartley matrices of length N, H N .They achieve a factorization of the matrix H N into a product of sparse matrices (at most, two nonzero entries per row and column) that allows an iterative calculation of H N x, for any x ∈ R N .Since the matrices are sparse and orthogonal, the factorization of H N generates a fast and low arithmetic cost DHT algorithms.The intraconnection of Hartley matrices of types (II), (III), and (IV) is expressed by means of other Hartley matrix of type (I), H N (I), is pursued by means of twiddle matrices T N and T N (that are direct sums of 1 and of rotationreflection matrices of order 2).Finally, factorization of H N (I) is achieved requiring permutations, scaling operations, butterfly operations, and plane rotations with small angles.The computational complexity is calculated for all types DHT-X, X = I, II, III, and IV but for comparison with our results we will consider the best result which is for X = I.The number of additions is denoted by α(DHT-I, N) and the number of multiplications by μ(DHT-I, N): As seen in the paper, the operation error follows the IEEE precision arithmetic u = 2 −24 or u = 2 −53 depending on the precision of the mantissa (24 or 53 bits, resp.).The roundoff algorithmic errors are related to the structure of the involved matrices and for direct calculation the round-off error is evaluated as a squared distance bounded by an expression ≈ k N u.The numerical stability is measured by k N that can be understood as the relative error on the output vector (of the mapping previously defined).For any X, a different expression for k N is obtained for the corresponding DHT-X(N).All k N expressions are similar and have linear dependence of log 2 N.For example, the normwise forward stability for DHT-I(N) is 4 3 As far as we can compare this very deep and strong theoretical approach with our method that is rather empirical, the results that can be taken into account are the computational cost and the stability of the algorithms.To make easier the comparison with our paper in what concerns the number of operations to be performed, a recursive formulation of (DHT-I, N) and μ(DHT-I, N) = for N = 2 n has been deduced from (34):  The initial values are for n =1, following (34): The comparison between ( 19) and (36) (WS(n) versus μ(n) and SS(n) versus α(n)) outlines that α(n) and μ(n) increase at a higher speed than WS(n) and SS(n), respectively, (i) for all n, α(n) > SS(n), (ii) for n > 6, μ(n) > WS(n).
The value of the normwise forward stability in the case of DHT-I (N) is (((4/3) In order to compare with our results The comparison between Tables 7 and 13 shows that for 16 bits (fragmentation lengths k = 2 and k = 4), for 32 bits data (k = 2, 4, and 8) and for 64 bits data (k = 2, 4, and 8) our algorithm behaves better.

DCT
The search for recursive algorithms with regular structure and less computation time remains an active research area.The recursive algorithms for computing 1D DCT are highly regular and modular [41][42][43][44][45][46][47].However, a great number of cycles are required to compute the 2D transformation by using 1D recursive structures.For computing the 2D DCT by row-column approaches, the row (column) transforms of the input 2D data are first determined.A transposition memory is required to store those temporal results.Finally, the 2D DCT results are obtained by the column (row) transforms of the transposed data.The RAM is usually adopted as the transposition memory.This approach has disadvantages such as higher-power consumption and long access time.Chen et al. develop in 2004 a new recursive structure with fast and regular recursion to achieve fewer recursive cycles without using any transposition memory structure [48].The 2D recursive DCT/IDCT algorithms are developed considering that the data with the same transform base can be pre-added such that the recursive cycles can be reduced.First, the 2D DCT/IDCT is decomposed into four portions which can be carried out either by 1D DCT or 1D DST (discrete sine transform).Based on the use of Chebyshev polynomials, efficient transform kernels are obtained for the 1D DCT and the DST.A reduction on the number or recursive cycles is achieved by a further folding on the inputs of the transform kernels.Considering other fast algorithms, the N × N DCT which maps the 2D index of the input sequence into the new 1D index is decomposed into N length-N 1D DCTs [49,50].Table 14 presents the number of multiplication and addition operations for these fast algorithms, for the case of 4 × 4 DCTs.Our proposal can be compared by assimilating the weighted sums and the multiplications (see (32)).
The number of operations required for our proposal is lower than those required for the existing methods.Table 15 shows the number of recursive cycles for different N ×N DCT recursive structures in five different algorithms [43,44,[46][47][48].In [48], a recursive cycle represents the time delay needed for computing the 2D DCT cosine transform for a pair of frequency indexes.The circuit involves two parallel identical block diagrams, both with a condensed 1D DCT/DST IIR filter which obtains the corresponding input data from a recursive input buffer in order to perform the partial calculation of the transform.In the last stage, the transform is recombined by a sum of the two partial results.So, the overall time delay for the 2D may be the same as for the 1D and the comparison with our proposal can be done as we assimilate the number of recursive cycles with the number of weighted sums to be performed following (32).
It can be outlined that our proposal has a better performance than the other ones, namely fast and recursive algorithms, in what concerns the number of recursive cycles.In [48] the chip area can be estimated as we depict the hardware recursive circuitry.Table 16 summarizes the hardware devices of the recursive architecture compared with our proposal for 4 × 4 DCT transform.
It can be observed that the devices for the implementation of the recursive architecture are numerous.Therefore, greater values for N × N may imply an increase of the chip area; the reason is the growth of the storing memory required for the buffers and for the number of outputs of the demultiplexer.Reference [48] does not offer any estimation of the time delay of the calculation.Our proposal implementation is very simple and has no variation related to the amount of devices when the number of calculated values varies.With respect to the time delay of the calculation in [48], as far as we can suppose, it can be estimated by analyzing the critical path of the depicted circuit.It seems to be higher than our proposal's one.

CONCLUSIONS
This paper has presented an approach to the scalability problem caused by the exploding requirements of computing resources in function calculation methods.The fundamentals of our proposal claim that the use of a more complete primitive, namely a weighted sum, converts the calculation of the function values into a recursive operation defined by a twoinput table.The strength of the method is concerned with the fact that the operation to be performed is the same for the evaluation of different functions (elementary or not).Therefore, only the table must be changed because it holds the features of the concrete evaluated function in the parameter values.This method provides a linear computational cost when some conditions are fulfilled.Image processing transforms that involve combined trigonometric functions provide an interesting application field.A generic calculation scheme has been developed for the DFT as paradigm.Other image transforms namely the DHT and the DCT/DST are analyzed under the scope of the DFT.When comparing with other well-known proposals, it has been confirmed that our approach provides a good trade-off between hardware resource and time delay saving as well as encouraging partial results in what concerns error contention.

Figure 2 :
Figure 2: Arithmetic processor for the crossed paired evaluation.

Figure 3 :
Figure 3: Growing rates s(n) and m(n) versus n.

Table 2 :
Definition of the operation ⊕

Table 3 :
Definition of the operation ⊕ for k = 2.

Table 5 :
Arithmetic processor estimations of area cost and time delay for 16 bits and one-bit fragmented data.

Table 6 :
Relationship between area, time delay, and fragment length k, for 16 bits data for processor 2.

Table 7 :
Data length evolution and error versus number of calculated values for n = p = 16, 32, and 64 bits.

Table 8 :
Critical path of the basic calculation module in the BDA architecture.

Table 9 :
Comparison between the hardware needed by BDA and our architecture implementations.

Table 10 :
Comparison between the BDA and our architecture implementations in terms of τ a and τ t .

Table 11 :
Comparison between our proposal and other ones.

Table 12 :
Lowest known operation counts (real multiplications + additions) for power-of-two DHT and corresponding DFT algorithms versus our proposal (weighted sums + simple sums).

Table 14 :
Number of multiplication and addition operations for different 4 × 4 DCTs.

Table 15 :
Number of recursive cycles for different N × N DCT recursive structures.

Table 16 :
Comparison between the hardware needed by the recursive architecture versus that of the implementation of our proposal for 4 × 4 DCT transform.

Table 7 ,
the previous formula has been calculated for the cases u = 2 −16 , 2 −32 , and 2 −64 bits and for different values of N.