Fast Recursive Computation of Sliding DHT with Arbitrary Step

Short-time (sliding) transform based on discrete Hartley transform (DHT) is often used to estimate the power spectrum of a quasi-stationary process such as speech, audio, radar, communication, and biomedical signals. Sliding transform calculates the transform coefficients of the signal in a fixed-size moving window. In order to speed up the spectral analysis of signals with slowly changing spectra, the window can slide along the signal with a step of more than one. A fast algorithm for computing the discrete Hartley transform in windows that are equidistant from each other is proposed. The algorithm is based on a second-order recursive relation between subsequent equidistant local transform spectra. The performance of the proposed algorithm with respect to computational complexity is compared with the performance of known fast Hartley transform and sliding algorithms.


Introduction
In today's world, the study of real-life phenomena begins with sensors that convert the description of phenomena into digital signals. Since real sensors are sensitive to internal thermal, electronic, and environmental noise, the sensors' output must be intelligently processed for subsequent use in a computer system. Short-time signal processing [1] is an appropriate technique for such a system, since it is capable of adaptively filtering long-term stationary (in wide-sense) and quasi-stationary signals. An example of quasi-stationary signal is speech. Frequency components of this signal change over time when the position of tongue and mouth changes, but on a short period of time they might be considered stable because tongue cannot move very quickly. Short-time processing is used in a variety of applications such as radar emitter recognition [2], heart sound classification [3], sensor noise removal [4], gearbox fault diagnosis [5], speech processing and spectral analysis [6], and adaptive linear filtering [7].
Basically, short-time transform is a time series of windowed signal transforms that can be used to perform time-frequency analysis. Local signal processing in a moving window (special case of short-time processing) can be carried out as follows: at each position of the running window, the coefficients of an orthogonal transform of the signal are modified to estimate only the central window element [8]. For time-varying signals, it is preferred that the window length be small enough so that the windowed signal would be approximately stationary. On the other hand, the window size should be large enough to provide a reasonable frequency resolution of the local signal. Typically, with local sliding processing, the window moves with a step of one along the signal. However, if it is necessary to process the signal at a high speed and the signal spectrum changes slowly over time, the window can move with a step of more than one. This approach has recently been proposed to design preview tools for multiresolution signal and image restoration [9].
Discrete Hartley transform (DHT) [10] is used in various practical applications such as speech spectral analysis [11], feature extraction and sea surface modeling [12], data compression [13], and signal interpolation [14]. Four types of DHT are classified [15,16]. Discrete Hartley and Fourier transforms of real signals and their relations with discrete sinusoidal transformations are discussed in detail [17,18]. Note that for real signals, all DHT coefficients can be used in the design of transform domain scalar filters. Representation of a signal in the DHT sliding domain is especially suitable for processing of time-varying, quasi-stationary data. Since calculating the DHT at each position of a running window is an intensive task, recursive algorithms can be used. First-order shift properties were obtained [19]. Since this approach is not very efficient in terms of computational complexity, four types of the sliding DHT were derived based on second-order recursive equations [20].
In this work, a recursive algorithm is proposed for fast computation of the DHT in a window sliding along the signal more than one sample at a time. The main contributions of the work are as follows: • Second-order recursive equation between three consecutive equidistant DHT spectra is obtained using the z-transform technique. • Efficient sliding DHT algorithm is proposed using the properties of discrete sinusoidal functions and the recursive equation.

•
Computational complexity and running time of the proposed sliding algorithm are compared with the known sliding and fast DHT algorithms.
The paper is organized as follows. In Section 2, a second-order recursive equation between three consecutive equidistant DHT spectra is derived. In Section 3, a fast forward algorithm for computing the sliding DHT is suggested. Section 3 also presents a fast inverse sliding DHT algorithm. The performance of the proposed and common DHT algorithms in terms of computational complexity and runtime is given and discussed in Section 4. Finally, our conclusions are presented in Section 5.

Second-Order Equation for Recursive Computation of Sliding DHT
We recall here the definition of DHT. The following notations are used: cs N (rs) ≡ cos 2π N rs , sn N (rs) ≡ sin 2π N rs , snc N (rs) ≡ sn N (rs) sn N (s) cas N (rs) ≡ cos 2π N rs + sin 2π N rs , cas N (rs) ≡ cos 2π N rs − sin 2π N rs , W s N ≡ e i 2π N s = cs N (s) + isn N (s), where the subscript N is the transform order, r is an integer, and s = 0, 1, . . . N − 1. For clarity, the normalization factor 1/ √ N is disregarded until the inverse transform.
Sliding DHT with a step of p can be defined as follows: where x(kp) ; k = . . . , −N 1 , −N 1 + 1, . . . , 0, 1, . . . , N 2 , N 2 + 1 . . . .} is a signal to be processed; y s (kp), s = 0, 1, . . . N − 1 are the transform coefficients around time kp; N = N 1 + N 2 + 1 is the size of the moving window; N 1 and N 2 are arbitrary integer values. A recursive relationship between three consecutive sliding spectra is given as [20]: where f s (k) = (x(k + N 2 This is a linear inhomogeneous difference equation defined on k. For fixed s, (2) can be considered as a linear difference equation with constant coefficients. Linear time-invariant systems defined by such equations can be analyzed using the unilateral z-transform [1]. Suppose that we deal with a causal linear system (region of convergence is outside a circle on the z-plane) and nonzero initial conditions y s (0) and y s (1) are given. Applying the unilateral z-transform to (2) and using its shift property, we get: where Y s (z) and F s (z) are the z-transforms of y s (k) and f s (k), respectively. We express Y s (z) as: For s 0, there are two different roots of the denominator; that is, q 1 (s) = W s N and q 2 (s) = W −s N . For s = 0, the roots are equal to q 1 (s) = q 2 (s) = 1. Let us carry out the partial fraction expansion of the equation: and obtain The inverse z-transform can be computed as follows [1]: where u(k) is the unit-step function defined as 1, for k ≥ 0 and 0, for k < 0. Using (5), (4) can be represented as follows: Using the shifting and convolution properties, the inverse transform is given as: Substituting (7) into (9) and taking into account that d s (k) = 0 for k < 2, q 1 (s)q 1 (s) = 1 and 2cs N (s) = q 1 (s) + q 2 (s), we arrive at: with k ≥ 2.
Given two initial values y s (0) and y s (p) (p > 1), one can obtain y s (1) from (10) as: Substituting (11) into (10), for arbitrary k we obtain: Finally, the DHT at the window position 2p is recursively computed using y s (0) and y s (p) as: This equation gives relationship between three consecutive equidistant DHT spectra, computed at times 0, p, and 2p.

Special Values of Discrete Sinusoidal Functions
Given below are special values (0, 1, − 1) of discrete sinusoidal functions used in the performance analysis of the proposed sliding algorithm. Let us consider the following functions: cs N (rs), snc N (rs), and cas N (s) with s = 0, 1, . . . , N − 1. The special values are shown in Table 1. Here, r, l, and N are arbitrary integer values, and the binary variable b takes the values {0, 1}. So, it is necessary to find such values of the variables l and b in order to obtain integer values of s in the range from 0 to N − 1. Table 1. Special values of the discrete sinusoidal functions.

Functions
Values The functions are used in the proposed algorithm in the operations of addition and multiplication; that is, when the functions are equal to 0, then the corresponding addition and multiplication operations can be discarded, and when they are equal to either −1 or 1, then the corresponding multiplication operations can be canceled. Therefore, the use of the special values can reduce the computational complexity of the algorithm.

Design of Forward Sliding Algorithm
Equation (13) for p > 1 can be rewritten as: where At each window, the number of additions required for calculation of B(r), A(r); r = 1, . . . p is 3p − 2, since the coefficients ∆(r); r = 0, . . . p − 1 have already been computed and stored at the position p of the window. Suppose that N is odd. Let g r , h + r , h − r be the greatest common factors of N and r, r +1, r −1, respectively. The property of symmetry of the discrete function snc N (rs) is given as follows: snc N (rs) = snc N (r(N − s)) ; s = 1, . . .  Table 1). So, the number of multiplications can be estimated as: (15) Let us denote SUM Q (s) ≡ Additional expenses are required for calculating initial p coefficients. Suppose that N is even. Equation (14) can be represented as follows: Let g r ,h + r , h − r be the greatest common factors of N 2 and r, r + 1, r −1, respectively. In this case, the property of symmetry of the discrete function snc N (rs) is given as follows: snc N (rs) = snc N (r(N − s)) = (−1) r−1 snc N r N 2 ± s , snc N r N 2 = (−1) r−1 r ; s = 0, . . .  Table 1). The number of zeros of cs N (ps) ; s = 1, . . . N − 1} is equal to 2g p (see the first row of Table 1). If N/8 is integer (see the third row of Table 1), then the number of zeros of the function cas N (s) ; s = 1, . . . N − 1} equals C cas = 2, otherwise C cas = 0. It means that the zeros of the function can cancel the operations of multiplication and addition in the first term of (17) (see the last line of the equation) at the corresponding frequencies s. The number of ±1 of the function cas N (s) ; s = 0, . . . N − 1} is equal to 2C cas . Finally, the number of multiplications can be estimated as follows: Using the symmetry property of the functions SUM A (s) = SUM A (N − s), SUM B (s) = SUM B (N − s); s = 1, . . . N 2 − 1 , the number of additions can be estimated as: Additional expenses are needed for calculating initial p coefficients. To more clearly show the design of the algorithm, we provide an example for computing the sliding DHT coefficients. Assume that N 1 = 7, N 2 = 8, N = 16, p = 2, the DHT is computed at the window position 2p y s (2p), s = 0, 1, ...15 . We borrow from the window position p two coefficients 0 = x 9 − x −7 ; 1 = x 10 − x −6 , and calculate the following auxiliary data: DHT coefficients are computed as follows: . The computational complexity of the algorithm is 25 multiplications and 61 additions.

Design of Inverse Sliding Algorithm
The inverse sliding DHT computes only the element x(kp) of the window. The inverse algorithm can be written as follows: where N = N 1 + N 2 + 1. If x(kp) is the central window element, that is, N = 2N 1 , then one can simplify the inverse transform to: Finally, for N 1 = 0 the inverse transform is given as: The proposed algorithms require only one multiplication and N − 1 additions.

Results and Discussion
In this section, using computer simulation, we analyze the performance of the proposed algorithm with respect to the computational complexity and execution time and compare it with that of fast DHT and conventional recursive algorithms. Among fast DHT algorithms, the most popular are fast radix-2 [21,22]. The recursive sliding (with a step of one) DHT algorithms [19,20] are carried out p times for computing the DHT spectra at equidistant positions kp of the window. Tables 2 and 3 show the computational complexity in terms of multiplications and additions, respectively, of the proposed, fast radix-2 DHT, and known sliding algorithms at a fixed window position for p = 2 when N varies.  It can be seen that the proposed algorithm for p > 2 outperforms the conventional sliding DHT algorithms. Note that the proposed algorithm is more efficient than the fast and conventional recursive algorithms when the window length increases.
In modern processors, the execution times of floating point multiplication and addition are comparable. So, further we will estimate the computational complexity with respect to flop counts (real additions and multiplications). Table 4 shows the computational complexity of the tested algorithms for N = 256 when p varies. Fast DHT [22] 2724 2724 2724 2724 2724 Ref. [19] 2556 3834 5112 6390 7668 Ref. [20] 1912 It can be seen that the proposed algorithm is more efficient than the fast DHT algorithm when the step p < 5. This boundary value of the step, when the proposed algorithm is still better than the fast DHT algorithm, increases with increasing the window length.
Obviously, the execution time of any algorithm depends on the characteristics of a computer used in a particular implementation of the algorithm. Now, with the help of computer simulation, we want to illustrate how the theoretical computational complexity of the tested algorithms relates to the execution time of the implemented algorithms. Computer simulation was performed on a laptop with an Intel Core i7-2630QM processor with 8 GB of RAM using the MATLAB R2012b. Experiments were carried out 100 times to ensure statistically correct results. Average time results for each algorithm were calculated. Figure 1 shows the performance of the tested algorithms in terms of runtime: DHT is the discrete Hartley transform given by definition, Fast DHT is implemented in Matlab, SDHT is sliding DHT [20], ALG is the proposed algorithm. It can be seen that the obtained results are in good agreement with the results in Table 4. Only the performance of the fast algorithm implemented in the MATLAB is slightly lower than that of the algorithm [22].
There are four types of DHT [16] that are suitable for the time varying processing of different signal models. In this paper, the second order recursive equation and fast recursive algorithm have been proposed for only one type of discrete Hartley transform (DHT-I). In the future, the same approach can be used to derive recursive equations and design fast recursive algorithms for other types of DHT, which will allow the system to work effectively with various signal models.

Conclusions
A second-order recursive equation between three consecutive equidistant DHT spectra was obtained utilizing of the unilateral z-transform technique. Using the properties of discrete sinusoidal functions and the recursive equation, a fast sliding DHT algorithm was proposed. A fast inverse sliding DHT transform was also presented. The computational complexity of the proposed sliding algorithm was compared with the known running and fast DHT algorithms. The proposed algorithm was implemented on a laptop, and it was shown that the theoretical computational complexity and execution time of the implemented algorithm are in good agreement.
Funding: This research received no external funding.

Conflicts of Interest:
The author declares no conflict of interest. It can be seen that the obtained results are in good agreement with the results in Table 4. Only the performance of the fast algorithm implemented in the MATLAB is slightly lower than that of the algorithm [22].
There are four types of DHT [16] that are suitable for the time varying processing of different signal models. In this paper, the second order recursive equation and fast recursive algorithm have been proposed for only one type of discrete Hartley transform (DHT-I). In the future, the same approach can be used to derive recursive equations and design fast recursive algorithms for other types of DHT, which will allow the system to work effectively with various signal models.

Conclusions
A second-order recursive equation between three consecutive equidistant DHT spectra was obtained utilizing of the unilateral z-transform technique. Using the properties of discrete sinusoidal functions and the recursive equation, a fast sliding DHT algorithm was proposed. A fast inverse sliding DHT transform was also presented. The computational complexity of the proposed sliding algorithm was compared with the known running and fast DHT algorithms. The proposed algorithm was implemented on a laptop, and it was shown that the theoretical computational complexity and execution time of the implemented algorithm are in good agreement.