A Factorization Scheme for Some Discrete Hartley Transform Matrices

Discrete transforms such as the discrete Fourier transform (DFT) and the discrete Hartley transform (DHT) are important tools in numerical analysis. The successful application of transform techniques relies on the existence of efficient fast transforms. In this paper some fast algorithms are derived. The theoretical lower bound on the multiplicative complexity for the DFT/DHT are achieved. The approach is based on the factorization of DHT matrices. Algorithms for short blocklengths such as $N \in \{3, 5, 6, 12, 24 \}$ are presented.


INTRODUCTION
Discrete transforms defined over finite or infinite fields have been playing a relevant role in numerical analysis. A striking example is the discrete Fourier Transform, which has found applications in several areas. Another relevant example concerns the discrete Hartley transform (DHT) [1], the discrete version of the integral transform introduced by Hartley in [2]. Besides its numerical side appropriateness, the DHT has proven over the years to be a powerful tool [3,4,5].
A decisive factor for applications of the DFT has been the existence of fast transforms (FT) for computing it [6]. Fast Hartley transforms also exist and are deeply connected to the DHT applications [7,8]. Recent promising applications of discrete transforms concern the use of finite field Hartley transforms [9] to design digital multiplex systems, efficient multiple access systems [10] and multilevel spread spectrum sequences [11].
Discrete transforms presenting a low multiplicative complexity have been an object of interest for a long time. Very efficient algorithms such as the Prime Factor Algorithm (PFA) or Winograd Fourier Transform Algorithm (WFTA) have also been used [12,13]. The minimal multiplicative complexity, µ(·), of the one-dimensional DFT for all possible sequence lengths, N, can be computed by converting the DFT into a set of multi-dimensional cyclic convolutions. A lower bound on the multiplicative complexity of a DFT is given in [14,Theorem 5.4,p. 98]. The values of µ DFT (N) for some short blocklengths are given in Table 1 (some local minima of µ).
In this paper, some FTs are presented, which meet the minimal multiplicative complexity. There is a simple relationship between the DHT and the DFT of a given real discrete signal v.
Therefore, a fast algorithm for the DHT is also a fast algorithm for the DFT and vice-versa [14, Corollary 6.9].
Besides being a real transform, the DHT is also involutionary, i.e., the kernel of the inverse transform is the same as the one of the direct transform (self-inverse transform). Since the DHT is a more symmetrical version of a discrete transform, this symmetry is exploited so as to derive a FT that requires the minimal number of real floating point multiplications.

COMPUTING THE 3-POINT DHT
Let Note that the irrational elements of H 3 have the same decimal part, i.e., except from their integer part, they have the same absolute value (see that ( . So let us make the following decomposition: Since the 2 nd and 3 rd columns of the new matrix H ′ 3 have the same elements (taken in absolute values), we can consider new variables v 1 + v 2 and v 1 − v 2 . Thus, this substitution yields the following matrix equation: So the transform can be expressed by: Observe that the new matrix H ′′ 3 can be splitted in two new matrices, as shown below: Joining the above equations in a single statement, we have that: where a = ( One can recognize the pre-addition matrix A, the multiplication matrix B and the post-addition matrix C. This algorithm introduces a new kind of matrices denoted by L. We will name them "layer matrix". As the notation was explained, we can now express the entire algorithm compactly by the following equation: This algorithms has only one nontrivial multiplication and 7 additions. In Figure 1, there is a schematic diagram for this transform.

COMPUTING THE 5-POINT DHT
Let v H ←→ V be a 5-point DHT pair. The corresponding matrix formulation is now V = H 5 v, where where a = 1/4 We can combine the 2 nd and the 4 th columns as well as the 3 rd and the 5 th ones using a Hadamard transform unit of length 2 (a butterfly). As a result of the process of combining columns, we achieve the following matrix factorization: where A 1 is the pre-addition matrix, B is the multiplication matrix and C i are the post-addition matrices. This matrices are detailed below: Now let us work in the matrixÃ 2 . Note that it contains multiplicative elements, namely e and f , and four additions. We can go further and factorize this matrix in such a way that purely multiplicative and additive matrices appear. This can be done by the following method: Thus, in Equation 10 one should replace the matrixÃ 2 by its factorization. The full decomposition of the original transform matrix H 5 is then achieved.
The arithmetic complexity of this algorithm is 3 multiplications and 17 additions. The schematic diagram is depicted in Figure 2.
Using the Hadamard transform to combine the 1 st and the 4 th columns, the 2 nd and the 5 th columns, and finally, the where A 1 is the first pre-addition matrix. The matrix A 1 is detailed below: The form of the matrix A 1 is the same for all transforms of even blocklength. The explanation to this fact is given by the lemma below.

Lemma 1 The pre-addition matrix A 1 of a Hadamard decomposition of a even blocklength DHT has the following construction:
where Had 2 is the Hadamard matrix, ⊗ is the direct product and I N  Using the same strategy described in the 3-point DHT algorithm, we can take aside the integer part of some elements of the matrix H ′ 6 . This procedure yields to a new more "balanced" matrix. These steps are represented by the following equation: Carrying out the procedure of combining columns which "agree", we will have the next pre-addition matrix A 2 : This makes the matrix H ′ 6 be written as H ′ 6 = H ′′ 6 · A 2 , as seen in this equation:  Now observe that the factorization of H ′′′ 6 yields the multiplication matrix, B, and the post-additions matrix, C.
where a = ( We have then completed the algorithm, and it can be represented by: This algorithm has two multiplications and 20 additions and is depicted in Figure 3.

COMPUTING THE 12-AND THE 24-POINT DHT
The procedure used to derive the 3-and 6-point DHT fast transforms can be extended to other blocklengths, such as 12 and 24. We derived these algorithms and achieved the arithmetic complexity showed in Table 2.
In Figure 4, we see a diagram of the 24-point DHT fast transform, where the shorter transforms (3-, 6-, 12-point) are embedded.
The algorithms proposed so far can be described in a general framework according to the following proposition.

Proposition 1
The DHT decomposition has the following general formulation V = (C n B n A n + L n−1 )C n−1 B n−1 A n−1 · · · + L 2 C 2 B 2 A 2 + L 1 C 1 where n is the number of "layers" in the decomposition.

CONCLUSIONS
Short blocklength DHT fast algorithms that achieve the lower bound on the multiplicative complexity were derived.
Low values for additive complexity were also found. These algorithms can be implemented in digital signal processors capable of low-power consumption.

ACKNOWLEDGMENTS
This work was partially supported by CNPq and CAPES.