Fast Algorithms for Short-Length Odd-Time and Odd-Frequency Discrete Hartley Transforms

Abstract

This paper proposes fast algorithms for odd-time and odd-frequency discrete Hartley transforms (DHTs) with normalization if the size of the input sequence ranges from 2 to 8. Fast algorithms for small-sized inputs can be included as modules in fast algorithms for discrete transforms of large-length input sequences. The existing fast odd-time DHT and odd-frequency DHT algorithms are primarily radix-type algorithms, specifically, radix-2 and prime factor algorithms. However, the algorithms in the literature do not normalize the initial transform. This means that after applying such algorithms, additional N multiplications are required. In this paper, the structural approach is exploited supposing that the starting point to design the fast algorithms is a matrix vector product expression. The factorization of the matrices of DHT coefficients is produced to reduce computational complexity taking into account the repetition and arranging of the matrix entries. The strict mathematical background proves the correctness of the obtained algorithmic solutions. Also, the MATLAB R2023b environment was applied for testing of the performance of the proposed algorithms. The obtained factorizations of the odd-time DHT and odd-frequency DHT matrices allow us to reduce the number of multiplications by 69%, while the amount of additions is decreased by about of 5% for odd-time DHTs and 8% for odd-frequency DHTs if the length of input data is in the range from 2 to 8. A comparison is provided with the direct calculation of the matrix vector product. Additionally, the computational complexity for each obtained solution is compared with the computational complexity of the existing fast algorithms for odd-time DHTs and odd-frequency DHTs. In addition, data flow graphs have been created for the proposed odd-time DHT and odd-frequency DHT algorithms. The modular space–time structure of the resulting data flow graphs is suitable for VLSI implementation.

1. Introduction

The DHT is widely applied in digital signal and image processing to obtain the frequency representation of data [1,2,3]. This transform is used in image compression [4], smart grid applications [5], wireless communication [6], audio watermarking [7], digital image watermarking [8,9,10], signal processing [11], and other applications. This transform can successfully supplement the discrete cosine transform (DCT) in different applications due to additional exploiting of the sine transform [1,2,3]. The latter allows for the processing of the signal details at high noise levels with better performance. For example, in [12], the DHT is applied in three-dimensional (3D) medical image compression. X-ray angiogram (XA) and magnetic resonance (MR) data with varied resolutions and with different numbers of frames/slices were processed. Each XA or MR data point was represented by a 3D array which was divided into blocks of size 8 × 8 × 8, and then a 3D DHT and a 3D zigzag scheme were employed. For comparison, the 3D DCT was applied with the corresponding zigzag pattern [13]. In the full 3D DHT-based codec implementation, only the first L zig-zag coefficients are preserved, encoded by entropy coding [14], and then transmitted and/or stored. However, to study the effect of the 3D DHT block, the bitrate in bit per voxel (bpv) must be calculated before entropy coding. As a result, both transforms have shown better performance at compressing XA data than MR data. Compared with the DCT on the XA data [15], approximate versions of the exact 3D DHT and the one proposed in [12] present better performance than the 3D DCT in terms of peak noise–signal ratio for the low bitrate ranging from 1 bpv to 6 bpv.

To increase robustness against various noise disturbances, similar advantages of the DHT are employed in image and audio watermarking [8,9,10]. Additionally, a significant reduction in the multiplicative complexity allows us to increase the computing efficiency and save resources [12].

To date, different DHT types have been exploited, specifically a generalized DHT [16,17], sliding DHT [18], two-dimensional DHT [19,20], 3D DHT [12], etc. However, the traditional DHT is still the most developed and researched transform compared with other types of DHT transforms [21,22,23,24]. In our paper, the odd-time DHT and odd-frequency DHT are considered. The odd-frequency DHT allows for data processing of frequencies that are centred in the range of those employed by the traditional DHT.

Direct computation of the different types of DHTs requires a huge number of arithmetic operations. Specifically, direct computation of the odd-time DHT or odd-frequency DHT of an N-point sequence requires $N^2$ multiplications which are especially time- and resource-consuming. To achieve a reduction in computational complexity, fast DHT algorithms have been designed [23,24,25]. But in articles devoted to the efficient implementation of the odd-time DHT and odd-frequency DHT, large-length inputs are considered [17,25,26,27,28,29,30,31]. It should be noted that more complex fast algorithms can be designed with modules containing odd-time DHT and odd-frequency DHT algorithms for short-length sequences [17]. In particular, in [21], a $2^{m-n}$-point odd-frequency DHT is used to develop a split algorithm for a $2^m$-point traditional DHT with reduced arithmetic complexity, where n is changed in the range from 1 to m. Once constructed, odd-time DHT and odd-frequency DHT algorithms for small-size inputs can be successfully exploited in different applications. As such, the development of fast odd-time DHT and odd-frequency DHT algorithms for short-length input sequences is a topical issue, and the appropriate approach for solving the determined parts of a general problem should be selected.

1.1. Related Papers

Most of the fast odd-time DHT and odd-frequency DHT algorithms are designed based on the divide-and-conquer strategy [17], which means obtaining a solution of a large problem via a set of solutions of smaller tasks. Bearing in mind this strategy, radix-type algorithms were developed including radix-2 and prime factor algorithms.

Thus, the radix-2 fast odd-time DHT and odd-frequency DHT algorithms were proposed in [17,26,27,28,29]. In [26,27] the proposed fast algorithms decompose the N-point odd-time DHT into two adjacent N/2-point odd-time DHTs. Such fast odd-time DHT algorithms save from 16% to 24% of the arithmetic operations in comparison with Hu’s algorithm if N varies from 16 to 64 [26].

In [17] the N-point odd-frequency DHT is decomposed into the N/2-point odd-frequency DHT and the N/2-point odd-time odd-frequency DHT. For odd-time DHT similar algorithms are also proposed. As a result, the computational complexity of the N-point odd-time DHT and odd-frequency DHT is reduced. In [28] the odd-frequency DHT of length N = $2^m$ was represented by several 16-point odd-frequency DHTs. Then the 16-point odd-frequency DHT is calculated by first-order moments. As a consequence, a needed number of additions and multiplications is lower than in other existing algorithms.

In [29,30] radix-2 fast algorithms are developed for the VLSI implementation of odd-time DHT and odd-frequency DHT of length N = $2^m$. The algorithms differ by parallelism as well as regular and simple computational structure with low arithmetic cost.

Developed in [26,27,28,29] the fast odd-time DHT and odd-frequency DHT algorithms apply the decimation-in-frequency approach. However, for real-time applications the fast algorithms based on the decimation-in-time approach are preferable. Presented in [31] the radix-2 fast odd-time DHT algorithm with decimation-in-time substantially saves the arithmetic operations relative to Hu’s algorithm. Also, in [31] a regular signal flow graph is constructed which can be used for odd-time DHT with different lengths of inputs.

In [32] a split-radix algorithm with decimation-in-time splits N-point odd-time DHT into N/2-point odd-time DHT for even indexed samples and two odd-time DHTs of length N/4 for odd-indexed samples. Thus the number of multiplications and additions was reduced up to 29% and 10% respectively in comparison with the radix-2 odd-time DHT algorithm [31].

Mentioned radix-2 odd-time DHT and odd-frequency DHT algorithms are cost-effective but when the input sequence length is not a power of two the computational overhead frequently arises. Thus the prime factor algorithms are introduced in [17] for sequence lengths N = p × q, where p and q are relatively prime. The first prime factor algorithm decomposes the N-point odd-time DHT into p q-point odd-time DHT and q p-point odd-time DCTs. A similar decomposition also takes place for odd-frequency DHT by replacing odd-time transforms with odd-frequency transforms. However, in this case, irregularity in the index mapping occurs. To implement the in-place computation, which is easily attained by the radix-2 algorithms, the second prime factor algorithm was developed [17]. The N-point odd-time DHT is decomposed into N/q q-point odd-frequency DCTs and q N/q-point odd-time DHTs. The N-point odd-frequency DHT computation is decomposed into q N/q-point odd-frequency DHTs and N/q q-point odd-time DCTs.

Along with the algorithms synthesized based on the divide-and-conquer strategy, the algebraic signal processing is taken into account to develop the fast odd-time DHT and odd-frequency DHT algorithms. Thus, in [33] the polynomial transforms, linear congruencies, and index permutations are applied to design the fast odd-time DHT and odd-frequency DHT algorithms. In [11] fast general-radix algorithms for a wide class of orthogonal trigonometric transforms are provided. This transform class includes the odd-time DHT and odd-frequency DHTs and related algorithms are obtained by manipulating the polynomial algebras instead of manipulating the transform expression.

The short review of existing fast odd-time DHT and odd-frequency DHT algorithms has allowed defining their drawbacks and formulation of unsolved parts of the general problem of reducing computational complexity.

1.2. The Main Contributions of the Paper

After review, we note that the existing fast odd-time DHT and odd-frequency DHT algorithms are developed using polynomial algebras or radix-type approach. However, several papers devoted to fast algorithms for orthogonal trigonometric transforms exploit a structural approach that compares the transform coefficient matrix with matrix patterns [34,35,36,37,38,39]. As an option, these matrix patterns are introduced in [34] based on Winograd’s approach to the fast algorithm development. The advantage of the structural approach is the factorization of the transform matrix on sparse, scaled orthogonal, rotational, or rotational-reflection matrices [38,39]. Moreover, based on the structural approach we can obtain the data flow graphs of odd-time DHT and odd-frequency DHT algorithms with only one multiplication along the critical path [34,35,36,37]. If there is more than one multiplication in the critical path of the algorithm, then this will create additional problems for the implementation of computations.

In [37] based on the structural approach efficient algorithms for small-sized odd-time DHT were developed that do not consider the normalizing coefficients. Many authors do this, naively thinking that the result can be multiplied by the normalizing constants at subsequent stages of data processing and combined with other calculations. However, this trick cannot always be implemented, and then the normalization still has to be introduced leading to the need to perform N additional multiplications. This means that to obtain true quantitative estimates of the computational complexity of the proposed solutions, it is necessary to include N additional multiplications for each N-point algorithm, otherwise such estimates will be incorrect. To avoid such shortcomings, it is necessary to consider the normalizing constants in the original transform matrix at the very beginning.

The article [34] describes a general approach to constructing reduced complexity algorithms for calculating matrix-vector products in those cases where the matrices have certain structures. In the paper, we apply this approach to specific transforms, namely, odd-time DHT and odd-frequency DHT. Taking into account the normalization of these transforms is influenced by the number of arithmetic operations, and in some cases on the matrix structure. As a result, the mentioned approach can also allow us to obtain new interesting solutions if normalization is kept in mind.

Therefore, in this article, we consider the case of factorization of the original transformation matrices, the entries of which are already initially multiplied by the corresponding normalizing constants. This paper aims to reduce the computational complexity of the odd-time DHT and odd-frequency DHT with normalization via developing fast algorithms. We base on the structural approach supposing that input sequence length N is changed from 2 to 8.

The article is organized as follows. The considered problem is introduced and the research subject is formulated in Section 1. Brief mathematical background is provided and the used notations are defined in Section 2. Section 3 is devoted to odd-time DHT and odd-frequency DHT matrix factorizations for different N from 2 to 8 and following obtaining of the fast odd-time DHT and odd-frequency DHT algorithms via the data flow graphs. In Section 4 and Section 5 the computational complexity of the proposed fast odd-time DHT and odd-frequency DHT algorithms is evaluated. In Section 6 the conclusions are provided.

2. Preliminary Background

The traditional DHT is expressed as follows [17]:

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

where casθ = sinθ + cosθ = $\sqrt{2}$cos(θ$-$π/4); $y_k$ is the output data sequence after the traditional DHT; $x_n$ is the sequence of input data; N is the number of signal samples.

The direct generalized DHT can be described by the following expression [17]:

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

where $y_k$ is the output sequence after the direct generalized DHT; $n_0$ is a time shift; $k_0$ is a frequency shift.

Using the expression (2) we have the traditional DHT if $n_0=0, k_0=0,$ the odd-time DHT if $n_0=1/2, k_0=0,$ the odd-frequency DHT if $n_0=0, k_0=1/2,$ the odd-time odd-frequency DHT if $n_0=1/2, k_0=1/2.$ When for the odd-time DHT with normalization the expression (2) is represented as follows [17]:

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

where $y_k$ is the output sequence after the odd-time DHT.

For the odd-frequency DHT with normalization we derive the expression (2) as follows [17]:

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

where $y_k$ is the odd-frequency DHT outputs, $k=0, 1, \dots , N-1$.

The odd-time DHT and odd-frequency DHT can be obtained as a matrix-vector product:

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

where ${\mathbf{C}}_N=\left[\begin{array}{ccc}\begin{array}{cc}{c}_{0,0}& c_{0,1}\end{array}& \cdots & c_{0,N-1}\\ \begin{array}{cc}{c}_{1,0}& c_{1,1}\end{array}& \cdots & c_{1,N-1}\\ ⋮& \ddots & ⋮\\ \begin{array}{cc}{c}_{N-1,0}& c_{N-1,1}\end{array}& \cdots & c_{N-1,N-1}\end{array}\right]$, $c_{kn}={\displaystyle \frac{1}{\sqrt{N}}}\mathrm{cas}{\displaystyle \frac{\pi \left(2n+1\right)k}{N}}$ for odd-time DHT, and $c_{kn}={\displaystyle \frac{1}{\sqrt{N}}}\mathrm{cas}{\displaystyle \frac{\pi \left(2k+1\right)n}{N}}$ for odd-frequency DHT, n, k = 0, …, N − 1; ${\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}}$.

Because the generalized DHT with normalization is orthonormal then the matrices of inverse transforms are transposed to the direct ones. That is, for odd-time DHT the inverse is odd-frequency DHT and vice versa. For traditional DHT the inverse transform is the same [40].

In the article the following notations and abbreviations are used [34]:

• ⊕ is the direct sum of two matrices;

• ${\mathbf{I}}_N$ is an order N identity matrix;

• ⊗ is a Kronecker product of two matrices;

• ${\mathbf{H}}_2$ is a 2 × 2 Hadamard matrix;

• $s_m^{(N)}$ is a multiplier;

• ${\mathbf{P}}_N$ is a permutation matrix;

• ${\mathbf{W}}_{M\times N}$ is a block-diagonal matrix of size M × N;

• ${\mathbf{W}}_N$, ${\mathbf{D}}_{\mathit{N}}$ are block-diagonal and diagonal matrices of size N × N, respectively;

• an empty cell in a matrix means that it contains zero;

• DHT is a discrete Hartley transform;

• VLSI is very large-scale integration;

• DCT is a discrete cosine transform;

• XA is a X-ray angiogram;

• MR is magnetic resonance;

• 3D is a three-dimensional data.

3. Reduced Complexity Algorithms for Short-Length Odd-Time DHT and Odd-Frequency DHT

3.1. Algorithms for 2-Point Odd-Time DHT and Odd-Frequency DHT

To develop the algorithm for two-point odd-time DHT with normalization we consider the expression:

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

where ${\mathbf{C}}_2=\left[\begin{array}{cc}{a}_2& a_2\\ a_2& -a_2\end{array}\right]$, $a_2={\displaystyle \frac{1}{\sqrt{2}}}\mathrm{cas}0={\displaystyle \frac{1}{\sqrt{2}}}\approx$ 0.7071, ${\mathit{b}}_3={\displaystyle \frac{1}{\sqrt{3}}}\mathrm{cas}{\displaystyle \frac{\pi }{3}}\approx$ 0.7887, $c_3={\displaystyle \frac{1}{\sqrt{3}}}\mathrm{cas}{\displaystyle \frac{2\pi }{3}}\approx$ 0.2113, ${\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}}$. The expression (6) for two-point odd-frequency DHT with normalization is the same because ${\mathbf{C}}_2={\mathbf{C}}_2^{\mathbf{T}}$.

Then the factorization of the odd-time and odd-frequency DHT matrix for N = 2 can be presented as follows:

(7)${\mathbf{Y}}_{2\times 1}={\mathbf{D}}_2{\mathbf{H}}_2{\mathbf{X}}_{2\times 1},$

where ${\mathbf{D}}_2=\mathrm{diag}\left(s_0^{(2)}, s_1^{(2)}\right),s_0^{(2)}=s_1^{(2)}=a_2$.

We have shown a data flow graph of the algorithm for the 2-point odd-time DHT and odd-frequency DHT in Figure 1. As can be seen, we have 2 additions as well as when using the direct method. The number of multiplication is also 2 against 4 multiplications for the direct matrix-vector product.

3.2. Algorithms for 3-Point Odd-Time DHT and Odd-Frequency DHT

The expression for 3-point odd-time DHT with normalization is as follows:

(8)${\mathbf{Y}}_{3\times 1}={\mathrm{C}}_3{\mathbf{X}}_{3\times 1},$

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}{a}_3& a_3& a_3\\ b_3& -a_3& -c_3\\ c_3& a_3& -b_3\end{array}\right]$ with $a_3={\displaystyle \frac{1}{\sqrt{3}}}\mathrm{cas}0={\displaystyle \frac{1}{\sqrt{3}}}\approx$ 0.5774, $b_3={\displaystyle \frac{1}{\sqrt{3}}}\mathrm{cas}{\displaystyle \frac{\pi }{3}}\approx$ 0.7887, $c_3={\displaystyle \frac{1}{\sqrt{3}}}\mathrm{cas}{\displaystyle \frac{2\pi }{3}}\approx$ 0.2113. For odd-frequency DHT with normalization we use ${\mathbf{C}}_3^{\mathrm{T}}$ instead of ${\mathbf{C}}_3$ in Equation (8).

To factorize the ${\mathbf{C}}_3$ matrix we swap the first and second rows. Then we alter the sign of the third row and third column of the resulted matrix ${\mathbf{C}}_3^{(a)}$ and decompose it into two components:

(9)${\mathbf{C}}_3^{(a)}={\mathbf{C}}_3^{(b)}+{\mathbf{C}}_3^{(c)},$

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

In the second column and second row of the matrix ${\mathbf{C}}_3^{(b)}$ we obtained identical entries that differ only in sign. Therefore, the number of operations can be reduced without additional transforms. Next, we eliminate the second row and second column containing only zero entries in matrix ${\mathbf{C}}_3^{(c)}$, and matrix ${\mathbf{C}}_2^{(a)}=\left[\begin{array}{cc}{b}_3& c_3\\ -c_3& -b_3\end{array}\right]$ was obtained.

We derive the factorization for the 3-point odd-time DHT based on the decomposition of matrix ${\mathbf{C}}_2^{(a)}$ [34,37]:

(10)${\mathbf{Y}}_{3\times 1}={\mathbf{P}}_3{\mathbf{W}}_{3\times 4}{\mathbf{W}}_4{\mathbf{D}}_4{\mathbf{W}}_{4\times 3}{\mathbf{W}}_3{\mathbf{X}}_{3\times 1},$

where ${\mathbf{D}}_4=\mathrm{diag}\left(s_0^{(3)}, s_1^{(3)}, s_2^{(3)}, s_3^{(3)}\right)$, $s_0^{(3)}=(b_3+c_3)/2$, $s_1^{(3)}=(b_3-c_3)/2$, $s_2^{(3)}=s_3^{(3)}=-a_3$,

${ \mathbf{W}}_4=({\stackrel{\mathrm{-}}{\mathbf{I}}}_2{\mathbf{H}}_2)\oplus {\mathbf{I}}_2,$

${\mathbf{P}}_3=\left[\begin{array}{ccc}& 1& \\ 1& & \\ & & -1\end{array}\right], {\mathbf{W}}_{3\times 4}=\left[\begin{array}{cccc}1& & 1& \\ & & -1& -1\\ & 1& 1& \end{array}\right], {\mathbf{W}}_3=\left[\begin{array}{ccc}1& & -1\\ & 1& \\ 1& & 1\end{array}\right], {\mathbf{W}}_{4\times 3}=\left[\begin{array}{ccc}1& & \\ & & 1\\ & 1& \\ & & 1\end{array}\right].$

Based on Equation (10) the factorization for the 3-point odd-frequency DHT is obtained as:

(11)${\mathbf{Y}}_{3\times 1}={\mathbf{W}}_3^{\mathrm{T}}{\mathbf{W}}_{4\times 3}^{\mathrm{T}}{\mathbf{D}}_4{\mathbf{W}}_4^{\mathrm{T}}{\mathbf{W}}_{3\times 4}^{\mathrm{T}}{\mathbf{P}}_3^{\mathrm{T}}{\mathbf{X}}_{3\times 1}.$

A data flow graph of the 3-point odd-time DHT algorithm is shown in Figure 2. Then the amount of multiplications can be reduced from 9 to 3 because $s_0^{(3)}=0.5$. However, the number of additions enlarges from 6 to 7.

A data flow graph of the 3-point odd-frequency DHT algorithm is shown in Figure 3. The number of arithmetic operations remains the same as for 3-point odd-time DHT.

3.3. Algorithms for 4-Point Odd-Time DHT and Odd-Frequency DHT

The four-point odd-time DHT with normalization is defined as follows:

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

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}{a}_4& a_4& a_4& a_4\\ b_4& 0& -b_4& 0\\ a_4& -a_4& a_4& -a_4\\ 0& b_4& 0& -b_4\end{array}\right]$ with $a_4={\displaystyle \frac{1}{2}}\mathrm{cas}0=$ 0.5, $b_4={\displaystyle \frac{1}{2}}\mathrm{cas}{\displaystyle \frac{\pi }{4}}\approx$ 0.7071. For odd-frequency DHT with normalization we use ${\mathbf{C}}_4^{\mathrm{T}}$ instead of ${\mathbf{C}}_4$ in Equation (12).

We permutate the rows of ${\mathbf{C}}_4$ according to ${\pi }_4=\left(\begin{array}{cc}1& 2\\ 1& 3\end{array}\begin{array}{cc}3& 4\\ 2& 4\end{array}\right)$. The permutation matrix is ${\mathbf{P}}_4=\left[\begin{array}{cccc}1& & & \\ & & 1& \\ & 1& & \\ & & & 1\end{array}\right]$. The obtained matrix ${\mathbf{C}}_4^{(a)}=\left[\begin{array}{cccc}{a}_4& a_4& a_4& a_4\\ a_4& -a_4& a_4& -a_4\\ b_4& 0& -b_4& 0\\ 0& b_4& 0& -b_4\end{array}\right]$ matches the pattern $\left[\begin{array}{cc}{\mathbf{A}}_2& {\mathbf{A}}_2\\ {\mathbf{B}}_2& -{\mathbf{B}}_2\end{array}\right]$ where ${\mathbf{A}}_2=\left[\begin{array}{cc}{a}_4& a_4\\ a_4& -a_4\end{array}\right]$ and ${\mathbf{B}}_2=\left[\begin{array}{cc}{b}_4& \\ & b_4\end{array}\right]$.

Applying the matrix templates from [34,37], we derive the expressions

(13)${\mathbf{C}}_4^{(a)}=({\mathbf{A}}_2\oplus {\mathbf{B}}_2)({\mathbf{H}}_2\otimes {\mathbf{I}}_2), {\mathbf{A}}_2=(a_4\oplus a_4){\mathbf{H}}_2, {\mathbf{B}}_2=(b_4\oplus b_4){\mathbf{I}}_2.$

Therefore the factorization of the 4-point odd-time DHT matrix is as follows:

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

where ${\mathbf{W}}_4^{(0)}={\mathbf{H}}_2\oplus {\mathbf{I}}_2$, ${\mathbf{D}}_4=\mathrm{diag}\left(s_0^{(4)}, s_1^{(4)}, s_2^{(4)}, s_3^{(4)}\right)$, $s_0^{(4)}=s_1^{(4)}=a_4, s_2^{(4)}=s_3^{(4)}=b_4,$ ${\mathbf{W}}_4^{(1)}={\mathbf{H}}_2\otimes {\mathbf{I}}_2$.

The data flow graph of the 4-point odd-time DHT algorithm is presented in Figure 4. As one can see, the number of multiplications is reduced from 4 to 2 because $s_0^{(4)}=s_1^{(4)}=a_4=0.5$, the number of additions is decreased from 8 to 6 compared to direct matrix-vector product.

Taking into consideration the Equation (14), the factorization for the four-point odd-frequency DHT matrix is expressed as

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

The algorithm of the 4-point odd-frequency DHT with normalization is provided as a data flow graph in Figure 5. The number of arithmetic operations is the same as for odd-time DHT.

3.4. Algorithms for 5-Point Odd-Time DHT and Odd-Frequency DHT

To develop the algorithm for five-point odd-time DHT with normalization we use the following expression

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

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}}$, $a_5$ = ${\displaystyle \frac{1}{\sqrt{5}}}\approx$ 0.4472, $b_5$ =${\displaystyle \frac{1}{\sqrt{5}}}\mathrm{cas}{\displaystyle \frac{\pi }{5}}\approx$0.6247, $c_5$= ${\displaystyle \frac{1}{\sqrt{5}}}\mathrm{cas}{\displaystyle \frac{3\pi }{5}}\approx$ 0.2871, $d_5$ = ${\displaystyle \frac{1}{\sqrt{5}}}\mathrm{cas}{\displaystyle \frac{2\pi }{5}}\approx$ 0.5635, $e_5$= $-{\displaystyle \frac{1}{\sqrt{5}}}\mathrm{cas}{\displaystyle \frac{4\pi }{5}}\approx$ 0.0989,

${\mathbf{C}}_5=\left[\begin{array}{ccccc}{a}_5& a_5& a_5& a_5& a_5\\ b_5& c_5& -a_5& -d_5& e_5\\ d_5& -b_5& a_5& -e_5& -c_5\\ c_5& e_5& -a_5& b_5& -d_5\\ -e_5& d_5& a_5& -c_5& -b_5\end{array}\right].$

For odd-frequency DHT with normalization in Equation (16) we use ${\mathbf{C}}_5^{\mathrm{T}}$ instead of ${\mathbf{C}}_5$.

We alter the order of rows and columns of ${\mathbf{C}}_5$ applying the permutations

${\pi }_5=\left(\begin{array}{cc}1& 2\\ 2& 3\end{array}\begin{array}{ccc}3& 4& 5\\ 1& 5& 4\end{array}\right), {\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 corresponding permutation matrices are as follows

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

Then we decompose the resulting matrix ${\mathbf{C}}_5^{(a)}$ into two submatrices:

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

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

In the third column and third row of the matrix ${\mathbf{C}}_5^{(b)}$ it is obtained identical entries that differ only in sign. Therefore, the number of operations can be reduced without additional transforms [36]. We eliminate in matrix ${\mathbf{C}}_5^{(c)}$ the row and column with only zero entries. The resulted matrix

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

matches the matrix template $\left[\begin{array}{cc}{\mathbf{A}}_2^{(a)}& {\mathbf{B}}_2^{(a)}\\ -{\mathbf{B}}_2^{(a)}& -{\mathbf{A}}_2^{(a)}\end{array}\right]$ where ${\mathbf{A}}_2^{(a)}=\left[\begin{array}{cc}{b}_5& c_5\\ d_5& -b_5\end{array}\right]$, ${\mathbf{B}}_2^{(a)}=\left[\begin{array}{cc}{e}_5& -d_5\\ -c_5& -e_5\end{array}\right]$.

The factorization of template $\left[\begin{array}{cc}{\mathbf{A}}_2^{(a)}& {\mathbf{B}}_2^{(a)}\\ -{\mathbf{B}}_2^{(a)}& -{\mathbf{A}}_2^{(a)}\end{array}\right]$ is known from [34], and we represent the matrix ${\mathbf{C}}_4^{(d)}$ as

(18)${\mathbf{C}}_4^{(d)}=({\stackrel{\mathrm{-}}{\mathbf{I}}}_2 \otimes {\mathbf{I}}_2)({\mathbf{H}}_2 \otimes {\mathbf{I}}_2)1/2[({\mathbf{A}}_2^{(a)}+{\mathbf{B}}_2^{(a)}) \oplus ({\mathbf{A}}_2^{(a)}-{\mathbf{B}}_2^{(a)})]({\mathbf{H}}_2 \otimes {\mathbf{I}}_2),$

where ${\stackrel{\mathrm{-}}{\mathbf{I}}}_2=\left[\begin{array}{cc}1& 0\\ 0& -1\end{array}\right]$. The resulting matrices ${\mathbf{A}}_2^{(a)}+{\mathbf{B}}_2^{(a)}=\left[\begin{array}{cc}{b}_5+e_5& c_5-d_5\\ d_5-c_5& -b_5-e_5\end{array}\right]$ and ${\mathbf{A}}_2^{(a)}-{\mathbf{B}}_2^{(a)}=\left[\begin{array}{cc}{b}_5-e_5& c_5+d_5\\ d_5+c_5& -b_5+e_5\end{array}\right]$ match the templates $\left[\begin{array}{cc}g& h\\ -h& -g\end{array}\right]$ and $\left[\begin{array}{cc}p& q\\ q& -p\end{array}\right]$, respectively. Therefore, we obtain the following factorizations from [34]:

(19)$$\begin{gathered} {\mathbf{A}}_2^{(a)}+{\mathbf{B}}_2^{(a)}={\stackrel{\mathrm{-}}{\mathbf{I}}}_2{\mathbf{H}}_2\times \mathrm{diag}((g+h)/2, (g- h)/2){\mathbf{H}}_2, \\ {\mathbf{A}}_2^{(a)}-{\mathbf{B}}_2^{(a)} ={\mathrm{T}}_{2\times 3}^{(4)}\times \mathrm{diag}(p- q, -p- q, q){\mathrm{T}}_{3\times 2}^{(3)}, \end{gathered}$$

where ${\mathrm{T}}_{2\times 3}^{(4)}=\left[\begin{array}{ccc}1& 0& 1\\ 0& 1& 1\end{array}\right]$, ${\mathrm{T}}_{3\times 2}^{(3)}=\left[\begin{array}{cc}1& 0\\ 0& 1\\ 1& 1\end{array}\right],$ g$=b_5+e_5,$ h$=c_5-d_5,$ p$=b_5-e_5,$ q$=c_5+d_5.$

Let us add matrices from decomposition (19) to the left of the diagonal matrix into the matrix ${ \mathbf{W}}_{6\times 7}=({\stackrel{\mathrm{-}}{\mathbf{I}}}_2{\mathbf{H}}_2)\oplus {\mathbf{I}}_2 \oplus {\mathrm{T}}_{2\times 3}^{(4)}$, and to the left of the diagonal matrix into the matrix ${ \mathbf{W}}_{7\times 6}={\mathbf{H}}_2 \oplus {\mathbf{I}}_2 \oplus {\mathrm{T}}_{3\times 2}^{(3)}$. Combining the expressions (17), (18), and (19) the matrix factorization for the five-point odd-time DHT with normalization is represented as

(20)${\mathbf{Y}}_{5\times 1}={\mathbf{P}}_5^{(1)}{\mathbf{W}}_{5\times 6}{\mathbf{W}}_6^{(1)}{\mathbf{W}}_{6\times 7}{\mathbf{D}}_7{\mathbf{W}}_{7\times 6}{\mathbf{W}}_6^{(0)}{\mathbf{W}}_{6\times 5}{\mathbf{P}}_5^{(0)}{\mathbf{X}}_{5\times 1}$

where

$$\begin{gathered} {\mathbf{D}}_7=\mathrm{diag}\left(s_0^{(5)}, s_1^{(5)}, s_2^{(5)}, s_3^{(5)},s_4^{(5)}, s_5^{(5)}, s_6^{(5)}\right), s_0^{(5)}=(b_5+e_5+c_5-d_5)/4; s_1^{(5)}=(b_5+e_5-c_5+d_5)/4; \\ {s}_2^{(5)}=a_5; s_3^{(5)}=a_5; s_4^{(5)}=(b_5-e_5-c_5-d_5)/2; s_5^{(5)}=(-b_5+e_5-c_5-d_5)/2; s_6^{(5)}=(c_5+d_5)/2; \\ {\mathbf{W}}_6^{(1)}=\left[\begin{array}{cccccc}1& & & & 1& \\ & 1& & & & 1\\ & & 1& & & \\ & & & 1& & \\ -1& & & & 1& \\ & -1& & & & 1\end{array}\right], {\mathbf{W}}_6^{(0)}=\left[\begin{array}{cccccc}1& & & & 1& \\ & 1& & & & 1\\ & & 1& & & \\ & & & 1& & \\ 1& & & & -1& \\ & 1& & & & -1\end{array}\right], \\ {\mathbf{W}}_{6\times 5}=\left[\begin{array}{ccccc}1& & & & \\ & 1& & & \\ & & 1& & \\ 1& 1& & 1& 1\\ & & & 1& \\ & & & & 1\end{array}\right], {\mathbf{W}}_{5\times 6}=\left[\begin{array}{cccccc}1& & -1& & & \\ & 1& 1& & & \\ & & 1& 1& & \\ & & 1& & 1& \\ & & -1& & & 1\end{array}\right]. \end{gathered}$$

We have added the matrix ${\mathbf{I}}_2$ into the crossing of third and fourth columns and rows of the matrices ${ \mathbf{W}}_{7\times 6},{ \mathbf{W}}_{6\times 7}$, ${\mathbf{W}}_6^{(1)}$ and ${\mathbf{W}}_6^{(0)}$. This allows the factorizing the matrix ${\mathbf{C}}_5^{(b)}$ using ${\mathbf{W}}_{6\times 5}$ and ${\mathbf{W}}_{5\times 6}$.

Next, we present the algorithm of the five-point odd-time DHT with normalization via the data flow graph in Figure 6. As one can see, the number of multiplications is reduced from 25 to 6 because $s_1^{(5)}=0.25$; the number of additions is enlarged from 20 to 23.

Based on Equation (20), the matrix factorization for the five-point odd-frequency DHT with normalization is represented as

(21)${\mathbf{Y}}_{5\times 1}={\left({\mathbf{P}}_5^{(0)}\right)}^{\mathrm{T}}{\mathbf{W}}_{6\times 5}^{\mathrm{T}}{\left({\mathbf{W}}_6^{(0)}\right)}^{\mathrm{T}}{\mathbf{W}}_{7\times 6}^{\mathrm{T}}{\mathbf{D}}_7{\mathbf{W}}_{6\times 7}^{\mathrm{T}}{\left({\mathbf{W}}_6^{(1)}\right)}^{\mathrm{T}}{\mathbf{W}}_{5\times 6}^{\mathrm{T}}{\left({\mathbf{P}}_5^{(1)}\right)}^{\mathrm{T}}{\mathbf{X}}_{5\times 1}$

Figure 7 shows a data flow graph of the algorithm for the five-point odd-frequency DHT with normalization. The number of arithmetic operations is the same as for 5-point odd-time DHT with normalization.

3.5. Algorithms for 6-Point Odd-Time DHT and Odd-Frequency DHT

Let us propose an algorithm for six-point odd-time DHT with normalization which is expressed as follows:

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

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}{a}_6& a_6& a_6& a_6& a_6& a_6\\ b_6& a_6& -c_6& -b_6& -a_6& c_6\\ b_6& -a_6& -c_6& b_6& -a_6& -c_6\\ a_6& -a_6& a_6& -a_6& a_6& {-a}_6\\ c_6& a_6& -b_6& c_6& a_6& -b_6\\ -c_6& a_6& b_6& c_6& -a_6& -b_6\end{array}\right]$

with $a_6={\displaystyle \frac{1}{\sqrt{6}}}\approx 0.4082,$ $b_6$=${\displaystyle \frac{1}{\sqrt{6}}}\mathrm{cas}{\displaystyle \frac{\pi }{6}}\approx$ 0.5577, $c_6$=${\displaystyle \frac{1}{\sqrt{6}}}\mathrm{cas}{\displaystyle \frac{2\pi }{3}}\approx$ 0.1494. For odd-frequency DHT with normalization in Equation (22) we use ${\mathbf{C}}_6^{\mathrm{T}}$ instead of ${\mathbf{C}}_6$.

Note that $c_6$ = $b_6-a_6$, thus we decompose ${\mathbf{C}}_6=a_6{\mathbf{C}}_6^{(a)}+{\mathbf{C}}_6^{(b)}$, where

${\mathbf{C}}_6^{(a)}=\left[\begin{array}{cccccc}1& 1& 1& 1& 1& 1\\ 0& 1& 1& 0& -1& -1\\ 0& -1& 1& 0& -1& 1\\ 1& -1& 1& -1& 1& -1\\ -1& 1& 0& -1& 1& 0\\ 1& 1& 0& -1& -1& 0\end{array}\right], {\mathbf{C}}_6^{(b)}=\left[\begin{array}{}\\ 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].$

Further the matrix ${\mathbf{C}}_6^{(a)}$ is added to factorization immediately as a submatrix of the matrix ${\mathbf{W}}_{8\times 10}$: ${ \mathbf{W}}_{8\times 10}={\mathbf{C}}_6^{(a)} \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]$. The matrix ${\mathbf{B}}_{2\times 4}$ is determined by the structure of the matrix ${\mathbf{C}}_6^{(b)}$ for which the second and sixth rows, the third and fifth rows are the same except on the sign. The matrices ${\mathbf{W}}_{10\times 6}$ and ${\mathbf{W}}_{6\times 8}$ will denote the order of the columns and the rows of the matrices ${\mathbf{C}}_6^{(a)}$ and ${\mathbf{B}}_{2\times 4}$ in resulted factorization. The identity submatrix ${\mathbf{I}}_6$ of ${\mathbf{W}}_{10\times 6}$ and ${\mathbf{W}}_{6\times 8}$ related to ${\mathbf{C}}_6^{(a)}$, and the rest submatrices of ${\mathbf{W}}_{10\times 6}$ and ${\mathbf{W}}_{6\times 8}$ related to ${\mathbf{B}}_{2\times 4}$:

${\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]$

Let us introduce $s_0^{\left(6\right) }=s_1^{\left(6\right)}=s_2^{\left(6\right)}=s_3^{\left(6\right)}=s_4^{\left(6\right) }=s_5^{\left(6\right)}=a_6$ to multiply the each row of matrix ${\mathbf{C}}_6^{(a)}$, and $s_6^{\left(6\right)}=s_7^{\left(6\right)}=b_6$ to multiply the each row of ${\mathbf{B}}_{2\times 4}$. Then the factorization of the matrix of the 6-point odd-time DHT with normalization can be expressed as [37]

(23)${\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},$

where ${\mathbf{D}}_8=\mathrm{diag}(s_0^{\left(6\right)}, s_1^{\left(6\right)}, s_2^{\left(6\right)}, s_3^{\left(6\right)}, s_4^{\left(6\right)}, s_5^{\left(6\right)}, s_6^{\left(6\right)}, s_7^{\left(6\right)})$.

Figure 8 shows the data flow graph of the algorithm of the 6-point odd-time DHT with normalization. The amount of multiplications and additions can be lessened from 36 to 8 and from 30 to 28, respectively.

Based on Equation (23), the matrix factorization for the 6-point odd-frequency DHT with normalization is expressed as

(24)${\mathbf{Y}}_{6\times 1}={\mathbf{W}}_{10\times 6}^{\mathrm{T}}{\mathbf{W}}_{8\times 10}^{\mathrm{T}}{\mathbf{D}}_8{\mathbf{W}}_{6\times 8}^{\mathrm{T}}{\mathbf{X}}_{6\times 1},$

Figure 9 shows the data flow graph of the algorithm implementing the 6-point odd-frequency DHT with normalization. As one can see, the amount of multiplications may be lessened from 36 to 8. Also we decreased the amount of additions from 30 to 24.

3.6. Algorithms for 7-Point Odd-Time DHT and Odd-Frequency DHT

Now we will design the algorithm for 7-point odd-time DHT with normalization using the expression

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

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}}$, $a_7$ = ${\displaystyle \frac{1}{\sqrt{7}}}\approx$ 0.3780, $b_7$ =${\displaystyle \frac{1}{\sqrt{7}}}\mathrm{cas}{\displaystyle \frac{\pi }{7}}\approx$ 0.5045, $c_7$ = ${\displaystyle \frac{1}{\sqrt{7}}}\mathrm{cas}{\displaystyle \frac{3\pi }{7}}\approx$ 0.4526, $d_7$ = ${\displaystyle \frac{1}{\sqrt{7}}}\mathrm{cas}{\displaystyle \frac{5\pi }{7}}\approx$ 0.0598, $e_7$= ${\displaystyle \frac{1}{\sqrt{7}}}\mathrm{cas}{\displaystyle \frac{2\pi }{7}}\approx$ 0.5312, $g_7$= $-{\displaystyle \frac{1}{\sqrt{7}}}\mathrm{cas}{\displaystyle \frac{6\pi }{7}}\approx$ 0.1765, $f_7$=${\displaystyle \frac{1}{\sqrt{7}}}\mathrm{cas}{\displaystyle \frac{4\pi }{7}}\approx$ 0.2844,

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

For odd-frequency DHT with normalization in Equation (25) the ${\mathbf{C}}_7^{\mathrm{T}}$ matrix must be used instead of ${\mathbf{C}}_7$.

Let us define the permutations ${\pi }_7$ and ${\pi }_8$:

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

After permutation of the rows of ${\mathbf{C}}_7$ with ${\pi }_7$ and permutation the columns of ${\mathbf{C}}_7$ with ${\pi }_8$, we alter the sign of the fifth, sixth, and seventh columns. The corresponding permutation matrices are introduced as

${\mathbf{P}}_7^{(1)}=\left[\begin{array}{ccccccc}& & & 1& & & \\ 1& & & & & & \\ & 1& & & & & \\ & & 1& & & & \\ & & & & & & 1\\ & & & & & 1& \\ & & & & 1& & \end{array}\right], {\mathbf{P}}_7^{(0)}=\left[\begin{array}{cccccc}1& & & & & \\ & 1& & & & \\ & & 1& & & \\ & & & 1& & \\ & & & & & -1\\ & & & & -1& \\ & & & -1& & \end{array}\right].$

The obtained matrix ${\mathbf{C}}_7^{(a)}$ was decomposed into two components:

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

where

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

We can reduce the number of operations by taking into account the similar entries in the fourth column and fourth row of the matrix ${\mathbf{C}}_7^{(b)}$ [36], which differ only in sign. Next, in matrix ${\mathbf{C}}_7^{(c)}$ we eliminate the rows and columns with only zero entries. The obtained matrix

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

acquires the structure $\left[\begin{array}{cc}{\mathbf{A}}_3& {\mathbf{B}}_3\\ {\mathbf{B}}_3& {\mathbf{A}}_3\end{array}\right]$ and can be decomposed as [34]:

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

Now we will decompose the matrices

${\mathbf{A}}_3+{\mathbf{B}}_3 =\left[\begin{array}{ccc}{b}_7{-g}_7& c_7+f_7& e_7+d_7\\ e_7+d_7& b_7{-g}_7& -(c_7+f_7)\\ c_7+f_7& -(e_7+d_7)& b_7{-g}_7\end{array}\right], {\mathbf{A}}_3-{\mathbf{B}}_3=\left[\begin{array}{ccc}{b}_7{+g}_7& c_7- f_7& d_7-e_7\\ e_7-d_7& -(b_7 {+g}_7)& -(c_7-f_7)\\ c_7-f_7& d_7-e_7& b_7{+g}_7\end{array}\right],$

representing them as 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]$[35,36]:

(28)${\mathbf{H}}_3={\mathrm{T}}_3^{(1)}{\mathrm{T}}_{3\times 4}\times \mathrm{diag}(s_0,s_1,s_2,s_3){ \times \mathrm{T}}_{4\times 3}{\mathrm{T}}_3^{(0)},$

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$; ${\mathrm{T}}_3^{(1)}=\left[\begin{array}{ccc}1& 1& \\ 1& -1& -1\\ 1& & 1\end{array}\right]$,

${\mathrm{T}}_{3\times 4}=\left[\begin{array}{cccc}1& & & \\ & 1& & -1\\ & & 1& -1\end{array}\right], {\mathrm{T}}_{4\times 3}=\left[\begin{array}{ccc}1& & \\ & 1& \\ & & 1\\ & 1& 1\end{array}\right], {\mathrm{T}}_3^{(0)}=\left[\begin{array}{ccc}1& 1& 1\\ 1& & -1\\ & 1& -1\end{array}\right].$

We need to modify the matrices $({\mathbf{A}}_3+{\mathbf{B}}_3)$, $({\mathbf{A}}_3-{\mathbf{B}}_3)$ to make them consistent with the expression of the circular convolution for N = 3. For that, in the matrix $({\mathbf{A}}_3+{\mathbf{B}}_3)$ the sign of entries in the first column and first row was altered by multiplying left and right on ${\mathbf{P}}_3^{(1)}=\left[\begin{array}{ccc}-1& & \\ & 1& \\ & & 1\end{array}\right]$. In the $({\mathbf{A}}_3-{\mathbf{B}}_3)$ matrix the sign of the second-row entries was changed by multiplying left on ${\mathbf{P}}_3^{(0)}=\left[\begin{array}{ccc}1& & \\ & -1& \\ & & 1\end{array}\right]$. As a result, we obtain the matrices:

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

Let us apply the expression (28) to $({\mathbf{A}}_3+{\mathbf{B}}_3{)}^{(0)}$, $({\mathbf{A}}_3 -{\mathbf{B}}_3{)}^{(0)}.$ For the matrix $({\mathbf{A}}_3+{\mathbf{B}}_3{)}^{(0)}$ the elements of circular convolution are $h_0=b_7{-g}_7, h_1=-e_7-d_7, h_2=-c_7-f_7$. For the matrix $({\mathbf{A}}_3-{\mathbf{B}}_3{)}^{(0)}$ we have $h_0=b_7{+g}_7, h_1=-e_7+d_7, h_2=c_7-f_7$. Further we form the matrices ${ \mathbf{W}}_{10\times 8}={\mathrm{T}}_{4\times 3}\oplus {\mathbf{I}}_2\oplus {\mathrm{T}}_{4\times 3}$,${ \mathbf{W}}_{8\times 10}={\mathrm{T}}_{3\times 4} \oplus {\mathbf{I}}_2 \oplus {\mathrm{T}}_{3\times 4}$ to factorize the circular convolutions. The matrices ${\mathrm{T}}_3^{(1)}$ and ${\mathrm{T}}_3^{(0)}$ were multiplied on ${\mathbf{P}}_3^{(1)}$ and ${\mathbf{P}}_3^{(0)}$. And we obtained ${\mathbf{W}}_a={\mathrm{T}}_3^{(0)}{\mathbf{P}}_3^{(1)}=\left[\begin{array}{ccc}-1& 1& 1\\ -1& & -1\\ & 1& -1\end{array}\right]$, ${\mathbf{W}}_c={\mathbf{P}}_3^{(1)}{\mathrm{T}}_3^{(1)}=\left[\begin{array}{ccc}-1& -1& \\ 1& -1& -1\\ 1& & 1\end{array}\right]$, ${\mathbf{W}}_d{=\mathbf{P}}_3^{(0)}{\mathrm{T}}_3^{(1)}=\left[\begin{array}{ccc}1& 1& \\ -1& 1& 1\\ 1& & 1\end{array}\right]$

We have added the matrix ${\mathbf{I}}_2$ into direct sum in the matrices ${ \mathbf{W}}_{10\times 8},{ \mathbf{W}}_{8\times 10}$, ${\mathbf{W}}_8^{(0)}$ and ${\mathbf{W}}_8^{(1)}$. This allows the factorizing the matrix ${\mathbf{C}}_7^{(b)}$ using ${\mathbf{W}}_{7\times 8}$ and ${\mathbf{W}}_{8\times 7}$.

In sum, we derive the factorization of the matrix of the seven-point odd-time DHT with normalization:

(29)${\mathbf{Y}}_{7\times 1}={\mathbf{P}}_7^{(1)}{\mathbf{W}}_{7\times 8}{\mathbf{W}}_8{\mathbf{W}}_8^{(1)}{\mathbf{W}}_{8\times 10}{\mathbf{D}}_{10}{\mathbf{W}}_{10\times 8}{\mathbf{W}}_8^{(0)}{\mathbf{W}}_8{\mathbf{W}}_{8\times 7}{\mathbf{P}}_7^{(0)}{\mathbf{X}}_{7\times 1},$

where ${\mathbf{D}}_{10}=\mathrm{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^{(7)}, s_8^{\left(7\right)}, s_9^{(7)})$,

$$\begin{gathered} {\mathbf{W}}_8^{(0)}={\mathbf{W}}_a\oplus {\mathbf{I}}_2\oplus {\mathrm{T}}_3^{(0)}, {\mathbf{W}}_8^{(1)}={\mathbf{W}}_c\oplus {\mathbf{I}}_2\oplus {\mathbf{W}}_d, \\ {s}_0^{(7)}=(b_7{-g}_7-e_7-d_7-c_7-f_7)/6; s_1^{(7)}=(b_7{-g}_7{+c}_7+f_7)/2; s_2^{(7)}=(-e_7-d_7+c_7+f_7)/2; \\ {s}_3^{(7)}=(b_7{-g}_7-d_7-e_7+2c_7+{2f}_7)/6; s_4^{(7)}=a_7; s_5^{(7)}=a_7;s_6^{(7)}=(b_7{+g}_7-e_7+d_7+c_7-f_7)/6; \\ {s}_7^{(7)}=(b_7{+g}_7{-c}_7+f_7)/2; s_8^{(7)}=(-e_7+d_7-c_7+f_7)/2; s_9^{(7)}=(b_7{+g}_7-e_7+d_7-2c_7+{2f}_7)/6; \\ {\mathbf{W}}_{7\times 8}=\left[\begin{array}{cccccccc}1& & & -1& & & & \\ & 1& & 1& & & & \\ & & 1& -1& & & & \\ & & & 1& 1& & & \\ & & & 1& & 1& & \\ & & & -1& & & 1& \\ & & & 1& & & & 1\end{array}\right]& , & {\mathbf{W}}_{8\times 7}=\left[\begin{array}{ccccccc}1& & & & & & \\ & 1& & & & & \\ & & 1& & & & \\ & & & 1& & & \\ 1& 1& 1& & -1& -1& -1\\ & & & & 1& & \\ & & & & & 1& \\ & & & & & & 1\end{array}\right], \\ {\mathbf{W}}_8=\left[\begin{array}{cccccccc}1& & & & & 1& & \\ & 1& & & & & 1& \\ & & 1& & & & & 1\\ & & & 1& & & & \\ & & & & 1& & & \\ 1& & & & & -1& & \\ & 1& & & & & -1& \\ & & 1& & & & & -1\end{array}\right]. \end{gathered}$$

To calculate the 7-point odd-time DHT with normalization the fast algorithm is developed and presented in Figure 10 by data flow graph. As a consequence we can reduce the amount of multiplications from 49 to 10, but the amount of additions increased from 42 to 46.

Based on Equation (29), the matrix factorization for the 7-point odd-frequency DHT with normalization is obtained as

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

To compute the 7-point odd-frequency DHT with normalization we propose the fast algorithm shown in Figure 11 by data flow graph. The number of arithmetic operations is the same as for 7-point odd-time DHT with normalization.

3.7. Algorithms for 8-Point Odd-Time DHT and Odd-Frequency DHT

Let us develop the algorithm for eight-point odd-time DHT with normalization which is expressed as follows:

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

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}{a}_8& a_8& a_8& a_8& a_8& a_8& a_8& a_8\\ b_8& b_8& c_8& -c_8& -b_8& -b_8& -c_8& c_8\\ d_8& 0& -d_8& 0& d_8& 0& -d_8& 0\\ b_8& -b_8& c_8& c_8& -b_8& b_8& -c_8& -c_8\\ a_8& -a_8& a_8& -a_8& a_8& -a_8& a_8& -a_8\\ c_8& c_8& -b_8& b_8& -c_8& -c_8& b_8& -b_8\\ 0& d_8& 0& -d_8& 0& d_8& 0& -d_8\\ -c_8& c_8& b_8& b_8& c_8& -c_8& -b_8& -b_8\end{array}\right]$

with $a_8$ = ${\displaystyle \frac{1}{\sqrt{8}}}\approx$ 0.3536; $b_8$= ${\displaystyle \frac{1}{\sqrt{8}}}\mathrm{cas}{\displaystyle \frac{\pi }{8}}\approx$ 0.4619; $c_8$ = ${\displaystyle \frac{1}{\sqrt{8}}}\mathrm{cas}{\displaystyle \frac{5\pi }{8}}\approx$ 0.1913; $d_8$= ${\displaystyle \frac{1}{\sqrt{8}}}\mathrm{cas}{\displaystyle \frac{\pi }{4}}=$ 0.5. For odd-frequency DHT with normalization in (31) the ${\mathbf{C}}_8^{\mathrm{T}}$ matrix must be used instead of ${\mathbf{C}}_8$.

Let us define the permutation ${\pi }_8=\left(\begin{array}{cc}1& 2\\ 1& 5\end{array}\begin{array}{cccccc}3& 4& 5& 6& 7& 8\\ 3& 7& 4& 2& 6& 8\end{array}\right)$ to change the order of rows of ${\mathbf{C}}_8$. The permutation matrix ${\mathbf{P}}_8$ and the resulted matrix ${\mathbf{C}}_8^{(a)}$ can be expressed as

${\mathbf{P}}_8=\left[\begin{array}{cccccccc}1& & & & & & & \\ & & & & & 1& & \\ & & 1& & & & & \\ & & & & 1& & & \\ & 1& & & & & & \\ & & & & & & 1& \\ & & & 1& & & & \\ & & & & & & & 1\end{array}\right], {\mathbf{C}}_8^{(a)}=\left[\begin{array}{cccccccc}{a}_8& a_8& a_8& a_8& a_8& a_8& a_8& a_8\\ a_8& -a_8& a_8& -a_8& a_8& -a_8& a_8& -a_8\\ d_8& 0& -d_8& 0& d_8& 0& -d_8& 0\\ 0& d_8& 0& -d_8& 0& d_8& 0& -d_8\\ b_8& -b_8& c_8& c_8& -b_8& b_8& -c_8& -c_8\\ 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\\ -c_8& c_8& b_8& b_8& c_8& -c_8& -b_8& -b_8\end{array}\right].$

The resulted matrix ${\mathbf{C}}_8^{(a)}$ acquires the structure:

$\left[\begin{array}{cc}{\mathbf{A}}_4& {\mathbf{A}}_4\\ {\mathbf{B}}_4& -{\mathbf{B}}_4\end{array}\right]\mathrm{with} {\mathbf{A}}_4=\left[\begin{array}{cccc}{a}_8& a_8& a_8& a_8\\ a_8& -a_8& a_8& -a_8\\ d_8& 0& -d_8& 0\\ 0& d_8& 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 using the template from [34] we obtain ${\mathbf{C}}_8^{(a)}=({\mathbf{A}}_4 \oplus {\mathbf{B}}_4)({\mathbf{H}}_2 \otimes {\mathbf{I}}_4)$. The matrix ${\mathbf{A}}_4$ has the following structure:

$\left[\begin{array}{cc}{\mathbf{A}}_2^{(b)}& {\mathbf{A}}_2^{(b)}\\ {\mathbf{B}}_2^{(b)}& -{\mathbf{B}}_2^{(b)}\end{array}\right]\mathrm{with} {\mathbf{A}}_2^{(b)}=\mathrm{diag}(a_8, a_8){\mathbf{H}}_2\mathrm{and} {\mathbf{B}}_2^{(b)}=d_8{\mathbf{I}}_2.$

Using the decomposition of the extracted submatrices [34], the matrix ${\mathbf{A}}_4$ is represented as

(32)${\mathbf{A}}_4=\mathrm{diag}(a_8, a_8, d_8, d_8)({\mathbf{H}}_2 \oplus {\mathbf{I}}_2)({\mathbf{H}}_2 \otimes {\mathbf{I}}_2).$

The entries of the diagonal matrix we denoted as $s_0^{(8)}=a_8,$ $s_1^{(8)}=a_8,$ $s_2^{(8)}=d_8,$ $s_3^{(8)}=d_8.$

Next we swap the first and second column of the ${\mathbf{B}}_4$ matrix with permutation matrix ${\mathbf{P}}_2 \oplus {\mathbf{I}}_2,{\mathbf{P}}_2=\left[\begin{array}{cc}& 1\\ 1& \end{array}\right]$. The resulted matrix ${\mathbf{B}}_4^{(a)}=\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]$ is similar to the template $\left[\begin{array}{cc}{\mathbf{A}}_2^{(c)}& {\mathbf{B}}_2^{(c)}\\ {\mathbf{B}}_2^{(c)}& {\mathbf{A}}_2^{(c)}\end{array}\right]$ and can be decomposed as

(33)${\mathbf{B}}_4^{(a)}=({\mathbf{H}}_2 \otimes {\mathbf{I}}_2)1/2[({\mathbf{A}}_2^{(c)}+{\mathbf{B}}_2^{(c)}) \oplus ({\mathbf{A}}_2^{(c)}-{\mathbf{B}}_2^{(c)})]({\mathbf{H}}_2 \otimes {\mathbf{I}}_2),$

where $1/2({\mathbf{A}}_2^{(c)}+{\mathbf{B}}_2^{(c)})=1/2\left[\begin{array}{cc}-b_8+c_8& b_8+c_8\\ b_8+c_8& b_8-c_8\end{array}\right]$ matches the pattern $1/2\left[\begin{array}{cc}e& f\\ f& -e\end{array}\right]$, $e=-b_8+c_8$, $f=b_8+c_8$. It is represented as

(34)$1/2({\mathbf{A}}_2^{(c)}+{\mathbf{B}}_2^{(c)}) ={\mathrm{T}}_{2\times 3}^{(4)}\times \mathrm{diag}(e-f, -e- f, f){\mathrm{T}}_{3\times 2}^{(3)},$

where ${\mathrm{T}}_{2\times 3}^{(4)}$, ${\mathrm{T}}_{3\times 2}^{(3)}$ are defined as in Equation (19). The entries of diagonal matrix we denoted as $s_4^{(8)}=-b_8, s_5^{(8)}=-c_8, s_6^{(8)}=(b_8+c_8)/2.$ The matrix $1/2({\mathbf{A}}_2^{(c)}-{\mathbf{B}}_2^{(c)})$ is similarly decomposed with diagonal matrix entries $s_7^{(8)}=-b_8, s_8^{(8)}=c_8, s_9^{(8)}=(b_8-c_8)/2.$

Let us add matrices from decomposition (34) to the left of the diagonal matrix into the matrix ${\mathbf{W}}_{8\times 10}={\mathbf{I}}_4 \oplus {\mathrm{T}}_{2\times 3}^{(4)} \oplus {\mathrm{T}}_{2\times 3}^{(4)}$, and to the right of the diagonal matrix into the matrix ${\mathbf{W}}_{10\times 8}={\mathbf{I}}_4 \oplus {\mathrm{T}}_{3\times 2}^{(3)} \oplus {\mathrm{T}}_{3\times 2}^{(3)}$. Also we added matrices from decompositions (32), (33) to the left of the diagonal matrices into the matrix ${\mathbf{W}}_8^{(a)}={\mathbf{I}}_4 \oplus ({\mathbf{H}}_2\otimes {\mathbf{I}}_2)$, and to the right of the diagonal matrices into ${\mathbf{W}}_8^{(b)}={\mathbf{H}}_2{ \oplus \mathbf{I}}_2 \oplus {\mathbf{P}}_2 \oplus {\mathbf{I}}_2,{\mathbf{W}}_8^{(c)}=({\mathbf{H}}_2 \otimes {\mathbf{I}}_2)\oplus {\mathbf{I}}_4,{\mathbf{W}}_8^{(d)}={\mathbf{H}}_2\otimes {\mathbf{I}}_4$.

In sum, the 8-point odd-time DHT with normalization is factorized as

(35)${\mathbf{Y}}_{8\times 1}={\mathbf{P}}_8{\mathbf{W}}_8^{(a)}{\mathbf{W}}_{8\times 10}{\mathbf{D}}_{10}{\mathbf{W}}_{10\times 8}{\mathbf{W}}_8^{(a)}{\mathbf{W}}_8^{(b)}{\mathbf{W}}_8^{(c)}{\mathbf{W}}_8^{(d)}{\mathbf{X}}_{8\times 1},$

where ${\mathbf{D}}_{10}=\mathrm{diag}(s_0^{(8)}, s_1^{(8)}, s_2^{(8)}, s_3^{(8)}{, s}_4^{(8)}, s_5^{(8)}, s_6^{(8)}{, s}_7^{(8)}, s_8^{(8)}, s_9^{(8)})$.

To compute the 8-point odd-time DHT with normalization the data flow graph of the fast algorithm is designed and presented in Figure 12. As a consequence we can reduce the amount of multiplications from 48 to 8, because $s_2^{(8)}=s_3^{(8)}=0.5$, and the amount of additions from 48 to 28.

Based on Equation (35), the matrix factorization for the 8-point odd-frequency DHT with normalization is given as

(36)${\mathbf{Y}}_{8\times 1}={\left({\mathbf{W}}_8^{(d)}\right)}^{\mathrm{T}}{\left({\mathbf{W}}_8^{(c)}\right)}^{\mathrm{T}}{\left({\mathbf{W}}_8^{(b)}\right)}^{\mathrm{T}}{\left({\mathbf{W}}_8^{(a)}\right)}^{\mathrm{T}}{\mathbf{W}}_{10\times 8}^{\mathrm{T}}{\mathbf{D}}_{10}{\mathbf{W}}_{8\times 10}^{\mathrm{T}}{\left({\mathbf{W}}_8^{(a)}\right)}^{\mathrm{T}}{\mathbf{P}}_8^{\mathrm{T}}{\mathbf{X}}_{8\times 1},$

Figure 13 shows the data flow graph of the algorithm implementing the 8-point odd-frequency DHT with normalization. The amount of arithmetic operations is the same as for 8-point odd-time DHT with normalization.

4. Results

We confirmed the correctness of the odd-time DHT and odd-frequency DHT matrix factorizations obtained in Section 3 by experiment. Using the MATLAB R2023b software it was evaluated the coincidence of the results of the multiplication of factorized matrices by Formulas (7), (10), (14), (20), (23), (29), and (35) with the original odd-time DHT matrices. Also, it was evaluated the similarity of the results of the computing of the matrices of the odd-frequency DHT with the expressions (11), (15), (21), (24), (30), (36) using the MATLAB R2023b environment. The last formulas define the factorizations of the odd-frequency DHT matrices. The obtained coincidence of the factorized matrices with the initial odd-time DHT and odd-frequency DHT matrices for each N indicates the correctness of the developed algorithms.

Next, we evaluated the number of arithmetic operations for the calculation of the odd-time DHT and odd-frequency DHT with normalization via matrix-vector product and for the proposed algorithms. In general, the multiplications on data flow graphs are corresponded to the graph vertices, denoted by circles with factors. The amount of additions was evaluated by the amount of entering of two graph edges in one vertex. One must notice that these operations for zero elements of odd-time DHT and odd-frequency DHT matrices were not counted. In the matrices ${\mathbf{C}}_4$ and ${\mathbf{C}}_8$ we have four and eight zeros, respectively. Also if an element of the odd-time DHT or odd-frequency DHT matrix can be expressed as a power of two then the multiplication by this element can be implemented by a shift and has not counted as a multiplication. The proposed algorithms use the same number of multiplications to both the odd-time DHT and the odd-frequency DHT with normalization. Applying these algorithms we can reduce the number of multiplications required for odd-time or odd-frequency DHTs by nearly 69%, and the amount of additions by nearly of 5% for odd-time DHT and 8% for odd-frequency DHT if the length of input data is in the range from 2 to 8.

Table 1 shows the obtained results. In parentheses, we indicate the difference in the amount of operations in percent. The plus means that the amount of operations has increased relative to the matrix-vector product, and the minus means it has lessened. In Table 2 and Table 3 the number of multiplication and additions is presented for the existing odd-time DHT and odd-frequency DHT algorithms.

5. Discussion of Computational Complexity

Let us compare the number of multiplication and additions for the existing algorithms and for the proposed algorithms of odd-time DHT and odd-frequency DHT with normalization. As for the odd-time DHT, if N = 4 then the proposed algorithm is no different from algorithms known from the literature by the multiplication number. However if the existing algorithms do not include the normalization they require the additional four multiplications. If N = 8 then the amount of multiplications may be lessened with the proposed solution versus Hamood’s algorithms [31,32] by 33%. Relative to the prime factor algorithm [17], the amount of multiplications for N = 8 is the same, and less than for the radix-2 algorithm [27] by 20%. However, the considered existing algorithms of odd-time DHT do not include normalization for N = 8 which requires 8 additional multiplications. The amount of additions for N = 8 has increased by 17% relative to Hamood’s algorithms [31,32] and is the same as for algorithms in [17,27].

As for the odd-frequency DHT, if N = 4 then the obtained algorithm is no different from algorithms known from the literature by the amount of multiplications. However the absence of normalization requires the additional four multiplications. Besides, for N = 4, the proposed algorithm reduces the number of additions significantly (by 40–45%). If N = 8 then the number of multiplications may be significantly lessened with the proposed solution relative to the radix-2 algorithm [28,29] but exceeds by 33% the number of multiplications for the algorithm of the prime factor (N = $2^m$) without normalization [17,41,42]. The normalization will require 8 additional multiplications.

Due to the evaluation that has been done in [34], we can estimate the reduction of computational complexity of the proposed algorithms by Big O notation. The list represents the number of arithmetic operations that the algorithm performs to obtain the odd-time or odd-frequency DHT of size N. We can use the estimation of template complexity for even N directly. For odd N we have decomposed the initial transform matrix as the sum of two matrices. Therefore we have to add two multiplications and 2(N − 1) additions to the Big O notation of the template complexity. The obtained estimations are summarised in Table 4. Analysing the resulting Big O estimates it should be noted the following. The number of multiplications is reduced two times, in general, although the number of additions is increased by 2(N − 1) if the proposed algorithms are compared to the direct matrix-vector product.

6. Conclusions

In the paper, the fast algorithms for odd-time DHT and odd-frequency DHT with normalization are proposed for input signals of short length. The obtained algorithms for short-length input signals can subsequently be used as readymade modules of algorithms for long-length sequences.

Including the normalization in the fast algorithms for odd-time DHT and odd-frequency DHT allows us to use them directly in applications that involve recovering the processed signals or images, specifically, in image compression, signal and image denoising, image encryption, etc. Normalization ensures the orthonormality of transform and, as a result, the inverse transform matrix is obtained by a transposing of the direct transform matrix. The last property enables the construction of fast algorithms for inverse transform based on the fast algorithms for direct transform easily.

The fast algorithms for odd-time DHT and odd-frequency DHT were developed in the article based on the structural properties of the matrices of the transform coefficients. Optimization of the parameters of the odd-time DHT and odd-frequency DHT is not performed. We have applied the partition of the initial matrix into the submatrices with previous rearranging of the matrix rows and/or columns. The partitioning of the original matrix and rearranging its structure aims to find the factorization which possibly reduces further the computational complexity. Based on the obtained factorizations, we have constructed the fast algorithms presented by the data flow graphs. The critical path of each designed data flow graph has only one multiplication. This is a significant advantage of the proposed solutions. If the critical path of the data flow graph has more than one multiplication, additional difficulties arise with the presentation and further processing of the result of each successive multiplication.

The comparison of the computational complexity of the direct matrix-vector products with the computational complexity of the proposed algorithms has shown a reducing in the number of multiplications by 69% for the input sequence lengths from 2 to 8. At the same time, the amount of additions is lowered by nearly 5% for odd-time DHT and by 8% for odd-frequency DHT for the same range of the input sequence lengths.

The comparison of the designed fast algorithms for odd-time DHT and odd-frequency DHT with the existing algorithms shows that bearing in mind the structure of transform matrices enabled the increase in the fast algorithms’ efficiency against the efficiency of radix-type algorithms. This result enables easier operation in real time because the amount of needed signal processor resources is significantly reduced. Reducing the number of multiplications significantly contributes to the resource- and time-saving of signal processing and enables easier operation in real-time.

Author Contributions

Conceptualization, A.C.; methodology, A.C. and M.P.; software, M.P.; validation, A.C., M.P. and J.S.; formal analysis, A.C., M.P. and J.S.; investigation, M.P. and J.S.; writing—original draft preparation, M.P. and A.C.; writing—review and editing, A.C., M.P. and J.S.; supervision, A.C. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

Data are contained within the article.

Conflicts of Interest

The authors declare no conflicts of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1. The data flow graph of the proposed algorithm for computation of two-point odd-time DHT and odd-frequency DHT.

Figure 2. The data flow graph of the algorithm for the 3-point odd-time DHT.

Figure 3. The data flow graph of the algorithm for the 3-point odd-frequency DHT.

Figure 4. The data flow graph of the algorithm for implementation of 4-point odd-time DHT.

Figure 5. The data flow graph of the algorithm for implementation of 4-point odd-frequency DHT.

Figure 6. The data flow graph of the algorithm for the 5-point odd-time DHT.

Figure 7. The data flow graph of the algorithm for the 5-point odd-frequency DHT.

Figure 8. The data flow graph of the algorithm for the computation of the 6-point odd-time DHT.

Figure 9. The data flow graph of the algorithm for the computation of the 6-point odd-frequency DHT.

Figure 10. The data flow graph implementing the 7-point odd-time DHT.

Figure 11. The data flow graph implementing the 7-point odd-frequency DHT.

Figure 12. The data flow graph for the computing of the 8-point odd-time DHT.

Figure 13. The data flow graph of the algorithm for the 8-point odd-time DHT.

Table 1

The computational complexity of the matrix-vector product and the proposed fast algorithms for odd-time and odd-frequency DHTs.

N	Matrix-Vector Product		Proposed Odd-Time DHT Algorithms		Proposed Odd-Frequency DHT Algorithms
	Adds.	Mults.	Adds.	Mults.	Adds.
2	2	4	2 (0%)	2 (−50%)	2 (0%)
3	6	9	7 (+16%)	3 (−67%)	6 (0%)
4	8	4	6 (−25%)	2 (−50%)	6 (−25%)
5	20	25	23 (+15%)	6 (−76%)	23 (+15%)
6	30	36	28 (−7%)	8 (−78%)	24 (−20%)
7	42	49	46 (+10%)	10 (−80%)	46 (+10%)
8	48	48	28 (−42%)	8 (−83%)	28 (−42%)

Table 2

The number of operations for the existing odd-time DHT algorithms.

Algorithm	Year of Publication, Reference	N = 4		N = 8
		Mults.	Adds.	Mults.	Adds.
Prime factor	2003, [17]	2 (0%)	10 (−40%)	8 (0%)	28 (0%)
Radix-2	2011, [27]	–	–	10 (−20%)	28 (0%)
Radix-2 with decimation-in-time	2016, [31]	2 (0%)	6 (0%)	12 (−33%)	24 (+17%)
Split-radix	2020, [32]	2 (0%)	6 (0%)	12 (−33%)	24 (+17%)
Proposed	–	2	6	8	28

Table 3

The number of operations of the existing odd-frequency DHT algorithms.

Algorithm	Year of Publication, Reference	N = 4		N = 8
		Mults.	Adds.	Mults.	Adds.
Radix-2	2011, [28,29]	8 (−75%)	11 (−45%)	16 (−50%)	37 (−24%)
Prime factor, Wang	1992, [41]	2 (0%)	10 (−40%)	8 (0%)	28 (0%)
Prime factor, Hu, Chang, Ersoy	1992, [42]	2 (0%)	10 (−40%)	8 (0%)	28 (0%)
Prime factor, Bi, Zeng	2003, [17]	2 (0%)	10 (−40%)	6 (+33%)	26 (+8%)
Proposed	–	2	6	8	28

Table 4

The computational complexity of the proposed algorithms and the direct matrix-vector product in the Big O notation.

N	Proposed Algorithms		Matrix-Vector Product
	Adds.	Mults.	Adds.	Mults.
2, 4, 8	$N^2/2$	$N^2/2$	$N(N- 1)$	$N^2$
3, 5, 7	$(N+2)(N-1)$	${(N-1)}^2/2+2$
6	$(N+2)(N- 1)$	$N+2$

References

1. Jones, K. The regularized fast Hartley transform. Low-Complexity Parallel Computation of the FHT in One and Multiple Dimensions; 2nd ed. Springer: New York, NY, USA, 2021.

2. Patra, N.; Nayak, S.S. Discrete Hartley transform and its applications—A review. J. Integr. Sci. Technol.; 2022; 10, pp. 173-179.

3. Moraga, C. Generalized discrete Hartley transforms. Proceedings of the 39th International Symposium on Multiple-Valued Logic; Naha, Japan, 21–23 May 2009.

4. Vimal Kumar, C.; Raghavendra, D.; Eshwarappa, M.N. Design and implementation of VLSI DHT algorithm for image compression. Int. J. Eng. Trends Technol.; 2015; 23, pp. 218-223. [DOI: https://dx.doi.org/10.14445/22315381/IJETT-V23P241]

5. de Oliveira Nascimento, F.A. Hartley transform signal compression and fast power quality measurements for smart grid application. IEEE Trans. Power Deliv.; 2023; 99, pp. 1-11. [DOI: https://dx.doi.org/10.1109/TPWRD.2023.3302188]

6. Anoh, K.; Tanriover, C.; Ribeiro, M.V.; Adebisi, B.; See, C.H. On the fast DHT precoding of OFDM signals over frequency-selective fading channels for wireless applications. Electronics; 2022; 11, 3099. [DOI: https://dx.doi.org/10.3390/electronics11193099]

7. Budiman, G.; Suksmono, A.B.; Danudirdjo, D.; Pawellang, S. QIM-based audio watermarking with combined techniques of SWT-DST-QR-CPT using SS-based synchronization. Proceedings of the 6th International Conference on Information and Communication Technology (ICoICT); Bandung, Indonesia, 3–4 May 2018.

8. Lipinski, P.; Puchala, D. Digital image watermarking using fast parametric transforms. Bull. Pol. Acad. Sci. Tech. Sci.; 2019; 67, pp. 463-477. [DOI: https://dx.doi.org/10.24425/bpasts.2019.129646]

9. Kasban, H. A spiral-based image watermarking scheme using Karhunen–Loeve and discrete Hartley transforms. Multidimens. Syst. Signal Process.; 2017; 28, pp. 573-595. [DOI: https://dx.doi.org/10.1007/s11045-015-0361-4]

10. Zhang, X.; Su, Q. A spatial domain-based colour image blind watermarking scheme integrating multilevel discrete Hartley transform. Int. J. Intell. Syst.; 2021; 36, pp. 4321-4345. [DOI: https://dx.doi.org/10.1002/int.22461]

11. Voronenko, Y.; Puschel, M. Algebraic Signal Processing Theory: Cooley-Tukey Type Algorithms for Real DFTs. IEEE Trans. Signal Process.; 2009; 57, pp. 205-222. [DOI: https://dx.doi.org/10.1109/TSP.2008.2006152]

12. Coutinho, V.A.; Cintra, R.J.; Bayer, F.M. Low-complexity three-dimensional discrete Hartley transform approximations for medical image compression. Comput. Biol. Med.; 2021; 139, 105018. [DOI: https://dx.doi.org/10.1016/j.compbiomed.2021.105018]

13. Bozinovic, N.; Konrad, J. Scan order and quantization for 3D-DCT coding. Proc. SPIE Vis. Commun. Image Process.; 2003; 5150, 1205.

14. Rao, K.R.; Yip, P. The Transform and Data Compression Handbook; 1st ed. CRC Press LLC: Boca Raton, FL, USA, 2001.

15. Rao, K.R.; Yip, P. Discrete Cosine Transform: Algorithms, Advantages, Applications; 2nd ed. Academic Press: Cambridge, UK, 2014.

16. Prots’ko, I. Algorithm of efficient computation of generalised discrete Hartley transform based on cyclic convolutions. IET Signal Process.; 2014; 8, pp. 301-308. [DOI: https://dx.doi.org/10.1049/iet-spr.2013.0204]

17. Bi, G.; Zeng, Y. Transforms and Fast Algorithms for Signal Analysis and Representations; 1st ed. Birkhäuser: Boston, FL, USA, 2004.

18. Kober, V. Fast recursive computation of sliding DHT with arbitrary step. Sensors; 2020; 20, 5556. [DOI: https://dx.doi.org/10.3390/s20195556]

19. Jiang, L.; Shu, H.; Wu, J.; Wang, L.; Senhadji, L. A novel split-radix fast algorithm for 2-D discrete Hartley transform. IEEE Trans. Circuits Syst. I Regul. Pap.; 2010; 57, pp. 911-924. [DOI: https://dx.doi.org/10.1109/TCSI.2009.2028639]

20. Wu, J.; Shu, H.; Senhadji, L.; Luo, L.M. Radix-3x3 algorithm for the 2-D discrete Hartley transform. IEEE Trans. Circuits Syst. Part 2 Analog Digit. Signal Process.; 2008; 55, pp. 566-570.

21. Grigoryan, A.M. Novel algorithm for computing the 1-D discrete Hartley transform. IEEE Signal Process. Lett.; 2004; 11, pp. 156-159. [DOI: https://dx.doi.org/10.1109/LSP.2003.819876]

22. Hamood, M. New efficient algorithm for the discrete Hartley transform. IOP Conf. Ser. Mater. Sci. Eng.; 2020; 978, 012013. [DOI: https://dx.doi.org/10.1088/1757-899X/978/1/012013]

23. Chiper, D.F. A novel VLSI DHT algorithm for a highly modular and parallel architecture. IEEE Trans. Circuits Syst. II Express Briefs; 2013; 60, pp. 282-286. [DOI: https://dx.doi.org/10.1109/TCSII.2013.2251974]

24. Pyrgas, L.; Kitsos, P.; Skodras, A.N. An FPGA design for the two-band fast discrete Hartley transform. Proceedings of the IEEE International Symposium on Signal Processing and Information Technology (ISSPIT); Limassol, Cyprus, 12–14 December 2016.

25. Yang, X.; Megson, G.M.; Xing, Y.; Evans, D.J. A novel radix-3/9 algorithm for type-III generalized discrete Hartley transform. J. Circuits Syst. Comput.; 2006; 15, pp. 301-312. [DOI: https://dx.doi.org/10.1142/S0218126606003064]

26. Shu, H.; Wang, Y.; Senhadji, L.; Luo, L. Direct computation of type-II discrete Hartley transform. IEEE Signal Process. Lett.; 2007; 14, pp. 329-332. [DOI: https://dx.doi.org/10.1109/LSP.2006.888362]

27. Chiper, D.F. Fast radix-2 algorithm for the discrete Hartley transform of type II. IEEE Signal Process. Lett.; 2011; 18, pp. 687-689. [DOI: https://dx.doi.org/10.1109/LSP.2011.2170166]

28. Hua, X.; Liu, J. A novel fast algorithm for discrete Hartley transform of type-III base on the split-radix and first-order moments. Proceedings of the SPIE Eighth International Symposium on Multispectral Image Processing and Pattern Recognition; Wuhan, China, 26–27 October 2013.

29. Chiper, D.F. Radix-2 fast algorithm for computing discrete Hartley transform of type III. IEEE Trans. Circuits Syst. II Express Briefs; 2012; 59, pp. 297-301. [DOI: https://dx.doi.org/10.1109/TCSII.2012.2190863]

30. Chiper, D.F. An improved algorithm for the VLSI implementation of type II generalized DHT that allows an efficient incorporation of obfuscation technique. Proceedings of the 45th International Conference on Telecommunications and Signal Processing (TSP); Prague, Czech Republic, 13–15 July 2022.

31. Hamood, M.T. Fast algorithm for computing the discrete Hartley transform of type-II. Indones. J. Electr. Eng. Inform.; 2016; 4, pp. 120-125. [DOI: https://dx.doi.org/10.11591/ijeei.v4i2.220]

32. Hamood, M.T. Direct split-radix algorithm for fast computation of type-II discrete Hartley transform. Telecommun. Comput. Electron. Control; 2020; 18, pp. 3067-3072. [DOI: https://dx.doi.org/10.12928/telkomnika.v18i6.16100]

33. Huang, Y.-M.; Wu, J.-L. Polynomial transform based algorithms for computing two-dimensional generalized DFT, generalized DHT, and skew circular convolution. Signal Process.; 2000; 80, pp. 2255-2260. [DOI: https://dx.doi.org/10.1016/S0165-1684(00)00110-9]

34. Cariow, A. Strategies for the synthesis of fast algorithms for the computation of the matrix-vector product. J. Signal Process. Theory Appl.; 2014; 3, pp. 1-19. [DOI: https://dx.doi.org/10.7726/jspta.2014.1001]

35. Cariow, A.; Papliński, J. Algorithmic structures for realizing short-length circular convolutions with reduced complexity. Electronics; 2021; 10, 2800. [DOI: https://dx.doi.org/10.3390/electronics10222800]

36. Polyakova, M.; Witenberg, A.; Cariow, A. The design of fast type-V discrete cosine transform algorithms for short-length input sequences. Electronics; 2024; 13, 4165. [DOI: https://dx.doi.org/10.3390/electronics13214165]

37. Polyakova, M.; Cariow, A. Fast Type-II Hartley Transform Algorithms for Short-Length Input Sequences. Appl. Sci.; 2024; 14, 10719. [DOI: https://dx.doi.org/10.3390/app142210719]

38. Perera, S.M. Signal flow graph approach to efficient and forward stable DST algorithms. Linear Algebra Its Appl.; 2018; 542, pp. 360-390. [DOI: https://dx.doi.org/10.1016/j.laa.2017.05.050]

39. Perera, S.M.; Lingsch, L. Sparse matrix-based low-complexity, recursive, radix-2 algorithms for discrete sine transforms. IEEE Access; 2021; 9, pp. 141181-141198. [DOI: https://dx.doi.org/10.1109/ACCESS.2021.3120051]

40. Paraskevas, I.; Barbarosou, M.; Chilton, E. Hartley transform and the use of the whitened Hartley spectrum as a tool for phase spectral processing. J. Eng.; 2015; 2015, pp. 95-101. [DOI: https://dx.doi.org/10.1049/joe.2014.0350]

41. Wang, Z. A prime factor fast W transform algorithm. IEEE Trans. Signal Process.; 1992; 40, pp. 2361-2368. [DOI: https://dx.doi.org/10.1109/78.157241] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/18425179]

42. Hu, N.C.; Chang, H.I.; Ersoy, O.K. Generalized discrete Hartley transforms. IEEE Trans. Signal Process.; 1992; 40, pp. 2931-2940.