Fast Type-II Hartley Transform Algorithms for Short-Length Input Sequences

Abstract

This paper presents the type-II fast discrete Hartley transform (DHT-II) algorithms for input data sequences of lengths from 2 to 8. The starting point for developing the eight algorithms is the representation of DHT-II as a matrix–vector product. The underlying matrices usually have a good block structure. These matrices must then be successfully factorized to obtain a computational procedure that reduces the number of operations in computing the matrix–vector product. In some cases, it is necessary to pre-decompose the original matrices into submatrices and rearrange the rows and/or columns of the resulting matrices to find the factorizations that would substantially save the arithmetic operations. As a result of applying the pointed transformations, we synthesized the final algorithms with reduced computational complexity. The correctness of the obtained algorithmic solutions was theoretically justified using the rigorous mathematical background of each of them. Then, the complex algorithms were further tested using the MATLAB R2023b software to confirm their performance. Finally, an evaluation of the computational complexity for each obtained solution was compared with the computational complexity of the direct calculation of the matrix–vector product and existing fast DHT-II algorithms. The obtained factorizations of the DHT-II transformation matrices on average reduce the number of additions by 5% and the number of multiplications by 73% compared with the direct calculation of the matrix–vector product.

1. Introduction

Among other orthogonal transforms, the discrete Hartley transform (DHT) is a well-known and well-tested tool in digital signal and image processing [1,2,3]. This transform is used in smart grid applications [4], signal processing [5], digital image watermarking [6,7,8], audio watermarking [9], image compression and encryption [10,11], image enhancement [12], OFDM communication systems [13,14,15], face recognition [16], audio classification [17,18], etc. Its advantage is primarily that it can successfully supplement the discrete Fourier transform in many applications [1,2,3]. For example, in audio classification applications, features extracted from the representation of signals in the frequency domain are typically focused on the magnitude spectral content, while the phase spectral content is ignored [19]. The conventional Fourier phase spectrum is a highly discontinuous function. To extract the features in classification applications, function continuity is required. Thus, the Hartley phase spectrum was introduced as an alternative to the Fourier phase spectrum. The former encapsulates the phase content of the signal more efficiently if the sources of phase discontinuities are detected and overcome [20]. The latter is because the Hartley phase spectrum does not suffer from the phase “wrapping ambiguities” introduced due to the inverse tangent function employed in the Fourier phase spectrum computation [19].

Compared to wavelet analysis, the Hartley transform provides only global frequency information, although the ability to analyze signals in terms of both time and frequency may be obtained with windowing applied. As a result of windowing, a short DHT is obtained [21]. The DHT allows for efficient signal denoising, being a fixed global transform. Moreover, the traditional DHT can be adapted to detect fine signal details by using windows of different sizes [8,22]. Although DHT [10,11] is exploited in image compression infrequently, it is widely applied in audio and image watermarking [6,7,8,22] to enhance robustness against different noise attacks. Also, a significant increase in computing efficiency by saving arithmetic operations and memory storage is achieved [22].

To date, four types of discrete Hartley transforms have been formulated [23]. However, the most popular and most developed transform is still DHT-I [1,10,13,14,15]. The other transforms from this family have been studied much less. One important, but not the main, advantage of discrete orthogonal transforms is the huge number of so-called fast algorithms that have been designed for them, allowing for the acceleration of their implementation [23]. The direct computation of the DHT-II of an N-point sequence requires N² arithmetic operations. Therefore, fast DHT-II algorithms have been introduced to reduce the computational complexity and implementation cost [24,25,26,27,28,29].

It should be noted that the papers dealing with the efficient implementation of the DHT-II transform considered the case of large lengths of input data sequences [24,25,26,27,28,29]. However, short-length DHT-II algorithms are of particular interest, since these transforms are considered typical modules in the synthesis of more complex algorithms [23]. Once constructed, short-length DHT-II algorithms can be successfully applied in various projects to unify the development process of the final product. There are many papers devoted to the efficient implementation of traditional short-length trigonometric transforms, but they mainly concern cosine or sine transforms [30,31,32,33]. Then, the development of fast DHT-II algorithms for short-length input sequences becomes relevant. Below, the current state of the problem is analyzed to select an approach to the development of such fast DHT-II algorithms.

1.1. Related Papers

To develop fast algorithms for discrete orthogonal transforms the strategy of “divide and conquer” is applied [23,34,35,36,37]. This means dividing a large problem into smaller ones. Then, the solution to the big problem can be obtained by collecting the solutions of smaller problems, which are organized according to the dividing procedures. Often, an obtained solution to the big problem is modular and systematic, which is desirable for implementation. Based on the divide and conquer strategy, fast algorithms for discrete orthogonal transforms have been developed as split-radix algorithms and prime factor algorithms [23].

Split-radix algorithms assume that the length of the input data sequence can be expressed as N = q × 2^m, where m is a positive integer and q is an odd integer. Various sequence length values can be represented by setting different values of m and q. For example, the radix-2 algorithm recursively decomposes the length N DHT-II into two N/2-length DHT-IIs. In contrast, the split-radix algorithm decomposes the N-length DHT-II using different radices. In the literature, split-radix algorithms are developed based on decimation-in-time (DIT) and decimation-in-frequency decomposition [23]. The advantages of split-radix algorithms are substantial, saving the number of arithmetic operations and supporting a wide range of sequence lengths.

Despite the advantages of the cost-effective radix algorithms, computational overheads usually occur when the input sequence length is not matched to the assumptions of the available algorithm. Hence, prime factor algorithms were introduced for various sequence lengths, in particular, for sequence lengths other than a power of two. However, prime factor algorithms generally have a more complex computational structure. Specifically, they require irregular index mapping to reorder the input data sequence while computing. Therefore additional computation time and memory space are generally needed. Also, the desirable property of in-place computation, which is easily achieved by split-radix algorithms, is difficult to implement.

In the literature, the fast algorithms of both categories are described for DHT-II [24,25,26,27,28,29]. Thus, in [23], the prime factor algorithm initially decomposed N-point DHT-II into q p-point discrete cosine transforms and p q-point DHTs. Then, to reduce the required number of arithmetic operations and eliminate the irregular index mapping, the resulting decomposition was combined with the split-radix algorithm. A regular computational structure was achieved, providing flexibility for composite sequence lengths and substantially saving the computational complexity.

In [24], an algorithm of the DHT-II computation was presented when two sequences of N/2-length DHT-II adjacent coefficients were given. In this way, two new radix-2 algorithms were developed for computing the N-point DHT-II. Next, obtaining the N-point DHT-II was described based on the two adjacent N/2-point DHT-IIs. Compared with the traditional Hu’s algorithm, the arithmetic operations could be saved from 16% to 24% for N varying from 16 to 64 [24]. If N = 8, the number of arithmetic operations slightly increased.

In [25], the algorithm from [24] was further developed. The obtained radix-2 recursive algorithm required a reduced number of arithmetic operations compared with existing algorithms and could be easily implemented.

The algorithms proposed in [24,25] were based on the decimation-in-frequency approach. But, if DHT-II fast algorithms need to be developed for real-time applications, then the DIT approach is preferable. That is why in [26], the radix-2 DIT fast algorithm for DHT-II calculation was presented. It was described by a regular signal flow graph, which provided flexibility to different transform lengths, substantially reducing the arithmetic complexity compared with Hu’s algorithms.

In [27], a split-radix DIT algorithm for the fast calculation of the DHT-II was proposed. The splitting of DHT-II of length N into one DHT-II of length N/2 was implemented for even-indexed samples. The splitting of DHT-II of length N into two DHTs-II of length N/4 was applied for odd-indexed samples. The proposed algorithm considerably reduced the arithmetic complexity compared to the existing DHT-II algorithms. In particular, the number of multiplications could be decreased by 29% compared to the radix-2 DHT-II algorithm [26], and the number of additions was reduced by 10%.

In [28], a modular VLSI algorithm for the fast calculation of DHT-II was introduced. The author efficiently incorporated the obfuscation technique with low overheads, which presented the additional advantage of hardware security. The proposed VLSI algorithm applied modular and regular computational structures named quasi-band correlations. In addition, the use of the systolic arrays architectural paradigm allowed for the efficient implementation of this algorithm on a VLSI chip.

The aforementioned split-radix algorithms were focused on input sequences, with length N being the power of two. Fast algorithms for calculating DHT-II with length N being the power of three were also developed. For example, in [29], a fast radix-3 algorithm for the computation of the DHT-II was developed for input sequences of length N, where N = 3^m, m ≥ 2. Then, this algorithm was applied to the direct computation of DHT-II when three adjacent sequences of length N/3 DHT-II were given. The arithmetic operations could be saved by up to 29%. Furthermore, the new algorithm could be easily implemented in VLSI.

1.2. The Main Contributions of This Paper

As a result of the analysis of articles devoted to developing the fast DHT-II algorithms, the following should be noted [5]. Firstly, the mathematical notations used by the authors to describe the proposed fast algorithms are often quite complicated. As a result, the understanding of the essence of the suggested solutions may be difficult for practitioners unfamiliar with the polynomial theory, number theory, Galois field theory, etc. Secondly, the mentioned theories and the developed algorithms concern specialized tasks. However, reducing the computing operations is desirable when a wide range of problems are solved. To lessen these disadvantages, in [5], a structural approach was proposed to develop fast algorithms for matrix–vector multiplication. It was based on a deep analysis of the structure of base transform matrices, identifying individual features of the arrangement of identical entries. Then, if necessary, the matrix structure was changed for subsequent use of the submatrix factorization [5]. The structural approach did not limit the length of the original data sequence, for example, to a power of two or three [5]. The number of multiplication and addition operations necessary for matrix–vector product computation could be reduced using this incredibly simple and effective approach. Since DHT-II coefficient matrices are structured, this research aimed to develop reduced-complexity DHT-II algorithms for input sequences of lengths N = 2, 3, 4, 5, 6, 7, and 8 based on a structural approach.

This paper is structured as follows. Section 1 introduces the problem and defines the subject of this study. Section 2 contains a short mathematical background and defines the system of used notations and symbols. Section 3 contains the derivation of algorithms for DHT-II transform sizes from 2 to 8. Section 4 and Section 5 contain a discussion of this study’s results and a comparative evaluation of the computational costs required to implement the developed and known algorithms. Conclusions are presented in Section 6.

2. Preliminary Remarks

The discrete Hartley transform (DHT) is one of the orthogonal transforms used, among others, for the analysis and processing of sounds and signals. The DHT advantages are primarily because it can successfully replace the discrete Fourier transform in many applications. If we apply the DHT to real-valued signals, then real-valued results are obtained in contrast to the Fourier transform. The use of the latter to real-valued signals in general gives complex-valued results. The real-data processing significantly simplifies the implementation of calculations [1,2,3].

For example, the direct generalized DHT can be represented by the following expression [23]:

$$X_k={{\displaystyle \sum }}_{n=0}^{N-1}{x}_n{ \mathit{cas}(n+n}_0{\left)\right(k+k}_0){\displaystyle \frac{2\mathrm{\pi }}{N}}, k=0, 1, \dots , N-1,$$

(1)

where

$$\mathit{cas}\mathsf{\theta }=\mathit{cos}\mathsf{\theta }+\mathit{sin}\mathsf{\theta }=\sqrt{2}\mathit{cos}(\mathsf{\theta }-\pi /4);$$

(2)

$X_k$ is the output sequence after the direct generalized DHT is performed, $x_n$ is the sequence of input data, and N is the number of signal samples. By expression (1), we have the DHT-I if $n_0=0, k_0=0;$ the DHT-II if $n_0=1/2, k_0=0;$ the DHT-III if $n_0=0, k_0=1/2;$ and the DHT-IV if $n_0=1/2, k_0=1/2.$

Therefore the DHT-II can be obtained as follows [23]:

$$y_k={{\displaystyle \sum }}_{n=0}^{N-1}{x}_n \mathit{cas}{\displaystyle \frac{\pi \left(2n+1\right)k}{N}}, k=0, 1, \dots , N-1,$$

(3)

where $y_k$ is the output sequence after the DHT-II is performed.

In matrix notation, DHT-II can be represented as follows:

$${\mathbf{Y}}_{N\times 1}{=\mathbf{C}}_N{\mathbf{X}}_{N\times 1},$$

(4)

where ${\mathbf{Y}}_{N\times 1}={\left[y_0{, y}_1, \dots , y_{N-1}\right]}^{\mathrm{T}}$, ${\mathbf{X}}_{N\times 1}={\left[x_0{, x}_1,\dots , x_{N-1}\right]}^{\mathrm{T}}$, $c_{\mathit{kl}}=\mathit{cas}{\displaystyle \frac{\mathrm{\pi }\left(2n+1\right)k}{N}}$, k, l = 0, …, N − 1.

In this paper, we use the following markings and signs [5]:

• ${\mathbf{I}}_N$ is an order N identity matrix;
• ${\mathbf{H}}_2$ is a 2 × 2 Hadamard matrix;
• ${\mathbf{1}}_{N\times M}$ is a N × M matrix of ones (a matrix where every element is equal to one);
• ⊗ is the Kronecker product of two matrices;
• ⊕ is the direct sum of two matrices.

DHT-II in matrix notation is as follows:

$$\left[\begin{array}{c}\begin{array}{c}{y}_0\\ y_1\end{array}\\ \begin{array}{c}⋮\\ y_{N-1}\end{array}\end{array}\right]=\left[\begin{array}{ccc}\begin{array}{cc}{c}_{0,0}& c_{0,1}\\ c_{1,0}& c_{1,1}\end{array}& \cdots & \begin{array}{c}{c}_{0,N-1}\\ c_{1,N-1}\end{array}\\ ⋮& \ddots & ⋮\\ \begin{array}{cc}{c}_{N-1, 0}& c_{N-1, 1} \end{array}& \cdots & c_{N-1,N-1}\end{array}\right]\left[\begin{array}{c}\begin{array}{c}{x}_0\\ x_1\end{array}\\ \begin{array}{c}⋮\\ x_{N-1}\end{array}\end{array}\right]$$

(5)

An empty cell in a matrix means that it contains zero. The multipliers are marked as $s_m^{\left(N\right)}$.

3. Fast Algorithms for Small-Size DHT-II

3.1. Algorithm for 2-Point DHT-II

Let us elaborate the algorithm for two-point DHT-II. The expression for two-point DHT-II is as follows:

$${\mathbf{Y}}_{2\times 1}{=\mathbf{C}}_2{\mathbf{X}}_{2\times 1},$$

(6)

where ${\mathbf{Y}}_{2\times 1}= {\left[y_0{, y}_1\right]}^{\mathrm{T}}$, ${\mathbf{X}}_{2\times 1}= {\left[x_0{, x}_1\right]}^{\mathrm{T}}$, ${\mathbf{C}}_2=\left[\begin{array}{cc}{a}_2& a_2\\ a_2& {-a}_2\end{array}\right]$, $a_2$ = 1.

Then, the expression for DHT-II for N = 2 can be presented as follows:

$${\mathbf{Y}}_{2\times 1}{=\mathbf{H}}_2{\mathbf{X}}_{2\times 1}.$$

(7)

Figure 1 shows a data flow graph of the synthesized algorithm for the two-point DHT-II. As can be seen, the number of addition operations is 2, which is the same as when using the direct method.

3.2. Algorithm for 3-Point DHT-II

Now, the algorithm for the three-point DHT-II can be developed. The expression for the three-point DHT-II is as follows:

$${\mathbf{Y}}_{3\times 1}{=\mathbf{C}}_3{\mathbf{X}}_{3\times 1},$$

(8)

where ${\mathbf{Y}}_{3\times 1}={\left[y_0{, y}_1{, y}_2\right]}^{\mathrm{T}}$, ${\mathbf{X}}_{3\times 1}= {\left[x_0{, x}_1, x_2\right]}^{\mathrm{T}}$, ${\mathbf{C}}_3=\left[\begin{array}{ccc}1& c_3& b_3\\ 1& -1& 1\\ 1& {-b}_3& {-c}_3\end{array}\right]$ with $b_3$ = 0.366, $c_3$ = 1.366.

We denote $c_3$ = 1 + $b_3$ and decompose the matrix ${\mathbf{C}}_3$ into two components:

$${\mathbf{C}}_3{=\mathbf{C}}_3^{\left(a\right)}{+\mathbf{C}}_3^{\left(b\right)},$$

(9)

where ${\mathbf{C}}_3^{\left(a\right)}=\left[\begin{array}{ccc}1& b_3& b_3\\ 1& -1& 1\\ 1& {-b}_3& {-b}_3\end{array}\right]$, ${\mathbf{C}}_3^{\left(b\right)}=\left[\begin{array}{ccc}& 1& \\ & & \\ & & -1\end{array}\right]$.

Taking into account properties of structural matrices [5,32], the computational procedure for the three-point DHT-II is represented by the expression

$${\mathbf{Y}}_{3\times 1}{=\mathbf{W}}_{3\times 5}{\mathbf{D}}_5{\mathbf{W}}_{5\times 3}{\mathbf{X}}_{3\times 1},$$

(10)

where ${\mathbf{W}}_{5\times 3}=\left[\begin{array}{ccc}& 1& \\ 1& & \\ 1& -1& 1\\ & 1& 1\\ & & 1\end{array}\right]$, ${\mathbf{D}}_5=\mathit{diag}\left({1, 1, 1, s}_0^{\left(3\right)}, 1\right)$, $s_0^{\left(3\right)}{=b}_3,$ ${\mathbf{W}}_{3\times 5}=\left[\begin{array}{ccccc}1& 1& & 1& \\ & & 1& & \\ & 1& & -1& -1\end{array}\right]$.

A data flow graph of the proposed algorithm for the three-point DHT-II is shown in Figure 2. As can be seen, we are able to reduce the number of multiplication operations from 4 to 1. However, the number of addition operations is increased from 6 to 7.

3.3. Algorithm for 4-Point DHT-II

Let us design the algorithm for the four-point DHT-II. The four-point DHT-II is expressed as follows:

$${\mathbf{Y}}_{4\times 1}{=\mathbf{C}}_4{\mathbf{X}}_{4\times 1},$$

(11)

where ${\mathbf{Y}}_{4\times 1}={\left[y_0{, y}_1{, y}_2{, y}_3\right]}^{\mathrm{T}}$, ${\mathbf{X}}_{4\times 1}={\left[x_0{, x}_1, x_2{, x}_3\right]}^{\mathrm{T}}$, ${\mathbf{C}}_4=\left[\begin{array}{cccc}1& a_4& 1& 0\\ 1& 0& -1& a_4\\ 1& {-a}_4& 1& 0\\ 1& 0& -1& a_4\end{array}\right]$ with $a_4$= 1.4142.

To change the order of the columns of ${\mathbf{C}}_4$, we define the permutation ${\mathrm{\pi }}_4$ in the following form: ${\mathrm{\pi }}_4=\left(\begin{array}{cc}1& 2\\ 1& 3\end{array}\begin{array}{cc}3& 4\\ 2& 4\end{array} \right)$.

After permutation of the columns of ${\mathrm{C}}_4$ according to ${\mathrm{\pi }}_4$, we obtain the matrix ${\mathbf{C}}_4^{\left(\mathbf{A}\right)}=\left[\begin{array}{cccc}1& 1& a_4& 0\\ 1& -1& 0& a_4\\ 1& 1& {-a}_4& 0\\ 1& -1& 0& a_4\end{array}\right]$ with the permutation matrix ${\mathbf{P}}_4=\left[\begin{array}{cccc}1& & & \\ & & 1& \\ & 1& & \\ & & & 1\end{array}\right]$.

Taking into account the properties of structural matrices [5,30], the computational procedure for the four-point DHT-II is represented by the expression

$${\mathbf{Y}}_{4\times 1}{=\mathbf{W}}_4^{\left(1\right)}{\mathbf{D}}_4{\mathbf{W}}_4^{\left(0\right)}{\mathbf{P}}_4{\mathbf{X}}_{4\times 1},$$

(12)

where ${\mathbf{W}}_4^{\left(0\right)}{=\mathbf{H}}_2\oplus {\mathbf{I}}_2$, ${\mathbf{D}}_4=\mathit{diag}\left({1, 1, s}_0^{\left(4\right)}{, s}_0^{\left(4\right)}\right)$, $s_0^{\left(4\right)}{=a}_4,$ ${\mathbf{W}}_4^{\left(1\right)}{=\mathbf{H}}_2\otimes {\mathbf{I}}_2$.

A data flow graph of the proposed algorithm of the four-point DHT-II is presented in Figure 3. In particular, the number of multiplication operations may be reduced from 4 to 2 and the number of addition operations is decreased from 8 to 6.

3.4. Algorithm for 5-Point DHT-II

Let us obtain the algorithm for the five-point DHT-II. The five-point DHT-II is expressed as follows:

$${\mathbf{Y}}_{5\times 1}{=\mathbf{C}}_5{\mathbf{X}}_{5\times 1},$$

(13)

where ${\mathbf{Y}}_{5\times 1}={\left[y_0{, y}_1{, y}_2{, y}_3{, y}_4\right]}^{\mathrm{T}}$, ${\mathbf{X}}_{5\times 1}= {\left[x_0{, x}_1, x_2{, x}_3{, x}_4\right]}^{\mathrm{T}}$, $b_5$ = 1.3968, $c_5$= 0.6420, $d_5$ = 1.2601, $e_5$= 0.2212,

$${\mathbf{C}}_5=\left[\begin{array}{ccccc}1& b_5& d_5& c_5& {-e}_5\\ 1& c_5& {-b}_5& e_5& d_5\\ 1& -1& 1& -1& 1\\ 1& {-d}_5& {-e}_5& b_5& {-c}_5\\ 1& e_5& {-c}_5& {-d}_5& {-b}_5\end{array}\right].$$

To change the order of rows and columns of ${\mathbf{C}}_5$, the permutations ${\mathrm{\pi }}_5$ and ${\mathrm{\pi }}_6$ are defined in the following form:

$${\mathrm{\pi }}_5=\left(\begin{array}{cc}1& 2\\ 3& 1\end{array} \begin{array}{ccc}3& 4& 5\\ 2& 5& 4\end{array} \right), {\mathrm{\pi }}_6=\left(\begin{array}{cc}1& 2\\ 1& 2\end{array} \begin{array}{ccc}3& 4& 5\\ 3& 5& 4\end{array} \right).$$

The rows of ${\mathbf{C}}_5$ are permutated according to ${\mathrm{\pi }}_5$ and the columns of ${\mathbf{C}}_5$ are permutated according to ${\mathrm{\pi }}_6$. After permutations, the matrix acquires the following structure:

$${\mathbf{C}}_5^{\left(a\right)}=\left[\begin{array}{ccccc}1& -1& 1& -1& 1\\ 1& b_5& d_5& {-e}_5& c_5\\ 1& c_5& {-b}_5& d_5& e_5\\ 1& e_5& {-c}_5& {-b}_5& {-d}_5\\ 1& {-d}_5& {-e}_5& {-c}_5& b_5\end{array}\right].$$

Then, the matrix ${\mathbf{C}}_5^{\left(a\right)}$ is decomposed into two components:

$${\mathbf{C}}_5^{\left(a\right)}{=\mathbf{C}}_5^{\left(b\right)}{+\mathbf{C}}_5^{\left(c\right)},$$

(14)

where ${\mathbf{C}}_5^{\left(b\right)}=\left[\begin{array}{ccccc}1& -1& 1& -1& 1\\ 1& & & & \\ 1& & & & \\ 1& & & & \\ 1& & & & \end{array}\right]$ and ${\mathbf{C}}_5^{\left(c\right)}=\left[\begin{array}{ccccc} & & & & \\ & b_5& d_5& {-e}_5& c_5\\ & c_5& {-b}_5& d_5& e_5\\ & e_5& {-c}_5& {-b}_5& {-d}_5\\ & {-d}_5& {-e}_5& {-c}_5& b_5\end{array}\right]$.

Matrix ${\mathbf{C}}_5^{\left(b\right)}$ has the same entries except on the sign in the first column and first row, which allows for reducing the number of operations without the need for further transformations [32]. After eliminating rows and columns containing only zero entries in matrix ${\mathbf{C}}_5^{\left(c\right)}$, we obtain matrix ${\mathbf{C}}_4^{\left(d\right)}$:

$${\mathbf{C}}_4^{\left(d\right)}=\left[\begin{array}{cccc}{b}_5& d_5& {-e}_5& c_5\\ c_5& {-b}_5& d_5& e_5\\ e_5& {-c}_5& {-b}_5& {-d}_5\\ {-d}_5& {-e}_5& {-c}_5& b_5\end{array}\right]$$

which matches the matrix pattern $\left[\begin{array}{cc}{\mathbf{A}}_2& {\mathbf{B}}_2\\ {-\mathbf{B}}_2& {-\mathbf{A}}_2\end{array}\right],$ where ${\mathbf{A}}_2=\left[\begin{array}{cc}{b}_5& d_5\\ c_5& {-b}_5\end{array}\right]$, ${\mathbf{B}}_2=\left[\begin{array}{cc}{-e}_5& c_5\\ d_5& e_5\end{array}\right]$.

Hence, the matrix ${\mathbf{C}}_4^{\left(d\right)}$ can be represented as [5,30]

$${\mathbf{C}}_4^{\left(d\right)}=({\overline{\mathbf{I}}}_2\otimes {\mathbf{I}}_2)({\mathbf{H}}_2\otimes {\mathbf{I}}_2{)1/2[(\mathbf{A}}_2{+\mathbf{B}}_2{)\oplus (\mathbf{A}}_2{-\mathbf{B}}_2{\left)\right](\mathbf{H}}_2\otimes {\mathbf{I}}_2),$$

(15)

where ${\overline{\mathrm{I}}}_2=\left[\begin{array}{cc}1& 0\\ 0& -1\end{array}\right]$. Considering the structures of the resulting matrices ${\mathbf{A}}_2{+\mathbf{B}}_2$ and (${\mathbf{A}}_2-{\mathbf{B}}_2$), we note that these matrices match the patterns $\left[\begin{array}{cc}a& b\\ b& -a\end{array}\right]$ and $\left[\begin{array}{cc}c& d\\ -d& -c\end{array}\right],$ respectively. Here, ${a=b}_5{ - e}_5{; b=c}_5{+d}_5{; c=b}_5{+e}_5; d=-{ c}_5{+d}_5$. Then [5,31],

$$\begin{array}{c}\hfill {\mathbf{A}}_2{+\mathbf{B}}_2={\mathbf{T}}_{2\times 3}^{\left(4\right)} \times \mathit{diag}(a-b, -a-b, b){\mathbf{T}}_{3\times 2}^{\left(3\right)}\\ \hfill {\mathbf{A}}_2{-\mathbf{B}}_2= {\overline{\mathbf{I}}}_2{\mathbf{H}}_2\times \mathit{diag}\left(\right(c+d)/2, (c-d)/2) {\mathbf{H}}_2\end{array}$$

(16)

where ${\mathbf{T}}_{2\times 3}^{\left(4\right)}=\left[\begin{array}{ccc}1& 0& 1\\ 0& 1& 1\end{array}\right]$, ${\mathbf{T}}_{3\times 2}^{\left(3\right)}=\left[\begin{array}{cc}1& 0\\ 0& 1\\ 1& 1\end{array}\right]$.

We define $s_0^{\left(5\right)}=(a-{b)/2; s}_1^{\left(5\right)}=-{ (a+b)/2; s}_2^{\left(5\right)}{=b/2; s}_3^{\left(5\right)}{=(c+d)/4; s}_4^{\left(5\right)}=(c-d)/4$. Then

$${\mathbf{C}}_4^{\left(d\right)}{\overline{\mathbf{I}}}_2\otimes {\mathbf{I}}_2{\left)\right(\mathbf{H}}_2\otimes {\mathbf{I}}_2{\left)\right(\mathbf{T}}_{2\times 3}^{\left(4\right)}\oplus ({\overline{\mathbf{I}}}_2{\mathbf{H}}_2\left)\right)\times \mathit{diag}(s_0^{\left(5\right)}{, s}_1^{\left(5\right)}{, s}_2^{\left(5\right)}{, s}_3^{\left(5\right)}{, s}_4^{\left(5\right)}) {(\mathbf{T}}_{3\times 2}^{\left(3\right)}\oplus {\mathbf{H}}_2{\left)\right(\mathbf{H}}_2\otimes {\mathbf{I}}_2).$$

(17)

Taking into account the properties of structural matrices [5,32], the computational procedure for the five-point DHT-II is represented by the expression

$${\mathbf{Y}}_{5\times 1}{=\mathbf{P}}_5^{\left(0\right)}{\mathbf{W}}_{5\times 7}{\mathbf{W}}_7^{\left(0\right)}{\mathbf{W}}_{7\times 8}{\mathbf{D}}_8{\mathbf{W}}_{8\times 7}{\mathbf{W}}_7^{\left(1\right)}{\mathbf{W}}_{7\times 10}{\mathbf{W}}_{10\times 5}{\mathbf{P}}_5^{\left(1\right)}{\mathbf{X}}_{5\times 1}$$

(18)

where

$${\mathbf{P}}_5^{\left(0\right)}=\left[\begin{array}{ccccc}& 1& & & \\ & & 1& & \\ 1& & & & \\ & & & & 1\\ & & & 1& \end{array}\right], {\mathbf{P}}_5^{\left(1\right)}=\left[\begin{array}{ccccc}1& & & & \\ & 1& & & \\ & & 1& & \\ & & & & 1\\ & & & 1& \end{array}\right]$$

$$\begin{array}{c}{ \mathbf{W}}_{8\times 7}{=1\oplus \mathbf{T}}_{3\times 2}^{\left(3\right)}\oplus {\mathbf{H}}_2\oplus {\mathbf{I}}_2, {\mathbf{W}}_7^{\left(1\right)}=1\oplus \left({\mathbf{H}}_2\oplus {\mathbf{I}}_2\right)\oplus {\mathbf{I}}_2, {\mathbf{W}}_{7\times 10}{=\mathbf{I}}_5\oplus {\mathbf{P}}_{2\times 5}^{\left(2\right)}, {\mathbf{W}}_{10\times 5}{=\mathbf{1}}_{2\times 1}\otimes {\mathbf{I}}_5, \\ {\mathbf{W}}_{7\times 8}{=1\oplus \mathbf{T}}_{2\times 3}^{\left(4\right)}\oplus ({\overline{\mathbf{I}}}_2{ \mathbf{H}}_2{)\oplus \mathbf{I}}_2, {\mathbf{D}}_8{=\mathit{diag}(1, s}_0^{\left(5\right)}{, s}_1^{\left(5\right)}{, s}_2^{\left(5\right)}{, s}_3^{\left(5\right)}{, s}_4^{\left(5\right)}, 1, 1),\\ {\mathbf{P}}_{2\times 5}^{\left(2\right)}=\left[\begin{array}{ccccc}& -1& 1& 1& -1\\ 1& & & & \end{array}\right],\\ {\mathbf{W}}_7^{\left(0\right)}=\left[\begin{array}{ccccccc}1& & & & & & \\ & 1& & 1& & & \\ & & 1& & 1& & \\ & -1& & 1& & & \\ & & -1& & 1& & \\ & & & & & 1& \\ & & & & & & 1\end{array}\right], {\mathbf{W}}_{5\times 7}=\left[\begin{array}{ccccccc}1& & & & & 1& \\ & 1& & & & & 1\\ & & 1& & & & 1\\ & & & 1& & & 1\\ & & & & 1& & 1\end{array}\right]\end{array}$$

A data flow graph of the proposed algorithm of the five-point DHT-II is presented in Figure 4. In particular, the number of multiplication operations may be reduced from 16 to 5, but the number of addition operations is increased from 20 to 23.

3.5. Algorithm for 6-Point DHT-II

We now propose an algorithm for the six-point DHT-II. The six-point DHT-II is expressed as follows:

$${\mathbf{Y}}_{6\times 1}{=\mathbf{C}}_6{\mathbf{X}}_{6\times 1},$$

(19)

where ${\mathbf{Y}}_{6\times 1}={\left[y_0{, y}_1{, y}_2{, y}_3{, y}_4{, y}_5\right]}^{\mathrm{T}}$, ${\mathbf{X}}_{6\times 1}={\left[x_0{, x}_1, x_2{, x}_3{, x}_4{, x}_5\right]}^{\mathrm{T}}$, ${\mathbf{C}}_6=\left[\begin{array}{cccccc}1& c_6& c_6& 1& b_6& {-b}_6\\ 1& 1& -1& -1& 1& 1\\ 1& {-b}_6& {-b}_6& 1& {-c}_6& c_6\\ 1& {-c}_6& c_6& -1& b_6& b_6\\ 1& -1& -1& 1& 1& -1\\ 1& b_6& {-b}_6& -1& {-c}_6& {-c}_6\end{array}\right]$ with $b_6$ = 0.366, $c_6$ = 1.366.

As can be seen, it is $c_6$ = 1 +$b_6$; then, we decompose

$${\mathbf{C}}_6{=\mathbf{C}}_6^{\left(a\right)}{+\mathbf{C}}_6^{\left(b\right)},$$

(20)

where

$${\mathbf{C}}_6^{\left(a\right)}=\left[\begin{array}{cccccc}1& 1& 1& 1& 0& 0\\ 1& 1& -1& -1& 1& 1\\ 1& 0& 0& 1& -1& 1\\ 1& -1& 1& -1& 0& 0\\ 1& -1& -1& 1& 1& -1\\ 1& 0& 0& -1& -1& -1\end{array}\right]\mathbf{,} {\mathbf{C}}_6^{\left(b\right)}=\left[\begin{array}{cccccc}& b_6& b_6& & b_6& {-b}_6\\ & & & & & \\ & {-b}_6& {-b}_6& & {-b}_6& b_6\\ & {-b}_6& b_6& & b_6& b_6\\ & & & & & \\ & b_6& {-b}_6& & {-b}_6& {-b}_6\end{array}\right].$$

Based on the properties of structural matrices [5,32] and denoted $s_0^{\left(6\right)}=b_6$, the computational procedure for the six-point DHT-II is represented by the expression

$${\mathbf{Y}}_{6\times 1}{=\mathbf{W}}_{6\times 8}{\mathbf{D}}_8{\mathbf{W}}_{8\times 10}{\mathbf{W}}_{10\times 6}{\mathbf{X}}_{6\times 1},$$

(21)

where

$$\begin{array}{c}{\mathbf{W}}_{10\times 6}=\left[\begin{array}{cccccc}1& & & & & \\ & 1& & & & \\ & & 1& & & \\ & & & 1& & \\ & & & & 1& \\ & & & & & 1\\ & 1& & & & \\ & & 1& & & \\ & & & & 1& \\ & & & & & 1\end{array}\right], {\mathbf{W}}_{6\times 8}=\left[\begin{array}{cccccccc}1& & & & & & 1& \\ & 1& & & & & & \\ & & 1& & & & -1& \\ & & & 1& & & & -1\\ & & & & 1& & & \\ & & & & & 1& & 1\end{array}\right],\\ { \mathbf{W}}_{8\times 10}{=C}_6^{\left(a\right)}\oplus {\mathbf{B}}_{2\times 4}, {\mathbf{B}}_{2\times 4}=\left[\begin{array}{cccc}1& 1& 1& -1\\ 1& -1& -1& -1\end{array}\right], {\mathbf{D}}_8{=\mathit{diag}(1, 1, 1, 1, 1, 1, s}_0^{\left(6\right)}{, s}_0^{\left(6\right)}).\end{array}$$

A data flow graph of the proposed algorithm of the six-point DHT-II is presented in Figure 5. We are able to reduce the number of multiplications from 16 to 2 and the number of addition operations from 30 to 26.

3.6. Algorithm for 7-Point DHT-II

To elaborate the algorithm for the seven-point DHT-II, the matrix notation of this transform is expressed as follows:

$${\mathbf{Y}}_{7\times 1}{=\mathbf{C}}_7{\mathbf{X}}_{7\times 1},$$

(22)

where ${\mathbf{Y}}_{7\times 1}={\left[y_0{, y}_1{, y}_2{, y}_3{, y}_4{, y}_5{, y}_6\right]}^{\mathrm{T}}$, ${\mathbf{X}}_{7\times 1}={\left[x_0{, x}_1, x_2{, x}_3{, x}_4{, x}_5{, x}_6\right]}^{\mathrm{T}}$, $b_7$=1.3349, $c_7$=1.1974, $d_7$=0.1583, $e_7$=1.4053, $g_7$=0.4671, $f_7$=0.7524,

$${\mathbf{C}}_7=\left[\begin{array}{ccccccc}1& b_7& e_7& c_7& f_7& d_7& -g_7\\ 1& c_7& -g_7& -e_7& -d_7& b_7& f_7\\ 1& d_7& -c_7& b_7& -g_7& -f_7& e_7\\ 1& -1& 1& -1& 1& -1& 1\\ 1& -e_7& f_7& g_7& -b_7& c_7& -d_7\\ 1& -f_7& b_7& d_7& e_7& g_7& -c_7\\ 1& g_7& -d_7& -f_7& -c_7& -e_7& -b_7\end{array}\right].$$

Now, we need to change the order of rows and columns. Let us define the permutations ${\mathrm{\pi }}_7,{ \mathrm{\pi }}_8$ in the following form:

$${\mathrm{\pi }}_7=\left(\begin{array}{cc}1& 2\\ 4& 1\end{array} \begin{array}{ccccc}3& 4& 5& 6& 7\\ 2& 3& 7& 6& 5\end{array} \right), {\mathrm{\pi }}_8=\left(\begin{array}{ccccccc}1& 2& 3& 4& 5& 6& 7\\ 1& 2& 3& 4& 7& 6& 5\end{array}\right).$$

Then, we permute rows and columns of ${\mathbf{C}}_7$ according to ${\mathrm{\pi }}_7,{ \mathrm{\pi }}_8$. After permutations, the matrix acquires the following structure:

$${\mathbf{C}}_7^{\left(a\right)}=\left[\begin{array}{ccccccc}1& -1& 1& -1& 1& -1& 1\\ 1& b_7& e_7& c_7& -g_7& d_7& f_7\\ 1& c_7& -g_7& -e_7& f_7& b_7& -d_7\\ 1& d_7& -c_7& b_7& e_7& -f_7& -g_7\\ 1& g_7& -d_7& -f_7& -b_7& -e_7& -c_7\\ 1& -f_7& b_7& d_7& -c_7& g_7& e_7\\ 1& -e_7& f_7& g_7& -d_7& c_7& -b_7\end{array}\right].$$

Further, the matrix ${\mathbf{C}}_7^{\left(a\right)}$ is decomposed into two components:

$${\mathbf{C}}_7^{\left(a\right)}{=\mathbf{C}}_7^{\left(b\right)}{+\mathbf{C}}_7^{\left(c\right)},$$

(23)

where

$$\begin{array}{c}{\mathbf{C}}_7^{\left(b\right)}=\left[\begin{array}{ccccccc}1& -1& 1& -1& 1& -1& 1\\ 1& & & & & & \\ 1& & & & & & \\ 1& & & & & & \\ 1& & & & & & \\ 1& & & & & & \\ 1& & & & & & \end{array}\right], \\ {\mathbf{C}}_7^{\left(c\right)}=\left[\begin{array}{ccccccc}& & & & & & \\ & b_7& e_7& c_7& {-g}_7& d_7& f_7\\ & c_7& {-g}_7& {-e}_7& f_7& b_7& {-d}_7\\ & d_7& {-c}_7& b_7& e_7& {-f}_7& {-g}_7\\ & g_7& {-d}_7& {-f}_7& {-b}_7& {-e}_7& {-c}_7\\ & {-f}_7& b_7& d_7& {-c}_7& g_7& e_7\\ & {-e}_7& f_7& g_7& {-d}_7& c_7& {-b}_7\end{array}\right].\end{array}$$

Matrix ${\mathbf{C}}_7^{\left(B\right)}$ has the same entries except on the sign in the first column and first row, which allows us to reduce the number of operations without the need for further transformations [31]. After eliminating the rows and columns containing only zero entries in matrix ${\mathbf{C}}_7^{\left(c\right)}$, we obtain matrix ${\mathbf{C}}_6^{\left(d\right)}$:

$${\mathbf{C}}_6^{\left(d\right)}=\left[\begin{array}{cccccc}{b}_7& e_7& c_7& {-g}_7& d_7& f_7\\ c_7& {-g}_7& {-e}_7& f_7& b_7& {-d}_7\\ d_7& {-c}_7& b_7& e_7& {-f}_7& {-g}_7\\ g_7& {-d}_7& {-f}_7& {-b}_7& {-e}_7& {-c}_7\\ {-f}_7& b_7& d_7& {-c}_7& g_7& e_7\\ {-e}_7& f_7& g_7& {-d}_7& c_7& {-b}_7\end{array}\right].$$

The obtained matrix acquires the structure $\left[\begin{array}{cc}{\mathbf{A}}_3& {\mathbf{B}}_3\\ {-\mathbf{B}}_3& {-\mathbf{A}}_3\end{array}\right],$ with ${\mathbf{A}}_3=\left[\begin{array}{ccc}{b}_7& e_7& c_7\\ c_7& {-g}_7& {-e}_7\\ d_7& {-c}_7& b_7\end{array}\right]$ and ${\mathbf{B}}_3=\left[\begin{array}{ccc}{-g}_7& d_7& f_7\\ f_7& b_7& {-d}_7\\ e_7& {-f}_7& {-g}_7\end{array}\right]$. Then [5,32],

$${\mathbf{C}}_6^{\left(d\right)}=({\overline{\mathrm{I}}}_2\otimes {\mathbf{I}}_3{\left)\right(\mathbf{H}}_2\otimes {\mathbf{I}}_3{)1/2[(\mathbf{A}}_3{+\mathbf{B}}_3{)\oplus (\mathbf{A}}_3{-\mathbf{B}}_3{\left)\right](\mathbf{H}}_2\otimes {\mathbf{I}}_3).$$

(24)

Submatrices ${(\mathbf{A}}_3{+\mathbf{B}}_3)$, ${(\mathbf{A}}_3{-\mathbf{B}}_3)$ of quasi-diagonal matrix ${(\mathbf{A}}_3{+\mathbf{B}}_3{)\oplus (\mathbf{A}}_3{-\mathbf{B}}_3)$ are represented as circular convolution matrices [32]:

$${\mathbf{A}}_3{+\mathbf{B}}_3=\left[\begin{array}{ccc}{ b}_7{-g}_7& { e}_7{+d}_7& c_7{+f}_7\\ { c}_7{+f}_7& b_7{-g}_7& {-(e}_7{+d}_7)\\ { e}_7{+d}_7& {-(c}_7{+f}_7)& b_7{-g}_7\end{array}\right], {\mathbf{A}}_3{-\mathbf{B}}_3=\left[\begin{array}{ccc}{ b}_7{+g}_7& { e}_7{-d}_7& c_7{-f}_7\\ { c}_7{-f}_7& {-(b}_7{+g}_7)& {-(e}_7{-d}_7)\\ {-(e}_7{-d}_7)& {-(c}_7{-f}_7)& b_7{+g}_7\end{array}\right].$$

Now, we will move onto dealing with matrices ${(\mathbf{A}}_3{+\mathbf{B}}_3)$ and ${(\mathbf{A}}_3{-\mathbf{B}}_3)$. In this case, the expressions for calculating the entries of a circular convolution matrix ${\mathbf{H}}_3=\left[\begin{array}{ccc}{ h}_0& { h}_2& { h}_1\\ { h}_1& { h}_0& { h}_2\\ { h}_2& { h}_1& { h}_0\end{array}\right]$ for N = 3 are as follows [32]:

$${\mathbf{H}}_3{=\mathbf{T}}_3^{\left(1\right)}{\mathbf{T}}_{3\times 4}{\times \mathit{diag}(s}_0{, s}_1{, s}_2{, s}_3{)\times \mathbf{T}}_{4\times 3}{\mathbf{T}}_3^{\left(0\right)},$$

(25)

where $s_0{=(h}_0{+h}_1{+h}_2{)/3, s}_1{=h}_0{-h}_2{, s}_2{=h}_1{-h}_2{, s}_3{=(h}_0{+h}_1{-2h}_2)/3$; ${\mathbf{T}}_3^{\left(1\right)}=\left[\begin{array}{ccc}1& 1& \\ 1& -1& -1\\ 1& & 1\end{array}\right]$, ${\mathbf{T}}_{3\times 4}=\left[\begin{array}{cccc}1& & & \\ & 1& & -1\\ & & 1& -1\end{array}\right]$, ${\mathbf{T}}_{4\times 3}=\left[\begin{array}{ccc}1& & \\ & 1& \\ & & 1\\ & 1& 1\end{array}\right]$, ${\mathbf{T}}_3^{\left(0\right)}=\left[\begin{array}{ccc}1& 1& 1\\ 1& & -1\\ & 1& -1\end{array}\right]$.

To make the ${(\mathbf{A}}_3{+\mathbf{B}}_3)$, ${(\mathbf{A}}_3{-\mathbf{B}}_3)$ matrices consistent with the circular convolution expression, we need to modify them. In the ${(\mathbf{A}}_3{+\mathbf{B}}_3)$ matrix, we change the sign of all terms in the first column and first row. In the ${(\mathbf{A}}_3{-\mathbf{B}}_3)$ matrix, we change the sign in the second column. In this way, we obtain the matrices:

$$\begin{array}{c}{{(\mathbf{A}}_3{+\mathbf{B}}_3)}^{\left(0\right)}=\left[\begin{array}{ccc}{ b}_7{-g}_7& { -(e}_7{+d}_7)& {-(c}_7{+f}_7)\\ {-(c}_7{+f}_7)& b_7{-g}_7& {-(e}_7{+d}_7)\\ {-( e}_7{+d}_7)& {-(c}_7{+f}_7)& b_7{-g}_7\end{array}\right],\\ {{(\mathbf{A}}_3{-\mathbf{B}}_3)}^{\left(0\right)}=\left[\begin{array}{ccc}{ b}_7{+g}_7& { -(e}_7{-d}_7)& c_7{-f}_7\\ { c}_7{-f}_7& b_7{+g}_7& {-(e}_7{-d}_7)\\ { -(e}_7{-d}_7)& c_7{-f}_7& b_7{+g}_7\end{array}\right].\end{array}$$

Using the three-point convolution algorithm, the values $s_i$, I = 0, 1, 2, 3 for matrices ${{(\mathbf{A}}_3{+\mathbf{B}}_3)}^{\left(0\right)}$ and ${{(\mathbf{A}}_3{-\mathbf{B}}_3)}^{\left(0\right)}$ take the following form:

$$\begin{array}{c}{s}_0^{\left(7\right)}{=(b}_7-{ g}_7-{ (e}_7{+d}_7) -{ (c}_7{+f}_7{\left)\right)/6; s}_1^{\left(7\right)}{=(b}_7-{ g}_7{+e}_7{+d}_7{)/2; s}_2^{\left(7\right)}{=(e}_7{+d}_7{ - c}_{7 }{- f}_7)/2;\\ s_3^{\left(7\right)}{=(b}_{7 }{- g}_7{- c}_{7 }{- f}_7{+2(e}_7{+d}_7\left)\right)/6; s_4^{\left(7\right)}{=(b}_{7 }{+g}_7{+c}_{7 }{- f}_7{-e}_7{+d}_7)/6; s_5^{\left(7\right)}{=(b}_{7 }{+g}_7{+e}_7-{ d}_7)/2;\\ s_6^{\left(7\right)}{=(c}_{7 }{- f}_7{+e}_7-{ d}_7)/2; s_7^{\left(7\right)}{=(b}_{7 }{+g}_7{+c}_{7 }{- f}_7{+2(e}_7-{ d}_7\left)\right)/6.\end{array}$$

Then,

$$\begin{array}{c}{\mathbf{C}}_6^{\left(d\right)}=({\overline{\mathbf{I}}}_2\otimes {\mathbf{I}}_3{\left)\right(\mathbf{H}}_2\otimes {\mathbf{I}}_3{\left)\right(\mathbf{T}}_3^{\left(1\right)}{\mathbf{T}}_{3\times 4}\oplus {\mathbf{T}}_3^{\left(1\right)}{\mathbf{T}}_{3\times 4}{)\times \mathit{diag}(s}_0^{\left(7\right)}{, s}_1^{\left(7\right)}{, s}_2^{\left(7\right)}{, s}_3^{\left(7\right)}{, s}_4^{\left(7\right)}{, s}_5^{\left(7\right)}{, s}_6^{\left(7\right)}{, s}_7^{\left(7\right)})\times \\ \times ({\mathbf{T}}_{4\times 3}{\mathbf{T}}_3^{\left(0\right)}\oplus {\mathbf{T}}_{4\times 3}{\mathbf{T}}_3^{\left(0\right)}{\left)\right(\mathbf{H}}_2\otimes {\mathbf{I}}_3).\end{array}$$

(26)

Based on the properties of structural matrices [5,32], the computational procedure for the seven-point DHT-II is represented by the expression

$${\mathbf{Y}}_{7\times 1}{=\mathbf{P}}_7^{\left(1\right)}{\mathbf{W}}_{7\times 9}{\mathbf{W}}_9{\mathbf{W}}_{9\times 11}{\mathbf{D}}_{11}{\mathbf{W}}_{11\times 9}{\mathbf{W}}_{9\times 14}{\mathbf{W}}_{14\times 7}{\mathbf{P}}_7^{\left(0\right)}{\mathbf{X}}_{7\times 1},$$

(27)

where

$$\begin{array}{c}{\mathbf{P}}_7^{\left(0\right)}=\left[\begin{array}{cccccc}1& & & & & \\ & 1& & & & \\ & & 1& & & \\ & & & 1& & \\ & & & & & 1\\ & & & & 1& \\ & & & 1& & \end{array}\right],\\ {\mathbf{P}}_7^{\left(1\right)}=\left[\begin{array}{ccccccc}& 1& & & & & \\ & & 1& & & & \\ & & & 1& & & \\ 1& & & & & & \\ & & & & & & 1\\ & & & & & 1& \\ & & & & 1& & \end{array}\right],\\ { \mathbf{W}}_{11\times 9}{=1\oplus (\mathbf{T}}_{4\times 3}{\mathbf{T}}_3^{\left(0\right)}{)\oplus (\mathbf{T}}_{4\times 3}{\mathbf{T}}_3^{\left(0\right)}{)\oplus \mathbf{I}}_2, { \mathbf{W}}_{9\times 11}{=1\oplus (\mathbf{T}}_3^{\left(1\right)}{\mathbf{T}}_{3\times 4}{)\oplus (\mathbf{T}}_3^{\left(1\right)}{\mathbf{T}}_{3\times 4}{)\oplus \mathbf{I}}_2,\\ {\mathbf{W}}_{9\times 14}=1\oplus \left({\mathbf{H}}_2\otimes {\mathbf{I}}_3\right)\oplus {\mathbf{P}}_{2\times 7}, {\mathbf{W}}_9=1\oplus ({\overline{\mathbf{I}}}_2\otimes {\mathbf{I}}_3{\left)\right(\mathbf{H}}_2\otimes {\mathbf{I}}_3)\oplus {\mathbf{I}}_2, {\mathbf{W}}_{14\times 7}{=1}_{2\times 1}\otimes {\mathbf{I}}_7,\\ {\mathbf{D}}_{11}{=\mathit{diag}(1, s}_0^{\left(7\right)}{, s}_1^{\left(7\right)}{, s}_2^{\left(7\right)}{, s}_3^{\left(7\right)}{, s}_4^{\left(7\right)}{, s}_5^{\left(7\right)}{, s}_6^{\left(7\right)}{, s}_7^{\left(7\right)}, 1, 1),\\ {\mathbf{P}}_{2\times 7}=\left[\begin{array}{ccccccc}& -1& 1& -1& 1& -1& 1\\ 1& & & & & & \end{array}\right],\\ {\mathbf{W}}_{7\times 9}=\left[\begin{array}{ccccccccc}1& & & & & & & 1& \\ & 1& & & & & & & 1\\ & & 1& & & & & & 1\\ & & & 1& & & & & 1\\ & & & & 1& & & & 1\\ & & & & & 1& & & 1\\ & & & & & & 1& & 1\end{array}\right].\end{array}$$

Figure 6 shows a data flow graph of the synthesized algorithm for the seven-point DHT-II. As can be seen, we are able to reduce the number of multiplications from 36 to 8, but the number of addition operations is increased from 42 to 46.

3.7. Algorithm for 8-Point DHT-II

Let us design the algorithm for the eight-point DHT-II. The eight-point DHT-II is expressed as follows:

$${\mathbf{Y}}_{8\times 1}{=\mathbf{C}}_8{\mathbf{X}}_{8\times 1},$$

(28)

where ${\mathbf{Y}}_{8\times 1}={\left[y_0{, y}_1{, y}_2{, y}_3{, y}_4{, y}_5{, y}_6{, y}_7\right]}^{\mathrm{T}}$, ${\mathbf{X}}_{8\times 1}={\left[x_0{, x}_1, x_2{, x}_3{, x}_4{, x}_5{, x}_6{, x}_7\right]}^{\mathrm{T}}$, ${\mathbf{C}}_8=\left[\begin{array}{cccccccc}1& b_8& d_8& b_8& 1& c_8& 0& -c_8\\ 1& b_8& 0& -b_8& -1& c_8& d_8& c_8\\ 1& c_8& -d_8& c_8& 1& -b_8& 0& b_8\\ 1& -c_8& 0& c_8& -1& b_8& -d_8& b_8\\ 1& -b_8& d_8& -b_8& 1& -c_8& 0& c_8\\ 1& -b_8& 0& b_8& -1& -c_8& d_8& -c_8\\ 1& -c_8& -d_8& -c_8& 1& b_8& 0& -b_8\\ 1& c_8& 0& -c_8& -1& -b_8& -d_8& -b_8\end{array}\right]$ with $b_8$ = 1.3066, $c_8$ = 0.5412, $d_8$ = 1.4142.

Let us define the permutation

$${\mathrm{\pi }}_8=\left(\begin{array}{cc}1& 2\\ 1& 5\end{array} \begin{array}{cccccc}3& 4& 5& 6& 7& 8\\ 3& 7& 2& 4& 6& 8\end{array} \right)$$

to change the order of columns of ${\mathbf{C}}_8$:

$${\mathbf{C}}_8^{\left(a\right)}=\left[\begin{array}{cccccccc}1& 1& d_8& 0& b_8& b_8& c_8& -c_8\\ 1& -1& 0& d_8& b_8& -b_8& c_8& c_8\\ 1& 1& -d_8& 0& c_8& c_8& -b_8& b_8\\ 1& -1& 0& -d_8& -c_8& c_8& b_8& b_8\\ 1& 1& d_8& 0& -b_8& b_8& -c_8& c_8\\ 1& -1& 0& d_8& -b_8& b_8& -c_8& -c_8\\ 1& 1& -d_8& 0& -c_8& -c_8& b_8& -b_8\\ 1& -1& 0& -d_8& c_8& -c_8& -b_8& -b_8\end{array}\right]$$

After permutations, the matrix ${\mathbf{C}}_8^{\left(a\right)}$ acquires the structure

$$\left[\begin{array}{cc}{\mathbf{A}}_4& {\mathbf{B}}_4\\ {\mathbf{A}}_4& {-\mathbf{B}}_4\end{array}\right] \mathrm{with} {\mathbf{A}}_4=\left[\begin{array}{cccc}1& 1& d_8& 0\\ 1& -1& 0& d_8\\ 1& 1& {-d}_8& 0\\ 1& -1& 0& {-d}_8\end{array}\right] \mathrm{and} {\mathbf{B}}_4=\left[\begin{array}{cccc}{b}_8& b_8& c_8& {-c}_8\\ b_8& {-b}_8& c_8& c_8\\ c_8& c_8& {-b}_8& b_8\\ {-c}_8& c_8& b_8& b_8\end{array}\right].$$

Then, ${\mathbf{C}}_8^{\left(a\right)}{=(\mathbf{H}}_2\otimes {\mathbf{I}}_4{\left)\right(\mathbf{A}}_4\oplus {\mathbf{B}}_4)$. Based on the properties of structural matrices [30,31], the matrix ${\mathbf{A}}_4$ has the following structure:

$$\left[\begin{array}{cc}{\mathbf{A}}_2^{\left(a\right)}& {\mathbf{B}}_2^{\left(\mathbf{A}\right)}\\ {\mathbf{A}}_2^{\left(\mathbf{A}\right)}& {-\mathbf{B}}_2^{\left(\mathbf{A}\right)}\end{array}\right] \mathrm{with} {\mathbf{A}}_2^{\left(\mathbf{A}\right)}{=\mathbf{H}}_2 \mathrm{and} {\mathbf{B}}_2^{\left(a\right)}{=d}_8{\mathbf{I}}_2.$$

Hence, the matrix ${\mathbf{A}}_4$ is represented as

$${\mathbf{A}}_4{=(\mathbf{H}}_2\otimes {\mathbf{I}}_2{\left)\right(\mathbf{H}}_2\oplus { \mathit{diag}(s}_0^{\left(8\right)}{, s}_1^{\left(8\right)}\left)\right)$$

where $s_0^{\left(8\right)}$=${ s}_1^{\left(8\right)}=$ $d_8$.

In the ${\mathbf{B}}_4$ matrix, we change the sign of all terms in the fourth column and fourth row. The obtained matrix ${\mathbf{B}}_4^{\left(a\right)}$ matches the pattern $\left[\begin{array}{cc}{\mathbf{A}}_2^{\left(b\right)}& {\mathbf{B}}_2^{\left(b\right)}\\ {\mathbf{B}}_2^{\left(b\right)}& {-\mathbf{A}}_2^{\left(b\right)}\end{array}\right]$. Therefore, we decompose the obtained matrix ${\mathbf{B}}_4^{\left(a\right)}$ as

$${\mathbf{B}}_4^{\left(a\right)}{=(\mathbf{T}}_{23}^{\left(4\right)}\otimes {\mathbf{I}}_2{\left)\right[(\mathbf{A}}_2^{\left(b\right)}-{ \mathbf{B}}_2^{\left(b\right)}{) \oplus (-\mathbf{A}}_2^{\left(b\right)}-{ \mathbf{B}}_2^{\left(b\right)}{) \oplus \mathbf{B}}_2^{\left(\mathbf{B}\right)}{\left]\right(\mathbf{T}}_{32}^{\left(3\right)}\otimes {\mathbf{I}}_2),$$

(29)

where ${\mathbf{T}}_{23}^{\left(4\right)}$ and ${\mathbf{T}}_{32}^{\left(3\right)}$ are defined as in (16), ${\mathbf{A}}_2^{\left(b\right)}-{ \mathbf{B}}_2^{\left(b\right)}=\left[\begin{array}{cc}{b}_8{-c}_8& b_8{-c}_8\\ b_8{-c}_8& {-b}_8{+c}_8\end{array}\right]$, ${-\mathbf{A}}_2^{\left(b\right)}{- \mathbf{B}}_2^{\left(b\right)}=\left[\begin{array}{cc}{-b}_8{-c}_8& {-b}_8{-c}_8\\ {-b}_8{-c}_8& {-b}_8{+c}_8\end{array}\right]$, ${\mathbf{B}}_2^{\left(b\right)}=\left[\begin{array}{cc}{c}_8& c_8\\ c_8& {-c}_8\end{array}\right]$.

Next, we present ${\mathbf{A}}_2^{\left(b\right)}-{ \mathbf{B}}_2^{\left(b\right)}$, ${-\mathbf{A}}_2^{\left(b\right)}{- \mathbf{B}}_2^{\left(b\right)}$, and ${\mathbf{B}}_2^{\left(b\right)}$ as [5]

$$\begin{array}{c}{\mathbf{A}}_2^{\left(b\right)}-{ \mathbf{B}}_2^{\left(b\right)}{=\mathbf{H}}_2{\times \mathit{diag}(s}_2^{\left(8\right)}{, s}_3^{\left(8\right)}{), s}_2^{\left(8\right)}{=s}_3^{\left(8\right)}{=b}_8{-c}_8,\\ {-\mathbf{A}}_2^{\left(b\right)}{- \mathbf{B}}_2^{\left(b\right)}{=\mathbf{H}}_2{\times \mathit{diag}(s}_4^{\left(8\right)}{, s}_5^{\left(8\right)}{), s}_4^{\left(8\right)}{=s}_5^{\left(8\right)}{=-b}_8{-c}_8,\\ {\mathbf{B}}_2^{\left(b\right)}{=\mathbf{H}}_2{\times \mathit{diag}(s}_6^{\left(8\right)}{, s}_7^{\left(8\right)}{), s}_6^{\left(8\right)}{=s}_7^{\left(8\right)}{=c}_8.\end{array}$$

Relying on the properties of structural matrices [5,32], the computational procedure for the eight-point DHT-II is represented by the expression

$${\mathbf{Y}}_{8\times 1}{=\mathbf{W}}_8{\mathbf{W}}_{8\times 10}{\mathbf{W}}_{10}{\mathbf{D}}_{10}{\mathbf{W}}_{10\times 8}{\mathbf{P}}_8{\mathbf{X}}_{8\times 1},$$

(30)

where

$$\begin{array}{c}{\mathbf{W}}_8{=\mathbf{H}}_2\otimes {\mathbf{I}}_4, {\mathbf{W}}_{8\times 10}{=(\mathbf{H}}_2\otimes {\mathbf{I}}_2{)\oplus (\mathbf{T}}_{23}^{\left(4\right)}\otimes {\mathbf{I}}_2), {\mathbf{W}}_{10\times 8}{=\mathbf{I}}_4\oplus {(\mathbf{T}}_{3\times 2}^{\left(3\right)}\otimes {\mathbf{I}}_2),\\ {\mathbf{W}}_{10}{=\mathbf{H}}_2\oplus {\mathbf{I}}_2\oplus {\mathbf{H}}_2\oplus {\mathbf{H}}_2\oplus {\mathbf{H}}_2, {\mathbf{D}}_{10}{=\mathit{diag}(1, 1, s}_0^{\left(8\right)}{, s}_1^{\left(8\right)}{, s}_2^{\left(8\right)}{, s}_3^{\left(8\right)}{, s}_4^{\left(8\right)}{, s}_5^{\left(8\right)}{, s}_7^{\left(8\right)}),\\ {\mathbf{W}}_8=\left[\begin{array}{cccccccc}1& & & & 1& & & \\ & 1& & & & 1& & \\ & & 1& & & & 1& \\ & & & 1& & & & -1\\ 1& & & & -1& & & \\ & 1& & & & -1& & \\ & & 1& & & & -1& \\ & & & 1& & & & 1\end{array}\right],\\ {\mathbf{P}}_8=\left[\begin{array}{cccccccc}1& & & & & & & \\ & & & & 1& & & \\ & & 1& & & & & \\ & & & & & & 1& \\ & 1& & & & & & \\ & & & 1& & & & \\ & & & & & 1& & \\ & & & & & & & -1\end{array}\right].\end{array}$$

A data flow graph of the proposed algorithm of the eight-point DHT-II is presented in Figure 7. Based on this algorithm, the number of multiplications can be reduced from 40 to 8 and the number of additions from 48 to 26.

4. Results

The experiment confirmed the correctness of the DHT-II matrix factorizations obtained in Section 3. Namely, the coinciding results of the multiplication of factorized matrices from Section 3 (Equations (7), (10), (12), (18), (21), (27), and (30)) with the original DHT-II matrix was evaluated using the MATLAB R2023b software. For this, the DHT-II matrices were directly calculated for different input sequence sizes N using the expression (3) from Section 3. Next, the factorized matrices were obtained using the expressions (7), (10), (12), (18), (21), (27), and (30) for N = 2, 3, 4, 5, 6, 7, and 8, respectively. Further, the results of the direct calculation of the DHT-II matrix were compared to the results of the multiplication of matrices included in the DHT-II factorizations. The results’ consistency was evaluated, indicating the correctness of the developed factorizations.

Based on the DHT-II matrix factorization, the proposed fast algorithms of the DHT-II calculation were represented in Section 3 by the data flow graphs in Figure 1, Figure 2, Figure 3, Figure 4, Figure 5, Figure 6 and Figure 7. To test the correctness of the constructed data flow graphs for the proposed algorithms, 210 sequences of random normally distributed numbers of length N from 2 to 8 were generated. Each sequence was transformed by a direct matrix–vector product using a data flow graph. The coinciding results obtained for each N indicated the correctness of the constructed data flow graphs for the proposed algorithms.

Next, the numbers of multiplications and additions for the direct DHT-II calculation method and the proposed algorithms were estimated and compared. To ensure the accuracy of the estimates, zero elements in the base DHT-II matrices were not taken into account. This meant that the estimate of the multiplicative complexity of the calculations was calculated using the formula N²-Q and the additive complexity was calculated using the formula N(N − 1)-Q, where Q is the number of positions with zero elements in the base transform matrix. So, in the matrix ${\mathbf{C}}_4$, four zeros appeared; in the matrix ${\mathbf{C}}_8,$ eight zeros were found. Overall, the number of multiplications on the data flow graphs was determined by the number of graph vertices, indicated by circles with factors. The number of additions was determined by the number of occurrences of two graph edges in one vertex. The proposed algorithms significantly reduced the number of multiplications in the DHT-II algorithms for a 1D signal, with the number of samples in the range of 3 to 8. At the same time, the number of additions was slightly reduced or even increased. The number of additions was reduced by an average of 5% and the number of multiplications was reduced by an average of 73% with the number of samples in the range of 3 to 8. The obtained results are shown in

Table 1. The number of operations of the direct method and the proposed algorithms.

N	Direct Method		Proposed Algorithms
	Additions	Multiplications	Additions	Multiplications
2	2	0	2 (0%)	0 (0%)
3	6	4	7 (+16%)	1 (−75%)
4	8	4	6 (−25%)	2 (−50%)
5	20	16	23 (+15%)	5 (−69%)
6	30	16	26 (−13%)	2 (−88%)
7	42	36	46 (+10%)	8 (−78%)
8	48	40	26 (−35%)	8 (−80%)

. The percentage difference in the number of operations is indicated in parentheses, a plus means that the number of operations increased compared to the direct method, and a minus means it was reduced.

The numbers of multiplications and additions for the algorithms known from the literature are shown in

Table 2. The number of operations for the algorithms known from the literature.

Algorithm	Reference, Year of Publication	N = 4		N = 8
		Multiplications	Additions	Multiplications	Additions
Hu’s radix-2 DIT algorithm	[24], 2007	2 (0%)	18 (−67%)	12 (−33%)	48 (−46%)
Prime factor algorithm	[23], 2004	2 (0%)	10 (−40%)	8 (0%)	28 (−7%)
Radix-2 algorithm	[24], 2007	4 (−50%)	18 (−67%)	14 (−43%)	42 (−38%)
Radix-2 algorithm	[25], 2011	–	–	10 (−20%)	28 (−7%)
Radix-2 DIT algorithm	[26], 2016	2 (0%)	6 (0%)	12 (−33%)	24 (+8%)
Split-radix algorithm	[27], 2020	2 (0%)	6 (0%)	12 (−33%)	24 (+8%)
Proposed algorithm	–	2	6	8	26

. Also, the percentage difference in the number of operations is indicated in parentheses.

5. Discussion of Computational Complexity

Analyzing the obtained results, we note that for N = 4, the proposed algorithm did not differ from the existing algorithms by the number of multiplications. If N = 8, the number of multiplications may be reduced with the proposed solution by 33% compared to Hu’s algorithm [24] and Hamood’s algorithms [26,27]. Compared to the prime factor algorithm [23], the number of multiplications for N = 8 was not reduced. The number of additions for N = 8 was reduced by 7–46% compared to Hu’s algorithm [24], the prime factor algorithm [23], and the radix-2 algorithms [24,25]. Compared to Hamood’s algorithms [26,27], the number of additions for N = 4 was not reduced and the number of additions for N = 8 was increased by 8%. For N = 4, the number of additions was reduced by 40% and 67% compared to the algorithms from [23,24], respectively.

The constructive difference between the proposed DHT-II algorithms and the existing algorithms is that the critical path in the signal flow graph of any of the presented algorithms contains only one multiplication. If there is more than one multiplication on the critical path of the algorithm, then the multiplication of two n-bit operands results in a 2n-bit product. Repeated multiplication requires more manipulation of the operands and, therefore, requires more time and effort compared to multiplication alone. In fixed-point devices, this may result in overflow–underflow processing. To maintain accuracy, double access to memory is required both when writing and when reading. Using floating point arithmetic requires exponent alignment, mantissa addition, etc. Thus, the presence of more than one multiplication on the critical path for the known algorithms is a noticeable drawback.

Through this research on the proposed DHT-II algorithms, the following shortcomings of the traditional algorithms were determined. The number of additions and multiplications for the algorithms described in the literature, in general, is greater compared to the proposed algorithms. In addition, traditional algorithms often limit the length of the input data sequence to powers of two or three. However, the application of the proposed algorithms may be limited by the following. First, the structural approach is more suitable for constructing fast algorithms for short data sequences. For long data sequences, the structure of the transformation matrices is more difficult to identify. Second, the efficiency of synthesizing fast algorithms based on the structural approach depends on the structural properties of the transformation matrices, namely, how efficiently these matrices can be transformed so that their block structure corresponds to the matrix patterns formulated in [5].

6. Conclusions

As a result of research, the problem of designing the fast algorithms of type-II DHT was solved for short-length input signals based on a structural approach. The algorithms for short-length input signals are of special interest because they are subsequently used as building blocks of algorithms for long-length sequences.

After a literature review, it was decided to use the structural properties of the type-II DHT coefficients matrix to develop fast algorithms for this transform. The proposed algorithms differ from existing solutions by applying pre-decomposition of the original matrices into submatrices, as well as rearranging the rows and/or columns of the resulting matrices to find the factorizations that would reduce the computational complexity. The relations between the values of the matrix entries were also considered. In an experiment, the computational complexity of the proposed algorithms was compared with the direct calculation of matrix–vector products. The resulting DHT-II matrix factorizations reduced the number of additions by 5% on average and the number of multiplications by 73% for DHT-II lengths in the range of 3 to 8.

As a result, the amount of resources used on the signal processor was significantly reduced while allowing for easier operation in real-time. A significant contribution to the speedup of signal processing came from the reduction in multiplications since multiplications are more labor-intensive and time-expensive to use than additions. To achieve computational efficiency, the structure of the transform matrix and not the parameters of DHT-II were optimized.

A comparison of the proposed algorithms based on the structural approach with algorithms known from the literature showed that the efficiency of the proposed algorithms significantly exceeded the efficiency of radix or split-radix algorithms and was close to the efficiency of prime factor algorithms.

As for prospects for further research, we note the following. In [38], as an example, the application of the Hartley transform in wavelet-like decomposition was shown. The application of known orthogonal transformations in multiscale decomposition and reconstruction significantly reduced computational costs and gave the same effect as wavelet decomposition and reconstruction [39,40]. Thus, we assume that the wavelet-like decomposition can be applied to increase the performance of texture image segmentation for skin disease diagnosis [41].