1348 IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 45, NO. 5, MAY 1997
[5] M. Honig, U. Madhow, and S. Verd´u, “Blind adaptive multiuser detec-tion,” IEEE Trans. Inform. Theory, vol. 41, pp. 944–960, July 1995. [6] P. Comon, “Independent component analysis, A new concept?,” Signal
Processing, vol. 36, pp. 287–314, Apr. 1994.
[7] C. Jutten and J. H´erault, “Blind separation of sources, Part I,” Signal Processing, vol. 24, pp. 1–10, July 1991.
[8] N. Delfosse and P. Loubaton, “Adaptive blind separation of independent sources: A deflation approach,” Signal Processing, vol. 45, pp. 59–83, 1995.
[9] A. J. Bell and T. J. Sejnowski, “An information-maximization approach to blind separation and blind deconvolution,” Neural Comput., vol. 7, pp. 1129–1159, Nov. 1995.
[10] E. Moreau and O. Macchi, “High-order contrasts for self-adaptive source separation,” Int. J. Adaptive Cont. Signal Processing., vol. 10, pp. 19–46, 1996.
[11] J. F. Cardoso and B. Laheld, “Equivariant adaptive source separation,” IEEE Trans. Signal Processing, vol. 44, pp. 3017–3030, Dec. 1996. [12] D. N. Godard, “Self-recovering equalization and carrier tracking in
two-dimensional data communication systems,” IEEE Trans. Commun., vol. COM-28, pp. 1867–1875/2523–2526, Nov. 1980.
[13] J. Lundell and B. Widrow, “Application of the constant modulus adaptive beamformer to constant and nonconstant modulus algorithms,” in Proc. 22nd Asilomar Conf. Signals, Syst. Comput., Pacific Grove, CA, Nov. 1987.
[14] J. G. Proakis, Digital Communications. Singapore: McGraw-Hill, 1995.
[15] C. B. Papadias and A. Paulraj, “A space-time constant modulus algo-rithm for SDMA systems,” in Proc. IEEE/VTS 46th Veh. Technol. Conf. (VTC-96), Atlanta, GA, Apr. 1996.
[16] G. Strang, Linear Algebra and its Applications. San Diego, CA: Harcourt Brace Jovanovich, 1988.
Normal Bases Expansion of the Discrete Cosine Transform Ming-Chwen Yang, Ja-Ling Wu, and Yuh-Ming Huang
Abstract— In this correspondence, we show that the discrete cosine transform (DCT) can be obtained by projecting the discrete Fourier trans-form from the extension field to the basefield. Applying the framework of projection operator, a fast fully recursive algorithm for computing the DCT is also presented.
I. INTRODUCTION
In their recent work [1], Hong et al. have drawn out the relationship between the discrete Hartley transform (DHT) and the discrete Fourier transform (DFT) from the viewpoints of field extension and projection. This approach not only demonstrates the intimate connection among various transforms but also presents a powerful framework for developing fast transformation algorithms based on the standard FFT algorithms.
In this correspondence, the result of [1] is extended to derive the discrete cosine transform (DCT), which is an important discrete
Manuscript received March 16, 1995; revised Decmber 9, 1996. The associate editor coordinating the review of this paper and approving it for publication was Dr. Nurgun Erdol.
M.-C. Yang and Y.-M. Huang are with the Department of Computer Science and Information Engineering, National Taiwan University, Taipei, Taiwan (e-mail: [email protected]; [email protected]).
J.-L. Wu is with the Department of Information Engineering, National Chi Nan University, Puli, Nantou, Taiwan and the Department of Computer Science and Information Engineering, National Taiwan University, Taipei, Taiwan (e-mail: [email protected]).
Publisher Item Identifier S 1053-587X(97)03335-7.
Fig. 1. Relationship amongx(n); y(n); Y (k); Y0(k), and X(k).
orthogonal transform in the area of data compression, from the normal bases expansion of the DFT. A new interpretation of the forward and the inverse DCT transforms is then presented. Based on the proposed technique, a fully recursive fast DCT algorithm can be derived more easily in a systematic way.
II. A NEWINTERPRETATION OF THEDISCRETECOSINETRANSFORM
The DCT of anN-point real sequence x(n) is given by [2], [3] X(k) =
N01 n=0
x(n) cos k(2n + 1)2N ; k = 0; 1; . . . ; N 0 1 (1) where, for convenience, the1=p2 normalization factor for X(0) is not included. For simplicity of expression, assume thatN is even. Define a new sequencey(n) by permuting x(n) as follows:
y(n) = x(2n);x(2N 0 2n 0 1); n = N=2; . . . ; N 0 1:n = 0; 1; . . . ; N=2 0 1 (2) The DFT of y(n) is defined as
Y (k) =N01
n=0
y(n)!nk
N ; k = 0; 1; . . . ; N 0 1 (3)
where !N = e , and j = p01. Equations (2) and (3) can,
respectively, be represented in a more compact form as
fy(n)g = P fx(n)g n = 0; 1; . . . ; N 0 1 (4) fY (k)g = Ffy(n)g; n; k = 0; 1; . . . ; N 0 1: (5) It follows thatP and F are invertible. In order to have a fixed basis, Y (k) must be multiplied by an index k dependent weight factor to formY0(k). That is, the k-variable dependent weight factor will be embedded inY0(K) (the so-called weighted Fourier coefficient). In other words,Y0(k) can be represented in a more compact form as
fY0(k)g = MfY (k)g = Y (k)1!k0
4N ; k = 0; 1; . . . ; N 01: (6)
In order to have a better understanding and for the ease of expla-nation, the relations among x(n); y(n); X(k); Y (k); and Y0(k) are illustrated in Fig. 1. In the figure, the map between Y0(k) and X(k) is what we are looking for. By composition of maps, we will show that the so-called basefield transformB describes exactly the relation presented in (1). Note thatx(n) and X(k) reside in the basefield K (real field), whereasY (k) and Y0(k) are in the finite extension field F (complex field) of K.
1053–587X/97$10.001997 IEEE
IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 45, NO. 5, MAY 1997 1349
Since is a linear functional on FK, then there exists an 2 F such that [4]
() = Tr( ); 8 2 F: (7)
By applying the relation of (7) to the weighted Fourier coefficients Y0(k), it follows that
(Y0(k)) = Tr( Y0(k)) = Y0(k) + 3Y03(k);
0 k N 0 1: (8) As it was shown in [4], is a linear functional, and the map
: Y0(k) 7! X(k) = Tr( Y0(k)) (9)
defines a one-to-one correspondence betweenY0(k) and X(k). Since Y (k) is the Fourier transform of a real sequence, it satisfies the conjugacy relation, i.e.,Y3(k) = Y (N 0 k). In addition, it is also satisfied withY03(k) = Y0(N 0 k) because
Y0(N 0 k) = Y (N 0 k) 1 !(N0k)0 4N
= Y3(k) 1 !k04N 3= Y03(k): (10) Consider the action of on a conjugacy class fY0(k); Y0(N 0 k)g. X(k)= (Y0(k))=Tr( Y0(k))= Y0(k)+ 3Y03(k): (11)
X(N 0 k) = (Y0(N 0 k)) = Tr( Y0(N 0 k))
= 3Y0(k) + Y03(k): (12)
Notice that (11) and (12) show that the DCT coefficients X(k) andX(N 0 k) can be interpreted as the linear combination of the conjugacy class of the weighted Fourier coefficients with respect to the normal basis generated by . The corresponding matrix M for the projection of conjugacy class can then be obtained
X(k) X(N 0 k) M 1 Y 0(k) Y0(N 0 k) = 3 3 1 Y 0(k) Y0(N 0 k) (13)
where is equal to 12!4N. To have an invertible basefield transform, it must guarantee thatM01exists. In addition, it is easy to derive that
M01 3
3 01
= 3 3 (14)
where = 2 3 = !4N0 . In other words, the weighted Fourier coefficients can be obtained from the DCT coefficients as follows:
Y0(k)
Y0(N 0 k) = 3
3 1 X(N 0 k) :X(k) (15)
Since ; M; F; and P are bijective maps, the basefield transform B can be obtained as X(k) = Bfx(n)g = M F P fx(n)g = M Ffy(n)g = fY0(k)g = Tr( Y0(k)) = Y0(k) + 3Y03(k) = Real(2 Y0(k)) = Real 2 !4Nk0 N01 n=0 y(n)!nk N = Real !k4N N=201 n=0 x(2n)!Nnk + N01 n=N=2 x(2N 0 2n 0 1)!nk N = Real !k 4N n:even x(n)!nk=2 N + n:odd x(n)!(N0(n+1)=2)kN =Real n:even x(n)!(n=2+1=4)k N + n:odd x(n)!0(n=2+1=4)k N = N01 n=0 x(n) cos (2n + 1)kN (16)
where “” denotes the operation of composition. Equation (16) shows that the DCT coefficients can be obtained as a projection of the weighted Fourier coefficients from the extension field to the basefield. In addition, the corresponding inverse relation can be obtained by applyingB01(= P01F01M01 01) to X(k), that is, applying F01 M01 01 toX(k) to derive y(n) and then applying P01
to y(n) to obtain x(n). y(n) = 1 N N01 k=0 Y (k)!0nk N = 1N N01 k=0 Y0(k)! 0k 4N !0nkN = 1 N N01 k=0 ( X(k) + 3X(N 0 k))! 0k 4N !0nkN = 1N X(0)+ N01 k=1 X(k)!N0n+ k+!N0k 4N X(N 0k)!N0nk = 1N X(0) +N01 k=1 X(k)!N0n+ k+ X(k)!Nn+ k = 1N X(0) + N01 k=1 2X(k) cos (4n + 1)kN : ApplyingP01 to y(n), it follows that
x(n) = 1N X(0) + N01 k=1 2X(k) cos (2n + 1)kN = cn N01 k=0 X(k) cos (2n + 1)kN (17) where cn= 1 N; for n = 0 2 N; otherwise. (18) Equation (17) shows that the inverse DCT can be obtained by applying the same technique as the forward DCT. In addition, this result coincides with the one given in [2].
III. A NEWFULLYRECURSIVEFASTDCT ALGORITHM
In this section, a fast algorithm for the DCT is derived from the viewpoints of field projection. If we assume the input is a real sequence, then only the Fourier kernel part and the twiddle factors contain complex numbers in Y0(k). Therefore, only the
above two factors need to be projected from Y0(k) to X(k).
Based on the idea of [1], the coefficients of the kernel part and the twiddle factors are expanded with respect to different bases f 0; 1g and f 0; 1g = f1; 0ig, respectively. As pointed out
previously, the linear functionals and 01 can be expanded by using the normal basesf 0; 1g = f ; 3g = f12!N=24N ;12!0N=24N g
1350 IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 45, NO. 5, MAY 1997
andf 0; 1g = f ; 3g = f!4N0N=2; !4NN=2g, respectively. Thus, the
fast DCT algorithm can be derived as follows: Y (k) = N01 n=0 y(n)!nk N = N=201 n=0 y(n)!Nnk+ N=201 n=0 y(N 0 n 0 1)!(N0n01)kN = N=201 n=0 x(2n)!nk N + N=201 n=0 x(2n + 1)!0(n+1)k N :
By applying theM operator, we can map Y (k) into Y0(k) as follows: Y0(k) = N=201 n=0 x(2n)!Nnk+ 0 + N=201 n=0 x(2n+1)!0(n+1)k+ 0N =!0k 4N N=201 n=0 x(2n)!nk+ 0N + N=201 n=0 x(2n+1)!N0nk0 0 = 1 l=0 !0k(l) l N=201 n=0 x(2n) 1 j=0 !(j)nk+ 0 j + N=201 n=0 x(2n + 1) 1 j=0 !(j) 0nk0 0 j = 1 l;j=0 N=201 n=0 x(2n)!(l) 0k!nk+ 0(j) l j + 1 l;j=0 N=201 n=0 x(2n + 1)!0k(l)!0nk0 0(j) l j:
Projecting the above equation onto the real axis via yields X(k) = 1 l;j=0 N=201 n=0 x(2n)!(l) 0k!(j)nk+ 0 Tr( l j) + 1 l;j=0 N=201 n=0 x(2n + 1)!0k(l)!0nk0 0(j) Tr( l j) = N=201 n=0 x(2n) cos k2Ncos (4n + 2)k2N + sin k2N sin (4n + 2)k2N + N=201 n=0 x(2n + 1) cos k2Ncos (04n 0 2)k2N + sin k2N sin (04n 0 2)k2N = N=201 n=0 (x(2n) + x(2n + 1)) cos (2n + 1)k2 1 N=2 cos k2N + N=201 n=0 (x(2n) 0 x(2n + 1)) sin (2n + 1)k2 1 N=2 sin k2N: (19) Equation (19) shows that theN-point DCT can be decomposed into a N=2-point type-II DCT and N=2-point type-II discrete sine transform (DST) [5]. As shown in [5], theN-point type-II DCT and type-II DST require onlyO(N log2N) multiplications and additions. Thus,
(19) can be obtained inO(N log2N) multiplications and additions as well. If we set G(k) = N=201 n=0 (x(2n) + x(2n + 1)) cos (2n + 1)k2 1 N=2 and H(k) = N=201 n=0 (x(2n) 0 x(2n + 1)) sin (2n + 1)k2 1 N=2 then X(N 0 k) = (Y0(N 0 k)) = 0G(k) sin k 2N + H(k) cos k2N; k = 1; 2; . . . ; N=2 0 1: (20) Thus, an O(N log2N) fully recursive fast DCT algorithm can be obtained as X(k) X(N 0 k) = cos k 2N sin2Nk 0 sink 2N cos2Nk G(k) H(k) ; k = 1; 2; . . . ; N=2 0 1 (21) where X(0) = H(0) and X(N=2) = cos 4 1 H(N=2): IV. CONCLUSION
Normal bases expansion provides a general framework for con-structing the isomorphic mapping between the basefield and the extension field transforms. Based on this technique, in this correspon-dence, the DCT can be treated as the projection of the weighted DFT. Moreover, a new fast DCT algorithm can be obtained in a systematic way, and the convolution property of the basefield transform can be derived accordingly.
REFERENCES
[1] J. Hong, M. Vetterli, and P. Duhamel, “Basefield transforms with the convolution property,” Proc. IEEE, vol. 82, pp. 400–412, Mar. 1994. [2] J. Makhoul, “A fast cosine transform in one and two dimensions,” IEEE
Trans. Acoust., Speech, Signal Processing, vol. ASSP-28, pp. 27–34, Feb. 1980.
[3] H. Malvar, “Fast computation of discrete cosine transform through fast Hartley transform,” Electron. Lett., vol. 22, no. 7, pp. 352–353, Mar. 1986.
[4] R. Lidl and H. Niederreiter, Introduction to Finite Fields and Their Applications. New York: Cambridge Univ. Press, 1986.
[5] K. R. Rao and P. Yip, Discrete Cosine Transform—Algorithms, Advan-tages Applications. New York: Academic, 1992.
