The Arithmetic Cosine Transform: Exact and Approximate Algorithms

The Arithmetic Cosine Transform: Exact and Approximate Algorithms

R. J. Cintra R. J. Cintra is with the Signal Processing Group, Departamento de Estatística, Universidade Federal de Pernambuco. This work was done while visiting the Department of Electrical and Computer Engineering, University of Calgary, Calgary, AB, Canada. E-mail: rjdsc@stat.ufpe.org,    V. S. Dimitrov V. S. Dimitrov is with the Advanced Technology Information Processing Systems (ATIPS) Laboratory, University of Calgary, Calgary, AB, Canada. This work was done during second author’s sabbatical leave at Nanyang Technological University, Singapore. Email: vdvsd103@gmail.com
Abstract

In this paper, we introduce a new class of transform method — the arithmetic cosine transform (ACT). We provide the central mathematical properties of the ACT, necessary in designing efficient and accurate implementations of the new transform method. The key mathematical tools used in the paper come from analytic number theory, in particular the properties of the Riemann zeta function. Additionally, we demonstrate that an exact signal interpolation is achievable for any block-length. Approximate calculations were also considered. The numerical examples provided show the potential of the ACT for various digital signal processing applications.

Keywords Discrete cosine transform, arithmetic transform algorithms, nonuniform sampling

1 Introduction

The arithmetic Fourier transform (AFT) has been emerged in 1988 as a signal processing tool [1]. In a broad sense, the AFT is a number-theoretic algorithm for spectrum evaluation. One of the main features of the AFT is its virtually multiplication-free nature. Although initial versions of the AFT procedure could only compute the real part of DFT [1, 2], additional improvements due to Reed et al. expanded the AFT methodology and the DFT could be fully evaluated [3, 4].

Arithmetic methods lack a larger degree of adoption mainly because data is required to be sampled nonuniformly. In fact, the AFT algorithm imposes a sampling scheme that can be related to Farey fractions [5]. This can be a sensitive issue since incoming signals are usually uniformly sampled. Therefore, in order to obtain the necessary nonuniform samples, literature essentially suggests two methods: (i) signal oversampling or (ii) sample interpolation. On the one hand, oversampling can be a major drawback because it requires an above Nyquist rate sampling. Thus, this is often discarded as an option [3]. On the other hand, to obtain nonuniform samples from uniformly sampled data, interpolation methods have been considered. Archived literature lacks an exact procedure to interpolate nonuniform samples obtained from arithmetic transform methods. Therefore, in a simplistic way, zero-order approximations and linear interpolation methods have been considered [6]. Although not furnishing exact computations, these crude interpolation methods could attain acceptable approximations, when large block-lengths were considered [4]. However, for small block-lengths, the implied interpolation errors could be large enough to totally preclude a meaningful computation.

Nevertheless, recent existing works have applied the AFT algorithm in a number of ways. In [7], the AFT is considered as an alternative to the Goertzel algorithm for the single spectral component evaluation. The AFT could also provide frequency domain testing units for built-in self-test routines with reduced hardware overhead [8]. Additionally, hardware considerations have been directed to enhance the associated nonuniform sampling of the AFT [9]. Finally, the AFT has been considered as a tool for DCT evaluation [10].

In all above mentioned applications, the AFT procedure as described in [1, 3] was considered. In particular, even when the DCT was required as in [10], it was obtained by means of the AFT [2]. Indeed, the DFT spectrum can be mapped into the DCT spectrum at the cost of extra computations.

The goal of this paper is twofold. First, we aim to introduce a class of purely arithmetic algorithms for the DCT, termed arithmetic cosine transforms (ACT). Such new methods are the main contribution of this work. To accomplish this objective, we examine arithmetic averages based on the existing arithmetic Fourier transform. Then, we advance new arithmetic averages specially tailored for the DCT computation.

Second, an exact interpolation method for the ACT is sought. Such exact interpolation could provide solid mathematical grounds for arithmetic transform methods, constituting a theoretical advance in the area. In fact, any arithmetic transform method could only furnish an exact computation if a precise interpolation method could be devised. In particular, we discuss an interpolation scheme that could provide an exact spectrum evaluation, even when short block-lengths are considered. Moreover, an asymptotic analysis for the proposed interpolation approach is elaborated and a heuristic approximation algorithm is devised.

The article unfolds as follows. In Section 2, we examine some of the existing tools employed in the AFT theory and propose a new algorithm for the DCT. Section 3 is devoted to the interpolation issues that arise from the introduced ACT. In Section 4, the discussed arithmetic transform is formatted in terms of matrix representation. Conclusions and final remarks are given in Section 5.

2 Arithmetic cosine transform

2.1 Preliminary concepts

Among the several types of DCTs, we separated the DCT-II, which can be regarded as the most employed form [11]. This transformation relates two -dimensional vectors, and , according to the following pair of relations [11]

 Vk =√2NαkN−1∑n=0vncos(πk(n+1/2)N),k=0,1,…,N−1, (1)
 vn =√2NN−1∑k=0αkVkcos(πk(n+1/2)N),n=0,1,…,N−1,

where the coefficients , are given by

 αk={1/√2,if {}k=0,1,otherwise.

Hereafter, we refer to this transformation simply as DCT.

Before defining and examining the ACT, some fundamental tools are necessary. In what follows, we introduce the generalized Möbius inversion formula tailored for finite series. Additionally, some useful identities for trigonometric sums are also presented. We state these preambulary results below.

Theorem 1 (Generalized Möbius Inversion Formula for Finite Series)

Let be a sequence (e.g., signal samples) such that it is nonnull for and null for . Admit another sequence defined as

 gn=⌊N/n⌋∑k=1akfkn,

where is a sequence of scalar coefficients and is the floor function. Then,

 fn=⌊N/n⌋∑m=1bmgmn,

where is the Dirichlet inverse sequence of , given that it exists [12].

Proof:  It follows by applying the arguments described in the proof of Theorem 3 from [4, p. 469] into the Möbius inversion formula as shown in [13, p. 556].

If we consider the unitary sequence , for , then the associated Dirichlet inverse sequence is the Möbius sequence , for The Möbius function  [14] is defined over the positive integers and is given by

 μ(n) =⎧⎨⎩1,if {}n=1,(−1)q,if {}ncan be factorized into {}qdistinct primes,0,if {}nis divisible by a square number.

In this particular case, arithmetic transform literature simply refers to the above theorem as the Möbius inversion formula for finite series [4, 3].

The following lemmas are pivotal for subsequent developments. They link usual trigonometric functions to number theoretic behavior functions, a connection that is not always trivial [15, 16].

Lemma 1 (Reed’s Lemma)

Let and be positive integers. Then,

and

 k−1∑m=0sin(2mk′kπ)=0.

We amplified Lemma 1 with the next two results.

Lemma 2

Let and be positive integers. Then,

 2k−1∑m=0cos(πmk′k)={2k,if {}2k|k′,0,otherwise.
Proposition 1

Let and be integers and be a real quantity. Then,

 k−1∑m=0cos(2πk′k(m+α))=⎧⎪⎨⎪⎩k,if {}k′=0,kcos(2πk′kα),if {}k|k′, {}k′≠0,0,otherwise.

Proof:  It follows from usual trigonometric manipulations and an application of Lemma 1.

2.2 Arithmetic averages

Based on existing definitions and concepts inherited from the AFT theory, our initial attempt to define an arithmetic transform procedure for the DCT is detailed below. Thus, we consider the following definition for AFT-like averages.

Definition 1 (kth AFT-like average)

Let the th average be defined as

 Sk≜12k2k−1∑m=0vmNk−12,k=1,2,…,N−1.

Except for the -shift, which is naturally imposed by the DCT kernel, the above concept was previously described in [17, 4, 3] as a framework for the AFT.

Let us investigate the consequences of this definition. Relaxing the integer index constraint of the DCT formulae and substituting the inverse DCT relation into the above described th average, we must have

 Sk =12k2k−1∑m=0vmNk−12 =12k2k−1∑m=0√2NN−1∑k′=0αk′Vk′cos⎛⎝πk′(mNk−12+12)N⎞⎠ =√2N12k2k−1∑m=0(γV0+N−1∑k′=0Vk′cos(πmk′k)) =√2NγV0+√2N12k2k−1∑m=0N−1∑k′=0Vk′cos(πmk′k),k=1,2,…,N−1,

where . By interchanging the summations, it holds that

 Sk =√2NγV0+√2N12kN−1∑k′=0Vk′2k−1∑m=0cos(πmk′k),k=1,2,…,N−1. (2)

Invoking Lemma 2, we have

 Sk =√2NγV0+√2N12kN−1∑k′=0Vk′{2k,if % {}2k|k′,0,otherwise,} =√2NγV0+√2NN−1∑k′=0Vk′{1,if {}2k|k′,0,otherwise,},k=1,2,…,N−1.

Let . Then,

 Sk =√2NγV0+√2N⌊N−12k⌋∑s=0V2sk,k=1,2,…,N−1.

Without any loss of generality, let us assume that input signals have a null mean value. Therefore, it follows that is zero. This assumption does not affect the value of the remaining components of the DCT spectrum. Thus, we find that

 Sk

In order to invert this last expression, we must obtain the Dirichlet inverse of the sequence , as required by Theorem 1. However, a sequence can only admit the Dirichlet inverse if and only if its first element is nonnull [12, p. 18]. This is clearly not the case for the particular sequence under examination.

We conclude that the derivation of an arithmetic method for the DCT requires more than simply reusing the AFT averages. For the DCT computation, specifically tailored averages should be considered. In the next section, such specific derivation is sought. However, the mathematical manipulations and arguments shown above prove to be a useful framework. We show next that the proposition of DCT specific averages can be greatly benefited from the above exposition.

2.3 ACT averages

In order to derive the ACT, we adjusted the AFT averages as suggested in the following definition.

Definition 2 (kth ACT average)

Let the th ACT average be defined as

 Sk≜1kk−1∑m=0v2(m+β)Nk−12,k=1,2,…,N−1,

where is a fixed real number.

Considering algebraic manipulations similar to the ones that has led us from Definition 1 to equation (2) and admitting that the input signal has null mean, mutatis mutandis, we obtain

 Sk =√2N1kN−1∑k′=1Vk′k−1∑m=0cos(2π(m+β)k′k),k=1,2,…,N−1.

Thus, invoking Proposition 1, we have

 Sk =√2NN−1∑k′=1Vk′{cos(2πk′kβ),if {}k|k′,0,otherwise,},k=1,2,…,N−1.

Performing the substitution , it follows that

 Sk

Observe that the above expression is suitable for the application of the generalized Möbius inversion formula for finite series. Under the notation of Theorem 1, we recognize ,

Incidentally, not always the Dirichlet inverse of is well-defined. Only when , the existence of the Dirichlet inverse can be considered [12]. Thus, we must impose as a necessary condition for the derivation of the ACT procedure. This issue is precisely the point that Definition 1 misses.

However, finding the Dirichlet inverse of , say for arbitrary values of can be a cumbersome maneuver. Therefore, we separated two particular useful cases: (i) and (ii) . Notice that leads to a non-invertible sequence, since this makes .

For , we have the unitary sequence and , for This is usually the situation addressed in standard AFT analysis. On the other hand, setting yields , In this case, the Dirichlet inverse is not immediately recognized, but is can be obtained analytically. In the Appendix, we derive the sought Dirichlet inverse, which is given by

 bn={−μ(n),if {}nis odd,−2m−1μ(2−mn),if {}n=2ms, where {}sis odd. (3)

The first 32 terms of are listed below:

 −1,−1,1,−2,1,1,1,−4,0,1,1,2,1,1,−1,−8,1,0,1,2,−1,1,1,4,0,1,0,2,1,−1,1,−16.

In the context of digital signal processing this is a potentially useful sequence, since multiplying a given number by a power-of-two can be implemented by simple bit shift operations, which possess a low computational complexity.

Regardless the choice of , to invert the following expression given by

 Sk

we may consider an auxiliary sequence given by for all . Thus, a direct application of the generalized Möbius inversion formula for finite series just tells us that

 Gk =⌊N−1k⌋∑l=1blSkl,k=1,2,…,N−1,

where is the Dirichlet inverse of . Finally, undoing the auxiliary substitution we can write

 Vk (4)

Depending on whether is the unitary sequence or the alternate sequence , the inverse sequence is the Möbius sequence or the sequence described in (3), respectively.

Recognizing that unitary and Möbius sequences constitute a simpler pair of Dirichlet inverse sequences, in the rest of this paper we focus our attention to the case . However, all ensuing developments encompass the proposed alternative formulation () without significant modifications. Nevertheless, the case can furnish a link between Definitions 1 and 2. Indeed, for , we have

 Sk =1kk−1∑m=0v2(m+12)Nk−12 =1kk−1∑m=0v(2m+1)Nk−12 =1k2k−1∑m=1{}moddvmNk−12,k=1,2,…,N−1.

This last summation corresponds to the average shown in Definition 1, when only odd values of the dummy index are considered. In other words, if the samples required by Definition 1 are submitted to a downsampling by a factor of 2, then Definition 1 collapses to Definition 2 for . In a sense, this observation establishes a connection between both formalisms.

As far as the computational complexity of the Möbius inversion formulae are concerned, we can provide the following probabilistic reasoning. The probability that a randomly chosen integer is not divisible by a perfect square is  [12, p. 4]. Therefore, 61% of the values of the Möbius function are zeros; meaning that the computation of (or ) requires additions/subtractions on average.

Up to this point, we assumed that input signals have null mean. Now, let us remove this restriction. Let be the mean value of the input signal. Thus, an arbitrary signal can be converted into a null mean signal by simply subtracting the quantity . Therefore, the th ACT averages associated to can be manipulated as follows:

 S′k ≜1kk−1∑m=0(v2mNk−12−¯v) =1k(k−1∑m=0v2mNk−12)−¯v =Sk−¯v,k=1,2,…,N−1.

Consequently, equation (4) can be rearranged, yielding

 Vk =√N2⎛⎜ ⎜⎝⌊N−1k⌋∑l=1μ(l)Skl−⌊N−1k⌋∑l=1μ(l)¯v⎞⎟ ⎟⎠,k=1,2,…,N−1.

Considering the Mertens function , which is defined as  [18, p. 272], we obtain a more compact expression as

 Vk =√N2⎛⎜ ⎜⎝⌊N−1k⌋∑l=1μ(l)Skl⎞⎟ ⎟⎠−√N2¯vM(⌊N−1k⌋),k=1,2,…,N−1. (5)

The zeroth component is expressed separately by .

3 ACT interpolation

3.1 Fractional indices manipulation

Now we must address the problem of handling fractional index samples. Usually, arithmetic transform literature addresses this issue by prescribing the utilization of zero- or first-order interpolation methods. We argue that such a naive approach is not the proper way of interpolating.

Again, let us relax the integer index constraint of the DCT formulae. Considering a possibly noninteger quantity , the following construction is derived given by

 vr =√2NN−1∑k=0αkVkcos(πk(r+1/2)N).

Taking into account the direct transformation formula, we obtain

 vr =√2NN−1∑k=0αk(√2NαkN−1∑n=0vncos(πk(n+1/2)N))cos(πk(r+1/2)N).

Inverting the summation order yields

 vr =2NN−1∑n=0vnN−1∑k=0α2kcos(πk(n+1/2)N)cos(πk(r+1/2)N).

Then, let us define the ACT weighting function as

 wn(r) ≜2NN−1∑k=0α2kcos(πk(n+1/2)N)cos(πk(r+1/2)N) =−1N+2NN−1∑k=0cos(πk(n+1/2)N)cos(πk(r+1/2)N),n=0,1,…,N−1.

Thus, the samples associated to fractional indices can be obtained after a linear combination of the available uniformly obtained samples. Hence,

 vr =N−1∑n=0wn(r)vn. (6)

Notice also that if is integer, we have that

 wn(r)={1,if {}n=r,0,otherwise.

This fact stems from the orthogonality properties of the transformation kernel.

We further investigate the suggested weighting function. The next proposition states that the discussed weighting functions are indeed inherently normalized.

Proposition 2

The ACT weighting function satisfies the following summation formula

 N−1∑n=0wn(r)=1.

Proof:  Usual trigonometric manipulations furnish the derivations below:

 N−1∑n=0wn(r) =N−1∑n=0(−1N+2NN−1∑k=0cos(πk(n+1/2)N)cos(πk(r+1/2)N)) =−1+2NN−1∑k=0cos(πk(r+1/2)N)N−1∑n=0cos(πk(n+1/2)N).

However, for each , the inner summation can be expanded as follows:

 N−1∑n=0cos(πk(n+1/2)N) =cos(πk2N)N−1∑n=0cos(2π(k/2)nN)−sin(πk2N)N−1∑n=0sin(2π(k/2)nN).

Since dummy index runs from to , we have that never divides . Therefore, applying Lemma 1, we obtain that

 N−1∑n=0cos(πk(n+1/2)N) =cos(πk2N){N,if {}k=0,0,otherwise,} ={N,if {}k=0,0,otherwise.

Finally, returning to the previous double summation, we establish that

 N−1∑n=0wn(r) =−1+2NN−1∑k=0cos(πk(r+1/2)N){N,if {}k=0,0,otherwise,} =1.

Regardless of the considered block-length, the use of the ACT interpolation results in an exact calculation of the DCT spectrum. Indeed, no approximation was considered in any of our arguments.

3.2 An example: the 8-point ACT

To illustrate the ACT structure and its interpolation scheme, we devised an example. Intentionally, we selected the short block-length transformation furnished by the -point DCT. This particular block-length is widely adopted in several image and video coding standards, such as JPEG, MPEG-1, MPEG-2, H.261, and H.263 [19]. The -point DCT is also subject to an extensive analysis in [11]. In the following, we set .

The first step of the ACT procedure consists of the identification of the necessary interpolation points. According to Definition 2, these points are given by for and . Therefore, we find the following fractional indices: . Fig. 1 depicts a block diagram of the full algorithm. The ACT interpolation block calculates the required samples according to the discussed interpolation procedure.

The fractional index samples are then employed to obtain the ACT averages. The block of ACT averages simply implements the following set of equations:

 S1 =v−12 2S2 =S1+v152 3S3 =S1+2v296 4S4 =S2+2v72 5S5 =S1+2v2710+2v5910 6S6 =S2+S3+2v136 7S7 =S1+2v2514+2v5714+2v8914.

Finally, the ACT averages are combined with respect to the Möbius function (cf. (4)). The resulting calculation involves no approximations and furnishes the exact DCT spectrum.

3.3 Asymptotic analysis

In this section, we closely examine the weighting function required by the interpolation procedure. We aim to propose simpler interpolation expressions allowing an efficient computation of the ACT.

First, we note that the ACT weighting function can be formulated in closed form as detailed below. In fact, invoking elementary trigonometric identities, we can establish the following relations:

 wn(r) =−1N+2NN−1∑k=0cos(πk(n+1/2)N)cos(πk(r+1/2)N) =−1N+1NN−1∑k=0(cos(πk(n+r+1)N)+cos(πk(n−r)N)) =−1N+1NN−1∑k=0cos(πk(n+r+1)N)+1NN−1∑k=0cos(πk(n−r)N),n=0,1,…,N−1.

The above trigonometric summations can be given in terms of the Dirichlet kernel [20, p. 312]. Therefore, it holds that

 wn(r) =−1N+1N(12+12DN−1(πN(n+r+1)))+1N(12+12DN−1(πN(n−r))) =12N(DN−1(πN(n+r+1))+DN−1(πN(n−r))),

where denotes the Dirichlet kernel.

The similarity between the Dirichlet kernel and the function is apparent. Indeed, as , these functions are interchangeable in several situations [21]. However, even for small values of such an approximation is good [21, p. 180]. For instance, taking and centered functions at the origin, it follows that the mean square error (MSE) of implied approximation is less than over the interval . Thus, the discussed weighting function assumes a limiting form given by the sum of two translated sampling functions as

 limN→∞wn(r)=sinc(n+r+1)+sinc(n−r),

where .

Notice that this asymptotic expression connects the ACT interpolation to the function interpolation. Signal interpolation according to the function can be efficiently implemented in the time domain [22, 23]. Additionally, fractional delay FIR filtering methods offer another computational approach [24, 25].

The limiting form of leads us to draw some additional conclusions on its asymptotic behavior. Let denote the nearest integer function as implemented in C or Matlab programming languages and be a fractional interpolation point. We observed that, when , the asymptotic weighting function is essentially governed by a single function. The role of could be neglect, since its argument would be large enough. Indeed, this function approaches to zero according to . Thus, in this case we say that

 limN→∞wn(r)≈sinc(n−r).

Fig. 2 shows some plots of for several values of . These plots intuitively indicate the use of a linear approximation for the main lobe of the function. Additionally, we assumed that the effect of the remaining lobes is negligible. We also observed that, for , both functions are relevant since they overlap each other. This last situation was treated separately.

In terms of the above discussion, we derived an empirical approximation procedure that considers at most two uniform samples to render an interpolated sample. In Fig. 3, algorithmic details of the suggested approximate interpolation are given. In the proposed heuristic algorithm, we admit an auxiliary quantity , an error tolerance , and a scaling factor , for . This last quantity was found according to a standard linear fitting procedure.

Considering a conservative choice of , we employed the proposed heuristics to calculate the approximate DCT spectra of 256 randomly generated signal vectors. A short block-length was deliberately selected. The elements of the input signal were chosen to be distributed according to a standard uniform distribution. The resulting average MSE due to the discussed approximation was as low as . This figure is comparable to the MSE associated to some integer approximation algorithms for the -point DCT (e.g., integer cosine transform, -matrix transform) as detailed in [11].

4 Matrix-vector representation

Further insight on the nature of the ACT can be obtained when the previous constructions are represented in matrix-vector form. Let the input signal and its associated DCT spectrum be denoted by column vectors and , respectively. Additionally, consider the DCT matrix , whose elements are defined according to (1): , for . The above quantities are related by and  [11, p. 41].

In order to render the ACT structure in matrix-vector form, we need to introduce some special matrices.

Definition 3 (Möbius matrix)

The element of the -order Möbius matrix is given by , whenever divides , for . Otherwise, it is zero.

By construction, the Möbius matrix is upper triangular with unity diagonal elements. Thus, it is always non-singular and its determinant is unity for any dimension. The inverse of the Möbius matrix can be directly obtained without calling an inversion procedure. The element of is 1, if divides ; otherwise, it is zero. This fact stems from the Möbius inversion formula relations.

Additionally, we consider the extended Möbius matrix defined by

 M′≜diag(1,MN−1)=⎡⎢ ⎢ ⎢ ⎢⎣10⋯00⋮MN−10⎤⎥ ⎥ ⎥ ⎥⎦,

where the operator returns a diagonal matrix.

We may also admit a vector of averages expressed by , where the zeroth average is defined separately as , and the remaining values are the th ACT averages for . Notice that DC component of the spectrum is related to the zeroth average .

The ACT framework can provide a new insight into the DCT spectrum. In fact, equation (5) indicates that the DCT spectrum can be separated into two parts: (i) one due to the Möbius combination of the ACT averages and (ii) another due to the Mertens function. Let these two parts be termed the Möbius and the Mertens parts of the DCT spectrum denoted by and , respectively. Thus, the DCT spectrum can be decomposed as .

Joining the above structures it follows that the vector is related to the Möbius part of the DCT spectrum via the Möbius extended matrix according to

 V1=√N2⋅M′⋅S. (7)

In the light of the proposed ACT interpolation, we can recast the ACT averages in terms of the weighting function. Thus, invoking (6), we have

 Sk =1kk−1∑m=0v2mNk−12

Inverting the order of the summations, we may also obtain

 Sk =N−1∑n=0(1kk−1∑m=0wn(2mNk−12))vn,k=1,…,N−1.

For a fixed , the inner expression in parenthesis is an average of particular weighting values at different interpolating points. This term depends only on and being independent of the input vector . Then, we can separate this expression for better understanding. Let us define the th weighting average as

 Wk,n≜1kk−1∑m=0wn(2mNk−12),k=1,…,N−1, n=0,…,N−1.

These averages give rise to the construction of the following matrix:

 W=[Wk,n],k=1,…,N−1,n=0,…,N−1.

Since the rows of are convex combinations of the weighting function, Proposition 2 indicates that the sum of elements across rows of is equal to 1, i.e., , for . By augmenting the matrix with the inclusion of an extra row and padding, as described below,

 W′=⎡⎢ ⎢ ⎢ ⎢⎣1/N1/N⋯1/N0⋮W0⎤⎥ ⎥ ⎥ ⎥⎦,

we may write the following relation:

 S=α⋅W′⋅v,

where .

An application of this last expression into (7) furnishes the final matrix-vector form for the Möbius part of the considered procedure. Hence,

 V1=√N2⋅M′⋅α⋅W′⋅v.

Effectively, the transformation matrix that relates and is given by

 C1≜√N2⋅M′⋅α⋅W′.

Now considering the spectral part due to the Mertens function in (5), we advance the following additional matrix:

 C2

where is a column vector consisting of unit elements. Of course, if the input signal has null mean, then the Mertens part of the DCT spectrum is always zero: .

In view of the above, we must have

 V =V1+V2 =C1⋅v+C2⋅v =(C1+C2)⋅v.

This manipulation suggests that the DCT matrix can be decomposed as . Notice that when the input signal has null mean, a new formulation for the DCT matrix arises. In this case, we can establish that

 C⋅v=C1⋅v.

This alternative transformation matrix for the DCT can be considered as a starting point to derive new algorithms under the constraint of null mean input signals.

5 Conclusion and Remarks

The main contributions of this paper are (i) a new arithmetic transform and (ii) the elucidation of the arithmetic transform interpolation issues. The introduced arithmetic cosine transform is a number-theoretic based algorithm devoted to the DCT computation. Therefore, DSP applications that require a DCT evaluation are potential candidates for the use of the ACT method.

Differently from the standard AFT analysis, we generalized existing inversion formula, allowing new frameworks for the arithmetic transform theory. Besides the Möbius sequence, we identified alternative adequate sequences for the suggested method.

Moreover, the introduced interpolation method allows the exact computation of the DCT spectrum, even for small block-lengths. This is particularly distinct from the existing AFT algorithms, which (i) inevitably introduce approximation errors and (ii) tend to excel only when large block-lengths are considered. Using arithmetic methods (e.g., AFT) for small block-length transform evaluation was a challenging issue, for which area literature could not furnish adequate solutions until now.

The complexity of the ACT procedure is mainly due to the interpolation step. If the sampling process could be adjusted in order to collect samples natively in an appropriate nonuniform fashion, then the ACT computation would not need any sort of interpolation. Indeed, the DCT would be exactly computed after some few additions/subtractions associated to the Möbius function. Other than that, the arithmetic transform philosophy can concentrate the transform computational complexity into the interpolation block. This is a different perspective for the design of transform procedures. For instance, fast algorithms for fractional delay filtering [24] now receive an additional motivation under the arithmetic cosine transform paradigm.

Finally, the existence of efficient frameworks for bidimensional AFTs suggests a venue for extending the proposed ACT into the bidimensional case [26, 27, 28, 10]. Without much effort, the proposed method could be used as the fundamental building block for a bidimensional arithmetic cosine transform. On the other hand, a dedicate analysis for the bidimensional case could prove to be more appropriate. In particular, an application Feig-Winograd direct approach could avoid row-column methods [29]. This could be a more attractive method for the bidimensional ACT. In any case, the characterization of the exact interpolation for bidimensional signals is an open topic.

Acknowledgment

The first author wishes to express his gratitude to Prof. Hélio M. de Oliveira who provided several stimulating scientific discussions. This work was partially supported by the Department of Foreign Affairs, Trade and Development, Canada.

References

• [1] D. W. Tufts and G. Sadasiv, “The arithmetic Fourier transform,” IEEE ASSP Magazine, vol. 5, no. 1, pp. 13–17, 1988.
• [2] D. W. Tufts, Z. Fan, and Z. Cao, “Image processing and the arithmetic Fourier transform,” in SPIE High Speed Computing II, vol. 1058, Los Angeles, CA, May 1989, pp. 46–53.
• [3] I. S. Reed, M. T. Shih, T. K. Truong, E. Hendon, and D. W. Tufts, “A VLSI architecture for simplified arithmetic Fourier transform algorithm,” IEEE Transactions on Signal Processing, vol. 40, no. 5, pp. 1122–1133, May 1992.
• [4] I. S. Reed, D. W. Tufts, X. Yu, T. K. Truong, M. T. Shih, and X. Yin, “Fourier analysis and signal processing by use of the Möbius inversion formula,” IEEE Transactions on Acoustics, Speech and Signal Processing, vol. ASSP-38, no. 3, pp. 458–470, Mar. 1990.
• [5] N. Tepedelenlioglu, “A note on the computational complexity of the arithmetic Fourier transform,” IEEE Transactions on Acoustics, Speech and Signal Processing, vol. ASSP-37, no. 7, pp. 1146–1147, Jul. 1989.
• [6] B. T. Kelley and V. K. Madisetti, “Efficient VLSI architectures for the arithmetic Fourier transform (AFT),” IEEE Transactions on Signal Processing, vol. 41, no. 1, pp. 365–384, Jan. 1993.
• [7] J. B. Lima, R. M. Campello de Souza, H. M. Oliveira, and M. M. Campello de Souza, Faster DTMF Decoding, ser. Lecture Notes in Computer Science.    Springer-Verlag Berlin / Heidelberg, 2004, vol. 3124/2004, ch. Telecommunications and Networking - ICT 2004, pp. 510–515.
• [8] N. Venkateswaran and K. Bharath, “Frequency domain testing of general purpose processors at the instruction execution level,” in Second IEEE International Workshop on Electronic Design, Test and Applications, Perth, Australia, Jan. 2004, pp. 15–20.
• [9] G. Ray and M. Chen, “Analog to feature conversion,” in IEEE International Conference on Acoustics, Speech and Signal Processing, vol. 2, Honolulu, HI, Apr. 2007, pp. 365–368.
• [10] E. J. Tan, Z. Ignjatovic, and M. F. Bocko, “A CMOS image sensor with focal plane discrete cosine transform computation,” in IEEE International Symposium on Circuits and Systems, New Orleans, LA, May 2007, pp. 2395–2398.
• [11] V. Britanak, P. Yip, and K. R. Rao, Discrete cosine and sine transforms.    Amsterdam: Academic Press, 2007.
• [12] P. T. Bateman and H. G. Diamond, Analytic Number Theory.    New Jersey: World Scientific, 2004.
• [13] R. R. Goldberg and R. S. Varga, “Moebius inversion of Fourier transforms,” Duke Math. Journal, vol. 24, no. 4, pp. 553–560, Dec. 1956.
• [14] T. M. Apostol, Introduction to Analytic Number Theory.    New York: Springer-Verlag, 1984.
• [15] S.-C. Pei and K.-W. Chang, “Odd Ramanujan sums of complex roots of unity,” IEEE Signal Processing Letters, vol. 14, no. 1, pp. 20–23, Jan. 2007.
• [16] V. Laohakosol, P. Ruengsinsub, and N. Pabhapote, “Ramanujan sums via generalized Möbius functions and applications,” International Journal of Mathematics and Mathematical Sciences, 2006.
• [17] G. Fisher, D. Tufts, and R. Unnikrishnan, “VLSI implementation of the arithmetic Fourier transform,” in Proceedings of the 32nd Midwest Symposium on Circuits and Systems, vol. 2, Champaign, IL, Aug. 1989, pp. 800–803.
• [18] M. Schroeder, Number Theory in Science and Communication: With Applications in Cryptography, Physics, Digital Information, Computing, and Self-similarity.    New York: Springer-Verlag, 2008.
• [19] N. Roma and L. Sousa, “Efficient hybrid DCT-domain algorithm for video spatial downscaling,” EURASIP Journal on Advances in Signal Processing, vol. 2007, no. 2, pp. 30–30, 2007.
• [20] S. G. Krantz, Real Analysis and Foundations.    Boca Raton, FL: Chapman & Hall/CRC, 2005.
• [21] G. Bachman, L. Narici, and E. Beckenstein, Fourier and wavelet analysis.    New York: Springer-Verlag, 2000.
• [22] T. Schanze, “Sinc interpolation of discrete periodic signals,” IEEE Transactions on Signal Processing, vol. 43, no. 6, pp. 1502–1503, Jun. 1995.
• [23] S. R. Dooley and A. K. Nandi, “Notes on the interpolation of discrete periodic signals using sinc function related approaches,” IEEE Transactions on Signal Processing, vol. 48, no. 4, pp. 1201–1203, Apr. 2000.
• [24] V. Välimäki and T. I. Laakso, Nonuniform sampling: theory and practice.    New York: Springer-Verlag, 2001, ch. Fractional Delay Filters—Design and Applications.
• [25] T. I. Laakso, V. Valimaki, M. Karjalainen, and U. K. Laine, “Splitting the unit delay [FIR/all pass filters design],” IEEE Signal Processing Magazine, vol. 13, no. 1, pp. 30–60, Jan. 1996.
• [26] X.-J. Ge, N.-X. Chen, and Z.-D. Chen, “Efficient algorithm for 2-D arithmetic Fourier transform,” IEEE Transactions on Signal Processing, vol. 45, no. 8, pp. 2136–2140, Aug. 1997.
• [27] V. G. Atlas, D. G. Atlas, and E. I. Bovbel, “2-D arithmetic Fourier transform using the Bruns method,” IEEE Transactions on Circuits and Systems I: Fundamental Theory and Applications, vol. 44, no. 6, pp. 546–551, Jun. 1997.
• [28] D. W. Tufts and G. Sadasiv, “Arithmetic Fourier transform and adaptive delta modulation: A symbiosis for high speed computation,” in SPIE High Speed Computing, vol. 880, Los Angeles, CA, Jan. 1988, pp. 168–178.
• [29] E. Feig and S. Winograd, “Fast algorithms for the discrete cosine transform,” IEEE Transactions on Signal Processing, vol. 40, no. 9, pp. 2174–2193, Sep. 1992.
• [30] G. H. Hardy, E. M. Wright, R. Heath-Brown, J. Silverman, and A. Wiles, An Introduction to the Theory of Numbers.    New York: Oxford University Press, 2008.

Appendix A The Dirichlet inverse of {(−1)n}

In order to derive the Dirichlet inverse of the sequence , for , let us examine its associated Dirichlet series. A result from the theory of functions [30, p. 337] states that

 ∞∑n=1(−1)nns=−(1−21−s)ζ(s),R(s)>0,

where is the Riemann zeta function. Therefore, we directly have that the closed form of Dirichlet series of , , is

 A(s)=−(