A Short Survey on Arithmetic Transforms and the Arithmetic Hartley Transform

Arithmetic complexity has a main role in the performance of algorithms for spectrum evaluation. Arithmetic transform theory offers a method for computing trigonometrical transforms with minimal number of multiplications. In this paper, the proposed algorithms for the arithmetic Fourier transform are surveyed. A new arithmetic transform for computing the discrete Hartley transform is introduced: the Arithmetic Hartley transform. The interpolation process is shown to be the key element of the arithmetic transform theory.


Introduction and Historical Background
Despite the existence of fast algorithms for discrete transforms (e.g., fast Fourier transform, FFT), it is well known that the number of multiplications can significantly increase their computational (arithmetic) complexity. Even today, the multiplication operation consumes much more time than addition or subtraction. Table 1 brings the clock count of some mathematical operations as implemented for the Intel Pentium TM processor. Observe that multiplications and divisions can be by far more time demanding than additions, for instance. Sine and cosine function costs are also shown.
This fact stimulated the research on discrete transform algorithms that minimize the number of multiplications. The Bhatnagar's algorithm [1a], which uses Ramanujan numbers to eliminate multiplications (however, the choice of the transform blocklength is rather limited), is an example. Parallel to this, approximation approaches, which perform a trade-off between accuracy and computational complexity, have been proposed [2a-4a].
Arithmetic transforms emerged in this framework as an algorithm for spectrum evaluation, aiming the elimination of multiplications. Thus, it may possess a lower computational complexity. The theory of arithmetic transform is essentially based on Möbius function theorems [5a], offering only trivial multiplications, i.e., multiplications by {−1, 0, 1}. Therefore, only addition operations (except for multiplications by scale factors) are left to computation. Beyond the computational attractiveness, arithmetic transforms turned out to be naturally suited for parallel process- * R. J. Cintra Table 1: Clock count for some arithmetic instructions carried on a Pentium TM processor [1] Operation Clock count add 1-3 sub 1-3 fadd 1-7 fsub 1-7 mul (unsigned) 10-11 mul (signed) 10-11 div (unsigned) 17-41 div (signed) 22-46 fdiv 39 sin, cos  ing and VLSI design [2,3]. The very beginning ofthe research on arithmetic transforms dates back to 1903 when the German mathematician Ernst Heinrich Bruns 1 , depicted in Figure 1, published the Grundlinien des wissenschaftlichnen Rechnens [4], the seminal work in this field. In spite of that, the technique remained unnoticed for a long time. Forty-two years later, in Baltimore, U.S.A., the Hungarian Aurel Freidrich Wintner 2 , privately published a monograph entitled An Arithmetical Approach to Ordinary Fourier Series. This monograph presented an arithmetic method using Möbius function to calculate the Fourier series of even periodic functions.
In the quest to implement it, two other researchers played an important role: Dr. Oved Shisha of the U.R.I. Department of Mathematics and Dr. Charles Rader of Lincoln Laboratories. They were aware of Wintner's monograph and helped Tufts in many discussions. In 1988 The Arithmetic Fourier Transform by Tufts and Sadasiv was published in IEEE Acoustic, Speech, and Signal Processing (ASSP) Magazine [2].
Another breakthrough came in early 1990s is due to Emeritus Professor Dr. Irving S. Reed-the coinventor of the widely used Reed-Muller (1954) and Reed-Solomon (1964) codes. Author of hundreds of publications, Dr. Reed made important contributions to the area of signal processing. Specifically on arithmetic transforms, in 1990 Reed, Tufts and co-workers provided two fundamental contributions [5,6]. In [5], a reformulated version of Tufts-Sadasiv approach, the arithmetic Fourier transform (AFT) algorithm was able to encompass a larger class of signals and to compute Fourier series of odd periodic functions as well as even periodic ones. The publication of the 1992 A VLSI Architecture for Simplified Arithmetic Fourier Transform Algorithm by Dr. Reed and collaborators in the IEEE Transactions on ASSP [6] was another crucial slash on the subject. Indeed, that paper was previously presented at the International Conference on Application Specific Array Processors held in Princeton. However, the 1992 publication reached a vastly larger public, since it was published in a major journal. The new method, an enhancement of the last proposed algorithm [5], was re-designed to have a more balanced and computationally efficient performance. As a matter of fact, Reed et al. proved that the newly proposed algorithm was identical to Bruns' original method.
When the AFT was introduced, some concerns on the feasibility of the AFT were pointed out [7]. The main issue dealt with the number of samples required by the algorithm. However, later studies showed that the use of interpolation techniques on a sub-sampled set (e.g., zero-and first-order interpolation) could overcome these difficulties [8].
The conversion of the standard 1-D AFT into 2-D versions was just a matter of time. Many variants were proposed following the same guidelines of the 1-D case [9][10][11][12][13][14]. Further research was carried out seeking different implementations of the AFT. An alternative method [15] proposed a "Möbius-function-free AFT". Iterative [16] and adaptative approaches [17] were also examined. In spite of that, the most popular presentations of the AFT are still those found in [5,6].
Although the main and original motivation of the arithmetic algorithm was the computation of the Fourier Transform, further generalizations were performed and the arithmetic approach was utilized to calculate other transforms. Dr. Luc Knockaert of Department of Information Technology at Ghent University, Belgium, amplified the Bruns procedure, defining a generalized Möbius transform [18,19]. Moreover, four versions of the cosine transform was shaped in the arithmetic transform formalism [20].
Further generalization came in early 2000s with the definition of the Arithmetic Hartley Transform (AHT) [21,22]. These works constituted an effort to make arithmetical procedure applicable for the computation of trigonometrical transforms, other than Fourier transform. In particular the AHT computes the discrete Hartley transform 3 : the real, symmetric, Fourier-like discrete transform defined in 1983 by Emeritus Professor Ronald Newbold Bracewell in The Discrete Hartley Transform, an article published in the Journal of Optical Society of America.
In 1988 and then the technological state-of-art was dramatically different from that Bruns and Wintner found. Computational facilities and digital signal processing integrated circuits made possible AFT to leave theoretical constructs and reach practical implementations. Since its inception in engineering, the AFT was recognized as tool to be implemented with VLSI techniques. Tufts himself had observed that AFT could be naturally implemented in VLSI architectures [2]. Implementations were proposed in [3,6,12,[23][24][25][26][27][28][29][30][31][32][33]. Initial applications of the AFT took place in several areas: pattern matching techniques [34], measurement and instrumentation [35,36], auxiliary tool for computation of z-transform [37,38], and imaging [39]. This paper is organized in two parts. In Section 2, the mathematical evolution of the Arithmetic Fourier Transform is outlined. In Section 3, a summary of the major results on the Arithmetic Hartley Transform is shown. Interpolation issues are addressed and many points of the AFT are clarified, particularly the zero-order approximation.

The Arithmetic Fourier Transform
In this section, the three major developments of the arithmetic Fourier transform technique are presented. With emphasis on the theoretical groundwork, the AFT algorithms devised by Tufts, Sadasiv, Reed et alli are briefly outlined. In this work, k 1 |k 2 denotes that k 1 is a divisor of k 2 ; ⌊·⌋ is the floor function and [·] is the nearest integer function.
Lemma 1 Let k, k ′ and m be integers.
Proof: Consider the expression When Otherwise, we have: Therefore, it follows that: By taking the real and imaginary parts, we conclude the proof.
Definition 1 (Möbius µ-function) For a positive integer n, if p 2 |n for some prime p. ( A relevant lemma related to µ µ µ-function is stated below.

Theorem 1 (Möbius Inversion Formula for Finite Series)
Let n be a positive integer and f n a non-null sequence for 1 ≤ n ≤ N and null for n > N . If This is the finite version of the Möbius inversion formula [5a]. A proof can be found in [5].

Tufts-Sadasiv Approach
Consider a real even periodic function expressed by its Fourier series, as seen below: The components v k (t) represent the harmonics of v(t), given by: where a k is the amplitude of the kth harmonic. It was assumed, without loss of generality, that v(t) had unitary period and null mean (a 0 = 0). Furthermore, consider the N first harmonics as the only significant ones, in such a way that v k (t) = 0, for k > N (bandlimited approximation). Thus the summation of Equation 7 might be constrained to N terms.

Definition 2
The nth average is defined by for n = 1, 2, . . . , N . S n (t) is null for n > N .
After an application of Equations 7 and 8 into 9, it yielded: Proceeding that way, the nth average can be written in terms of the harmonics of v(t), instead of its samples (Definition 2). Since we assumed v n (t) = 0, n > N , only the first ⌊N/n⌋ terms of Equation 10 might possibly be nonnull. As a consequence the task was to invert Equation 10. Doing so, the harmonics can be expressed in terms of the averages, S n (t), which are derived from the samples of the signal v(t). The inversion is accomplished by invoking the Möbius inversion formula.
Proof: Some manipulation is needed. Substituting Equation 10 in Equation 11, it yields Now it is the tricky part of the proof.
According to Lemma 2, the inner summation can only be null if j/k = 1. In other words, the term v k (t) is the only survivor of the outer summation and the proof is completed.
The following aspects of the Reed-Tufts algorithm can be highlighted [2]: • This initial version of the AFT had a strong constraint: it can only handle even signals; • All computations are performed using only additions (except for few multiplications due to scaling); • The algorithm architecture is suitable for parallel processing, since each average is computed independently from the others; • The arithmetic transform theory is based on Fourier series, instead of the discrete Fourier transform.

Reed-Tufts Approach
Presented by Reed et al. in 1990 [5], this algorithm is a generalization of Tuft-Sadasiv method. The main constraint of the latter procedure (handling only with even signals) was removed, opening path for the computation of all Fourier series coefficient of periodic functions.
Let v(t) be a real T -periodic function, whose N -term finite Fourier series is given by where a 0 is the mean value of v(t). The even and odd coefficients of the Fourier series are a n and b n , respectively. Letv(t) denote the signal v(t) removed of its mean value a 0 . Consequently, A delay (shift) of αT inv(t) furnishes to the following expression:v where −1 < α < 1 and c n (α) = a n cos(2πnα) + b n sin(2πnα), (17) d n (α) = −a n sin(2πnα) + b n cos(2πnα).
In the sequel, the computation of the Fourier coefficients a n and b n based on c n (α) is outlined. Meanwhile, the formula for the nth average (Tufts-Sadasiv) is updated by the next definition.
Now the quantities c n (α) can be related to the averages, according to the following theorem.
Theorem 3 The coefficients c n (α) are computed via Möbius inversion formula for finite series and are expressed by Proof: Substituting the result of Equation 16 into Equation 19 furnishes the following expression: A direct application of Lemma 1 yields Invoking the Möbius inversion formula for finite series, the theorem is proved. Finally, the main result can be derived.
Theorem 4 (Reed-Tufts) The Fourier series coefficients a n and b n are computed by a n = c n (0), where k and m are determined by the factorization n = 2 k (2m + 1).
Joining these two sub-cases, it is easy to verify that The number of real multiplications and additions of this algorithm is given by [5] and respectively, where N is the blocklength of the transform.

Reed-Shih (Simplified AFT)
Introduced by Reed et al. [3], this algorithm is an evolution of that one developed by Reed and Tufts. Surprisingly, in this new method, the averages are re-defined in accordance to the theory created by H. Bruns [4] in 1903.

Definition 4 (Bruns Alternating Average) The 2nth
Bruns alternating average, B 2n (α), is defined by Invoking the definition of c n , applying Theorem 3 and Definition 3, the following theorem was derived.
Theorem 5 The coefficients c n (α) are given by the Möbius inversion formula for finite series as Proof: See [6]. Since a relation between the signal samples and the Bruns alternating averages was obtained, as well as an expression connecting the Bruns alternating averages to the c n coefficients, was available, few points are missing to compute the Fourier series coefficients. Actually, an expression relating the Fourier series coefficients (a n and b n ) to the coefficients c n is sought. Examining Equation 17, two conditions are distinguishable: • a n = c n (0); • b n = c n 1 4n . Calling Theorem 5, the next result was obtained.
Theorem 6 (Reed-Shih) The Fourier series coefficients a n and b n are computed by for n = 1, . . . , N .
Proof: The proof is similar to the proof of Theorem 4. For a blocklength N , the multiplicative and additive complexities are given by and respectively. The AFT algorithm proposed by Reed-Shih presents some improvements over previous algorithms: • The computation of both a n and b n has roughly the same computational effort. The algorithm is more "balanced" than Reed-Tufts algorithm; • The algorithm is naturally suited to a parallel processing implementation; • It is computationally less complex that Reed-Tufts algorithm.

An Example
In this subsection, we draw some comments in connection to an example originally proposed by Reed et al. [3]. Let v(t) be a signal with period T = 1 s. Consider the computation of the Fourier series coefficients up to the 5th harmonic.
According to the Reed-Shih algorithm, the coefficients a n and b n of the Fourier series of v(t) are expressed by Comparing these formulations with the ones of Reed-Tufts algorithm, one may note the balance in the computation of a n and b n . Both coefficients are obtained through similar matrices. • This algorithm is not naturally suited for uniform sampling.
• A uniform sampler utilized to obtain all the required samples would need a very high sampling rate. In the example illustrated here, a 120 Hz clock should be required to sample the necessary points for the computation of the Fourier series of a 1 Hz bandlimited signal.
Certainly these observations appear to be disturbing and seems to jeopardize the feasibility of the whole procedure. However, it is important to stress that this procedure furnishes the exact computation of the Fourier series coefficients.
An empirical solution to circumvent this problem is to interpolate. An interpolation based on uniformly sampled  Table 2 shows, for example, that the computation of B 4 (0) requires -among other samples -v 1 4 , which is clearly not available. To overcome this difficulty, a rounding operation can be introduced. Thus, the sample v 3 10 can be used whenever the algorithm called v 1 4 ( 10 1 4 /10 = 3/10). This rounding operation is also known as zero-order interpolation.
The accuracy of the AFT algorithm is deeply associated with the sampling period T 0 . If more precision is required, then one should expect to increase sampling rate, resulting in the introduction of smaller errors due to the interpolation scheme. Higher order of interpolation (e.g. first-order interpolation) can also be used to obtain more accurate estimations of the Fourier series coefficients. The following trade-off is quite clear accuracy versus order of interpolation.
However, for signals sampled at Nyquist rate (or close to), zero-order interpolation already leads to good results [5]. A detailed error analysis of interpolation schemes can be found in [5,20,37,40]. Further comments can be found in [3].
be an important tool with several applications, such as biomedical image compression, OFDM/CDMA systems, and ADSL transceivers. Searching the literature, no mention about a possible "Arithmetic Hartley Transform" to compute the DHT was found.
In this section, a condensation of the main results of the Arithmetic Hartley Transform is outlined. The method used to define the AHT turned out to furnish a new insight into the arithmetic transform. In particular, the role of interpolation is clarified. Additionally, it is mathematically shown that interpolation is a pivotal issue in arithmetic transforms. Indeed it determines the transform.
A new approach to arithmetic transform is adopted. Instead of considering uniformly sampled points extracted from a continuous signal v(t), the AHT is based on the purely discrete signal. Thus, the starting point of the development is the discrete transform definition, not the series expansion, as it was done in the development of the AFT algorithm. This approach is philosophically appealing, since in a final analysis a discrete transform relates two set of points, not continuous functions.
Let v be an N -dimensional vector with real elements. The DHT establishes a pair denoted by where cas x cos x + sin x is Hartley's "cosine and sine" kernel. The inverse discrete Hartley transform is then [6a] Proof: It follows directly from Lemma 1.
Similarly to the AFT theory, it is necessary to define averages S k , calculated from the time-domain elements. The averages are computed by This definition requires fractional index sampling. Analogously to the AFT methods, this fact seems to make further considerations impracticable, since only integer index samples are available. This subtle question is to be treated in the sequel. For now, we assume that the fractional index sample are somehow available.
An application of the inverse Hartley transform (Equation 42) in Equation 44 offered: Rearranging the summation order, simplifying, and calling Lemma 3, it yielded: For simplicity and without loss of generality, consider a signal v with zero mean value, i.e., 1 Clearly, this consideration has no influence on the values of V k , k = 0. An application of the modified Möbius inversion formula for finite series [5] is sufficient to obtain the final theorem to derive the Arithmetic Hartley Transform. According to Theorem 1, the result below follows.
To illustrate the application of above theorem, consider an 8-point DHT. Using Theorem 7, Equation 48, we obtain: The component V 0 = V 8 can be computed directly from the given samples, since it represents the mean value of the Figure 2 shows a diagram for this computation.
The above theorem and equations completely specified how to compute the discrete Hartley spectrum. Additionally, the inverse transformation can also be established. The following result is straightforward.

Corollary 1 Inverse discrete Hartley transform components can be computed by
where σ i ). A question arises: since the equations are the same, which spectrum is actually being evaluated? Fourier or Hartley spectrum? A clear understanding of underlying arithmetic transform mechanisms will be possible in the next section. Once more the reader is asked to put this question aside for a while, allowing further developments to be derived.
To sum it up, at this point two major questions are accumulated: (i) How to handle with fractional indexes? and (ii) How does same formulae result in different spectra? Interestingly, both questions had the same answer.
The arithmetic transform algorithm can be summarized in four major steps: In the rest of this paper, the step two is addressed. In the sequel, a mathematical method, explaining the importance of the interpolation process in the arithmetic algorithms, is derived.

Interpolation
Arithmetic transform theory usually prescribes zero-or first-order interpolation for the computation of spectral approximations [5,6,38]. In this section, it is shown that an interpolation process based on the known components (integer index samples) characterizes the definition of the fractional index components, v r , r ∈ N. This analysis allows a more encompassing perception of the interpolation mechanisms and gives mathematical tools for establishing validation constraints to such interpolation process. In addition, brief comments on the trade-off between accuracy and computational cost required by interpolation process close the section.

Ideal Interpolation
What does a fractional index discrete signal component really mean? The value of v r for a noninteger value r, r ∈ N, can be computed by Defining the Hartley weighting function by the value of the signal at fractional indexes can be found utilizing an N -order interpolation expressed by: It is clear that each transform kernel can be associated to a different weighting function. Consequently, a different interpolation process for each weighting function is required. In the arithmetic transform formalism, the difference from one transform to another resides in its interpolation process.
The weighting functions satisfy N −1 i=0 w i (r) = 1. If r is an integer number, then the orthogonality properties of cas(·) function [6a] make w r (r) = 1 and w i (r) = 0 (∀i = r). Therefore, no interpolation is needed.
After some trigonometrical manipulation, the interpolation weights for several kernels is expressed by closed formulae. Let Sa(·) be the sampling function,

Hartley kernel
To exemplify, Figure 3 shows two weighting functions used to compute v 10.1 and v 10.5 in the arithmetic Hartley transform. Figure 4 shows the weighting profile for the Fourier cosine and Hartley kernels. These functions are calculated by closed formulae.
With this proposition, the mathematical description of the AHT algorithm is completed. The derived formulae furnish the exact value of the spectral components. On the other hand, the computational complexity of the ideal interpolation implementation is similar to the direct implementation, i.e., computing the transform by its plain definition: To address this issue, a non-ideal interpolation scheme is proposed.

Non-Ideal Interpolation
According to the index generation (m N k ), the number R of points that require interpolation is upper bounded by R ≤ d∤N d − 1. This sum represents the number of samples with fractional index. Consequently, this approach can be attractive for large non-prime blocklength N with great number of factors, because it would require a smaller number of interpolations.  The maximum values are achieved at i = 0 or i = N/2 = 8 (central peak) and they corresponds to the unity. The value of the local maxima are exactly 1/2, occurring when r is integer. Note that each curve is almost null everywhere, except at the vicinities of i ≈ r or i ≈ N − r.
The next task is to find simpler formulae for the weighting functions, assuming large blocklength. Instead of using the exact weighting functions, the limit N→ ∞ is examined and used to derive asymptotic approximations of the weighting function.
In terms of the Hilbert transform and the Sa(·) function, the asymptotic weighting function for Hartley kernel is given bŷ or alternatively, where Hil denotes the Hilbert transform, Ca(x) cos x x , for x = 0, is the co-sampling function and δ(x) is the Dirac impulse.
Zero-order Interpolation. Rounding the fractional index provides the zero-order interpolation. The estimated (interpolated) signalv j is then expressed byv j = v [j] , where [·] is a function which rounds off its argument to its nearest integer. Examining the asymptotic behavior of the weighting function for Fourier cosine kernel, we derive the following results: Under the above assumptions, the Equation 52 furnishes: Thus, for even signal (v k = v N −k ) , the approximated value of the interpolated sample is roughly given byv r ≈ v [r] . It is straightforward to see that the influence of odd part of the signal vanishes in the zero-order interpolation. Besides zero-order interpolation is "blind" to the odd component of a signal. In fact, as show by the set of Equations 56, zero-order interpolation is an (indeed good) approximation to the cosine asymptotic weighting function. This puts some light on the role of the zero-order interpolations and its relation to the earliest versions of the arithmetic Fourier transform, which can not analyze odd periodic signals.
Zero-order interpolation has been intuitively employed in previous work by Tufts, Reed et alli [2,5,6]. Hsu, in his Ph.D. dissertation, derives an analysis of first-order interpolation effect [38]. Interpolation Order. Let M m be a set with the m (< N ) most significant coefficients w i (r). For zero-order interpolation, m = 1. Increasing the value of m, the interpolation process would gradually be improved, because more coefficients would be retained. Proceeding in this way, the following calculation performs a non-ideal interpolation: where η j∈Mm w j (r) is a normalization factor. Figure 5 presents a 32-point discrete Hartley transform of a particular signal computed by the plain definition of the DHT and by the arithmetic transform method using m = 2.

Conclusions
This paper supplied a short survey on arithmetic transforms. The arithmetic Fourier transform is explained and its three main 1-D versions were described. Furthermore, some comments on implementation, challenging points, and advantages of the arithmetic transforms were discussed via simple examples.
In this paper, the introduction of the AHT emphasized the key point of the arithmetic transforms: the interpolation process. It was shown that the fundamental equations of the arithmetic transform algorithms were essentially the same (regardless the kernel). This property could open path to the implementation of "universal transformers". In this type of construction, the circuitry would be the same for several transforms, except the interpolation module. A different interpolation module would reflect different transform (Fourier, Hartley, Fourier cosine). Finally, this paper could be taken as a starting point to those who want to investigate arithmetic transforms.

Acknowledgments
The first author would like to thank Emeritus Professor Dr. Irving S. Reed, University of Southern California, for kindly sending him a copy of Chin-Chi Hsu's Ph.D. thesis [38].