- Email: [email protected]
Fast algorithm for computing the forward and inverse MDCT in MPEG audio coding
ARTICLE IN PRESS
Signal Processing 86 (2006) 1055–1060 www.elsevier.com/locate/sigpro
Fast algorithm for computing the forward and inverse MDCT in MPEG audio coding T.K. Truonga,, P.D. Chena, T.C. Chengb a
Department of Information Engineering, College of Electrical and Information Engineering, I-Shou University, Ta-Hsu Hsiang, Koahsiung County, 84008, Taiwan, ROC b Department of Electrical Engineering-Electrophysics, University of Southern California, Los Angeles, USA Received 2 November 2004; received in revised form 15 May 2005; accepted 28 July 2005 Available online 31 August 2005
Abstract The modified discrete cosine transform (MDCT) is always employed in transform-coding schemes as the analysis/ synthesis filter bank. In this paper, an efficient algorithm for MDCT and inverse MDCT (IMDCT) computation for MPEG-1 audio layer III and MPEG-2 international audio-coding standards is proposed, using only the type-II DCT. Finally, the proposed algorithm is compared to the similar algorithms in this paper. r 2005 Elsevier B.V. All rights reserved. Keywords: MDCT/IMDCT; Audio coding; Type-II DCT; MPEG
1. Introduction The modified discrete cosine transform (MDCT)/inverse MDCT (IMDCT) is extensively used to realize the analysis/synthesis filter bank of time domain aliasing cancellation scheme for subband coding [1,2]. This transform has been adopted in several international standards and commercial products such as MPEG-1 [3], MPEG2 [4], and AC-3 [5] in audio coding. However, the computation of the MDCT and IMDCT in the Corresponding author. Tel.: +886 7 657 7711x6006;
fax: +886 7 657 7293. E-mail address: [email protected] (T.K. Truong).
MPEG audio coder can involve an extensive computation. Therefore, there is a need for a fast method to compute the forward and IMDCT so that it is naturally suitable for hardware or firmware implementations. During the past decade, many fast algorithms for computing the type-II DCT (DCT-II) and inverse DCT (IDCT-II) have been introduced [6–11]. Among them, Lee and Hou [6,7] developed the most attention fast algorithms for the 2m-point DCT. Also, Kok [10] suggested a new DCT-II algorithm which can be used to compute the evenlength DCT-II. In his algorithm, the DCT-II with even-length N ¼ 2M is decomposed into two DCTs of length M. If M is odd, it follows from
0165-1684/$ - see front matter r 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.sigpro.2005.07.027
ARTICLE IN PRESS 1056
T.K. Truong et al. / Signal Processing 86 (2006) 1055–1060
[11] that such a DCT of length M can be computed by the use of the identical length discrete Fourier transform (DFT) algorithm developed by Winograd [12]. Recently, using the Kok algorithm [10] together with the Heideman algorithm [11], Truong et al. [13] showed that a novel algorithm called the Kok–Heideman algorithm can be developed to compute an even-length DCT-II. In what follows, one will see that all of fast DCT-II algorithms described above can be used to compute the MDCT and IMDCT. In 1990, the modulated lapped transform equivalent to the MDCT was first proposed by Malver [14]. Typically, Duhamel et al. [15] proposed the fast implementation algorithm for the N-point MDCT. In other words, the N-point MDCT can be expressed in terms of a complex modified DFT computation of length N/4, where N ¼ 4M, which is implemented by the conventional split radix FFT or Winograd’s DFT algorithm. Recently, more efficient implementations of the forward and IMDCT for MPEG layer III were proposed in [16–22]. In 1994, Konstantinides [23] showed in detail that if N is a power of two, then the N-point MDCT and IMDCT can be converted first into the standard DCT-II and IDCT-II of half size, respectively. Such a DCT-II or IDCT-II can be then computed rapidly by the use of the Lee or Hou fast DCT algorithm. On the other hand, if N is a multiple of 4, i.e., N ¼ 4M, then a new algorithm for the MDCT/IMDCT was developed by the authors in [13]. In this algorithm, the MDCT and IMDCT can be converted first into the DCT-II of half length with pre- and postpermutations of data. Again, the Kok–Heideman algorithm can be utilized to compute the DCT-II of length N/2. The advantage of this algorithm is that it employs the same DCT-II to compute both the MDCT and IMDCT. However, this algorithm requires additional operations for computing the initial condition for the IMDCT decoder. In 2001, Lee [22] proposed a fast algorithm similar to Truong’s algorithm [13] for computing the MDCT and IMDCT. The advantage of this algorithm is that the computation of the initial condition needed in the IMDCT decoder given in [13] is completely avoided. Consequently, Lee’s algo-
rithm used to compute the IMDCT requires less operation than Truong’s algorithm. However, it requires slightly more operation than Truong’s algorithm for computing the MDCT. More recently, the authors [21] presented an algorithm, which can be used to realize the MDCT/IMDCT in a matrix approach. Actually, the ideas presented this algorithm are due to Lee’s algorithm. In this paper, by the use of the DCT-II, an efficient algorithm for computing the MDCT/IMDCT is presented. It uses Truong’s algorithm for computing the MDCT and Lee’s algorithm for computing the IMDCT. Finally, such a fast algorithm for the forward and IMDCT are compared with the known algorithms in terms of the number of arithmetic operations.
2. Some properties of the MDCT/IMDCT Let fyðnÞg be the input data sequence. The MDCT and its IMDCT [1,2] are given by N 1 X p N Y ðkÞ ¼ 2n þ 1 þ yðnÞ cos ð2k þ 1Þ , 2N 2 n¼0 0pkp
N 1 2
ð1Þ
and ^ ¼ yðnÞ
ðN=2Þ1 X
Y ðkÞ cos
k¼0
p N 2n þ 1 þ ð2k þ 1Þ , 2N 2
0pnpN 1,
ð2Þ
where N is the windows length and N/2 is the number of transform coefficients. For simplicity, assume that Y ðkÞ in (1) are MDCT coefficients after a proper windowing and time shifting of the ^ original sequence. The notation yðnÞ in (2) is the recovered data sequence obtained by the inverse transform of Y ðkÞ. Generally, it does not correspond to the original data sequence fyðnÞg. The MDCT of length N is usual to be a power of 2. In this paper, N is assumed to be a multiple of 4 for MPEG layer III. The substitution of N 1 n into n in (1) yields Y ðN 1 nÞ ¼ Y ðnÞ;
0pnp
N 1. 2
(3)
ARTICLE IN PRESS T.K. Truong et al. / Signal Processing 86 (2006) 1055–1060
From (3), the MDCT sequence in (1) has the ^ even antisymmetry property. yðnÞ has only N/2 independent variables and has the symmetry properties given by 3N 3N N þ n ¼ y^ 1 n ; 0pnp 1, y^ 4 4 4 (4) N N ^ y^ 1 n ¼ yðnÞ; 0pnp 1. (5) 2 4 ^ Therefore, if yðnÞ; n ¼ 0; 1; . . . ; ðN=4Þ 1 and ^ yðnÞ; n ¼ ð3N=4Þ 1; ð3N=4Þ 2; . . . ; N=2, are gi^ can be directly deduced ven, then other values yðnÞ from (4) and (5).
From (6), one observes that Y ð0Þ ¼ Y ð1Þ. Using this result and the recursive formula in (7), one has Y ð0Þ ¼ Y 0 ð0Þ=2,
N 1. (11) 2 Finally, using the sequence Y 0 ðkÞ obtained by (7), the sequence Y(K) for 0pkpðN=2Þ 1 can be obtained by the use of (10) and (11). &
Lemma 2. . Let ( 0 y n þ N4 ; 00 y ðnÞ ¼ y0 n 3N 4 ; Then
Instead of computing (1), an alternating technique is employed. To see this, replacing n by N 1 n in (1) yields
Y 0 ðkÞ ¼
N 1 X
0pkp
N 1 2
0pnp 3N 4 1; 3N 4 pnpN 1;
hp N
N 1. 2
(12)
i ð2n þ 1Þk ; ð13Þ
Proof. First, let n ¼ n0 þ N=4 (8) can now be rewritten in the form 1 hp i X N 0 0 0 Y ðkÞ ¼ y n þ cos ð2n0 þ 1Þk 4 N n0 ¼N=4
Lemma 1. If Y 0 ðkÞ ¼ Y ðkÞ þ Y ðk 1Þ for 0pkp
y00 ðnÞ cos
n¼0
N 1 X
ð6Þ
(10)
Y ðkÞ ¼ Y 0 ðkÞ Y ðk 1Þ; 1pkp
3. New MDCT algorithm
hp Y ðkÞ ¼ yðN 1 nÞ cos 2N n¼0 N N 2n þ 1 ð2k þ 1Þ ; 0pkp 1. 2 2
1057
þ (7)
3N=41 X n0 ¼0
hp i N y0 n0 þ cos ð2n0 þ 1Þk . 4 N ð14Þ
then
0
p N 0 0 Y ðkÞ ¼ 2n þ 1 y ðnÞ cos k , N 2 n¼0 N1 X
where y0 ðnÞ ¼ 2yðN 1 nÞ cos 0pnpN 1.
(8)
Also, let n þ N=4 ¼ n 3N=4 in the first summation and let n0 ¼ n in the second summation of (14). Then (14) is equivalent to Y 0 ðkÞ ¼
p N ð2n þ 1 Þ ; 2N 2 ð9Þ
hp i 3N cos ð2n þ 1Þk 2kp y0 n 4 N n¼3N=4 N 1 X
þ
ð3N=4Þ1 X n¼0
hp i N y nþ cos ð2n þ 1Þk . 4 N 0
ð15Þ Proof. Substituting (6) into (7) and using the trigonometric identities, one yield (8).
Since cos½ðp=NÞð2n þ 1Þk 2kp ¼ cos½ðp=NÞ ð2n þ 1Þk, the substitution this into (15) yields Eq. (13), thus completing the proof. &
ARTICLE IN PRESS T.K. Truong et al. / Signal Processing 86 (2006) 1055–1060
1058
Theorem 3. . For the known y00 ðnÞ as defined in (12), one yields ðN=2Þ1 X p ð2n þ 1Þk , X ðkÞ ¼ Y 0 ðkÞ ¼ xðnÞ cos 2ðN=2Þ n¼0 0pkp
N 1, 2
ð16Þ
where xðnÞ ¼ y00 ðnÞ þ y00 ðN 1 nÞ;
0pnp
N 1. 2
therefore composed of three successive steps of computation illustrated in the following algorithm. 1. Compute 8 3N y 4 1 n y 3N > 4 þ n bn ; > > > < 0pnp N 1; 4 3N xðnÞ ¼ N y 4 1 n bn ; y n > > 4 > > : N pnp N 1: 4 2
(17) Proof. The proof of this theorem is similar to that used in Lemma 2 in [23]. & From Theorem 3, it is obvious to see that (16) is an ðN=2Þ-point DCT-II which can be efficiently computed by the use of the Kok–Heideman algorithm. Suppose that yðnÞ is a sequence of length N. For 0pnpðN=4Þ 1, this implies that 3N=4pN 1 npN 1. Thus, from (12) and (9), (17) becomes N 00 00 0 xðnÞ ¼ y ðnÞ þ y ðN 1 nÞ ¼ y n þ 4 N 1n þ y0 4 hp i 3N ¼ 2y 1 n cos ð2n þ 1Þ 4 2N h i 3N p 2y þ n cos ð2n þ 1Þ 4 2N 3N 3N ¼ y 1n y þn 4 4 hp i 2 cos ð2n þ 1Þ . 2N Similarly, for N=4pnpN=2 1, (17) becomes N 0 0 5N 1n xðnÞ ¼ y n þ þy 4 4 N 3N ¼ yðn Þ yð 1 nÞ 4 4 hp i 2 cos ð2n þ 1Þ . 2N Let bn ¼ 2 cos½ðp=2NÞð2n þ 1Þ for 0pnpN=2 1 be the pre-computed constants. In general, the fast algorithm for computing the MDCT of yðnÞ in (1) is
2. Use the Kok–Heideman algorithm to compute the ðN=2Þ-point DCT-II of xðnÞ, i.e., Y 0 ðkÞ defined in (16). That is 0
Y ðkÞ ¼
ðN=2Þ1 X
xðnÞ cos
n¼0
for 0pkp
hp N
ð2n þ 1Þk
i
N 1. 2
3. Let Y ð0Þ ¼ Y 0 ð0Þ=2. Then compute Y ðkÞ ¼ Y 0 ðkÞ Y ðk 1Þ for 1pkpN=2 1. Note that the algorithm described above is a generalization of the one given in [13]. In particular, it can be seen that for the fast MDCT computation, only efficient 3-point and 9point DCT-II modules are needed. These short odd-length DCT-II prime factor modules can be computed by the Kok–Heideman algorithm.
4. New IMDCT algorithm Before developing a new IMDCT algorithm, let us briefly summarize the algorithms given in [13,22]. This will be helpful in the following developments. It is shown in [13] that the IMDCT in (2) can be obtained by using the (N/2)-point DCT-II as the forward MDCT. To see this, let ck ¼ cos½ðp=2NÞðN=2 þ 1Þð2k þ 1Þ be the pre-computed constants. For a given output sequence of (1), i.e., Y(k) for 0pkpðN=2Þ 1, the general algorithm for computing (2) can be shown to be
ARTICLE IN PRESS T.K. Truong et al. / Signal Processing 86 (2006) 1055–1060
composed of 4 successive steps of computation illustrated in the following algorithm.
1059
Table 1 Comparison of number of multiplications and additions for N ¼ 12, 36
1. Compute Y 0 ðkÞ ¼ Y ðkÞbk for
Length N ¼ 12
Length N ¼ 36
New algorithm
11/27-MDCT 11/23-IMDCT
43/129-MDCT 43/115-IMDCT
Truong’s algorithm [13]
11/27-MDCT 17/26-IMDCT
43/129-MDCT 61/128-IMDCT
Britanak’s algorithm [20]
14/36-MDCT 14/24-IMDCT
52/156-MDCT 56/124-IMDCT
Lee’s algorithm [22]
11/29-MDCT 11/23-IMDCT
43/133-MDCT 43/115-IMDCT
0pkpðN=2Þ 1. 2. Use the Kok–Heideman algorithm to compute the ðN=2Þ-point DCT-II of Y 0 ðkÞ; that is, ðN=2Þ1 X p ð2k þ 1Þn , y^ 00 ðnÞ ¼ Y 0 ðkÞ cos 2ðN=2Þ k¼0 0pnp
N 1. 2
3. Compute 8 00 y^ ðn þ N=4Þ; > > > < 0; y^ 0 ðnÞ ¼ > y^ 00 ð3N=4 nÞ; > > : 00 y^ ðn 3N=4Þ;
^ ¼ 4. Let yð0Þ
ðN=2Þ1 P
0pnpN=4 1; n ¼ N=4;
5. Summary
ðN=4Þ þ 1pnpð3N=4Þ 1; 3N=4pnpN 1;
^ ¼ Y ðkÞck . Then compute yðnÞ
k¼0
^ 1Þ for 1pnpN 1. y^ 0 ðnÞ yðn Note that the above algorithm is also a generalization of the one given in [13]. In the above algorithm, the computation of the ^ initial condition yð0Þ in step 4 is needed. Obviously, the number of multiplications and addi^ tions needed to realized yð0Þ are N=2 and ðN=2Þ 1, respectively. In order to overcome this problem, a modification of the above algorithm was developed by Lee, which is now so important to audio compress in what is called Lee’s algorithm. That is, the N-point IMDCT can be converted first into the standard DCT-IV with even-length N=2. Then such a DCT-IV can be computed by the use of the DCT-II. Finally, the DCT-II can be further decomposed into ðN=4Þpoint DCT-IIs. As mentioned earlier, the advantage of Lee’s algorithm [22] over Truong’s algorithm for computing the IMDCT is that it is no need to compute the initial condition, namely, ^ yð0Þ. As a consequence, Lee’s algorithm is used to compute the N-point IMDCT for N ¼ 12; 36. Actually, it is based on the DCT-II.
In summary, the new algorithm is composed of the MDCT algorithm given in Section 3 and the IMDCT given in Section 4. A comparison of this new algorithm with Truong’s algorithm, Britanak’s algorithm [20], and Lee’s algorithm in terms of the complexity of the MDCT and IMDCT for N ¼ 12; 36 is made in Table 1. From this table, one observes that the new algorithm requires fewer operations than those of known algorithm.
References [1] J.P. Princen, A.B. Bradley, Analysis/synthesis filter bank designs based on time domain aliasing cancellation, IEEE Trans. Acoust. Speech Signal Process. ASSP 34 (October 1986) 1153–1161. [2] J.P. Princen, A.W. Johnson, A.B. Bradley, Subband/ transform coding using filter bank designs based on time domain aliasing cancellation, Proceedings of IEEE International Conference on Acoustics, Speech, and Signal Processing, Dallas, TX, April 1987, pp. 2161–2164. [3] Information technology-coding of moving pictures and associated audio for digital storage media at up to about 1.5 Mbit/s-Part 3: Audio, ISO/IEC IS 11172-3, (MPEG1), 1992. [4] Information technology-generic coding of moving pictures and associated audio-Part 3: Audio, ISO/IEC IS 13818-3, (MPEG2), 1994. [5] Digital audio compression (AC-3) standard, Document of Advanced Television Systems Committee (ATSC): Audio Specialist Group T3/S7, December 1995.
ARTICLE IN PRESS 1060
T.K. Truong et al. / Signal Processing 86 (2006) 1055–1060
[6] B.G. Lee, A new algorithm to compute the discrete cosine transform, IEEE Trans. Acoust. Speech Signal Process. ASSP-2 (6) (December 1984) 1243–1245. [7] H.S. Hou, A fast recursive algorithm for computing the discrete cosine transform, IEEE Trans. Acoust. Speech Signal Process. ASSP-35 (10) (October 1987) 1455–1461. [8] M. Vetterli, H. Nussbaumer, Simple FFT and DCT algorithm with reduced number of operations, Signal Processing 6 (4) (August 1984) 267–278. [9] J. Makhoul, A fast cosine transform in one and two dimensions, IEEE Trans. Acoust. Speech Signal Process. ASSP-28 (February 1980) 27–34. [10] C.W. Kok, Fast algorithm for computing discrete cosine transform, IEEE Trans. Signal Process. 45 (3) (March 1997) 757–760. [11] M.T. Heideman, Computation of an odd-length DCT from a real-value DFT of the same length, IEEE Trans. Signal Process. 40 (1) (January 1992) 54–61. [12] S. Winogard, On computing the discrete Fourier transform, Math. Comp. 32 (January 1978) 175–199. [13] T.K. Truong, T.C. Cheng, J.H. Jeng, J.M. Kuo, C.L. Fu, Fast Algorithm for the Forward and Inverse MDCT in Audio Coding, US Pending Patent, February 2000. [14] H.S. Malvar, Lapped transforms for efficient transform/ subband coding, IEEE Trans. Acoust. Speech Signal Process. 38 (6) (June 1990) 969–978. [15] P. Duhamel, Y. Mahieux, J.P. Petit, A fast algorithm for the implementation of filter banks based on time domain aliasing cancellation, Proceedings of the IEEE Interna-
[16]
[17]
[18]
[19]
[20]
[21]
[22]
[23]
tional Conference on Acoustics, Speech, Signal Processing 1991, Toronto, ON, Canada, May 1991, pp. 2209–2212. Y.H. Fan, V.K. Madisetti, R.M. Mersereau, On fast algorithms for computing the inverse modified discrete cosine transform, IEEE Signal Process. Lett. 6 (3) (March 1999) 61–64. D.Y. Chan, J.F. Yang, S.Y. Chen, Regular implementation algorithm of time domain aliasing cancellation, IEE Proc. Vision, Image Signal Process. 143 (6) (December 1996) 387–392. C.H. Chen, B.D. Liu, J.F. Yang, Recursive architectures for realizing modified discrete cosine transform and its inverse, IEEE Trans. Circuits Systems-II: Analog Digital Signal Process. 50 (1) (January 2003) 38–45. C.M. Liu, W.C. Lee, A unified fast algorithm for cosine modulated filter banks in current audio coding standards, J. Audio Eng. Soc. 47 (12) (December 1999) 1061–1075. V.B. Britanak, K.R. Rao, An efficient implementation of the forward and inverse MDCT in MPEG audio coding, IEEE Signal Process. Lett. 8 (2) (February 2001) 48–51. M.H. Cheng, Y.H. Hsu, Fast IMDCT and MDCT algorithms—a matrix approach, IEEE Trans. Signal Process. 51 (1) (January 2003) 221–229. S.W. Lee, Improved algorithm for efficient computation of the forward and backward MDCT in MPEG audio coder, IEEE Trans. Circuits Syst.-II: Analog Digital Signal Process. 48 (10) (October 2001) 990–994. K. Konstantinides, Fast subband filtering in MPEG audio coding, IEEE Signal Process. Lett. 1 (2) (February 1994) 26–28.
Signal Processing 86 (2006) 1055–1060 www.elsevier.com/locate/sigpro
Fast algorithm for computing the forward and inverse MDCT in MPEG audio coding T.K. Truonga,, P.D. Chena, T.C. Chengb a
Department of Information Engineering, College of Electrical and Information Engineering, I-Shou University, Ta-Hsu Hsiang, Koahsiung County, 84008, Taiwan, ROC b Department of Electrical Engineering-Electrophysics, University of Southern California, Los Angeles, USA Received 2 November 2004; received in revised form 15 May 2005; accepted 28 July 2005 Available online 31 August 2005
Abstract The modified discrete cosine transform (MDCT) is always employed in transform-coding schemes as the analysis/ synthesis filter bank. In this paper, an efficient algorithm for MDCT and inverse MDCT (IMDCT) computation for MPEG-1 audio layer III and MPEG-2 international audio-coding standards is proposed, using only the type-II DCT. Finally, the proposed algorithm is compared to the similar algorithms in this paper. r 2005 Elsevier B.V. All rights reserved. Keywords: MDCT/IMDCT; Audio coding; Type-II DCT; MPEG
1. Introduction The modified discrete cosine transform (MDCT)/inverse MDCT (IMDCT) is extensively used to realize the analysis/synthesis filter bank of time domain aliasing cancellation scheme for subband coding [1,2]. This transform has been adopted in several international standards and commercial products such as MPEG-1 [3], MPEG2 [4], and AC-3 [5] in audio coding. However, the computation of the MDCT and IMDCT in the Corresponding author. Tel.: +886 7 657 7711x6006;
fax: +886 7 657 7293. E-mail address: [email protected] (T.K. Truong).
MPEG audio coder can involve an extensive computation. Therefore, there is a need for a fast method to compute the forward and IMDCT so that it is naturally suitable for hardware or firmware implementations. During the past decade, many fast algorithms for computing the type-II DCT (DCT-II) and inverse DCT (IDCT-II) have been introduced [6–11]. Among them, Lee and Hou [6,7] developed the most attention fast algorithms for the 2m-point DCT. Also, Kok [10] suggested a new DCT-II algorithm which can be used to compute the evenlength DCT-II. In his algorithm, the DCT-II with even-length N ¼ 2M is decomposed into two DCTs of length M. If M is odd, it follows from
0165-1684/$ - see front matter r 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.sigpro.2005.07.027
ARTICLE IN PRESS 1056
T.K. Truong et al. / Signal Processing 86 (2006) 1055–1060
[11] that such a DCT of length M can be computed by the use of the identical length discrete Fourier transform (DFT) algorithm developed by Winograd [12]. Recently, using the Kok algorithm [10] together with the Heideman algorithm [11], Truong et al. [13] showed that a novel algorithm called the Kok–Heideman algorithm can be developed to compute an even-length DCT-II. In what follows, one will see that all of fast DCT-II algorithms described above can be used to compute the MDCT and IMDCT. In 1990, the modulated lapped transform equivalent to the MDCT was first proposed by Malver [14]. Typically, Duhamel et al. [15] proposed the fast implementation algorithm for the N-point MDCT. In other words, the N-point MDCT can be expressed in terms of a complex modified DFT computation of length N/4, where N ¼ 4M, which is implemented by the conventional split radix FFT or Winograd’s DFT algorithm. Recently, more efficient implementations of the forward and IMDCT for MPEG layer III were proposed in [16–22]. In 1994, Konstantinides [23] showed in detail that if N is a power of two, then the N-point MDCT and IMDCT can be converted first into the standard DCT-II and IDCT-II of half size, respectively. Such a DCT-II or IDCT-II can be then computed rapidly by the use of the Lee or Hou fast DCT algorithm. On the other hand, if N is a multiple of 4, i.e., N ¼ 4M, then a new algorithm for the MDCT/IMDCT was developed by the authors in [13]. In this algorithm, the MDCT and IMDCT can be converted first into the DCT-II of half length with pre- and postpermutations of data. Again, the Kok–Heideman algorithm can be utilized to compute the DCT-II of length N/2. The advantage of this algorithm is that it employs the same DCT-II to compute both the MDCT and IMDCT. However, this algorithm requires additional operations for computing the initial condition for the IMDCT decoder. In 2001, Lee [22] proposed a fast algorithm similar to Truong’s algorithm [13] for computing the MDCT and IMDCT. The advantage of this algorithm is that the computation of the initial condition needed in the IMDCT decoder given in [13] is completely avoided. Consequently, Lee’s algo-
rithm used to compute the IMDCT requires less operation than Truong’s algorithm. However, it requires slightly more operation than Truong’s algorithm for computing the MDCT. More recently, the authors [21] presented an algorithm, which can be used to realize the MDCT/IMDCT in a matrix approach. Actually, the ideas presented this algorithm are due to Lee’s algorithm. In this paper, by the use of the DCT-II, an efficient algorithm for computing the MDCT/IMDCT is presented. It uses Truong’s algorithm for computing the MDCT and Lee’s algorithm for computing the IMDCT. Finally, such a fast algorithm for the forward and IMDCT are compared with the known algorithms in terms of the number of arithmetic operations.
2. Some properties of the MDCT/IMDCT Let fyðnÞg be the input data sequence. The MDCT and its IMDCT [1,2] are given by N 1 X p N Y ðkÞ ¼ 2n þ 1 þ yðnÞ cos ð2k þ 1Þ , 2N 2 n¼0 0pkp
N 1 2
ð1Þ
and ^ ¼ yðnÞ
ðN=2Þ1 X
Y ðkÞ cos
k¼0
p N 2n þ 1 þ ð2k þ 1Þ , 2N 2
0pnpN 1,
ð2Þ
where N is the windows length and N/2 is the number of transform coefficients. For simplicity, assume that Y ðkÞ in (1) are MDCT coefficients after a proper windowing and time shifting of the ^ original sequence. The notation yðnÞ in (2) is the recovered data sequence obtained by the inverse transform of Y ðkÞ. Generally, it does not correspond to the original data sequence fyðnÞg. The MDCT of length N is usual to be a power of 2. In this paper, N is assumed to be a multiple of 4 for MPEG layer III. The substitution of N 1 n into n in (1) yields Y ðN 1 nÞ ¼ Y ðnÞ;
0pnp
N 1. 2
(3)
ARTICLE IN PRESS T.K. Truong et al. / Signal Processing 86 (2006) 1055–1060
From (3), the MDCT sequence in (1) has the ^ even antisymmetry property. yðnÞ has only N/2 independent variables and has the symmetry properties given by 3N 3N N þ n ¼ y^ 1 n ; 0pnp 1, y^ 4 4 4 (4) N N ^ y^ 1 n ¼ yðnÞ; 0pnp 1. (5) 2 4 ^ Therefore, if yðnÞ; n ¼ 0; 1; . . . ; ðN=4Þ 1 and ^ yðnÞ; n ¼ ð3N=4Þ 1; ð3N=4Þ 2; . . . ; N=2, are gi^ can be directly deduced ven, then other values yðnÞ from (4) and (5).
From (6), one observes that Y ð0Þ ¼ Y ð1Þ. Using this result and the recursive formula in (7), one has Y ð0Þ ¼ Y 0 ð0Þ=2,
N 1. (11) 2 Finally, using the sequence Y 0 ðkÞ obtained by (7), the sequence Y(K) for 0pkpðN=2Þ 1 can be obtained by the use of (10) and (11). &
Lemma 2. . Let ( 0 y n þ N4 ; 00 y ðnÞ ¼ y0 n 3N 4 ; Then
Instead of computing (1), an alternating technique is employed. To see this, replacing n by N 1 n in (1) yields
Y 0 ðkÞ ¼
N 1 X
0pkp
N 1 2
0pnp 3N 4 1; 3N 4 pnpN 1;
hp N
N 1. 2
(12)
i ð2n þ 1Þk ; ð13Þ
Proof. First, let n ¼ n0 þ N=4 (8) can now be rewritten in the form 1 hp i X N 0 0 0 Y ðkÞ ¼ y n þ cos ð2n0 þ 1Þk 4 N n0 ¼N=4
Lemma 1. If Y 0 ðkÞ ¼ Y ðkÞ þ Y ðk 1Þ for 0pkp
y00 ðnÞ cos
n¼0
N 1 X
ð6Þ
(10)
Y ðkÞ ¼ Y 0 ðkÞ Y ðk 1Þ; 1pkp
3. New MDCT algorithm
hp Y ðkÞ ¼ yðN 1 nÞ cos 2N n¼0 N N 2n þ 1 ð2k þ 1Þ ; 0pkp 1. 2 2
1057
þ (7)
3N=41 X n0 ¼0
hp i N y0 n0 þ cos ð2n0 þ 1Þk . 4 N ð14Þ
then
0
p N 0 0 Y ðkÞ ¼ 2n þ 1 y ðnÞ cos k , N 2 n¼0 N1 X
where y0 ðnÞ ¼ 2yðN 1 nÞ cos 0pnpN 1.
(8)
Also, let n þ N=4 ¼ n 3N=4 in the first summation and let n0 ¼ n in the second summation of (14). Then (14) is equivalent to Y 0 ðkÞ ¼
p N ð2n þ 1 Þ ; 2N 2 ð9Þ
hp i 3N cos ð2n þ 1Þk 2kp y0 n 4 N n¼3N=4 N 1 X
þ
ð3N=4Þ1 X n¼0
hp i N y nþ cos ð2n þ 1Þk . 4 N 0
ð15Þ Proof. Substituting (6) into (7) and using the trigonometric identities, one yield (8).
Since cos½ðp=NÞð2n þ 1Þk 2kp ¼ cos½ðp=NÞ ð2n þ 1Þk, the substitution this into (15) yields Eq. (13), thus completing the proof. &
ARTICLE IN PRESS T.K. Truong et al. / Signal Processing 86 (2006) 1055–1060
1058
Theorem 3. . For the known y00 ðnÞ as defined in (12), one yields ðN=2Þ1 X p ð2n þ 1Þk , X ðkÞ ¼ Y 0 ðkÞ ¼ xðnÞ cos 2ðN=2Þ n¼0 0pkp
N 1, 2
ð16Þ
where xðnÞ ¼ y00 ðnÞ þ y00 ðN 1 nÞ;
0pnp
N 1. 2
therefore composed of three successive steps of computation illustrated in the following algorithm. 1. Compute 8 3N y 4 1 n y 3N > 4 þ n bn ; > > > < 0pnp N 1; 4 3N xðnÞ ¼ N y 4 1 n bn ; y n > > 4 > > : N pnp N 1: 4 2
(17) Proof. The proof of this theorem is similar to that used in Lemma 2 in [23]. & From Theorem 3, it is obvious to see that (16) is an ðN=2Þ-point DCT-II which can be efficiently computed by the use of the Kok–Heideman algorithm. Suppose that yðnÞ is a sequence of length N. For 0pnpðN=4Þ 1, this implies that 3N=4pN 1 npN 1. Thus, from (12) and (9), (17) becomes N 00 00 0 xðnÞ ¼ y ðnÞ þ y ðN 1 nÞ ¼ y n þ 4 N 1n þ y0 4 hp i 3N ¼ 2y 1 n cos ð2n þ 1Þ 4 2N h i 3N p 2y þ n cos ð2n þ 1Þ 4 2N 3N 3N ¼ y 1n y þn 4 4 hp i 2 cos ð2n þ 1Þ . 2N Similarly, for N=4pnpN=2 1, (17) becomes N 0 0 5N 1n xðnÞ ¼ y n þ þy 4 4 N 3N ¼ yðn Þ yð 1 nÞ 4 4 hp i 2 cos ð2n þ 1Þ . 2N Let bn ¼ 2 cos½ðp=2NÞð2n þ 1Þ for 0pnpN=2 1 be the pre-computed constants. In general, the fast algorithm for computing the MDCT of yðnÞ in (1) is
2. Use the Kok–Heideman algorithm to compute the ðN=2Þ-point DCT-II of xðnÞ, i.e., Y 0 ðkÞ defined in (16). That is 0
Y ðkÞ ¼
ðN=2Þ1 X
xðnÞ cos
n¼0
for 0pkp
hp N
ð2n þ 1Þk
i
N 1. 2
3. Let Y ð0Þ ¼ Y 0 ð0Þ=2. Then compute Y ðkÞ ¼ Y 0 ðkÞ Y ðk 1Þ for 1pkpN=2 1. Note that the algorithm described above is a generalization of the one given in [13]. In particular, it can be seen that for the fast MDCT computation, only efficient 3-point and 9point DCT-II modules are needed. These short odd-length DCT-II prime factor modules can be computed by the Kok–Heideman algorithm.
4. New IMDCT algorithm Before developing a new IMDCT algorithm, let us briefly summarize the algorithms given in [13,22]. This will be helpful in the following developments. It is shown in [13] that the IMDCT in (2) can be obtained by using the (N/2)-point DCT-II as the forward MDCT. To see this, let ck ¼ cos½ðp=2NÞðN=2 þ 1Þð2k þ 1Þ be the pre-computed constants. For a given output sequence of (1), i.e., Y(k) for 0pkpðN=2Þ 1, the general algorithm for computing (2) can be shown to be
ARTICLE IN PRESS T.K. Truong et al. / Signal Processing 86 (2006) 1055–1060
composed of 4 successive steps of computation illustrated in the following algorithm.
1059
Table 1 Comparison of number of multiplications and additions for N ¼ 12, 36
1. Compute Y 0 ðkÞ ¼ Y ðkÞbk for
Length N ¼ 12
Length N ¼ 36
New algorithm
11/27-MDCT 11/23-IMDCT
43/129-MDCT 43/115-IMDCT
Truong’s algorithm [13]
11/27-MDCT 17/26-IMDCT
43/129-MDCT 61/128-IMDCT
Britanak’s algorithm [20]
14/36-MDCT 14/24-IMDCT
52/156-MDCT 56/124-IMDCT
Lee’s algorithm [22]
11/29-MDCT 11/23-IMDCT
43/133-MDCT 43/115-IMDCT
0pkpðN=2Þ 1. 2. Use the Kok–Heideman algorithm to compute the ðN=2Þ-point DCT-II of Y 0 ðkÞ; that is, ðN=2Þ1 X p ð2k þ 1Þn , y^ 00 ðnÞ ¼ Y 0 ðkÞ cos 2ðN=2Þ k¼0 0pnp
N 1. 2
3. Compute 8 00 y^ ðn þ N=4Þ; > > > < 0; y^ 0 ðnÞ ¼ > y^ 00 ð3N=4 nÞ; > > : 00 y^ ðn 3N=4Þ;
^ ¼ 4. Let yð0Þ
ðN=2Þ1 P
0pnpN=4 1; n ¼ N=4;
5. Summary
ðN=4Þ þ 1pnpð3N=4Þ 1; 3N=4pnpN 1;
^ ¼ Y ðkÞck . Then compute yðnÞ
k¼0
^ 1Þ for 1pnpN 1. y^ 0 ðnÞ yðn Note that the above algorithm is also a generalization of the one given in [13]. In the above algorithm, the computation of the ^ initial condition yð0Þ in step 4 is needed. Obviously, the number of multiplications and addi^ tions needed to realized yð0Þ are N=2 and ðN=2Þ 1, respectively. In order to overcome this problem, a modification of the above algorithm was developed by Lee, which is now so important to audio compress in what is called Lee’s algorithm. That is, the N-point IMDCT can be converted first into the standard DCT-IV with even-length N=2. Then such a DCT-IV can be computed by the use of the DCT-II. Finally, the DCT-II can be further decomposed into ðN=4Þpoint DCT-IIs. As mentioned earlier, the advantage of Lee’s algorithm [22] over Truong’s algorithm for computing the IMDCT is that it is no need to compute the initial condition, namely, ^ yð0Þ. As a consequence, Lee’s algorithm is used to compute the N-point IMDCT for N ¼ 12; 36. Actually, it is based on the DCT-II.
In summary, the new algorithm is composed of the MDCT algorithm given in Section 3 and the IMDCT given in Section 4. A comparison of this new algorithm with Truong’s algorithm, Britanak’s algorithm [20], and Lee’s algorithm in terms of the complexity of the MDCT and IMDCT for N ¼ 12; 36 is made in Table 1. From this table, one observes that the new algorithm requires fewer operations than those of known algorithm.
References [1] J.P. Princen, A.B. Bradley, Analysis/synthesis filter bank designs based on time domain aliasing cancellation, IEEE Trans. Acoust. Speech Signal Process. ASSP 34 (October 1986) 1153–1161. [2] J.P. Princen, A.W. Johnson, A.B. Bradley, Subband/ transform coding using filter bank designs based on time domain aliasing cancellation, Proceedings of IEEE International Conference on Acoustics, Speech, and Signal Processing, Dallas, TX, April 1987, pp. 2161–2164. [3] Information technology-coding of moving pictures and associated audio for digital storage media at up to about 1.5 Mbit/s-Part 3: Audio, ISO/IEC IS 11172-3, (MPEG1), 1992. [4] Information technology-generic coding of moving pictures and associated audio-Part 3: Audio, ISO/IEC IS 13818-3, (MPEG2), 1994. [5] Digital audio compression (AC-3) standard, Document of Advanced Television Systems Committee (ATSC): Audio Specialist Group T3/S7, December 1995.
ARTICLE IN PRESS 1060
T.K. Truong et al. / Signal Processing 86 (2006) 1055–1060
[6] B.G. Lee, A new algorithm to compute the discrete cosine transform, IEEE Trans. Acoust. Speech Signal Process. ASSP-2 (6) (December 1984) 1243–1245. [7] H.S. Hou, A fast recursive algorithm for computing the discrete cosine transform, IEEE Trans. Acoust. Speech Signal Process. ASSP-35 (10) (October 1987) 1455–1461. [8] M. Vetterli, H. Nussbaumer, Simple FFT and DCT algorithm with reduced number of operations, Signal Processing 6 (4) (August 1984) 267–278. [9] J. Makhoul, A fast cosine transform in one and two dimensions, IEEE Trans. Acoust. Speech Signal Process. ASSP-28 (February 1980) 27–34. [10] C.W. Kok, Fast algorithm for computing discrete cosine transform, IEEE Trans. Signal Process. 45 (3) (March 1997) 757–760. [11] M.T. Heideman, Computation of an odd-length DCT from a real-value DFT of the same length, IEEE Trans. Signal Process. 40 (1) (January 1992) 54–61. [12] S. Winogard, On computing the discrete Fourier transform, Math. Comp. 32 (January 1978) 175–199. [13] T.K. Truong, T.C. Cheng, J.H. Jeng, J.M. Kuo, C.L. Fu, Fast Algorithm for the Forward and Inverse MDCT in Audio Coding, US Pending Patent, February 2000. [14] H.S. Malvar, Lapped transforms for efficient transform/ subband coding, IEEE Trans. Acoust. Speech Signal Process. 38 (6) (June 1990) 969–978. [15] P. Duhamel, Y. Mahieux, J.P. Petit, A fast algorithm for the implementation of filter banks based on time domain aliasing cancellation, Proceedings of the IEEE Interna-
[16]
[17]
[18]
[19]
[20]
[21]
[22]
[23]
tional Conference on Acoustics, Speech, Signal Processing 1991, Toronto, ON, Canada, May 1991, pp. 2209–2212. Y.H. Fan, V.K. Madisetti, R.M. Mersereau, On fast algorithms for computing the inverse modified discrete cosine transform, IEEE Signal Process. Lett. 6 (3) (March 1999) 61–64. D.Y. Chan, J.F. Yang, S.Y. Chen, Regular implementation algorithm of time domain aliasing cancellation, IEE Proc. Vision, Image Signal Process. 143 (6) (December 1996) 387–392. C.H. Chen, B.D. Liu, J.F. Yang, Recursive architectures for realizing modified discrete cosine transform and its inverse, IEEE Trans. Circuits Systems-II: Analog Digital Signal Process. 50 (1) (January 2003) 38–45. C.M. Liu, W.C. Lee, A unified fast algorithm for cosine modulated filter banks in current audio coding standards, J. Audio Eng. Soc. 47 (12) (December 1999) 1061–1075. V.B. Britanak, K.R. Rao, An efficient implementation of the forward and inverse MDCT in MPEG audio coding, IEEE Signal Process. Lett. 8 (2) (February 2001) 48–51. M.H. Cheng, Y.H. Hsu, Fast IMDCT and MDCT algorithms—a matrix approach, IEEE Trans. Signal Process. 51 (1) (January 2003) 221–229. S.W. Lee, Improved algorithm for efficient computation of the forward and backward MDCT in MPEG audio coder, IEEE Trans. Circuits Syst.-II: Analog Digital Signal Process. 48 (10) (October 2001) 990–994. K. Konstantinides, Fast subband filtering in MPEG audio coding, IEEE Signal Process. Lett. 1 (2) (February 1994) 26–28.