1

Regularity-constrained Fast Sine Transforms
Taizo Suzuki, Senior Member, IEEE, Seisuke Kyochi, Member, IEEE, and Yuichi Tanaka, Senior Member, IEEE

Abstract—This letter proposes a fast implementation of theregularity-constrained discrete sine transform (R-DST). The
original DST leaks the lowest frequency (DC: direct current)
components of signals into high frequency (AC: alternating
current) subbands. This property is not desired in many applications, particularly image processing, since most of the frequency
components in natural images concentrate in DC subband. The
characteristic of filter banks whereby they do not leak DC
components into the AC subbands is called regularity. Whilean R-DST has been proposed, it has no fast implementation
because of the singular value decomposition (SVD) in its internal
algorithm. In contrast, the proposed regularity-constrained fastsine transform (R-FST) is obtained by just appending a regularity
constraint matrix as a postprocessing of the original DST. When
the DST size is M × M (M = 2ℓ , ℓ ∈ N≥1 ), the regularity
constraint matrix is constructed from only M/2 − 1 rotation
matrices with the angles derived from the output of the DST
for the constant-valued signal (i.e., the DC signal). Since itdoes not require SVD, the computation is simpler and fasterthan the R-DST while keeping all of its beneficial properties.
An image processing example shows that the R-FST has finefrequency selectivity with no DC leakage and higher coding gain
than the original DST. Also, in the case of M = 8, the R-FST
saved approximately 0.126 seconds in a 2-D transformation of
512 × 512 signals compared with the R-DST because of fewer
extra operations.
Index Terms—DC leakage, discrete sine transform, fast implementation, regularity, rotation matrix.

I. I NTRODUCTION

THE discrete cosine transform (DCT) [1] is one of themost common tools in various signal processing applications. Among the eight types (I to VIII) of DCT, type II is
known to be especially effective for image processing and has
been used in many practical applications, such as image/videocoding (JPEG, AVC, HEVC, and VVC [2–5]) and restoration(denoising, deblurring, and so on [6–8]). In this letter, we focus
on the type-II DCT; hereafter, we term it simply “the DCT.”
The key properties of the DCT for image processing are as
follows. Firstly, the DCT has fine frequency selectivity, which
is suitable for natural images. Secondly, the DCT satisfies theregularity condition: if a linear transform satisfies the regularity condition, it does not leak the lowest frequency (DC: direct
current) components into the high frequency (AC: alternating
current) subbands. These two properties contribute to the DCT
producing a sparse image representation. Thirdly, the DCT hasvarious fast algorithms [9–13]. The traditional rotation matrix
This work was supported by JSPS Grant-in-Aid for Scientific Research (C)under Grant 22K04084.
T. Suzuki is with Faculty of Engineering, Information and Systems, University of Tsukuba, Ibaraki 305-8573, Japan (e-mail: taizo@cs.tsukuba.ac.jp).
S. Kyochi is with Department of Computer Science, Kogakuin University,
Tokyo 163-8677, Japan (e-mail: kyochi@cc.kogakuin.ac.jp).
Y. Tanaka is with Department of Electrical Engineering and Computer
Science, Tokyo University of Agriculture and Technology, Tokyo 184-8588,
Japan (e-mail: ytnk@cc.tuat.ac.jp).

factorizations [14–16] are the underlying methods of many fast
DCT algorithms.
On the other hand, the discrete sine transform (DST) [1] isalso one of the most useful frequency transforms. Although the
DST has several types, we focus on the type-II DST, which is acounterpart of the DCT in terms of sines and cosines; hereafter,
we call it simply “the DST.” The DST shares many properties
with the DCT, such as fine frequency selectivity and fast
implementations. However, it lacks the regularity condition.
This is not desired in many signal processing applications.
Our previous work [17] presented a regularity-constrained
DST (R-DST) that has fine frequency selectivity and, unlike
the original DST, satisfies the regularity condition; however,it has no fast implementation because of the singular value
decomposition (SVD) in its internal algorithm. In practice,
the resulting R-DST can be factorized into a multiplication of
the DST and an extra orthogonal matrix that is half the sizeof DST. However, more details have not been discussed.
This letter introduces a fast implementation of the R-DST,
called the regularity-constrained fast sine transform (R-FST).
The R-FST is the first DST variant to have the three desired
properties, i.e., fine frequency selectivity, regularity, and a fast
implementation. It is obtained by just appending a regularity
constraint matrix of a postprocessing of the original DST.
When the DST size is M × M (M = 2ℓ , ℓ ∈ N≥1 ), theregularity constraint matrix is constructed from only M/2 − 1
rotation matrices with the angles derived from the output of
the DST for the constant-valued signal (i.e., the DC signal).
The R-FST preserves the beneficial properties of the R-DST,
while its computation is much simpler. An image processing
example shows that the R-FST has fine frequency selectivity
with no DC leakage and higher coding gain [18] than that of
the original DST. Also, the R-FST achieves more time savings
than the R-DST because of fewer extra operations.
Notation: Boldface letters represent vectors and matrices.
Superscript ⊤ denotes transpose of a matrix/vector. Note thatthe size and dynamic range of the images in the figures havebeen adjusted for ease of visualization.
II. R EVIEW AND D EFINITIONS
A. Discrete Cosine and Sine Transforms
The (m, n)-elements (m, n ∈ N[0,M −1] ) of the M ×M DCT
C and DST S are defined asq 1(m = 0)M,(1)[C]m,n = q 2 cos π m n + 1
(m
=
̸
0)
M
M
2
q
 1 (−1)n(m = M − 1)
M[S]m,n = q, 2 sin π (m + 1) n + 1(m ̸= M − 1)
M
M
2
(2)

2

respectively. In fact, S is obtained from C by slightly manipulating C with order inversion and sign flip, as follows:
S = JCD,(3)

where the (m, n)-elements of J and D are expressed as(
1 (m = M − n − 1),
[J]m,n =0 (m ̸= M − n − 1)((−1)m (m = n)[D]m,n =,0(m ̸= n)(4)(5)

respectively.

The (m, n)-elements of an M × M rotation matrix Θ(i,j)
with a rotation angle θi,j (i, j ∈ N[0,M −1] ) for the i and jth
signals are represented as follows:Θ(i,j) m,n1(m = n ̸= i) or (m = n ̸= j)





cos
θ
(m = n = i)i,j.
= − cos θi,j (m = n = j)sin θi,j
(m
=
i
&
n
=
j)
or
(m
=
j
&
n
=
i)



0
(otherwise)(6)
Generally, an M × M orthogonal matrix Φ is constructed bycascading M (M −1)/2 rotation matrices with suitable rotation
angles,Φ=

i+1Yrank(S) = M − 1.

Θ(i,j) .

(7)

Note that the order of the rotation matrices is interchangeable.
Also, in some cases, the number of rotation matrices can beless than M (M − 1)/2; e.g., an 8 × 8 DST S [8] consists of
13 (< 8(8 − 1)/2 = 28) rotation matrices.
C. Regularity Condition of Orthogonal Matrix
If Φ satisfies the regularity condition, all the energy of theDC signal is concentrated in one transformed coefficient. Thisis represented as follows:
√
⊤

⊤
(8)
M 0 ··· 0 = Φ 1 1 ··· 1 ,
|
{z
}
|
{z
}
1

(12)

−1
Let us modify S. From (8), {sm }M
m=1 should be orthogonal
to s0 in order to impose the regularity condition on S. We canorthogonalize the odd rows of S in the following way.
Step 2: Set
e (0) = s0 0 s2 · · · sM −1 ⊤ .
S(13)e (0) also satisfiesSe (0) = M − 1.
rank S(14)

Thus, there is only one zero singular value and its corresponding right-singular vector (denoted as v(0) ) belongs to the nulle (0) . This implies thatspace of Se (0) v(0) = 0,Si=M −2 j=M −1

y with no DC leakage

This subsection reviews the design procedure of the original
R-DST [17]. Since the even rows of S do not satisfy the regularity condition, we replace these rows with ones approximatedby SVD. The procedure is as follows.
Step 1: Define an M × M modified DST S as⊤
S = s0 s1 · · · sM −1 ,
(10)
q
 1
(m = 0)M.
(11)[S]m,n = q
 2 sin π m n + 1  (m ̸= 0)
M
M
2
In short, it is constructed by replacing the (M − 1)th row of
S with that of the 0th row of C and then permuting the rows.
S satisfies

B. Rotation Matrix and Orthogonal Matrix

1
Y

D. Regularity-constrained Discrete Sine Transforms(15)e (0) is updated byi.e., v(0) satisfies the regularity condition. S(0)replacing 0 with v .
Step 3: Set⊤S(0) = s0 v(0) s2 · · · sM −1 .
(16)Note that v(0) can be explicitly represented as
r
h
i
1
(0)
(−1)n = [S]M −1,n
v
=
M
n

(17)

because the row of [S]M −1,n corresponding to the highest
frequency subband is orthogonal to {s0 , s2 , · · · , sM −1 }. It
clearly follows that
rank S(0) = M.
(18)

where, without loss of generality, we assume the DC subband
corresponds to the 0th element in the transformed coefficients.
Note that all rows except for the 0th row of C satisfy theregularity condition, while the even rows of S do not:
⊤× 0 × 0 · · · × 0 = S1,(9)Consequently, by repeating Steps 2 and 3, we can obtain an
orthogonal matrix S(M/2−1) whose odd rows are replaced bydifferent ones from the initial S(0) . Algorithm 1 is a summaryof the algorithm. Hereafter,

where × is a nonzero constant.

is termed R-DST.

S := S(M/2−1)

(19)

3

Algorithm 1 Design procedure for RDST
1: Set S is as in (11).
2: for k = 0 to M/2 − 1 doe (k) = · · · s2k 0 s2k+2 · · · ⊤ .
3:
Set S
4:
Find the right-singular vector v(k) corresponding to zero singularvalue.
⊤5:
Set S(k) = · · · s2k v(k) s2k+2 · · · .
6: end for
7: Output S = S(M/2−1) .
Fig. 1.

III. R EGULARITY- CONSTRAINED FAST S INE T RANSFORMS
This section describes the R-FST. First, how to impose theregularity condition on a general orthogonal matrix is given;
then it is applied to define the R-FST in Section III-A. After
that, design examples are given in Section III-B. Finally, othertypes are discussed in Section III-C.
A. Derivation and Definition
To impose the regularity condition on an M × M generalorthogonal matrix Φ in (7) including the M ×M DST S in (2)and (3), we multiply an M × M orthogonal matrix R, calleda regularity constraint matrix, after the original Φ, i.e., RΦ.
According to the procedure used to investigate the regularity
condition in (8), we define a(1) as the DC signal transformed
by Φ:
h
i⊤
(1)
(1)
a(1) = a(1)
= Φ1.
(20)
a
·
·
·
a
0
1
M −1
(k)

b
When the kth (k ∈ N≥1 ) rotation matrix Θ(i,j) , whose the(m, n)-elements are defined as
i
h
b (k)
Θ
(i,j)
 m,n

1
(m = n ̸= i) or (m = n ̸= j)



(k)

b

(m = n = i)cos θi,j(k)b= − cos θi,j (m = n = j)
,


(k)

b
sin
θ
(m
=
i
&
n
=
j)
or
(m
=
j
&
n
=
i)

i,j



0
(otherwise)(21)
in R and the related signals a(k) and a(k+1) are expressed as
(k)

(k)
b
a(k+1) = Θ,(i,j) a(k)it means that the jth element aj
(k+1)

aj

(k)

= ai

(22)

in a(k) changes into

(k)

(k)

(k)

sin θbi,j − aj cos θbi,j .

(23)

To make the jth rows of RΦ satisfy the regularity condition,(k)
let the angle θbi,j be(k)(k)θbi,j = arctan
(k+1)

aj

(24)

(k)

ai

(k)

(k)

for aj= 0, arbitrary ai (̸= 0), and arbitrary aj . To
completely impose the regularity condition on RΦ, let R be
the following matrix:R=1Y
j=M −1b (j) .
Θ
(0,j)

(25)

Implementation of M × M R-FST.

(k)b
Although there are various combinations of Θ(i,j) , we use(25), which is one of the simplest.
To impose the regularity condition on the DST, the resulting
R in (25) can be directly applied to S. However, it has
redundancy because the odd rows of S already satisfy theregularity condition, as shown in (9). Consequently, the M ×M
b is defined by just appending an M × M regularity
R-FST Se consisting of only M/2 − 1 rotation
constraint matrix R
matrices as a postprocessing of S, as follows (Fig. 1):1Yb (j)  JCD,b := 
(26)
ΘS(0,2j)| {z }
j=M/2−1

|

{z

S

}

e
R
(k)
(k)
where a0 and aj are redefined as

h
(k)
(k)
a(k) = a(k)0 a20 · · · aM −2
0(
S1(k = 1)= b (k−1)
Θ(0,2(k−1)) a(k−1) (k ≥ 2)

i⊤
0
(27)

(k)
b can be efficiently definedfor a2(k−1) = 0 (k ≥ 2). Hence, S
in a sequential manner.

B. Design Examples
This subsection shows some design examples of the M ×M
b [M ] and its identity except for the order and sign,
R-FST S
which do not affect the properties, with the M × M R-DST
S[M ] in [17].
Case of M = 2: From (2), we easily find that S [2] alreadysatisfies the regularity condition and can be directly definedb [2] as follows:
as S


[2]
1 1 1
[2]
b
S =S = √= S[2] .
(28)
1
−1
2
Case of M = 4: From (2), S [4] is expressed by
√
√
√
√

2 sin π82 cos π82 cos π82 sin π8

1
√ 1 π
√1 π
√−1 π √ −1 π  .
S [4] = 
2 cos 8 − 2 sin 8 − 2 sin 8
2 cos 8 
2
1
−1
1
−1
(29)
From (9), S [4] causes DC leakage only in the 2nd row,
√ 


⊤
0 cos π8 − sin π8
0 = S [4] 1.
2 cos π8 + sin π8
(30)

4

TABLE IC ODING GAINS [ D B].

DST
R-DST [17]
R-FST
(HT)

Fig. 2. Frequency responses of DSTs (blue dashed lines), R-FSTs (red solidlines), and HTs (green dotted lines): (left) 4 × 4 and (right) 8 × 8.

Barbara transformed with 8 × 8 (left) DST and (right) R-FST.

According to Section III-A, a 4 × 4 regularity constraint
e [4] is constructed from only a rotation matrix Θb (1) .
matrix R(0,2)
bConsequently, S[4]

is as follows:

1
1
1
[4]
1
1
(1)
[4]
b =Θ
b
S
= 
(0,2) S
−1 1
2| {z }
1 −1
e [4]
R

1
−1
1
1


1
−1
 = S[4] .
−1
−1
(31)

Cases of M ≥ 8: We can also calculate them in the sameb [8] with an
way of the case of M = 4. For example, we find Se [8] constructed from three
8 × 8 regularity constraint matrix R
rotation matrices, as follows:[8]b [8] = Θb (3) Θb (2) b (1)S= S[8] .
(0,6) (0,4) Θ(0,2) S
|
{z
}
e [8]
R

Note that the cases of M = 2 and 4 are coincidentally special
ones that are identical to the Hadamard transforms (HTs),
whereas the larger cases are not so.

M =4
4.73
7.17
7.17
(7.17)

M =8
5.09
7.72
7.72
(7.95)

M = 16
6.02
7.85
7.85
(8.19)

M = 32
7.24
8.09
8.09
(8.27)TABLE IIN UMBER OF EXTRA MULTIPLICATIONS (MUL) AND ADDITIONS (ADD).

R-DST [17]R-FST

Fig. 3.

M =2
5.05
5.05
5.05
(5.05)M =8
MUL ADD1612126

M = 16
MUL ADD
64
56
28
14

M = 32
MUL
ADD
256
240
60
30

R-DSTs in [17] are not shown because they were identical tothose of the R-FSTs. We find that all rows except for the 0th
row of the R-FST cause no DC leakage.
In addition, Table I shows that the coding gains [18],whereby a transform with a higher gain compacts more energy
into fewer coefficients and it is not depending on images, ofthe R-FSTs are greater than those of the DSTs and are identicalto those of the R-DST. Also, we can see that the R-FSTs inthe cases of M = 2 and 4 are identical to the HTs and onesin the larger cases are not so.
Moreover, it is clear that the R-FST is simpler and fasterthan the R-DST and it is easy to utilize the existing fastalgorithms for the following reasons:1) There are many fast algorithms for the DCT.
2) The DST is obtained by slightly manipulating the DCT.
3) The R-FST is implemented by just appending M/2 − 1
rotation matrices as a postprocessing of the DST.
The difference between the R-DST and R-FST is in 3). The
R-DST requires extra M 2 /4 multiplications and (M − 2)M/4
additions compared to the DST. In contrast, the R-FST requires
extra 2(M − 2) multiplications and M − 2 additions, that aremuch fewer than the R-DST. Table II shows the number ofextra multiplications and additions. To clarify this point, weevaluated the implementation time in the case of M = 8
as an example using MATLAB with Intel Core i9-11900K
CPU; the R-FST saved approximately 0.126 seconds in a 2-D
transformation of 512×512 signals compared with the R-DST.

C. Remark on Other Types
While this letter focuses on the type-II R-FST designs due
to limitation of space, we can easily extend the method to theother types of the DCT/DST in accordance with Section III-A
because the original DCT/DST matrices are orthogonal. We
e like S
b in (26). Moreover,can do so for type I by using R
the type-II DCT already satisfies the regularity condition, as
described in the Introduction.
IV. E XPERIMENTS
Figures 2 and 3 show the frequency responses of the DST
and R-FST and the Barbara (subband modes) transformed
with the 8×8 DST and R-FST. The results of the conventional

V. C ONCLUSION
This letter introduced the regularity-constrained fast sinetransform (R-FST). The R-FST is obtained by just appending
an M × M regularity constraint matrix, which is constructedfrom only M/2 − 1 rotation matrices with the angles derived
from the output of the M × M DST for the DC signal, as
a postprocessing of the original DST. An image processing
example showed that it has fine frequency selectivity with noDC leakage and higher coding gain than that of the original
DST. Also, the R-FST achieved more time savings than the
R-DST because of fewer extra operations.

5

R EFERENCES[1] V. Britanak, P. Yip, and K. R. Rao, Discrete Cosine and Sine Transforms: General Properties, Fast Algorithms and Integer Approximations,
Academic Press, 2014.
[2] G. K. Wallace, “The JPEG still picture compression standard,” IEEE
Trans. Consum. Electr., vol. 38, no. 1, pp. xviii–xxxiv, Feb. 1992.
[3] T. Wiegand, G. J. Sullivan, G. Bjøntegaard, and A. Luthra, “Overview
of the H.264/AVC video coding standard,” IEEE Trans. Circuits Syst.
Video Technol., vol. 13, no. 7, pp. 560–576, Jul. 2003.
[4] G. J. Sullivan, J.-R. Ohm, W.-J. Han, and T. Wiegand, “Overview of thehigh efficiency video coding (HEVC) standard,” IEEE Trans. Circuits
Syst. Video Technol., vol. 22, no. 12, pp. 1649–1668, Dec. 2012.
[5] B. Bross, Y.-K. Wang, Y. Ye, S. Liu, J. Chen, G. J. Sullivan, and J.-R.
Ohm, “Overview of the versatile video coding (VVC) standard and its
applications,” IEEE Trans. Circuits Syst. Video Technol., vol. 31, no.
10, pp. 3736–3764, Oct. 2021.
[6] M. Aharon, M. Elad, and A. Bruckstein, “K-SVD: An algorithm fordesigning overcomplete dictionaries for sparse representation,” IEEE
Trans. Signal Process., vol. 54, no. 11, pp. 4311–4322, Oct. 2006.
[7] A. Danielyan, V. Katkovnik, and K. Egiazarian, “BM3D frames andvariational image deblurring,” IEEE Trans. Image Process., vol. 21, no.
4, pp. 1715–1728, Apr. 2012.
[8] Y. Wang, C. Xu, S. You, C. Xu, and D. Tao, “DCT regularized extreme
visual recovery,” IEEE Trans. Image Process., vol. 26, no. 7, pp. 3360–
3371, Mar. 2017.
[9] T. D. Tran, “The BinDCT: Fast multiplierless approximation of the
DCT,” IEEE Signal Process. Lett., vol. 7, no. 6, pp. 141–144, Jun.
2000.
[10] J. Liang and T. D. Tran, “Fast multiplierless approximations of the DCT
with the lifting scheme,” IEEE Trans. Signal Process., vol. 49, no. 12,
pp. 3032–3044, Dec. 2001.
[11] P. Kumar, S. Y. Park, B. K. Mohanty, K. S. Lim, and C. Yeo, “Efficient
integer DCT architectures for HEVC,” IEEE Trans. Circuits Syst. Video
Technol., vol. 24, no. 1, pp. 168–178, Jan. 2014.
[12] T. Suzuki, S. Kyochi, Y. Tanaka, and M. Ikehara, “Multiplierless
lifting-based fast X transforms derived from fast Hartley transform
factorization,” Springer Mult. Syst. Signal Process., vol. 29, no. 1, pp.
99–118, Jan. 2018.
[13] X. Zhao, S.-H. Kim, Y. Zhao, H. E. Egilmez, M. Koo, S. Liu, J. Lainema,
and M. Karczewicz, “Transform coding in the VVC standard,” IEEE
Trans. Circuits Syst. Video Technol., vol. 31, no. 10, pp. 3878–3890,
Oct. 2021.
[14] W. H. Chen, C. H. Smith, and S. C. Fralick, “A fast computationalalgorithm for the discrete cosine transform,” IEEE Trans. Commun.,vol. 25, no. 9, pp. 1004–1009, Sep. 1977.
[15] Z. Wang, “Fast algorithms for the discrete W transform and forthe discrete Fourier transform,” IEEE Trans. Acoust. Speech. Signal
Process., vol. ASSP-32, no. 4, pp. 803–816, Apr. 1984.
[16] B. Lee, “A new algorithm to compute the discrete cosine transform,”
IEEE Trans. Acoust. Speech. Signal Process., vol. 32, no. 6, pp. 1243–
1245, Dec. 1984.
[17] S. Kyochi, T. Suzuki, and Y. Tanaka, “Directional analytic discretecosine frames,” IEEE Access, vol. 10, pp. 20964–20977, Feb. 2022.
[18] P. P. Vaidyanathan, Multirate Systems and Filter Banks, Englewood
Cliffs, NJ: Prentice Hall, 1992.