Sparse Matrix Based Low-Complexity, Recursive, and Radix-2 Algorithms for Discrete Sine Transforms

This paper presents factorizations of each discrete sine transform (DST) matrix of types I, II, III, and IV into a product of sparse, diagonal, bidiagonal, and scaled orthogonal matrices. Based on the proposed matrix factorization formulas, reduced multiplication complexity, recursive, and radix-2 DST I-IV algorithms are presented. We will present the lowest multiplication complexity DST-IV algorithm in the literature. The paper fills a gap in the self-recursive, exact, and radix-2 DST I-IIII algorithms executed via diagonal, bidiagonal, scaled orthogonal, and simple matrix factors for any input $n=2^{t} \; (t \geq 1)$ . The paper establishes a novel relationship between DST-II and DST-IV matrices using diagonal and bidiagonal matrices. Similarly, a novel relationship between DST-I and DST-III matrices is proposed using sparse and diagonal matrices. These interweaving relationships among DST matrices enable us to bridge the existing factorizations of the DST matrices with the proposed factorization formulas. We present signal flow graphs to provide a layout for realizing the proposed algorithms in DST-based integrated circuit designs. Additionally, we describe an implementation of algorithms based on the proposed DST-II and DST-III factorizations within a double random phase encoding (DRPE) image encryption scheme.


I. INTRODUCTION
Discrete sine transforms are real-valued transform matrices and a subclass of the discrete Fourier transform (DFT) matrix. The DFT matrix can be described as a Vandermonde matrix having nodes as the primitive n th roots of unity, where n is a positive integer. The fast Fourier transform (FFT) is an algorithm that can be used to compute the DFT and its inverse efficiently.
For a given vector x = [x 0 , x 1 , · · · , x n−1 ] T ∈ R n , its DST can be expressed as y = Sx, where S is the DST matrix. The DST matrix has four main variants, i.e., I-IV as defined below , n (j) = 1 for j ∈ {1, 2, · · · , n−1} and n ≥ 2 is an even integer. Here, we use the superscript to denote the type of a given DST matrix and the subscript to denote the order of a DST matrix. Among the DST I-IV matrices, S I n−1 and S IV n were introduced in [1], [2], and S II n and its inverse S III n were introduced in [3] to digital signal processing. The complete set of even discrete cosine transform (DCT) and DST variants was originally presented in [4]. DST-II is a complementary or alternative transform to DCT-II, which is used in transform coding. Like the DFT and DCT, these DST matrices hold linearity, convolution-multiplication, orthogonality, and shift properties. Moreover, the DSTs of types II and III are related via S III n = S II n T . In this paper, we will establish relationships VOLUME 9, 2021 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ between S II n and S IV n as well as S I n−1 and S III n . The matrix factorization for DST-I in [5] used the results from [6] to decompose DST-I into the DCT and DST. Though one can find orthogonal matrix factorizations for the DCT and DST in [7], the resulting algorithms in [7] are not completely recursive and hence do not lead to simple recursive algorithms. An alternative factorization for DCT-I in [8] and DST-I in [9] can be seen in [10] and [11], but the factorizations in the latter papers do not solely depend on DCT I-IV or DST I-IV. Moreover, [11] has used the same factorization for DST-II and IV as in [7]. The authors in [9], [12], however, produce radix-2, completely recursive, and stable DST I-IV algorithms. Furthermore, [12] established a connection between algebraic operations used in the sparse and scaled orthogonal factorization of DST I-IV matrices and signal flow graphs. We refer to many more historical remarks on DST factorizations and the corresponding signal flow graphs in [12]. We note here that the DCT algorithms in [8], [9], and [13] and the relationship between cosine and sine transform matrices in [10] can be used to obtain the DST algorithms in [12]. On the other hand, one can use an elegant factorization for Chebyshev-like Vandermonde matrices to factor the DCT and DST matrices as described in [14]. However, there are no simplest sparse factors, no explicit algorithms, and no complexity analysis to compute DCT and DST matrices by a vector in [14]. In this paper, we obtain simple factors, i.e, diagonal, bidiagonal, sparse, and scaled orthogonal factors, to compute DST matrices by a vector. Additionally, we present self-recursive/completely-recursive and radix-2 DST algorithms with the lowest multiplication complexity in the literature. In recent years, developments in calculating the DST have widely broadened the algorithm's applications. Methods based upon the recursive nature of Chebyshev-polynomials have led to a drastic decrease in hardware complexity for applications that only require the transmission of a few key transform coefficients, creating a DST algorithm that is suitable for very large scale integrated circuits (VLSI) [15]- [17]. Similarly, a first-order moment-based approach has been developed and can be used to calculate the 1-dimensional DST-II without the need for multiplications by using special accumulator arrays and without having explicit factors for the DST-II to decrease computation time for certain applications [18]. For complex inputs, sparse factorization-based methods offer an extremely reduced flop count when compared to the brute-force method. The DST-II algorithms proposed in [19] take this a step further by developing a radix-2 based even-odd output vector at each stage, but without having explicit factorization for the DST-II matrix for any n.
The development of the convolution-multiplication property for the families of DSTs and DCTs led way for the DST to become an alternative to the DFT in many filtering applications [20]. Similarly, the duality theorem proposed in [21] and [22] now allows for the avoidance of calculating the inverse DST (IDST) in applications where similar signals are commonly encountered. Recently, the DST has expanded to applications including underwater signal transmission [23], [24], analysis and noise estimation of speech [25], finite impulse response (FIR) frequency filtering [26], and feature extraction of modulated signals in orthogonal frequency division multiplexing (OFDM) systems [27]. Similarly, the noise filtering capabilities have led to a new DST-based approach to measuring the volatility of stock prices [28], [29]. In hardware applications, the DST provides advantages over DCT-based approaches for statistical processes like the first-order Markov sequences with specific boundary conditions, and it achieves a lower bit error rate (BER) than the DCT for signals with low correlation [30]. It has also been shown that the DST may be implemented in the design of fractional order differentiators [31]. With the ability to develop DST algorithms with regularity, modularity, pipelining capability, and local connectivity, the DST will continue to expand its role in VLSI implementation [15]- [19], [26], [30], [31].
Although the DCT has traditionally played the pivotal role of frequency decomposition for data compression in methods like JPEG and MPEG, recent works have suggested that the DST may be applied to develop new methods of compression with increased performance for image and audio analysis as well as high-efficiency video encoding (HEVC) [22], [32]- [35]. The DST has proven to be an effective tool in the medical field as well, combining multiple images of different modes into a singular, highly detailed image for more accurate diagnoses [36], [37]. Expansion upon the fractional DST has also broadened the transform's role in the field of cyber security, particularly for applications in image encryption [38]- [40]. It was shown in [38] that the DST increases the entropy of encryption schemes that rely on processing within transform domains. This paper explores this concept further in Section VIII.
We have observed that the lowest multiplicative complexity, radix-2, and recursive DST algorithms having diagonal, bidiagonal, sparse, and scaled orthogonal matrices are missing elements in the literature. Thus, in this paper, we follow a similar approach in deriving DCT II/III algorithms in [41] and the DCT I/IV algorithms in [42] to obtain novel lowest multiplication complexity, recursive, and radix-2 DST algorithms execute via diagonal, bidiagonal, sparse, and scaled orthogonal matrices. We also obtain novel relationships between DST-II and DST-IV matrices, and DST-I and DST-III matrices in terms of diagonal, bidiagonal, and sparse matrices.
We first introduce frequently used notation in Section II. In Section III, we obtain simple, diagonal, bidiagonal, sparse, and scaled orthogonal factors for the DST matrices. In Section IV, we derive relationships between DST-I and DST-III matrices, and also between DST-II and DST-IV matrices, using bidiagonal, diagonal, and sparse matrices. Next, in Section V, we state recursive and radix-2 DST algorithms to compute the DST matrices by a vector based on the factorization of the DST matrices described in Section III. In Section VI, we establish arithmetic complexities for computing the proposed DST algorithms. Within the same section, we compare the arithmetic complexities in computing the proposed DST algorithms with some known DST algorithms. In Section VII, we show the connection between the algebraic operations used to compute the proposed DST algorithms and the signal flow graph building blocks to provide an architecture for the DST-based integrated circuit design. Finally, we describe in Section VIII an application of the proposed DST algorithms within a DRPE image encryption scheme.

II. FREQUENTLY USED NOTATIONS
First, we provide notations in [14] to distinguish and show the novelty of the proposed DST matrix factorization in comparison to the results in [14]. Next, we will introduce additional notations that we use throughout the paper.
A. FACTORIZATION FOR DST-I AND DST-II IN [14] Here, we provide formulas in computing DST-I/II matrices in [14]. The relationship F n = W Q · V Q , where F n is the DCT or DST matrix of types I-IV, W Q is the weight matrix, and V Q is the Chebyshev-like polynomial Vandermode matrix for the DCT or DST matrices, was first introduced in [43]. Later in [14], the polynomial Vandermonde matrix V Q = Q k (x j ) n,n−1 j,k=1,0 was further split into self contained factors to obtain factorization for the DCT and DST matrices.
The weight matrix for DST-I in [14] and [43] for odd N is given by The weight matrix for DST-II in [14] and [43] for even n is given by The self-contained factorization for the Chebyshev-like polynomial Vandermode matrix for DST-I (for odd N ) and DST-II (for even n) was derived in [14] and stated as follows where  for k = 0, 1, · · · , N − 1 or n − 1, respectively for DST-I and DST-II.
The factorization for V odd was derived for DST-I (for odd N ) and has the form where and with respect to the nodes and the classical Chebyshev polynomials T k (d) = cos(k arccos d).
The factorization for V odd was derived for DST-II (for even n) and has the form where and , · · · , 1 cos (n−1)π 2n   , By using the above factorization formulas for the Chebyshev-like polynomial Vandermonde matrices and the relationship F n = W Q · V Q in [43], the following results were obtained in [14] to factor DST-I (for odd N ) and DST-II (for even n) matrices and The factorization formulas (14) and (15) were obtained by using the factorization formula (3) with the replacement of (from [43]), respectively. Although it says in [14] that the first matrix in the factorization (3) is the odd-even permutation matrix, it must actually be the transpose (inverse) of the odd-even permutation matrix. We also note here that the factorization formulas to compute DST-III/IV matrices in [14] are not included because the proposed DST-III/IV matrix factorization formulas are completely different from the factorization formulas in [14].

B. OTHER NOTATION
We state here diagonal, bidiagonal, sparse, and orthogonal matrices which are frequently used in this paper. For a given vector x ∈ R n and an integer n ≥ 3, let us introduce an evenodd permutation matrix P n by an odd-even permutation matrixP n bỹ For even n s.t. n ≥ 4, we introduce a sparse matrix a bidiagonal matrix diagonal matrices and where , and a matrix with three liftings (we will show later that the following sparse and orthogonal matrix,Q n , is the same as the matrix Q n with three lifting steps performed) .

III. SPARSE, ORTHOGONAL, BIDIAGONAL, DIAGONAL, AND SELF-CONTAINED FACTORS FOR DST MATRICES
In this section, we will propose formulas to factor DST I-IV matrices into the product of sparse, diagonal, bidiagonal, and orthogonal matrices. Among the factorization formulas, the self-contained factorization of the DST-II matrix will be used to establish a factorization formula for the DST-III matrix. Also, an orthogonal matrix with three lifting steps and the factorization formula of the DST-II matrix are utilized to obtain a factorization formula for the DCT-IV matrix.
A. SPARSE, ORTHOGONAL, DIAGONAL, AND SELF-CONTAINED FACTORS TO COMPUTE DST-I MATRIX Let us state a novel self-enclosed factorization to compute the DST-I matrix using simple, diagonal, sparse, and orthogonal matrices. Lemma 1 (Self-enclosed DST-I factorization): For an even integer n ≥ 4, the matrix S I n−1 can be factored into the form Proof: We factor the weight matrix W Q S I defined in the equation (1) into the product of matrices s.t.
where W even = diag sin( 2kπ n ) . Now, we use the trigonometric identities to factor G S I in equation (7) and G S I in equation (8) into the product of diagonal and sparse matrices s.t.
Thus, we substitute this relationship into the equation (30) and modify the last entry in the (2,2) block of the middle matrix in (30). This is modified to fix the scaling because the n 2 , n 2 entry of G S I is 1 √ 2 and scaling factor of S I n 2 −1 is 4 n . Then, we add a scaling factor of 1 √ 2 toH n−1 (i.e. modify the n − 1, n 2 entry ofH n−1 ) to obtain the factorization proposed in (26).
Remark 2: By comparing the factorizations of the DST-I matrix in (14) and (26) (i.e. distinguishing the proposed selfcontained factorization of the DST-I matrix from the DST-I factorization in [14]) one can observe that (a). the weight matrix is congregated in (26) as opposed to scattered weight matrices in (14), (b). the weight matrix W Q S I in (2) was factored into the product of permutation and odd/even weight matrices in (27), (c). the matrix G S I in (29) has been simplified and factored into the product of sparse and diagonal matrices as opposed to the matrix G S I in (7), (d). the matrix G S I in (29) has been simplified and factored into the product of diagonal matrices as opposed to G S I in (8) (26) is in the simplest form as opposed to the factorization of the DST-I matrix in (14).

, (f). the matrix H n−1 in (26) is orthogonal as opposed to the nonorthogonal matrix H S I in (4). Thus, the proposed self-contained factorization of the DST-I matrix in
Remark 3: The factorization of the DST-I matrix established in Lemma 1 also yields to the following alternative factorization formula through a permutation where P n−1 is an even-odd permutation matrix and H n−1 is defined in (22).

B. BIDIAGONAL, DIAGONAL, ORTHOGONAL, AND SELF-CONTAINED FACTORS TO COMPUTE DST-II MATRIX
Let us state a novel self-contained factorization to compute the DST-II matrix using sparse, bidiagonal, diagonal, and orthogonal matrices. Lemma 4 (Self-enclosed DST-II factorization): For an even integer n ≥ 4, the matrix S II n can be factored into the form Proof: We factor the weight matrix W Q S II defined in the equation (2) into the product of matrices s.t. where . Now, we can use the inverse, reciprocal, and cofunction trigonometric identities to factor G S II in the equation (11) and G S II in the equation (12) into the product of bidiagonal, diagonal, and sparse matrices s.t. (17). While factoring the matrix G S II , the n 2 , n 2 entry of B n 2 is modified into √ 2 because the n 2 , n 2 entry of W even is defined as 1 √ 2 in (33). Hence by (33), (34),

and using
the V odd defined in (10), we can obtain whereH n is given in (20). By following [43], we can obtain Thus, from the equations (35) and (36) Finally, by using the above equation and the trigonometric identities to simplify the formulaĨ n 2 W −1 oddĨ n 2 , we can obtain the factorization formula proposed in (32). (15) and (32) (i.e. distinguishing the proposed self-contained factorization of the DST-II matrix from the DST-II factorization in [14]) one can observe that (a). there are not many weight matrices spread in (32) like in (15), (b). a factorization for W Q S II is presented via a permutation and odd/even weight matrices in (33) as opposed to W Q S II in (2), (c). a factorization for G S II is presented via a simple bidiagonal and diagonal matrices in (34) as opposed to G S II in (11), (d). a factorization for G S II is presented via antidiagonal and diagonal matrices in (34) as opposed to G S II in (12), (e). an orthogonal matrixH n is given in (20) as opposed to the scaled orthogonal matrix H S II in (15), and (f). identified the exact permutation matrix for V Q S II within the Lemma 4 as opposed to V Q S II in (3). Thus, we have expressed an exact, simplified, and self-contained factorization formula to compute the DST-II matrix using bidiagonal, diagonal, sparse, and orthogonal factors in (32) as opposed to (15). We propose a novel factorization formula to compute the DST-IV matrix using bidiagonal, diagonal, sparse, and orthogonal matrices. The proposed factorization of the DST-IV matrix is based on the self-enclosed DST-II factorization in Lemma 4. Before we proceed, let us recall the DST-IV matrix factorization formula in [9], [12] s.t.

Remark 5: By comparing the factorizations of the DST-II matrix in
To derive a novel factorization formula to compute the DST-IV matrix, we will use the above factorization formula and apply three lifting steps [44] for the rotation matrix cos θ sin θ − sin θ cos θ and the rotation/reflection matrix cos θ sin θ sin θ − cos θ . Let us first derive a factorization formula for the rotational/rotational-reflection butterfly matrix Q n defined in (24). Lemma 8 (Butterfly matrix Q n is the same as the matrix with three lifting stepsQ n ): Let n ≥ 2 be an even integer. The matrix with three lifting steps,Q n , is the same as the rotational/rotational-reflection matrix Q n and can be factored in the form: where T n 2 = tan (2k+1)π 8n n 2 −1
Proof: Let us consider the block products inQ n explicitly: To verify the equivalency ofQ n and Q n , we will consider the corresponding blocks in both matrices.
Next, we use the immediate result to state a novel recursive factorization formula to compute the DST-IV matrix using the DST-II, sparse, and orthogonal matrices. We recall here that the DST-II matrix can further be factored into the product of sparse, bidiagonal, diagonal, and orthogonal matrices as shown in Lemma 4 followed by Remark 6.
where V n is defined in (23) andQ n is defined in (41). Proof: Lemma 8 shows thatQ n in (41) is the same as Q n in (24). This result together with the explicit derivation for the factorization of the DST-IV matrix in [9] and [12] gives the result (47).
Remark 10: The DCT I-IV matrix factorizations proposed in [41], [42] and the interweaving relationships among cosine and sine matrices in [10] can be utilized to obtain variants of the proposed DST I-IV matrix factorization formulas.

IV. RELATIONSHIP AMONG DST MATRICES AND CONNECTION AMONG DST FACTORIZATIONS IN THE LITERATURE
We first present a relationship between DST-I and DST-III matrices and next between the DST-II and DST-IV matrices using bidiagonal, diagonal, and sparse matrices. After these relationships are established, we will observe connections among the proposed DST matrix factorizations with the existing factorization of the DST matrices.   . . .
Let us evaluate the elements of A. First, we consider elements in the first row, but the first to the second-to-last column of (49), i.e., j = 0 and k = 0, 1, . . . , n 2 − 2.
The immediate result is true for any even integer n ≥ 2 satisfying the following relationship.
Corollary 12 (Relationship between S I n−1 and S III n ): For an even integer n ≥ 2, the matrix S III n is related to the matrix S I n−1 in the following way Proof: This is trivial from the Lemma 11.  · · · sin ( n 2 −1)(n−1)π n Let us evaluate the elements of M . First, we consider elements in the first row of (56), i.e., j = 0 and k = 0, 1, · · · , n/2 − 1.
Thus, from (57), (58), and (59), we prove the claim. The immediate result is true for any even integer n ≥ 2 satisfying the following relationship.
Corollary 14 (Relationship between S II n and S IV n ): For an even integer n ≥ 2, the matrix S IV n is related to the matrix S II n in the following way S IV n = B n S II n [W ] n .
Proof: This is trivial from the Lemma 13.

A. CONNECTION BETWEEN THE PROPOSED DST FACTORIZATION AND THE EXISTING DST FACTORIZATIONS
In this section, we utilize the relationship among DST matrices established in Section IV to identify a relationship among the proposed DST factorization and the existing factorization of DST matrices. The relationship between DST-II and DST-IV will be utilized to establish a connection between the proposed DST-II factorization with the existing DST-II factorization in [7], [9], [11], [12], [45]- [47]. We will also show the DST-III factorization in [7], [9]- [12], [46] can be seen as the proposed self-contained DST-III factorization. The relationship between DST-I and DST-III will be utilized to establish a connection between the proposed DST-I factorization with the existing DST-I factorization in [7], [9]- [12], [46], [47].

Corollary 15 (Connection between traditional DST-I and the proposed DST-I): For an even integer n ≥ 4, the matrix S I n−1 can be factorized into the form
Proof: This follows immediately from Lemma 1, Remark 3, and Lemma 11.
Corollary 16 (Connection between traditional DST-II and the proposed DST-II): For an even integer n ≥ 4, the matrix S II n can be factored into the form Proof: This follows immediately from Lemma 4, Remark 6, and Lemma 13.
The following result shows that the proposed factorization of the DST-III matrix can be converted into the existing DST-III factorization in [7], [9]- [12], [46].  (26), (32), and (39) in this paper, are simpler than the existing interdependent factorizations for the DST-I/II/III matrices. Moreover, the proposed DST factorizations for the types I/II/III can be seen as a real variant of the self-contained factorization of the DFT matrices in [48], [49]. Also, the proposed DST-I/II/III factorizations are easy to implement and establish a bridge to the existing DST factorizations in [7], [9]- [12], [45]- [47]. We will see in Section V that the proposed algorithms are self-recursive as opposed to completely recursive algorithms in the literature. Moreover, the new algorithms require only simple, sparse, scaled orthogonal, diagonal, and bidiagonal matrices to compute the DST matrices by a vector as shown in the next section.

V. SELF/COMPLETELY RECURSIVE AND RADIX-2 DST-I/II/III/IV ALGORITHMS
The factorizations of DST I-IV matrices established in Section III lead to self/completely recursive and radix-2 algorithms. These algorithms can be used to compute DST matrices by a vector using simple, sparse, diagonal, bidiagonal, VOLUME 9, 2021 and orthogonal matrices. In this section, we present lowcomplexity, recursive, and radix-2 DST I-IV algorithms. Similar to [12], we remove the scaling factor 1 √ 2 from H n , H n−1 , V n , S II 2 , and S III 2 for further reduction of multiplication complexity. Compared to most of the existing algorithms, the proposed algorithms are easy to implement and attain the lowest multiplication complexity.
Before stating the algorithms, let us define the following notation to denote diagonal and bidiagonal matrices which will be used hereafter for n ≥ 4.
Also, from now onwards we call √ 2V n =Ṽ n . We recall here from the notations that  Proof: The number of additions and multiplications required to compute the DST-IV algorithm can be found by using the number of additions and multiplications required to compute the DST-II algorithm at n 2 . By using the DST-IV algorithm, we get The structures ofQ n andṼ n gives us Following the arithmetic complexity in computing DST-II algorithm in Lemma 24 and using the equations (72) and (71), we get Simplifying the above yields the multiplication complexity of the modified DST-IV algorithm as in (70).
Remark 28: (a) The paper [54] presents the total lowest arithmetic operations in calculating DST-IV using the DFT and DCT-IV while leaving an open question: whether the DST-IV algorithm in [54] leads to practical gains in performance on real computers. To compute DST-IV, the authors of [54] used twiddle factors but these factors can be overlapped with other operations in the CPU [55]. Also to compute DST-IV, authors in [54] have used scaled factors of the form ±1 ± i tan(2πk/n) and ± cot(2πk/n) ± i, where i 2 = −1 but these factors contain singular functions. Hence, the practical gains of the lowest flop count for DST-IV calculation in [54] may lead to a numerically ill-posed problem. We will show in Section VI-B, that the proposed DST-IV algorithm attains the lowest multiplication complexity algorithm to compute the DST-IV matrix by a vector. (b) The arithmetic complexity of the proposed DST-IV is the same as that of lowest multiplication complexity DCT-IV in [42] because DST-IV is exactly equivalent to a DCT-IV in which the outputs are reversed and every other input is multiplied by -1 (or vice-versa) [50], [54]. Although the total flop count of the DST-I algorithm in [47] is lower than all other DST-I algorithms in [7], [10], [12] and the proposed algorithm, the sin1(x, n) algorithm is presented using simple, sparse, diagonal, and scaled orthogonal matrices as opposed to the DST-I factorization in [47].
Based on the counts in computing DST-II/III algorithms, the proposed DST-II/III algorithms have the lowest multiplication counts. We recall here that the authors in [53] calculated the lowest flop count split-radix DST-II/III algorithms. By considering the fact that the authors in [53] have reduced the flop count only by cutting off multiplication counts while preserving the same addition counts as in [56], the explicit multiplication count in computing splitradix DST-II/III algorithms in [53] is given via (7/18)nt + (10/27)n − (1/9)(−1) t t + (7/54)(−1) t + 1/2. In this situation, the split-radix DST-II/III algorithms attain the lowest multiplication counts for n ≥ 32 but without using simple,    sparse, diagonal, bidiagonal, and scaled orthogonal matrices. Furthermore, the split-radix DST-II/III algorithms leave a question based on the practical gains in performance on real computers as the algorithms in [53] use singular functions and hence may lead to ill-posed problems. By considering the fact that the authors in [53] have reduced the flop count in computing DST-IV matrix by a vector only by cutting off multiplication counts, the explicit multiplication count in computing the split-radix DST-IV algorithm in [53] is given via (5/9)nt + (37/27)n + (2/9)(−1) t t − (10/27)(−1) t − 2. Even though the total flop count of the split-radix DST-IV algorithm in [54] is lower than all the other DST-IV algorithms in the literature, the multiplication counts of the proposed sin4(x, n) algorithm and the DST-IV algorithm in [47] is the lowest in the literature. Moreover, the proposed DST-IV algorithm is presented using simple, sparse, diagonal, bidiagonal, and scaled orthogonal matrices as opposed to the DST-IV factorization in [47] and any other existing DST-IV algorithms in the literature. Note that if the factorization of the DCT or DST does not preserve orthogonality, the resulting DCT or DST algorithms lead to inferior numerical stability [61].
To sum up, the proposed self/completely recursive and radix-2 sin2(x, n), sin3(x, n), and sin4(x, n) algorithms can be computed via simple, sparse, diagonal, bidiagonal, and scaled orthogonal matrices having the lowest multiplication counts in the literature.

VII. SIGNAL FLOW GRAPHS FOR DST-I/II/III/IV ALGORITHMS
Signal flow graphs can be used to represent a physical system and its controllers with the help of a system of linear equations. We present signal flow graphs of the proposed DST algorithms while illustrating the connection between system components, associated equations, and the simplicity of the proposed factorizations of the DST I-IV matrices. The signal flow graphs are drawn using decimation-in-frequency. It is possible to convert these signal flow graphs into decimationin-time flow graphs.

VIII. ENCRYPTION USING DST ALGORITHMS
In [62], a 4f encryption scheme was proposed which introduces a chaotic mapping to the established double random phase encoding (DPRE) scheme originally proposed in [63]- [65]. The DPRE process applies to a phase image a random phase mask followed by a Fourier transform twice in series. While the DPRE scheme maximizes the entropy of the encrypted image, the linearity of the scheme creates vulnerabilities [66], [67]. Thus, [62] proposed an encryption scheme that applies a chaotic pixel mapping based on the famous 3D Lorenz System after the first Fourier transform. The Lorenz system, also commonly known as the Lorenz path or Lorenz attractor, was first described in [68]. The path itself is the solution to the set of ordinary differential equations where σ , ρ, and β are system parameters and the initial conditions are x(0) = x 0 , y(0) = y 0 , z(0) = z 0 . Given time to diverge, this set of equations will produce completely different paths with even the slightest difference in the system parameters. This chaotic nature, along with the relative speed by which the set of equations may be solved, makes the Lorenz system a valuable tool for mapping and permutation techniques.
In the following, we describe image encryption based on the DPRE encryption technique in [62] and the proposed DST algorithms. f (x, y), where x and y are locations in the spatial domain, is converted into a phase image p 0 (x, y) = e iπf (x,y) . 2) A random phase mask is applied to each pixel as defined by p 1 (x, y) = p 0 (x, y)e 2πir(x,y) , where r is a random matrix having the same size as the original image.

1) An image with pixel values
3) The image is transformed into the Fourier domain using the 2D DFT and the proposed sin2 algorithm in 2D. 4) A chaotic mapping technique is applied based on the 3D Lorenz system. The Lorenz system is solved using Matlab's ode45 function. First, the image is decomposed into n × n blocks, where n is a multiple of the total length of the image. The red, green, and blue color channels of each block are associated with the position vectors of the Lorenz path in time. These blocks are then sorted in ascending, row-major order based on the value of the associated Lorenz system position. This produces a new image p L . 5) A random phase mask is applied to each pixel as defined by p 2 (u, v) = p L (u, v)e 2πis(u,v) , where s is a random matrix of the same size as the original image, u and v are positions in the Fourier domain. 6) The proposed sin2 algorithm and the DFT are applied to p 2 to produce the fully encrypted image. The final result is an image of high entropy with a Gaussian distribution of pixel values across all channels. To decrypt the image, the process is simply performed in reverse using the same random phase mask keys, Lorenz system parameters, and the inverse of the respective transform. For the decryption of the image, the proposed 2D sin3 algorithm and the 2D inverse-DFT are used. For a block length n less than or equal to 1 8 th of the total image length (i.e. n = 32 for a 256 × 256 image), attempts to decrypt the image using Lorenz system initial conditions with even the most minute difference will fail. Blocks which are 1 4 th the total image length (i.e. n = 64 for a 256 × 256 image) will reveal some information if transformed into the Fourier domain using the DFT; however, we will show below that images transformed using the proposed DST algorithms will withstand these attacks. Figure 5 shows the original image which was used for testing purposes and four attempts to decrypt using an incorrect key for the Lorenz path. The corresponding image histograms are also displayed. The Lorenz parameter key used to encrypt the image is  employing larger mapping blocks when the proposed sin2 algorithm is used, as shown by Figures 5.g and 5.h. When the encryption scheme uses these block sizes with the DFT, however, the scheme starts to become liable to attack. As illustrated by Figure 5.j, the decryption attack does not yield an even distribution of pixel intensities across the channels. Likewise, Figure 5.i starts to reveal the characteristic blue and red nose of the baboon in the original image.

IX. CONCLUSION
We have proposed novel factorization formulas to compute DST I-IV matrices using the product of simple, sparse, diagonal, bidiagonal, and scaled orthogonal matrices. The proposed DST matrix factorizations are used to establish self/completely recursive and radix-2 DST algorithms. We have also established a novel algebraic relationship between DST-II and DST-IV matrices using diagonal and bidiagonal matrices and described connections between traditional DST-II/III algorithms and the proposed DST-II/III algorithms. Moreover, we have established another novel algebraic relationship between DST-I and DST-III matrices using sparse and diagonal matrices. We have also described connections between traditional DST-I algorithms and the proposed DST-I algorithm. The arithmetic complexity of the proposed algorithms is analyzed to show that the proposed algorithms are much more efficient than the existing algorithms. Specifically, the DST-IV algorithm attains the lowest multiplication complexity for all transform matrix sizes n ≥ 8 in the literature. The DST-II/III algorithms are the lowest multiplication complexity, radix-2, and self-recursive algorithms and execute using simple, sparse, diagonal, bidiagonal, and scaled orthogonal matrices for the transform matrix sizes n ≥ 8. The presented signal flow graphs illustrate the simplicity of the proposed DST algorithms and provide a foundational framework for the DST-based integrated circuits. The robust performance of the proposed algorithms within software applications is explored through implementation within a DPRE encryption scheme with a chaotic Lorenz path mapping. The proposed DST algorithms in 2D can increase encryption performance when compared to an implementation using the 2D DFT.