- Email: [email protected]
DST-IV computation via the MDCT
Signal Processing 83 (2003) 1803 – 1813 www.elsevier.com/locate/sigpro
Short communication
The fast DCT-IV/DST-IV computation via the MDCT Vladimir Britanak Institute of Informatics, Slovak Academy of Sciences, Dubravska cesta 9, 842 37 Bratislava, Slovak Republic Received 19 April 2002; received in revised form 6 February 2003
Abstract The discrete cosine transform of type IV (DCT-IV) and corresponding discrete sine transform of type IV (DST-IV) have played key role in the e1cient implementation of orthogonal lapped transforms and perfect reconstruction cosine-modulated 2lter banks such as the oddly stacked modi2ed discrete cosine transform (MDCT) or equivalently, the modulated lapped transform (MLT). However, the DCT-IV and DST-IV of double sizes are related to two variants of 2lter banks de2ned by Dolby Labs AC-3 digital audio compression algorithm. Since these two variants of 2lter banks are e1ciently computed by recently proposed new fast algorithm for the oddly stacked MDCT (Signal Processing 82 (2002) 433), it is shown that the e1cient DCT-IV and DST-IV computation can be realized via the MDCT of double size. The careful analysis of regular structure of the new fast MDCT algorithm allows to extract a new DCT-IV/DST-IV computational structure and to suggest a new sparse matrix factorization of the DCT-IV matrix. Finally, the new DCT-IV/DST-IV computational structure provides an alternative e1cient implementation of the forward and inverse MDCT in layer III of MPEG (MP3) audio coding. ? 2003 Elsevier Science B.V. All rights reserved. Keywords: Discrete cosine transform; Discrete sine transform; Modi2ed discrete cosine transform; Modi2ed discrete sine transform; Modulated lapped transform; MPEG audio coding
1. Introduction The discrete cosine transform of type IV (DCT-IV) and corresponding discrete sine transform of type IV (DST-IV) have been introduced into digital signal processing by Jain [13] as alternative transforms for spectral analysis. Jain has shown that basis functions of the DCT-IV and DST-IV, originally referred to the even discrete cosine-2 and even discrete sine-3 transforms, respectively, are eigenvectors of the parametrized symmetric tridiagonal Jacobi matrix. Subsequently, the DCT-IV and DST-IV have become the members of the complete set of DCTs and DSTs known as the class of discrete trigonometric transforms [40]. The DCT-IV and DST-IV are intrinsically related to the corresponding type of the generalized discrete Fourier transform (GDFT), called the odd-time odd-frequency DFT (O2 DFT) [1], and generalized discrete Hartley transform (GDHT) [39] or equivalently generalized discrete W transform (GDWT) [40] for real data sequences.
E-mail address: [email protected] (V. Britanak). 0165-1684/03/$ - see front matter ? 2003 Elsevier Science B.V. All rights reserved. doi:10.1016/S0165-1684(03)00109-9
1804
V. Britanak / Signal Processing 83 (2003) 1803 – 1813
Because of frequency shift of the DCT-IV and DST-IV basis functions, the DCT-IV and DST-IV are not useful for signal coding, and they have found applications in spectrum estimation and adaptive 2ltering. However, both the DCT-IV and DST-IV have played key role in the e1cient implementation of orthogonal lapped transforms [18,22] and perfect reconstruction cosine-modulated 2lter banks [7,12,14,17–25,28,32,41], such as the oddly stacked modi2ed discrete cosine transform (MDCT) [7,12,17] and corresponding modi2ed discrete sine transform (MDST) [25], the modulated lapped transform (MLT) [14,18,19,22,24,28,32], extended lapped transform (ELT) [20–22], modulated lapped biorthogonal transform (MLBT) [23,24], and recently introduced the modulated complex lapped transform (MCLT) [25], and nonuniform MCLT (NMCLT) [41]. For subband decomposition of a signal are also used the so-called real-valued polyphase 2lter banks, where all channels are shifted versions of the same prototype low-pass 2lter, and these 2lter banks are derived directly from the DCT-IV and DST-IV transform kernels [12]. As an example, Dolby Labs AC-3 digital audio compression algorithm [9] besides the MDCT for oddly stacked system [26] has adopted additional two variants of such polyphase 2lter banks, and they are referred to the 2rst and second short MDCTs. Over two decades many fast algorithms for the e1cient DCT-IV and DST-IV computation have been developed. In general, the DCT-IV or DST-IV can be computed by two approaches: indirectly or directly. In the indirect approach other unitary or orthogonal transforms are used, such as discrete Fourier transform (DFT) and its fast implementation, complex-valued FFT [12,13,22] or real-valued FFT [10,11], discrete Hartley transform (DHT) [16] or DCT of type II (DCT-II) [5,6,15,28,29]. Actually, the DCT-IV and DCT-II have an intimate relationship. Since N -point DCT-II can be decomposed into N=2-point DCT-II (even-indexed coe1cients) and N=2-point DCT-IV (odd-indexed coe1cients), this fact has led to widely accepted conclusion [14]: fast N -point DCT-IV algorithms can be derived indirectly from existing fast 2N -point DCT-II algorithms [30]. On the other hand, the DCT-IV of size N may be converted by a simple method to the DCT-II of the same size including additional multiplications and additions [5,6,15]. The DCT-IV and DST-IV are intrinsically related to O2 DFT of real data sequences. The fast algorithms based on O2 DFT [3,12,27,29] have a speci2c place in the class of indirect algorithms, because by those algorithms the simultaneous (on-line) DCT-IV/DST-IV computation is achieved. Direct algorithms for the e1cient DCT-IV and DST-IV computation are based on sparse matrix factorization of the DCT-IV matrix [5,8,31,33–38,42,43]. The 2rst attempt to derive a direct DCT-IV matrix factorization was presented by Chen et al. [8]. Although the DCT-II matrix was decomposed into DCT-II and DCT-IV matrices of half sizes, they also discussed and proposed the factorization of the DCT-IV matrix, which at that time, it was not recursive and very e1cient. Wang reconsidered this result [35] and he proposed a decomposition of the DCT-IV matrix consisting of a product od computationally simple sine/cosine butterGy and binary matrices [36]. Subsequently, Suehiro and Hatori [33,34] improved Wang’s decomposition in computational e1ciency identifying its more e1cient basic computational part. Later, Wang [37] developed the DCT-IV matrix factorization based on the DCT-III (being inverse of DCT-II) and corresponding DST-III of half sizes. Since there exists a simple relation between DCT-IV and DST-IV matrices [36], the fast DST-IV computation can be realized by existing DCT-IV algorithms. In this paper, a new fast algorithm for the DCT-IV and DST-IV computation is presented as a direct consequence of the recently proposed new fast MDCT algorithm [4]. The DCT-IV and DST-IV of double sizes are related to two variants of 2lter banks de2ned by Dolby Labs AC-3 digital audio compression algorithm. Since these two variants of 2lter banks are e1ciently computed by the new fast MDCT algorithm, it is shown that the e1cient DCT-IV and DST-IV computation can be realized via the MDCT of double size. The careful analysis of regular structure of the new fast MDCT algorithm allows to extract a new DCT-IV/DST-IV computational structure representing the uni2ed fast DCT-IV/DST-IV computation and their inverses for any N being an even integer. Subsequently, a new sparse matrix factorization of the DCT-IV matrix is suggested as the direct consequence of the sparse matrix factorization derived for the MDCT matrix. Finally, the new DCT-IV/DST-IV computational structure provides an alternative e1cient implementation of the forward and inverse MDCT in layer III of MPEG (MP3) audio coding.
V. Britanak / Signal Processing 83 (2003) 1803 – 1813
1805
2. Denitions and notations Let {x n }, n=0; 1; : : : ; N −1 represents an input data sequence. The DCT-IV and its inverse are, respectively, de2ned as [13,40] ckIV =
2 x n cos (2n + 1)(2k + 1) ; N 4N N −1
k = 0; 1; : : : ; N − 1;
(1)
n=0
xn =
N −1
ckIV cos
k=0
(2n + 1)(2k + 1) ; 4N
n = 0; 1; : : : ; N − 1:
(2)
The corresponding DST-IV and its inverse are, respectively, de2ned as [13,40] skIV =
2 x n sin (2n + 1)(2k + 1) ; N 4N N −1
k = 0; 1; : : : ; N − 1;
(3)
n=0
xn =
N −1
skIV sin
k=0
(2n + 1)(2k + 1) ; 4N
n = 0; 1; : : : ; N − 1:
(4)
The DCT-IV and DST-IV are involutory, i.e., they are orthogonal and symmetric. It means that fast algorithms for the forward and inverse DCT-IV/DST-IV are the same except for a normalization factor. AC-3 digital audio compression algorithm developed originally by Dolby Labs de2nes the forward and inverse MDCT, respectively, as [9] ck =
2 (2n + 1)(2k + 1) + (2k + 1)(1 + ) ; x n cos N 2N 4 N −1 n=0
N=2−1
xˆn =
k=0
where
−1 = 0 +1
ck cos
(2n + 1)(2k + 1) + (2k + 1)(1 + ) ; 2N 4
k = 0; 1; : : : ;
N − 1; 2
n = 0; 1; : : : ; N − 1;
(5)
(6)
for the 2rst short transform (N = 256); for the long transform (N = 512); for the second short transform (N = 256);
{x n } in (5) represents the windowed input data sequence, whereas {xˆn } in (6) represents recovered timedomain aliased data sequence. Hence, besides the long MDCT being oddly stacked MDCT [26], AC-3 de2nes additional two variants of cosine-modulated 2lter banks called the 2rst and second short MDCTs. Denote them as {ck } and {sk }. It is clear that they are related to the DCT-IV and DST-IV of double sizes, respectively. We recall that the new fast MDCT algorithm [4] can be also adopted for the alternate computation of short MDCTs in AC-3, and therefore, it is assumed that a reader is familiar with the new fast MDCT algorithm and sparse matrix factorization of the MDCT matrix. The MDCT-modi2ed generalized signal Gow graph for the alternate computation of short MDCTs is shown in Fig. 1 for N = 16. Compared to the original one [4], the N=4-point discrete sine transform of type II (DST-II) in its lower half has been replaced by N=4-point DCT-II with proper preceding sign changes.
1806
V. Britanak / Signal Processing 83 (2003) 1803 – 1813
Fig. 1. MDCT-modi2ed generalized signal Gow graph adopted for the computation of short MDCTs in AC-3 for N = 16.
The speci2c de2nition of the MDCT in AC-3, the relation between two variants of 2lter banks and DCT-IV/DST-IV, regular structure of the new MDCT algorithm, and the sparse matrix factorization of the MDCT matrix provide the possibility to derive a new fast algorithm for the alternate DCT-IV/DST-IV computation and their inverses via the MDCT.
3. New fast algorithm for the DCT-IV/DST-IV computation via MDCT The input data sequence {x n }, n = 0; 1; : : : ; N − 1 is assumed to be not windowed necessarily. At 2rst, let us extend the data sequence {x n } by zero padding to the length M = 2N as x n ; n = 0; 1; : : : ; N − 1; yn = (7) 0; n = N; N + 1; : : : ; 2N − 1: Then, the DCT-IV and DST-IV of {yn } are, respectively, de2ned as ckIV =
M −1 2 (2n + 1)(2k + 1) ; yn cos M 2M n=0
k = 0; 1; : : : ; M − 1;
(8)
V. Britanak / Signal Processing 83 (2003) 1803 – 1813
skIV =
M −1 2 yn sin (2n + 1)(2k + 1) ; M 2M
k = 0; 1; : : : ; M − 1
1807
(9)
n=0
and they are related to two variants of cosine-modulated 2lter banks in AC-3 with the length M =2N . Consider the new fast MDCT algorithm and its modi2ed generalized signal Gow graph for the computation of short MDCTs shown in Fig. 1. We can observe that after the second butterGy stage, the data sequence {yn } with the length M = 2N is decimated to the length N . Knowing a priori that the second half of extended data sequence is equal to zero, the 2rst and second butterGy stages can be implemented without any addition, and hence, without no increase of arithmetic operations. Furthermore, the DCT-IV given by (8) is related to the 2rst short MDCT in AC-3 by ckIV = ck ;
k = 0; 1; : : : ; N − 1;
(10)
whereas the DST-IV given by (9) is related to the second short MDCT in AC-3 by skIV = (−1)k+1 sk ;
k = 0; 1; : : : ; N − 1:
(11)
Thus, for a given N , the DCT-IV given by (8) is exactly the 2rst short MDCT, and the DST-IV given by (9) corresponds to the second short MDCT except for sign changes. We note that relations between two variants of 2lter banks and DCT-IV/DST-IV is quite similar to that of derived between the DCT-II and DCT-IV/DST-IV [28]. As a result, the new fast MDCT algorithm with minor modi2cations can be also used for the fast DCT-IV and DST-IV computation as follows: 1. Extend N -point input data sequence {x n } by zero padding to the length M = 2N according to (7). 2. For the DCT-IV of {yn } compute the 2rst short MDCT in AC-3. The DCT-IV coe1cients are given by (10). For the DST-IV of {yn } compute the second short MDCT in AC-3 and change signs of the DST-IV coe1cients according to (11). The inverse DCT-IV and DST-IV are realized by the new fast algorithm for inverse MDCT computation as follows: 1. Extend N -point DCT-IV/DST-IV sequence {ckIV }={skIV } by zero padding to the length M = 2N according to (7). 2. For the inverse DCT-IV of {ckIV } compute the 2rst short inverse MDCT in AC-3. For the inverse DST-IV of {skIV } change signs of the coe1cients according to (11) and compute the second short inverse MDCT in AC-3. Alternatively, since the DCT-IV and DST-IV matrices are involutory, the inverse DCT-IV and DST-IV computations can be realized only by the corresponding forward short MDCTs. However, in the case of inverse DST-IV computation there is minor modi2cation. Before the second short MDCT computation, the sign changed DST-IV coe1cients have to be reversed in order. The recovered time-domain data sequence will be also in reverse order. The exact number of arithmetic operations required for the DCT-IV/DST-IV computation and their inverses via the MDCT for data sequences with length N = 2n is given by (N=2)(n + 2) multiplications and (3N=2)n additions. The careful analysis of the MDCT modi2ed generalized signal Gow graph in Fig. 1 reveals that the new computational structure for the uni2ed DCT-IV/DST-IV computation and their inverses can be directly extracted from it. Let {x n }, n = 0; 1; : : : ; N − 1 be the input data sequence extended by zero padding to the length 2N . The analysis of the 2rst and second butterGy stages applied to the extended data sequence shows that
1808
V. Britanak / Signal Processing 83 (2003) 1803 – 1813
butterGy stages generate just a permutation of the original data sequence {x n }, n = 0; 1; : : : ; N − 1. For clarity, an example for N = 8 is shown x0 x1 x2 x3
x0 = x0 + x15 x1 = x1 + x14 x2 = x2 + x13 x3 = x3 + x12
x4 x5 x6 x7
x4 = x4 + x11 x5 = x5 + x10 x6 = x6 + x9 x7 = x7 + x8
x8 = 0 x9 = 0 x10 = 0 x11 = 0
x7 = x7 − x8 x6 = x6 − x9 x5 = x5 − x10 x4 = x4 − x11
x12 = 0 x13 = 0 x14 = 0 x15 = 0
x3 = x3 − x12 x2 = x2 − x13 x1 = x1 − x14 x0 = x0 − x15
DCT-IV xn = x7−n =0
DST-IV xn = x7−n =0
x3 − x4 x2 − x5 x1 − x6 x0 − x7
x3 x2 x1 x0
−x4 −x5 −x6 −x7
x0 − x7 x1 − x6 x2 − x5 x3 − x4
−x7 −x6 −x5 −x4
x0 x1 x2 x3
To be more speci2c, for the DCT-IV computation the following permutation is performed: yn = xN=2−1−n ;
yn+N=2 = −xN −1−n ;
n = 0; 1; : : : ;
N − 1; 2
(12)
N − 1: 2
(13)
and for the DST-IV computation the permutation yn = −x n+N=2 ;
yN −1−n = xN=2−1−n ;
n = 0; 1; : : : ;
From above example it can be seen that permuted data sequences for the DCT-IV and DST-IV computation are actually in reverse order of each other. Consequently, we can eliminate the 2rst and second butterGy stages and replace them by appropriate permutation of the original data sequence. Now for a given N we are able to extract from the MDCT-modi2ed generalized signal Gow graph the new DCT-IV/DST-IV computational structure which will depend on two identical DCT-IIs of half sizes. It is shown in Fig. 2 for N = 16. The symbols in brackets correspond to the DST-IV computation. With respect to Fig. 1 and Eqs. (10) and (11), the corresponding output DCT-IV and DST-IV coe1cients are, respectively, obtained as IV = z2k ; c2k
IV c2k+1 = −zN −2k−2 ;
k = 0; 1; : : : ;
N −1 2
(14)
and IV = −z2k ; s2k
IV s2k+1 = −zN −2k−2 ;
k = 0; 1; : : : ;
N − 1: 2
(15)
The new DCT-IV/DST-IV computational structure represents the uni2ed DCT-IV/DST-IV computation and their inverses for any N being an even integer. In speci2c case of N = 2n , the exact number of arithmetic operations required for the DCT-IV/DST-IV computation is (N=2)(n+2) multiplications and (3N=2)n additions. This is the lowest achievable complexity for the DCT-IV [22].
V. Britanak / Signal Processing 83 (2003) 1803 – 1813
1809
Fig. 2. The new DCT-IV/DST-IV computational structure for N = 16.
4. New factorization of the DCT-IV matrix In fact, the new DCT-IV/DST-IV computational structure de2nes also a new sparse matrix factorization of the DCT-IV matrix. Let CNIV be the orthogonal matrix of order N . Then CNIV can be decomposed into the following sparse matrix product: CNIV
= QN
1 IN=2−1 −1 −JN=2−1
II −JN=2−1 CN=2 JN=2 0 −IN=2−1
0 II DN=2 CN=2
GN
JN=2
0
0
−JN=2
;
(16)
II where the matrix product CN=2 JN=2 denotes N=2 × N=2 DCT-II matrix with reversed order of its columns. The II matrix product CN=2 DN=2 denotes the N=2 × N=2 DCT-II matrix with sign changed of its odd-indexed rows, whereby DN=2 = diag{1; −1; 1; : : : ; (−1)N=2+1 } is the diagonal matrix. GN is the rotation matrix (Given’s plane
1810
V. Britanak / Signal Processing 83 (2003) 1803 – 1813
rotations) de2ned by (N − 1) cos 4N GN = (N − 1) sin 4N
cos
(N − 3) 4N
−sin ·
· cos 4N sin 4N
−sin 4N cos 4N
· sin
(N − 3) 4N
(N − 3) 4N
· cos
(N − 3) 4N
(N − 1) −sin 4N ; (N − 1) cos 4N (17)
and QN is the permutation matrix for reordering the coe1cients according to (14). The new factorization of the DCT-IV matrix is the direct consequence of the sparse matrix factorization derived for the MDCT matrix [4]. Finally, considering the MLT and its fast implementation based on DCT-IV [22], the alternative e1cient implementation of the forward and inverse MDCT computation in MP3 audio coding can be derived using the new DCT-IV/DST-IV computational structure. After the permutation applied to the input data sequence {x n } [22] yn+N=4 = x n − xN=2−1−n ;
yN=4−1−n = −x n+N=2 − xN −1−n ;
n = 0; 1; : : : ;
N − 1; 4
(18)
the MDCT is converted to the DCT-IV of {yn }, n = 0; 1; : : : ; N=2 − 1. Signal Gow graphs for the forward and inverse MDCT in MP3 are shown for N = 12 and 36 in Figs. 3(a) and (b), respectively. In order to minimize total arithmetic complexity we can use optimized odd-length DCT-II modules [2]. Three-point DCT-II module requires one multiplication, four additions and one shift, whereas 9-point DCT-II module requires eight multiplications, 34 additions and two shifts. Then arithmetic complexity of the forward MDCT computation in MP3 for N = 12 is 11 multiplications, 27 additions and two shifts, whereas for N = 36 is 43 multiplications, 129 additions and four shifts. The inverse MDCT requires exactly N=2 less additions than that of the forward MDCT. 5. Conclusions Although, both the DCT-IV and DST-IV play key role in the e1cient implementation of orthogonal lapped transforms and perfect reconstruction cosine-modulated 2lter banks such as the oddly stacked MDCT or equivalently MLT, it is shown that the DCT-IV and DST-IV can be e1ciently computed by the new fast MDCT algorithm for oddly stacked system [4], and hence via the MDCT. It follows from fact that the DCT-IV and DST-IV of double sizes are related to two variants of 2lter banks de2ned by Dolby Labs AC-3 digital audio compression algorithm, and these two variants of 2lter banks can be e1ciently realized by the new fast MDCT algorithm. Moreover, by careful analysis of regular structure of the new fast MDCT algorithm, the
V. Britanak / Signal Processing 83 (2003) 1803 – 1813
1811
(a)
(b) Fig. 3. Signal Gow graphs for the forward and inverse MDCT computation in MP3 (a) for N = 12 and (b) for N = 36.
new DCT-IV/DST-IV computational structure is extracted and new sparse matrix factorization of the DCT-IV matrix is suggested. Finally, the new DCT-IV/DST-IV computational structure provides the alternative e1cient implementation of the forward and inverse MDCT in MP3 audio coding. Acknowledgements The author would like to thank the anonymous reviewers for their helpful comments and suggestions, which signi2cantly improved the clarity of presentation in this paper. This work was partially supported by Slovak Scienti2c Grant Agency VEGA, project No. 2/1100/21.
1812
V. Britanak / Signal Processing 83 (2003) 1803 – 1813
References [1] G. Bonnerot, M. Bellanger, Odd-time odd-frequency discrete Fourier transform for symmetric real-valued series, Proc. IEEE 64 (March 1976) 392–393. [2] V. Britanak, The re2ned e1cient implementation of the MDCT in MP3 and comparison with other methods, Technical Report II SAS-2002-02, Institute of Informatics, Slovak Academy of Sciences, Bratislava, September 2002. [3] V. Britanak, K.R. Rao, The fast generalized discrete Fourier transforms: a uni2ed approach to the discrete sinusoidal transforms computation, Signal Process. 79 (December 1999) 135–150. [4] V. Britanak, K.R. Rao, A new fast algorithm for the uni2ed forward and inverse MDCT/MDST computation, Signal Process. 82 (March 2002) 433–459. [5] S.C. Chan, K.L. Ho, Direct methods for computing discrete sinusoidal transforms, IEE Proc.: Radar Signal Process. Part F 137 (December 1990) 433–442. [6] S.C. Chan, K.L. Ho, Fast algorithms for computing the discrete cosine transform, IEEE Trans. Circuits Syst. II: Analog Digital Signal Process. 39 (March 1992) 185–190. [7] D.Y. Chan, J.F. Yang, S.Y. Chen, Regular implementation algorithms of time domain aliasing cancellation, IEE Proc. Vision Image Signal Process. 143 (December 1996) 387–392. [8] W.H. Chen, C.H. Smith, S.C. Fralick, A fast computational algorithm for the discrete cosine transform, IEEE Trans. Commun. COM-25 (September 1977) 1004–1009. [9] Digital Audio Compression (AC-3) Standard, Document of Advanced Television Systems Committee (ATSC), Audio Specialist Group T3/S7, December 20, 1995. [10] O.K. Ersoy, Fourier-related transforms, fast algorithms and applications, Prentice-Hall, Upper Saddle River, NJ, 1997, pp. 388– 452 (Chapter 9). [11] O.K. Ersoy, N.C. Hu, A uni2ed approach to the fast computation of all discrete trigonometric transforms, Proceedings of the IEEE ICASSP’87, Dallas, TX, April 1987, pp. 1843–1846. [12] R. Gluth, Regular FFT-related transform kernels for DCT/DST-based polyphase 2lter banks, Proceedings of the IEEE ICASSP’91, Toronto, Canada, May 1991, pp. 2205 –2208. [13] A.K. Jain, A sinusoidal family of unitary transforms, IEEE Trans. Pattern Anal. Machine Intell. PAMI-1 (October 1979) 356–365. [14] Z. Jiang, A.N. Willson, Fast cosine-modulated 2lter banks, Proceedings of the International Conference on Digital Signal Processing (DSP’95), Limassol, Cyprus, June 1995, pp. 188–193. [15] C.W. Kok, Fast algorithm for computing discrete cosine transform, IEEE Trans. Signal Process. 45 (March 1997) 757–760. [16] P.S. Kumar, K.M.M. Prabhu, A set of new algorithms for DCTs and DSTs, Int. J. of Electron. 76 (April 1994) 553–586. [17] Ch.M. Liu, W.Ch. Lee, A uni2ed fast algorithm for cosine modulated 2lter banks in current audio coding standards, J. Audio Eng. Soc. 47 (12) (December 1999) 1061–1075. [18] H.S. Malvar, Lapped transforms for e1cient transform/subband coding, IEEE Trans. Acoust. Speech Signal Process. 38 (June 1990) 969–978. [19] H.S. Malvar, Fast algorithm for modulated lapped transform, Electron. Lett. 27 (April 1991) 775–776. [20] H.S. Malvar, Extended lapped transforms: fast algorithms and applications, Proceedings of the IEEE ICASSP’91, Toronto, Canada, May 1991, pp. 1797–1800. [21] H.S. Malvar, Extended lapped transforms: properties, applications, and fast algorithms, IEEE Trans. Signal Process. 40 (November 1992) 2703–2714. [22] H.S. Malvar, Signal Processing with Lapped Transforms, Artech House, Norwood, MA, 1992. [23] H.S. Malvar, Biorthogonal and nonuniform lapped transforms for transform coding with reduced blocking and ringing artifacts, IEEE Trans. on Signal Process. 46 (April 1998) 1043–1053. [24] H.S. Malvar, Fast algorithms for orthogonal and biorthogonal modulated lapped transforms, Proceedings of the IEEE Symposium on Advances in Digital Filtering and Signal Processing, Victoria, Canada, June 1998, pp. 159 –163. [25] H.S. Malvar, A modulated complex lapped transform and its application to audio processing, Proceedings of the IEEE ICASSP’99, Phoenix, AR, May 1999, pp. 1421–1424. [26] J.P. Princen, A.W. Johnson, A.B. Bradley, Subband/transform coding using 2lter bank designs based on time domain aliasing cancellation, Proceedings of the IEEE ICASSP’87, Dallas, TX, April 1987, pp. 2161–2164. [27] N. Rama Murthy, M.N.S. Swamy, On the algorithms for the computation of even discrete cosine transform-2 (EDCT-2) of real sequences, IEEE Trans. on Circuits Syst. 37 (May 1990) 625–627. [28] N. Rama Murthy, M.N.S. Swamy, On a novel decomposition of the DCT and its application, IEEE Trans. Signal Process. 41 (January 1993) 480–485. [29] N. Rama Murthy, M.N.S. Swamy, On the on-line computation of DCT-IV and DST-IV transforms, IEEE Trans. Signal Process. 43 (May 1995) 1249–1251. [30] K.R. Rao, P. Yip, Discrete cosine transform: algorithms, advantages applications, Academic Press, Boston, 1990. [31] K.R. Rao, P.C. Yip (Eds.), The transform and data compression handbook, CRC Press, Boca Raton, FL, 2001, pp. 117–195 (Chapter 4).
V. Britanak / Signal Processing 83 (2003) 1803 – 1813
1813
[32] S. Shlien, The modulated lapped transform, its varying forms, and application to audio coding, IEEE Trans. Speech Audio Process. 5 (July 1997) 359–366. [33] N. Suehiro, M. Hatori, Fast algorithms for the discrete Fourier transform and others transforms, Proceedings of the IEEE ICASSP’86, Tokyo, Japan, April 1986, pp. 217–220. [34] N. Suehiro, M. Hatori, Fast algorithms for the discrete Fourier transform and others transforms, IEEE Trans. Acoust. Speech Signal Process. ASSP-34 (June 1986) 642–644. [35] Z. Wang, Reconsideration of a fast computational algorithm for the discrete cosine transform, IEEE Trans. Commun. COM-31 (January 1983) 121–123. [36] Z. Wang, Fast algorithms for the discrete W transform and for the discrete Fourier transform, IEEE Trans. Acoust. Speech Signal Process. ASSP-32 (August 1984) 803–816. [37] Z. Wang, On computing the discrete Fourier and cosine transforms, IEEE Trans. Acoust. Speech Signal Process. ASSP-33 (October 1985) 1341–1344. [38] Z. Wang, Fast discrete sine transform algorithms, Signal Process. 19 (February 1990) 91–102. [39] Z. Wang, Comments on generalized discrete Hartley transforms, IEEE Trans. Signal Process. 43 (July 1995) 1711–1712. [40] Z. Wang, B.R. Hunt, The discrete W transform, Appl. Math. Comput. 16 (1985) 19–48. [41] Z. Xiong, H.S. Malvar, A nonuniform modulated complex lapped transform, IEEE Signal Process. Lett. 8 (September 2001) 257–260. [42] P. Yip, K.R. Rao, Fast decimation-in-time algorithms for a family of discrete sine and cosine transforms, Circuits Syst. Signal Process. 3 (4) (1984) 387–408. [43] P. Yip, K.R. Rao, Fast decimation-in-frequency algorithms for a family of discrete sine and cosine transforms, Circuits Syst. Signal Process. 7 (1) (1988) 3–19.
Short communication
The fast DCT-IV/DST-IV computation via the MDCT Vladimir Britanak Institute of Informatics, Slovak Academy of Sciences, Dubravska cesta 9, 842 37 Bratislava, Slovak Republic Received 19 April 2002; received in revised form 6 February 2003
Abstract The discrete cosine transform of type IV (DCT-IV) and corresponding discrete sine transform of type IV (DST-IV) have played key role in the e1cient implementation of orthogonal lapped transforms and perfect reconstruction cosine-modulated 2lter banks such as the oddly stacked modi2ed discrete cosine transform (MDCT) or equivalently, the modulated lapped transform (MLT). However, the DCT-IV and DST-IV of double sizes are related to two variants of 2lter banks de2ned by Dolby Labs AC-3 digital audio compression algorithm. Since these two variants of 2lter banks are e1ciently computed by recently proposed new fast algorithm for the oddly stacked MDCT (Signal Processing 82 (2002) 433), it is shown that the e1cient DCT-IV and DST-IV computation can be realized via the MDCT of double size. The careful analysis of regular structure of the new fast MDCT algorithm allows to extract a new DCT-IV/DST-IV computational structure and to suggest a new sparse matrix factorization of the DCT-IV matrix. Finally, the new DCT-IV/DST-IV computational structure provides an alternative e1cient implementation of the forward and inverse MDCT in layer III of MPEG (MP3) audio coding. ? 2003 Elsevier Science B.V. All rights reserved. Keywords: Discrete cosine transform; Discrete sine transform; Modi2ed discrete cosine transform; Modi2ed discrete sine transform; Modulated lapped transform; MPEG audio coding
1. Introduction The discrete cosine transform of type IV (DCT-IV) and corresponding discrete sine transform of type IV (DST-IV) have been introduced into digital signal processing by Jain [13] as alternative transforms for spectral analysis. Jain has shown that basis functions of the DCT-IV and DST-IV, originally referred to the even discrete cosine-2 and even discrete sine-3 transforms, respectively, are eigenvectors of the parametrized symmetric tridiagonal Jacobi matrix. Subsequently, the DCT-IV and DST-IV have become the members of the complete set of DCTs and DSTs known as the class of discrete trigonometric transforms [40]. The DCT-IV and DST-IV are intrinsically related to the corresponding type of the generalized discrete Fourier transform (GDFT), called the odd-time odd-frequency DFT (O2 DFT) [1], and generalized discrete Hartley transform (GDHT) [39] or equivalently generalized discrete W transform (GDWT) [40] for real data sequences.
E-mail address: [email protected] (V. Britanak). 0165-1684/03/$ - see front matter ? 2003 Elsevier Science B.V. All rights reserved. doi:10.1016/S0165-1684(03)00109-9
1804
V. Britanak / Signal Processing 83 (2003) 1803 – 1813
Because of frequency shift of the DCT-IV and DST-IV basis functions, the DCT-IV and DST-IV are not useful for signal coding, and they have found applications in spectrum estimation and adaptive 2ltering. However, both the DCT-IV and DST-IV have played key role in the e1cient implementation of orthogonal lapped transforms [18,22] and perfect reconstruction cosine-modulated 2lter banks [7,12,14,17–25,28,32,41], such as the oddly stacked modi2ed discrete cosine transform (MDCT) [7,12,17] and corresponding modi2ed discrete sine transform (MDST) [25], the modulated lapped transform (MLT) [14,18,19,22,24,28,32], extended lapped transform (ELT) [20–22], modulated lapped biorthogonal transform (MLBT) [23,24], and recently introduced the modulated complex lapped transform (MCLT) [25], and nonuniform MCLT (NMCLT) [41]. For subband decomposition of a signal are also used the so-called real-valued polyphase 2lter banks, where all channels are shifted versions of the same prototype low-pass 2lter, and these 2lter banks are derived directly from the DCT-IV and DST-IV transform kernels [12]. As an example, Dolby Labs AC-3 digital audio compression algorithm [9] besides the MDCT for oddly stacked system [26] has adopted additional two variants of such polyphase 2lter banks, and they are referred to the 2rst and second short MDCTs. Over two decades many fast algorithms for the e1cient DCT-IV and DST-IV computation have been developed. In general, the DCT-IV or DST-IV can be computed by two approaches: indirectly or directly. In the indirect approach other unitary or orthogonal transforms are used, such as discrete Fourier transform (DFT) and its fast implementation, complex-valued FFT [12,13,22] or real-valued FFT [10,11], discrete Hartley transform (DHT) [16] or DCT of type II (DCT-II) [5,6,15,28,29]. Actually, the DCT-IV and DCT-II have an intimate relationship. Since N -point DCT-II can be decomposed into N=2-point DCT-II (even-indexed coe1cients) and N=2-point DCT-IV (odd-indexed coe1cients), this fact has led to widely accepted conclusion [14]: fast N -point DCT-IV algorithms can be derived indirectly from existing fast 2N -point DCT-II algorithms [30]. On the other hand, the DCT-IV of size N may be converted by a simple method to the DCT-II of the same size including additional multiplications and additions [5,6,15]. The DCT-IV and DST-IV are intrinsically related to O2 DFT of real data sequences. The fast algorithms based on O2 DFT [3,12,27,29] have a speci2c place in the class of indirect algorithms, because by those algorithms the simultaneous (on-line) DCT-IV/DST-IV computation is achieved. Direct algorithms for the e1cient DCT-IV and DST-IV computation are based on sparse matrix factorization of the DCT-IV matrix [5,8,31,33–38,42,43]. The 2rst attempt to derive a direct DCT-IV matrix factorization was presented by Chen et al. [8]. Although the DCT-II matrix was decomposed into DCT-II and DCT-IV matrices of half sizes, they also discussed and proposed the factorization of the DCT-IV matrix, which at that time, it was not recursive and very e1cient. Wang reconsidered this result [35] and he proposed a decomposition of the DCT-IV matrix consisting of a product od computationally simple sine/cosine butterGy and binary matrices [36]. Subsequently, Suehiro and Hatori [33,34] improved Wang’s decomposition in computational e1ciency identifying its more e1cient basic computational part. Later, Wang [37] developed the DCT-IV matrix factorization based on the DCT-III (being inverse of DCT-II) and corresponding DST-III of half sizes. Since there exists a simple relation between DCT-IV and DST-IV matrices [36], the fast DST-IV computation can be realized by existing DCT-IV algorithms. In this paper, a new fast algorithm for the DCT-IV and DST-IV computation is presented as a direct consequence of the recently proposed new fast MDCT algorithm [4]. The DCT-IV and DST-IV of double sizes are related to two variants of 2lter banks de2ned by Dolby Labs AC-3 digital audio compression algorithm. Since these two variants of 2lter banks are e1ciently computed by the new fast MDCT algorithm, it is shown that the e1cient DCT-IV and DST-IV computation can be realized via the MDCT of double size. The careful analysis of regular structure of the new fast MDCT algorithm allows to extract a new DCT-IV/DST-IV computational structure representing the uni2ed fast DCT-IV/DST-IV computation and their inverses for any N being an even integer. Subsequently, a new sparse matrix factorization of the DCT-IV matrix is suggested as the direct consequence of the sparse matrix factorization derived for the MDCT matrix. Finally, the new DCT-IV/DST-IV computational structure provides an alternative e1cient implementation of the forward and inverse MDCT in layer III of MPEG (MP3) audio coding.
V. Britanak / Signal Processing 83 (2003) 1803 – 1813
1805
2. Denitions and notations Let {x n }, n=0; 1; : : : ; N −1 represents an input data sequence. The DCT-IV and its inverse are, respectively, de2ned as [13,40] ckIV =
2 x n cos (2n + 1)(2k + 1) ; N 4N N −1
k = 0; 1; : : : ; N − 1;
(1)
n=0
xn =
N −1
ckIV cos
k=0
(2n + 1)(2k + 1) ; 4N
n = 0; 1; : : : ; N − 1:
(2)
The corresponding DST-IV and its inverse are, respectively, de2ned as [13,40] skIV =
2 x n sin (2n + 1)(2k + 1) ; N 4N N −1
k = 0; 1; : : : ; N − 1;
(3)
n=0
xn =
N −1
skIV sin
k=0
(2n + 1)(2k + 1) ; 4N
n = 0; 1; : : : ; N − 1:
(4)
The DCT-IV and DST-IV are involutory, i.e., they are orthogonal and symmetric. It means that fast algorithms for the forward and inverse DCT-IV/DST-IV are the same except for a normalization factor. AC-3 digital audio compression algorithm developed originally by Dolby Labs de2nes the forward and inverse MDCT, respectively, as [9] ck =
2 (2n + 1)(2k + 1) + (2k + 1)(1 + ) ; x n cos N 2N 4 N −1 n=0
N=2−1
xˆn =
k=0
where
−1 = 0 +1
ck cos
(2n + 1)(2k + 1) + (2k + 1)(1 + ) ; 2N 4
k = 0; 1; : : : ;
N − 1; 2
n = 0; 1; : : : ; N − 1;
(5)
(6)
for the 2rst short transform (N = 256); for the long transform (N = 512); for the second short transform (N = 256);
{x n } in (5) represents the windowed input data sequence, whereas {xˆn } in (6) represents recovered timedomain aliased data sequence. Hence, besides the long MDCT being oddly stacked MDCT [26], AC-3 de2nes additional two variants of cosine-modulated 2lter banks called the 2rst and second short MDCTs. Denote them as {ck } and {sk }. It is clear that they are related to the DCT-IV and DST-IV of double sizes, respectively. We recall that the new fast MDCT algorithm [4] can be also adopted for the alternate computation of short MDCTs in AC-3, and therefore, it is assumed that a reader is familiar with the new fast MDCT algorithm and sparse matrix factorization of the MDCT matrix. The MDCT-modi2ed generalized signal Gow graph for the alternate computation of short MDCTs is shown in Fig. 1 for N = 16. Compared to the original one [4], the N=4-point discrete sine transform of type II (DST-II) in its lower half has been replaced by N=4-point DCT-II with proper preceding sign changes.
1806
V. Britanak / Signal Processing 83 (2003) 1803 – 1813
Fig. 1. MDCT-modi2ed generalized signal Gow graph adopted for the computation of short MDCTs in AC-3 for N = 16.
The speci2c de2nition of the MDCT in AC-3, the relation between two variants of 2lter banks and DCT-IV/DST-IV, regular structure of the new MDCT algorithm, and the sparse matrix factorization of the MDCT matrix provide the possibility to derive a new fast algorithm for the alternate DCT-IV/DST-IV computation and their inverses via the MDCT.
3. New fast algorithm for the DCT-IV/DST-IV computation via MDCT The input data sequence {x n }, n = 0; 1; : : : ; N − 1 is assumed to be not windowed necessarily. At 2rst, let us extend the data sequence {x n } by zero padding to the length M = 2N as x n ; n = 0; 1; : : : ; N − 1; yn = (7) 0; n = N; N + 1; : : : ; 2N − 1: Then, the DCT-IV and DST-IV of {yn } are, respectively, de2ned as ckIV =
M −1 2 (2n + 1)(2k + 1) ; yn cos M 2M n=0
k = 0; 1; : : : ; M − 1;
(8)
V. Britanak / Signal Processing 83 (2003) 1803 – 1813
skIV =
M −1 2 yn sin (2n + 1)(2k + 1) ; M 2M
k = 0; 1; : : : ; M − 1
1807
(9)
n=0
and they are related to two variants of cosine-modulated 2lter banks in AC-3 with the length M =2N . Consider the new fast MDCT algorithm and its modi2ed generalized signal Gow graph for the computation of short MDCTs shown in Fig. 1. We can observe that after the second butterGy stage, the data sequence {yn } with the length M = 2N is decimated to the length N . Knowing a priori that the second half of extended data sequence is equal to zero, the 2rst and second butterGy stages can be implemented without any addition, and hence, without no increase of arithmetic operations. Furthermore, the DCT-IV given by (8) is related to the 2rst short MDCT in AC-3 by ckIV = ck ;
k = 0; 1; : : : ; N − 1;
(10)
whereas the DST-IV given by (9) is related to the second short MDCT in AC-3 by skIV = (−1)k+1 sk ;
k = 0; 1; : : : ; N − 1:
(11)
Thus, for a given N , the DCT-IV given by (8) is exactly the 2rst short MDCT, and the DST-IV given by (9) corresponds to the second short MDCT except for sign changes. We note that relations between two variants of 2lter banks and DCT-IV/DST-IV is quite similar to that of derived between the DCT-II and DCT-IV/DST-IV [28]. As a result, the new fast MDCT algorithm with minor modi2cations can be also used for the fast DCT-IV and DST-IV computation as follows: 1. Extend N -point input data sequence {x n } by zero padding to the length M = 2N according to (7). 2. For the DCT-IV of {yn } compute the 2rst short MDCT in AC-3. The DCT-IV coe1cients are given by (10). For the DST-IV of {yn } compute the second short MDCT in AC-3 and change signs of the DST-IV coe1cients according to (11). The inverse DCT-IV and DST-IV are realized by the new fast algorithm for inverse MDCT computation as follows: 1. Extend N -point DCT-IV/DST-IV sequence {ckIV }={skIV } by zero padding to the length M = 2N according to (7). 2. For the inverse DCT-IV of {ckIV } compute the 2rst short inverse MDCT in AC-3. For the inverse DST-IV of {skIV } change signs of the coe1cients according to (11) and compute the second short inverse MDCT in AC-3. Alternatively, since the DCT-IV and DST-IV matrices are involutory, the inverse DCT-IV and DST-IV computations can be realized only by the corresponding forward short MDCTs. However, in the case of inverse DST-IV computation there is minor modi2cation. Before the second short MDCT computation, the sign changed DST-IV coe1cients have to be reversed in order. The recovered time-domain data sequence will be also in reverse order. The exact number of arithmetic operations required for the DCT-IV/DST-IV computation and their inverses via the MDCT for data sequences with length N = 2n is given by (N=2)(n + 2) multiplications and (3N=2)n additions. The careful analysis of the MDCT modi2ed generalized signal Gow graph in Fig. 1 reveals that the new computational structure for the uni2ed DCT-IV/DST-IV computation and their inverses can be directly extracted from it. Let {x n }, n = 0; 1; : : : ; N − 1 be the input data sequence extended by zero padding to the length 2N . The analysis of the 2rst and second butterGy stages applied to the extended data sequence shows that
1808
V. Britanak / Signal Processing 83 (2003) 1803 – 1813
butterGy stages generate just a permutation of the original data sequence {x n }, n = 0; 1; : : : ; N − 1. For clarity, an example for N = 8 is shown x0 x1 x2 x3
x0 = x0 + x15 x1 = x1 + x14 x2 = x2 + x13 x3 = x3 + x12
x4 x5 x6 x7
x4 = x4 + x11 x5 = x5 + x10 x6 = x6 + x9 x7 = x7 + x8
x8 = 0 x9 = 0 x10 = 0 x11 = 0
x7 = x7 − x8 x6 = x6 − x9 x5 = x5 − x10 x4 = x4 − x11
x12 = 0 x13 = 0 x14 = 0 x15 = 0
x3 = x3 − x12 x2 = x2 − x13 x1 = x1 − x14 x0 = x0 − x15
DCT-IV xn = x7−n =0
DST-IV xn = x7−n =0
x3 − x4 x2 − x5 x1 − x6 x0 − x7
x3 x2 x1 x0
−x4 −x5 −x6 −x7
x0 − x7 x1 − x6 x2 − x5 x3 − x4
−x7 −x6 −x5 −x4
x0 x1 x2 x3
To be more speci2c, for the DCT-IV computation the following permutation is performed: yn = xN=2−1−n ;
yn+N=2 = −xN −1−n ;
n = 0; 1; : : : ;
N − 1; 2
(12)
N − 1: 2
(13)
and for the DST-IV computation the permutation yn = −x n+N=2 ;
yN −1−n = xN=2−1−n ;
n = 0; 1; : : : ;
From above example it can be seen that permuted data sequences for the DCT-IV and DST-IV computation are actually in reverse order of each other. Consequently, we can eliminate the 2rst and second butterGy stages and replace them by appropriate permutation of the original data sequence. Now for a given N we are able to extract from the MDCT-modi2ed generalized signal Gow graph the new DCT-IV/DST-IV computational structure which will depend on two identical DCT-IIs of half sizes. It is shown in Fig. 2 for N = 16. The symbols in brackets correspond to the DST-IV computation. With respect to Fig. 1 and Eqs. (10) and (11), the corresponding output DCT-IV and DST-IV coe1cients are, respectively, obtained as IV = z2k ; c2k
IV c2k+1 = −zN −2k−2 ;
k = 0; 1; : : : ;
N −1 2
(14)
and IV = −z2k ; s2k
IV s2k+1 = −zN −2k−2 ;
k = 0; 1; : : : ;
N − 1: 2
(15)
The new DCT-IV/DST-IV computational structure represents the uni2ed DCT-IV/DST-IV computation and their inverses for any N being an even integer. In speci2c case of N = 2n , the exact number of arithmetic operations required for the DCT-IV/DST-IV computation is (N=2)(n+2) multiplications and (3N=2)n additions. This is the lowest achievable complexity for the DCT-IV [22].
V. Britanak / Signal Processing 83 (2003) 1803 – 1813
1809
Fig. 2. The new DCT-IV/DST-IV computational structure for N = 16.
4. New factorization of the DCT-IV matrix In fact, the new DCT-IV/DST-IV computational structure de2nes also a new sparse matrix factorization of the DCT-IV matrix. Let CNIV be the orthogonal matrix of order N . Then CNIV can be decomposed into the following sparse matrix product: CNIV
= QN
1 IN=2−1 −1 −JN=2−1
II −JN=2−1 CN=2 JN=2 0 −IN=2−1
0 II DN=2 CN=2
GN
JN=2
0
0
−JN=2
;
(16)
II where the matrix product CN=2 JN=2 denotes N=2 × N=2 DCT-II matrix with reversed order of its columns. The II matrix product CN=2 DN=2 denotes the N=2 × N=2 DCT-II matrix with sign changed of its odd-indexed rows, whereby DN=2 = diag{1; −1; 1; : : : ; (−1)N=2+1 } is the diagonal matrix. GN is the rotation matrix (Given’s plane
1810
V. Britanak / Signal Processing 83 (2003) 1803 – 1813
rotations) de2ned by (N − 1) cos 4N GN = (N − 1) sin 4N
cos
(N − 3) 4N
−sin ·
· cos 4N sin 4N
−sin 4N cos 4N
· sin
(N − 3) 4N
(N − 3) 4N
· cos
(N − 3) 4N
(N − 1) −sin 4N ; (N − 1) cos 4N (17)
and QN is the permutation matrix for reordering the coe1cients according to (14). The new factorization of the DCT-IV matrix is the direct consequence of the sparse matrix factorization derived for the MDCT matrix [4]. Finally, considering the MLT and its fast implementation based on DCT-IV [22], the alternative e1cient implementation of the forward and inverse MDCT computation in MP3 audio coding can be derived using the new DCT-IV/DST-IV computational structure. After the permutation applied to the input data sequence {x n } [22] yn+N=4 = x n − xN=2−1−n ;
yN=4−1−n = −x n+N=2 − xN −1−n ;
n = 0; 1; : : : ;
N − 1; 4
(18)
the MDCT is converted to the DCT-IV of {yn }, n = 0; 1; : : : ; N=2 − 1. Signal Gow graphs for the forward and inverse MDCT in MP3 are shown for N = 12 and 36 in Figs. 3(a) and (b), respectively. In order to minimize total arithmetic complexity we can use optimized odd-length DCT-II modules [2]. Three-point DCT-II module requires one multiplication, four additions and one shift, whereas 9-point DCT-II module requires eight multiplications, 34 additions and two shifts. Then arithmetic complexity of the forward MDCT computation in MP3 for N = 12 is 11 multiplications, 27 additions and two shifts, whereas for N = 36 is 43 multiplications, 129 additions and four shifts. The inverse MDCT requires exactly N=2 less additions than that of the forward MDCT. 5. Conclusions Although, both the DCT-IV and DST-IV play key role in the e1cient implementation of orthogonal lapped transforms and perfect reconstruction cosine-modulated 2lter banks such as the oddly stacked MDCT or equivalently MLT, it is shown that the DCT-IV and DST-IV can be e1ciently computed by the new fast MDCT algorithm for oddly stacked system [4], and hence via the MDCT. It follows from fact that the DCT-IV and DST-IV of double sizes are related to two variants of 2lter banks de2ned by Dolby Labs AC-3 digital audio compression algorithm, and these two variants of 2lter banks can be e1ciently realized by the new fast MDCT algorithm. Moreover, by careful analysis of regular structure of the new fast MDCT algorithm, the
V. Britanak / Signal Processing 83 (2003) 1803 – 1813
1811
(a)
(b) Fig. 3. Signal Gow graphs for the forward and inverse MDCT computation in MP3 (a) for N = 12 and (b) for N = 36.
new DCT-IV/DST-IV computational structure is extracted and new sparse matrix factorization of the DCT-IV matrix is suggested. Finally, the new DCT-IV/DST-IV computational structure provides the alternative e1cient implementation of the forward and inverse MDCT in MP3 audio coding. Acknowledgements The author would like to thank the anonymous reviewers for their helpful comments and suggestions, which signi2cantly improved the clarity of presentation in this paper. This work was partially supported by Slovak Scienti2c Grant Agency VEGA, project No. 2/1100/21.
1812
V. Britanak / Signal Processing 83 (2003) 1803 – 1813
References [1] G. Bonnerot, M. Bellanger, Odd-time odd-frequency discrete Fourier transform for symmetric real-valued series, Proc. IEEE 64 (March 1976) 392–393. [2] V. Britanak, The re2ned e1cient implementation of the MDCT in MP3 and comparison with other methods, Technical Report II SAS-2002-02, Institute of Informatics, Slovak Academy of Sciences, Bratislava, September 2002. [3] V. Britanak, K.R. Rao, The fast generalized discrete Fourier transforms: a uni2ed approach to the discrete sinusoidal transforms computation, Signal Process. 79 (December 1999) 135–150. [4] V. Britanak, K.R. Rao, A new fast algorithm for the uni2ed forward and inverse MDCT/MDST computation, Signal Process. 82 (March 2002) 433–459. [5] S.C. Chan, K.L. Ho, Direct methods for computing discrete sinusoidal transforms, IEE Proc.: Radar Signal Process. Part F 137 (December 1990) 433–442. [6] S.C. Chan, K.L. Ho, Fast algorithms for computing the discrete cosine transform, IEEE Trans. Circuits Syst. II: Analog Digital Signal Process. 39 (March 1992) 185–190. [7] D.Y. Chan, J.F. Yang, S.Y. Chen, Regular implementation algorithms of time domain aliasing cancellation, IEE Proc. Vision Image Signal Process. 143 (December 1996) 387–392. [8] W.H. Chen, C.H. Smith, S.C. Fralick, A fast computational algorithm for the discrete cosine transform, IEEE Trans. Commun. COM-25 (September 1977) 1004–1009. [9] Digital Audio Compression (AC-3) Standard, Document of Advanced Television Systems Committee (ATSC), Audio Specialist Group T3/S7, December 20, 1995. [10] O.K. Ersoy, Fourier-related transforms, fast algorithms and applications, Prentice-Hall, Upper Saddle River, NJ, 1997, pp. 388– 452 (Chapter 9). [11] O.K. Ersoy, N.C. Hu, A uni2ed approach to the fast computation of all discrete trigonometric transforms, Proceedings of the IEEE ICASSP’87, Dallas, TX, April 1987, pp. 1843–1846. [12] R. Gluth, Regular FFT-related transform kernels for DCT/DST-based polyphase 2lter banks, Proceedings of the IEEE ICASSP’91, Toronto, Canada, May 1991, pp. 2205 –2208. [13] A.K. Jain, A sinusoidal family of unitary transforms, IEEE Trans. Pattern Anal. Machine Intell. PAMI-1 (October 1979) 356–365. [14] Z. Jiang, A.N. Willson, Fast cosine-modulated 2lter banks, Proceedings of the International Conference on Digital Signal Processing (DSP’95), Limassol, Cyprus, June 1995, pp. 188–193. [15] C.W. Kok, Fast algorithm for computing discrete cosine transform, IEEE Trans. Signal Process. 45 (March 1997) 757–760. [16] P.S. Kumar, K.M.M. Prabhu, A set of new algorithms for DCTs and DSTs, Int. J. of Electron. 76 (April 1994) 553–586. [17] Ch.M. Liu, W.Ch. Lee, A uni2ed fast algorithm for cosine modulated 2lter banks in current audio coding standards, J. Audio Eng. Soc. 47 (12) (December 1999) 1061–1075. [18] H.S. Malvar, Lapped transforms for e1cient transform/subband coding, IEEE Trans. Acoust. Speech Signal Process. 38 (June 1990) 969–978. [19] H.S. Malvar, Fast algorithm for modulated lapped transform, Electron. Lett. 27 (April 1991) 775–776. [20] H.S. Malvar, Extended lapped transforms: fast algorithms and applications, Proceedings of the IEEE ICASSP’91, Toronto, Canada, May 1991, pp. 1797–1800. [21] H.S. Malvar, Extended lapped transforms: properties, applications, and fast algorithms, IEEE Trans. Signal Process. 40 (November 1992) 2703–2714. [22] H.S. Malvar, Signal Processing with Lapped Transforms, Artech House, Norwood, MA, 1992. [23] H.S. Malvar, Biorthogonal and nonuniform lapped transforms for transform coding with reduced blocking and ringing artifacts, IEEE Trans. on Signal Process. 46 (April 1998) 1043–1053. [24] H.S. Malvar, Fast algorithms for orthogonal and biorthogonal modulated lapped transforms, Proceedings of the IEEE Symposium on Advances in Digital Filtering and Signal Processing, Victoria, Canada, June 1998, pp. 159 –163. [25] H.S. Malvar, A modulated complex lapped transform and its application to audio processing, Proceedings of the IEEE ICASSP’99, Phoenix, AR, May 1999, pp. 1421–1424. [26] J.P. Princen, A.W. Johnson, A.B. Bradley, Subband/transform coding using 2lter bank designs based on time domain aliasing cancellation, Proceedings of the IEEE ICASSP’87, Dallas, TX, April 1987, pp. 2161–2164. [27] N. Rama Murthy, M.N.S. Swamy, On the algorithms for the computation of even discrete cosine transform-2 (EDCT-2) of real sequences, IEEE Trans. on Circuits Syst. 37 (May 1990) 625–627. [28] N. Rama Murthy, M.N.S. Swamy, On a novel decomposition of the DCT and its application, IEEE Trans. Signal Process. 41 (January 1993) 480–485. [29] N. Rama Murthy, M.N.S. Swamy, On the on-line computation of DCT-IV and DST-IV transforms, IEEE Trans. Signal Process. 43 (May 1995) 1249–1251. [30] K.R. Rao, P. Yip, Discrete cosine transform: algorithms, advantages applications, Academic Press, Boston, 1990. [31] K.R. Rao, P.C. Yip (Eds.), The transform and data compression handbook, CRC Press, Boca Raton, FL, 2001, pp. 117–195 (Chapter 4).
V. Britanak / Signal Processing 83 (2003) 1803 – 1813
1813
[32] S. Shlien, The modulated lapped transform, its varying forms, and application to audio coding, IEEE Trans. Speech Audio Process. 5 (July 1997) 359–366. [33] N. Suehiro, M. Hatori, Fast algorithms for the discrete Fourier transform and others transforms, Proceedings of the IEEE ICASSP’86, Tokyo, Japan, April 1986, pp. 217–220. [34] N. Suehiro, M. Hatori, Fast algorithms for the discrete Fourier transform and others transforms, IEEE Trans. Acoust. Speech Signal Process. ASSP-34 (June 1986) 642–644. [35] Z. Wang, Reconsideration of a fast computational algorithm for the discrete cosine transform, IEEE Trans. Commun. COM-31 (January 1983) 121–123. [36] Z. Wang, Fast algorithms for the discrete W transform and for the discrete Fourier transform, IEEE Trans. Acoust. Speech Signal Process. ASSP-32 (August 1984) 803–816. [37] Z. Wang, On computing the discrete Fourier and cosine transforms, IEEE Trans. Acoust. Speech Signal Process. ASSP-33 (October 1985) 1341–1344. [38] Z. Wang, Fast discrete sine transform algorithms, Signal Process. 19 (February 1990) 91–102. [39] Z. Wang, Comments on generalized discrete Hartley transforms, IEEE Trans. Signal Process. 43 (July 1995) 1711–1712. [40] Z. Wang, B.R. Hunt, The discrete W transform, Appl. Math. Comput. 16 (1985) 19–48. [41] Z. Xiong, H.S. Malvar, A nonuniform modulated complex lapped transform, IEEE Signal Process. Lett. 8 (September 2001) 257–260. [42] P. Yip, K.R. Rao, Fast decimation-in-time algorithms for a family of discrete sine and cosine transforms, Circuits Syst. Signal Process. 3 (4) (1984) 387–408. [43] P. Yip, K.R. Rao, Fast decimation-in-frequency algorithms for a family of discrete sine and cosine transforms, Circuits Syst. Signal Process. 7 (1) (1988) 3–19.