- Email: [email protected]
Transform domain technique for windowing the DCT and DST
Journal of the Franklin Institute 339 (2002) 111–120
Brief Communication
Transform domain technique for windowing the DCT and DST$ B.G. Sherlock*, Y.P. Kakad Department of Electrical and Computer Engineering, University of North Carolina at Charlotte, 9201 University City Boulevard, Charlotte, NC 28223, USA Received 5 April 2000; received in revised form 4 October 2000; accepted 17 October 2001
Abstract In this paper, an algorithm is developed to apply Hann, Hamming, Blackman and related windows directly in the transform domain for the discrete cosine transform and discrete sine transform. These algorithms are useful in applications where windowing is required in order to minimize edge effects caused by implicit symmetries in the transform domain that are not replicated in the real-world data. Examples of such applications include data communication, adaptive system identification and filtering, real-time analysis of financial market data, etc. Software implementations in C language are also given. r 2002 The Franklin Institute. Published by Elsevier Science Ltd. All rights reserved. Keywords: Discrete cosine transform; Discrete sine transform; Fast transform; Windowing; Real-time signal processing
1. Introduction Since its introduction in 1974 [1] the discrete cosine transform (DCT) has found widespread application in digital signal and image processing for the purposes of pattern recognition, Wiener filtering, data compression, communication and several other areas. The reason for this success is the ability of the DCT to decorrelate signal data; specifically, the DCT has been shown to provide a close approximation to the ideal Karhunen–-Loe" ve transform (KLT) for data that exhibits Markov-1 statistics $
This research was supported in part by the Office of Naval Research Grant Number N00014-99-10091. The material in this paper appears in preliminary form in conference CSCC 2000, Athens, July 2000. *Corresponding author. E-mail addresses: [email protected] (B.G. Sherlock), [email protected] (Y.P. Kakad). 0016-0032/02/$22.00 r 2002 The Franklin Institute. Published by Elsevier Science Ltd. All rights reserved. PII: S 0 0 1 6 - 0 0 3 2 ( 0 1 ) 0 0 0 5 8 - 8
B.G. Sherlock, Y.P. Kakad / Journal of the Franklin Institute 339 (2002) 111–120
112
with a correlation coefficient r > 0:9 [1,6]. The KLT is known to be optimal with respect to the following performance measures: variance distribution [2], estimation using the mean-square error criterion [3,4], and the rate distortion function [5]. Although the KLT is optimal, it is a signal-dependent transform and there is no general algorithm for its fast computation. This paper investigates the computational properties of the DCT and discrete sine transform (DST) as applied to signal data which are already in the respective transform domains and has to be windowed. In numerous applications, particularly for real-time data processing, it is impossible or inconvenient to form the DCT of the entire signal. Rather, a suitable transform length N is chosen, and the first N data points are transformed. In practice, it is usually customary to multiply the data by a window function [7,8] before the DCT and DST are taken. Windowing serves to reduce edge effects which occur because the DCT treats its input data as having symmetries and periodicities that do not actually occur in the data. In this paper, algorithms enabling the application of Hann, Hamming, Blackman and related windows in the DCT and DST transform domains are developed. This is a computationally efficient technique whenever a variety of windows must be applied. Furthermore, the algorithms can be used to update the transforms of shifting data in real time in the presence of windowing, when used in conjunction with DCT/DST updating algorithms developed earlier by the authors [9–12].
2. Definitions of DCT and DST A non-windowed signal may be regarded as windowed by the ‘‘boxcar’’ window wðxÞ ¼ 1 since, multiplication by this window does not alter the signal data. The DCT of a signal f ðxÞ of length N is defined by rffiffiffiffiffi N1 X 2 ð2x þ 1Þ up CðuÞ ¼ ; f ðxÞ cos aðuÞ N 2N x¼0
ð1Þ
where x ¼ 0; 1; 2; :::; N 1; u ¼ 0; 1; 2; :::; N 1; ( pffiffiffi 1= 2 if u mod N ¼ 0; aðuÞ ¼ 1 otherwisel: The DST was introduced in 1978 [13,14] with the following definition for a signal f ðxÞ of length N: rffiffiffiffiffi N 1 X 2 ð2x þ 1Þðu þ 1Þp SA ðuÞ ¼ ; f ðxÞ sin bðuÞ N 2N x¼0
B.G. Sherlock, Y.P. Kakad / Journal of the Franklin Institute 339 (2002) 111–120
113
where x ¼ 0; 1; 2; :::; N 1; u ¼ 0; 1; 2; :::; N 1; ( pffiffiffi 1= 2 if u ¼ N 1; bðuÞ ¼ 1 otherwise: In some of the more recent literature [6,9], the definition of the DST is written in the following equivalent form: rffiffiffiffiffi N X 2 ð2x 1Þ up ; SB ðuÞ ¼ f ðxÞ sin aðuÞ N 2N x¼1 where x ¼ 1; 2; :::; N; u ¼ 1; 2; :::; N; and aðuÞ has the same definition as for the DCT. This paper will make use of a third equivalent form of the DST, expressed as: rffiffiffiffiffi N 1 X 2 ð2x þ 1Þ up ; ð2Þ SðuÞ ¼ f ðxÞ sin aðuÞ N 2N x¼0 where x ¼ 0; 1; 2; ; :::; N 1; u ¼ 1; 2; 3; :::; N and aðuÞ has the same definition as for the DCT. Note that the above three alternate forms represent precisely the same transform, namely the DST of type II [6,14]. The differences between SA ðuÞ; SB ðuÞ and SðuÞ simply indicate different choices of range for the spatial and frequency indices. We have chosen to use the form expressed in Eq. (2) because the argument of the trigonometrical function is the same as the argument found in the DCT of Eq. (1). This results in a notational simplification both of the derivations and of the final expressions. In the following section, an algorithm is developed for the update of DCT and DST sequence to reflect r additional data points and the removal of r old points from the signal.
3. Windowing in the transform domain The N point DCT treats data as if it had the following reflective symmetries and periodicities: f ðxÞ ¼ f ðxÞ ¼ f ð2N 1 xÞ ¼ f ðx þ 2NÞ: The corresponding symmetries for the DST are: f ðxÞ ¼ f ðxÞ ¼ f ð2N 1 xÞ ¼ f ðx þ 2NÞ: In practice, most real-world signals do not exhibit the above symmetries. If the ‘‘boxcar’’ window is applied to such a signal, the DST will treat it as if there were discontinuities at its edges. (For the DCT there will not be spurious discontinuities, but less serious edge effects will still occur.) Ringing effects near the edges of filtered signals may occur as a result of such spurious discontinuities [15]. Such effects can be
114
B.G. Sherlock, Y.P. Kakad / Journal of the Franklin Institute 339 (2002) 111–120
reduced by applying a more appropriate window. In addition to selecting a portion of the signal, the window modifies this portion to make or tend to make it continuous at the edges when regarded as periodically repeated, or smoother, under the implicit symmetries. Several types of windows have been described in the literature [16,8].
3.1. Algorithms for Hamming, Hann and Blackman windows These windows are defined for x ¼ 0; ::::; N 1 as follows [8]: Hann :
wðxÞ ¼ 0:5ð1 cosð2px=NÞÞ
Hamming : wðxÞ ¼ 0:54 0:46 cosð2px=NÞ
Blackman : wðxÞ ¼ 0:42 0:5 cosð2px=NÞ þ 0:08 cosð4px=NÞ: A general algorithm is developed for all these three windows using the following generic window: wðxÞ ¼ a0 þ a1 cosð2px=NÞ þ a2 cosð4px=NÞ: Thus, given the signal f ðxÞ; x ¼ 0; ::::; N 1; let the windowed signal be fw ðxÞ ¼ wðxÞf ðxÞ 2px 4px þ a2 cos f ðxÞ: ¼ a0 þ a1 cos N N Taking the DCT we have rffiffiffiffiffi N 1 X 2 2px ð2x þ 1Þup cos f ðxÞ cos Cw ðuÞ ¼ a0 CðuÞ þ a1 aðuÞ N N 2N x¼0 rffiffiffiffiffi N 1 X 2 4px ð2x þ 1Þup cos : þ a2 f ðxÞ cos aðuÞ N N 2N x¼0
ð3Þ
B.G. Sherlock, Y.P. Kakad / Journal of the Franklin Institute 339 (2002) 111–120
115
The second term on the right-hand side of Eq. (3), using the definitions of CðuÞ and SðuÞ and trigonometric identities, can be simplified as rffiffiffiffiffi N 1 X 2 2px ð2x þ 1Þup cos f ðxÞ cos a1 aðuÞ N N 2N x¼0 rffiffiffiffiffi
N 1 X 2 1 2px ð2x þ 1Þup ¼ a1 f ðxÞ cos aðuÞ N 2 N 2N x¼0 2px ð2x þ 1Þup þ þcos N 2N rffiffiffiffiffi
N 1 X 1 2 ð2x þ 1Þðu 2Þp þ 2p f ðxÞ cos ¼ a1 aðuÞ 2 N 2N x¼0 ð2x þ 1Þðu þ 2Þp 2p þcos 2N " rffiffiffiffiffi 1 X 1 p 2N ð2x þ 1Þðu 2Þp f ðxÞ cos ¼ a1 aðuÞ cos 2 N N x¼0 2N rffiffiffiffiffi N 1 p 2X ð2x þ 1Þðu 2Þp f ðxÞ sin sin N N x¼0 2N rffiffiffiffiffi N 1 p 2X ð2x þ 1Þðu þ 2Þp þ cos f ðxÞ cos N N x¼0 2N # rffiffiffiffiffi N 1 p 2X ð2x þ 1Þðu þ 2Þp þsin f ðxÞ sin N N x¼0 2N h 1 p Cðu 2Þ p Sðu 2Þ sin ¼ a1 aðuÞ cos 2 N aðu 2Þ N aðu 2Þ p Cðu þ 2Þ p Sðu þ 2Þ þ sin þ cos : N aðu þ 2Þ N aðu þ 2Þ Similarly, the third term on the right-hand side is simplified as rffiffiffiffiffi N 1 X 2 4px xð2x þ 1Þ up cos f ðxÞ cos aðuÞ a2 N N 2N x¼0 1 2p Cðu þ 4Þ 2p Sðu þ 4Þ þ sin ¼ a2 aðuÞ cos 2 N aðu þ 4Þ N aðu þ 4Þ 2p Cðu 4Þ 2p Sðu 4Þ þ cos sin : N aðu 4Þ N aðu 4Þ
ð4Þ
ð5Þ
116
B.G. Sherlock, Y.P. Kakad / Journal of the Franklin Institute 339 (2002) 111–120
Combining Eqs. (4) and (5) yields the final algorithm: 8 9 p Cðu þ 2Þ Cðu 2Þ > > > > > = < cos N aðu þ 2Þ þ aðu 2Þ > a1 Cw ðuÞ ¼ a0 CðuÞ þ aðuÞ > 2 p Sðu þ 2Þ Sðu 2Þ > > > > > þsin ; : N aðu þ 2Þ aðu 2Þ 8 9 2p Cðu þ 4Þ Cðu 4Þ > > > > > > þ cos < N aðu þ 4Þ aðu 4Þ = a2 ; u ¼ 0; 1; ::::; N 1: ð6Þ þ aðuÞ > 2 2p Sðu þ 4Þ Sðu 4Þ > > > > > þsin ; : N aðu þ 4Þ aðu 4Þ Analogously, we may obtain the following result for the DST: 8 9 p Sðu þ 2Þ Sðu 2Þ > > > > > = < cos N aðu þ 2Þ þ aðu 2Þ > a1 Sw ðuÞ ¼ a0 SðuÞ þ aðuÞ > 2 p Cðu 2Þ Cðu þ 2Þ > > > > > ; : þsin N aðu 2Þ aðu þ 2Þ 9 8 2p Sðu þ 4Þ Sðu 4Þ > > > > > > þ cos < N aðu þ 4Þ aðu 4Þ = a2 ; u ¼ 1; 2; ::::; N: þ aðuÞ > 2 2p Cðu 4Þ Cðu þ 4Þ > > > > > þsin ; : N aðu 4Þ aðu þ 4Þ
ð7Þ
Therefore, the DCT and DST of the non-windowed signal defined by Eqs. (1) and (2), respectively are used in applying the above equations to give the DCT and DST of the windowed signal. In applying Eqs. (6) and (7), terms such as Sðu þ 4Þ will result in arguments out of the range of u specified in the definitions of SðuÞ and CðuÞ: These may be resolved by applying the symmetries of the DCT and DST which are listed below for all integer u: CðuÞ ¼ CðuÞ ¼ Cð2N þ uÞ ¼ Cð2N uÞ ¼ Cð4N uÞ ¼ Cð4N þ uÞ; CðN uÞ ¼ CðN þ uÞ;
CðNÞ ¼ 0;
SðuÞ ¼ SðuÞ ¼ Sð2N uÞ ¼ Sð2N þ uÞ ¼ Sð4N uÞ ¼ Sð4N þ uÞ; SðN þ uÞ ¼ SðN uÞ;
Sð0Þ ¼ 0:
3.2. Implementation C language procedure implementing the update as described in Eqs. (6) and (7) is given in Figs. 1 and 2. This implementation is intended only to illustrate the algorithm and is not optimized for processing speed. A more efficient implementation would extend the DCT and DST once and for all before entering the main loop, rather than calculating the extended coefficients repeatedly for each value of u:
B.G. Sherlock, Y.P. Kakad / Journal of the Franklin Institute 339 (2002) 111–120
117
double extendDCT2(double C[], int u, int N){ // Given the DCT2 in array C[0..N-1], this function uses // the symmetries of the definition to extend the values // to all integer u, not just 0..N-1. int twoN,fourN, sign = 1; twoN=N+N; fourN= twoN+twoN; if(u<0)u = -u; u = u % fourN; if (u>twoN) u = fourN-u; if(u>N){ u = twoN-u; sign = -sign; } if (u==N){return(0.0);}else{return(sign*C[u]); } } double extendDST2(double S[], int u, int N){ // Given the DST2 in array S[1..N], this function uses // the symmetries of the definition to extend the values // to all integer u, not just 1..N. int twoN,fourN,sign = 1; twoN=N+N; fourN= twoN+twoN; if(u<0){ u = -u; sign = -sign;} u = u % fourN; if (u>=twoN){u -= twoN;sign = -sign;} if(u>N) u=twoN-u; if(u==0){return(0.0);}else{return(sign*S[u]);} } Fig. 1. C language procedures to extend the DCT2 and DST2 beyond the usual index ranges. These are used by the windowing procedure of Fig. 2.
4. Computational complexity The algorithms proposed have computational complexity of order N: The number of arithmetic operations required to window both the DCT and DST in the transform domain using the proposed algorithms is:
Hann :
16N þ 8 operations ð8N additions; 8N þ 8 multiplicationsÞ;
Hamming or Blackman :
30N þ 16 operations ð16N additions; 14N þ 8 multiplicationsÞ:
118
B.G. Sherlock, Y.P. Kakad / Journal of the Franklin Institute 339 (2002) 111–120
window_transform(double C[],double S[],double Cw[],double Sw[],int N, double a0, double a1, double a2){ // converts the DCT2 and DST2 of a length N signal into the // DCT2 and DST2 of the same signal with a window applied. // The window is w(x) = a0 + a1 cos 2 pi x/N + a2 cos 4 pi x/N // // C[0..N-1] DCT to be windowed (input) // S[1..N] DST to be windowed (input) // Cw[0..N-1] DCT windowed (output) // Sw[1..N] DST windowed (output) // N signal and transform length // a0,a1,a2 coefficients of window // // For Blackman window, a0=0.42, a1=-0.5, a2=0.08 // For Hamming window, a0=0.54, a1=-0.46, a2=0 // For Hann window, a0=0.5, a1=-0.5, a2=0 // int u; double double double double
Cum2,Cup2,Cum4,Cup4,Sum2,Sup2,Sum4,Sup4; cos_pivN,sin_pivN, root2, invroot2,cos_2pivN,sin_2pivN; Cterm1,Cterm2,Cterm3,Cterm4,Call4; Sterm1,Sterm2,Sterm3,Sterm4,Sall4;
cos_pivN = cos(M_PI/N); cos_2pivN = cos(2.0*M_PI/N); sin_pivN = sin(M_PI/N); sin_2pivN = sin(2.0*M_PI/N); root2=sqrt(2.0); invroot2=1.0/sqrt(2.0); for (u=0; u<=N;u++){ Cum2= extendDCT2(C,u-2,N); Cup2= extendDCT2(C,u+2,N); Sum2= extendDST2(S,u-2,N); Sup2= extendDST2(S,u+2,N); Cum4= extendDCT2(C,u-4,N); Cup4= extendDCT2(C,u+4,N); Sum4= extendDST2(S,u-4,N); Sup4= extendDST2(S,u+4,N); Cterm1 = cos_pivN*Cup2 + sin_pivN*Sup2; Cterm2 = cos_pivN*Cum2 - sin_pivN*Sum2; Cterm3 = cos_2pivN*Cup4 + sin_2pivN*Sup4; Cterm4 = cos_2pivN*Cum4 - sin_2pivN*Sum4; Sterm1 = cos_pivN*Sup2 - sin_pivN*Cup2; Sterm2 = cos_pivN*Sum2 + sin_pivN*Cum2; Sterm3 = cos_2pivN*Sup4 - sin_2pivN*Cup4; Sterm4 = cos_2pivN*Sum4 + sin_2pivN*Cum4; if(u== N-2){Cterm1 *= root2; Sterm1 *= root2;} if(u== 2){ Cterm2 *= root2; Sterm2 *= root2;} if(u== N-4){Cterm3 *= root2; Sterm3 *= root2;} if(u== 4){ Cterm4 *= root2; Sterm4 *= root2;} Call4 = (a1*(Cterm1+Cterm2) + a2*(Cterm3+Cterm4))/2.0; Sall4 = (a1*(Sterm1+Sterm2) + a2*(Sterm3+Sterm4))/2.0; if ((u==0)||(u==N)){Call4 *= invroot2; Sall4 *= invroot2;} Cw[u]= a0*C[u] + Call4; Sw[u]= a0*S[u] + Sall4; } } Fig. 2. C language implementation of the algorithm to perform windowing in the discrete cosine and sine transform domains, for the Blackman, Hamming, or Hann window. This procedure calls the subroutines of Fig. 1 to extend the DCT2 and DST2 beyond the usual index range.
B.G. Sherlock, Y.P. Kakad / Journal of the Franklin Institute 339 (2002) 111–120
119
The conventional approach requires N multiplications to apply the window function, and 2N log2 N operations for each transform if the fastest known DCT/ DST algorithms are used, where N is a power of 2. This yields a total of N þ 4N log2 N operations, for an overall computational complexity of order N log2 N: It is easily seen that our proposed approach is more efficient than the conventional approach when the transform length N is at least 4 for the Hann window, or 256 for the Hamming/Blackman windows. Here we have quoted the smallest powers of 2 for which our approach gives superior performance. For values of N that are not powers of 2, the computational advantage offered by our approach will be considerably greater, because of the additional computation required to calculate DCT/DST for such signal lengths. These algorithms would be of particular interest in applications where frequent updates are required in the presence of a sliding window, e.g. where new data points arrive at frequent intervals, and it is desired to maintain the DCT/DST of the N most recent data points.
5. Conclusion General, computationally efficient algorithms for applying Hann, Hamming, Blackman, and related windows to the DCT and DST directly in the transform domain data were analytically derived. These can be used in applications where the transforms must be updated in the presence of a window function, as new data is received in real time. Software implementations which demonstrate the correctness of developed algorithms in C language are also provided.
Acknowledgements This research was supported in part by the Office of Naval Research Grant Number N00014-99-1-0091. The authors would like to express their gratitude for the useful and detailed evaluation of the manuscript by the anonymous reviewers. The clarity and completeness of this paper has been substantially improved as a result of their suggestions.
References [1] N. Ahmed, T. Natarajan, K.R. Rao, Discrete cosine transform, IEEE Trans. Comput. C-23 (1974) 90–93. [2] H.C. Andrews, Multidimensional rotations in feature selection, IEEE Trans. Comput. C-20 (1971) 1045–1051. [3] W.K. Pratt, Generalized Wiener filtering computation techniques, IEEE Trans. Comput. C-21 (1972) 636–641.
120
B.G. Sherlock, Y.P. Kakad / Journal of the Franklin Institute 339 (2002) 111–120
[4] J. Pearl, Walsh processing of random signals, IEEE Trans. Electromagn. Components EMC-13 (1971) 137–141. [5] J. Pearl, H.C. Andrews, W.K. Pratt, Performance measures for transform data coding, IEEE Trans. Commun. Technol. COM-20 (1972) 411–415. [6] K.R. Rao, P. Yip, Discrete Cosine Transform: Algorithms, Advantages, Applications, Academic Press, Boston, 1990. [7] H. Kwakernaak, R. Sivan, Modern Signals and Systems, Prentice-Hall, Englewood Cliffs, NJ, 1991. [8] F.J. Harris, On the use of windows for harmonic analysis with the discrete fourier transform, Proc. IEEE 66 (1978). [9] P. Yip, K.R. Rao, On the shift property of DCT’s and DST’s, IEEE Trans. Acoust. Speech Signal Process. ASSP-35 (1987) 404–406. [10] B.G. Sherlock, D.M. Monro, Moving fast fourier transform, IEE Proc. F 139 (4) (1992) 279–282. [11] B.G. Sherlock, Windowed discrete fourier transform for shifting data, Signal Process. 74 (1999) 169– 177. [12] B.G. Sherlock, Y.P. Kakad, Windowing the discrete cosine transform in the transform domain, Fourth World Multiconference on Circuits, Systems, Communications and Computers (CSCC 2000), Athens, Greece, July 10–15, 2000, pp. 320–324. [13] H.B. Kekre, J.K. Solanki, Comparative performance of various trigonometric unitary transforms for transform image coding, Int. J. Electron. 44(3) (1978) 305–315. [14] Z. Wang, A fast algorithm for the discrete sine transform implemented by the fast cosine transform, IEEE Trans. Acoust. Speech Signal Process. ASSP-30 (1982). [15] R.C. Gonzalez, R.E. Woods, Digital Image Processing, Addison-Wesley, Reading, MA, 1992. [16] N.C. Geckinli, D. Yavuz, Some novel windows and a concise tutorial comparison of window families, IEEE Trans. Acoust. Speech Signal Process. ASSP-26 (1978) 501–507.
Brief Communication
Transform domain technique for windowing the DCT and DST$ B.G. Sherlock*, Y.P. Kakad Department of Electrical and Computer Engineering, University of North Carolina at Charlotte, 9201 University City Boulevard, Charlotte, NC 28223, USA Received 5 April 2000; received in revised form 4 October 2000; accepted 17 October 2001
Abstract In this paper, an algorithm is developed to apply Hann, Hamming, Blackman and related windows directly in the transform domain for the discrete cosine transform and discrete sine transform. These algorithms are useful in applications where windowing is required in order to minimize edge effects caused by implicit symmetries in the transform domain that are not replicated in the real-world data. Examples of such applications include data communication, adaptive system identification and filtering, real-time analysis of financial market data, etc. Software implementations in C language are also given. r 2002 The Franklin Institute. Published by Elsevier Science Ltd. All rights reserved. Keywords: Discrete cosine transform; Discrete sine transform; Fast transform; Windowing; Real-time signal processing
1. Introduction Since its introduction in 1974 [1] the discrete cosine transform (DCT) has found widespread application in digital signal and image processing for the purposes of pattern recognition, Wiener filtering, data compression, communication and several other areas. The reason for this success is the ability of the DCT to decorrelate signal data; specifically, the DCT has been shown to provide a close approximation to the ideal Karhunen–-Loe" ve transform (KLT) for data that exhibits Markov-1 statistics $
This research was supported in part by the Office of Naval Research Grant Number N00014-99-10091. The material in this paper appears in preliminary form in conference CSCC 2000, Athens, July 2000. *Corresponding author. E-mail addresses: [email protected] (B.G. Sherlock), [email protected] (Y.P. Kakad). 0016-0032/02/$22.00 r 2002 The Franklin Institute. Published by Elsevier Science Ltd. All rights reserved. PII: S 0 0 1 6 - 0 0 3 2 ( 0 1 ) 0 0 0 5 8 - 8
B.G. Sherlock, Y.P. Kakad / Journal of the Franklin Institute 339 (2002) 111–120
112
with a correlation coefficient r > 0:9 [1,6]. The KLT is known to be optimal with respect to the following performance measures: variance distribution [2], estimation using the mean-square error criterion [3,4], and the rate distortion function [5]. Although the KLT is optimal, it is a signal-dependent transform and there is no general algorithm for its fast computation. This paper investigates the computational properties of the DCT and discrete sine transform (DST) as applied to signal data which are already in the respective transform domains and has to be windowed. In numerous applications, particularly for real-time data processing, it is impossible or inconvenient to form the DCT of the entire signal. Rather, a suitable transform length N is chosen, and the first N data points are transformed. In practice, it is usually customary to multiply the data by a window function [7,8] before the DCT and DST are taken. Windowing serves to reduce edge effects which occur because the DCT treats its input data as having symmetries and periodicities that do not actually occur in the data. In this paper, algorithms enabling the application of Hann, Hamming, Blackman and related windows in the DCT and DST transform domains are developed. This is a computationally efficient technique whenever a variety of windows must be applied. Furthermore, the algorithms can be used to update the transforms of shifting data in real time in the presence of windowing, when used in conjunction with DCT/DST updating algorithms developed earlier by the authors [9–12].
2. Definitions of DCT and DST A non-windowed signal may be regarded as windowed by the ‘‘boxcar’’ window wðxÞ ¼ 1 since, multiplication by this window does not alter the signal data. The DCT of a signal f ðxÞ of length N is defined by rffiffiffiffiffi N1 X 2 ð2x þ 1Þ up CðuÞ ¼ ; f ðxÞ cos aðuÞ N 2N x¼0
ð1Þ
where x ¼ 0; 1; 2; :::; N 1; u ¼ 0; 1; 2; :::; N 1; ( pffiffiffi 1= 2 if u mod N ¼ 0; aðuÞ ¼ 1 otherwisel: The DST was introduced in 1978 [13,14] with the following definition for a signal f ðxÞ of length N: rffiffiffiffiffi N 1 X 2 ð2x þ 1Þðu þ 1Þp SA ðuÞ ¼ ; f ðxÞ sin bðuÞ N 2N x¼0
B.G. Sherlock, Y.P. Kakad / Journal of the Franklin Institute 339 (2002) 111–120
113
where x ¼ 0; 1; 2; :::; N 1; u ¼ 0; 1; 2; :::; N 1; ( pffiffiffi 1= 2 if u ¼ N 1; bðuÞ ¼ 1 otherwise: In some of the more recent literature [6,9], the definition of the DST is written in the following equivalent form: rffiffiffiffiffi N X 2 ð2x 1Þ up ; SB ðuÞ ¼ f ðxÞ sin aðuÞ N 2N x¼1 where x ¼ 1; 2; :::; N; u ¼ 1; 2; :::; N; and aðuÞ has the same definition as for the DCT. This paper will make use of a third equivalent form of the DST, expressed as: rffiffiffiffiffi N 1 X 2 ð2x þ 1Þ up ; ð2Þ SðuÞ ¼ f ðxÞ sin aðuÞ N 2N x¼0 where x ¼ 0; 1; 2; ; :::; N 1; u ¼ 1; 2; 3; :::; N and aðuÞ has the same definition as for the DCT. Note that the above three alternate forms represent precisely the same transform, namely the DST of type II [6,14]. The differences between SA ðuÞ; SB ðuÞ and SðuÞ simply indicate different choices of range for the spatial and frequency indices. We have chosen to use the form expressed in Eq. (2) because the argument of the trigonometrical function is the same as the argument found in the DCT of Eq. (1). This results in a notational simplification both of the derivations and of the final expressions. In the following section, an algorithm is developed for the update of DCT and DST sequence to reflect r additional data points and the removal of r old points from the signal.
3. Windowing in the transform domain The N point DCT treats data as if it had the following reflective symmetries and periodicities: f ðxÞ ¼ f ðxÞ ¼ f ð2N 1 xÞ ¼ f ðx þ 2NÞ: The corresponding symmetries for the DST are: f ðxÞ ¼ f ðxÞ ¼ f ð2N 1 xÞ ¼ f ðx þ 2NÞ: In practice, most real-world signals do not exhibit the above symmetries. If the ‘‘boxcar’’ window is applied to such a signal, the DST will treat it as if there were discontinuities at its edges. (For the DCT there will not be spurious discontinuities, but less serious edge effects will still occur.) Ringing effects near the edges of filtered signals may occur as a result of such spurious discontinuities [15]. Such effects can be
114
B.G. Sherlock, Y.P. Kakad / Journal of the Franklin Institute 339 (2002) 111–120
reduced by applying a more appropriate window. In addition to selecting a portion of the signal, the window modifies this portion to make or tend to make it continuous at the edges when regarded as periodically repeated, or smoother, under the implicit symmetries. Several types of windows have been described in the literature [16,8].
3.1. Algorithms for Hamming, Hann and Blackman windows These windows are defined for x ¼ 0; ::::; N 1 as follows [8]: Hann :
wðxÞ ¼ 0:5ð1 cosð2px=NÞÞ
Hamming : wðxÞ ¼ 0:54 0:46 cosð2px=NÞ
Blackman : wðxÞ ¼ 0:42 0:5 cosð2px=NÞ þ 0:08 cosð4px=NÞ: A general algorithm is developed for all these three windows using the following generic window: wðxÞ ¼ a0 þ a1 cosð2px=NÞ þ a2 cosð4px=NÞ: Thus, given the signal f ðxÞ; x ¼ 0; ::::; N 1; let the windowed signal be fw ðxÞ ¼ wðxÞf ðxÞ 2px 4px þ a2 cos f ðxÞ: ¼ a0 þ a1 cos N N Taking the DCT we have rffiffiffiffiffi N 1 X 2 2px ð2x þ 1Þup cos f ðxÞ cos Cw ðuÞ ¼ a0 CðuÞ þ a1 aðuÞ N N 2N x¼0 rffiffiffiffiffi N 1 X 2 4px ð2x þ 1Þup cos : þ a2 f ðxÞ cos aðuÞ N N 2N x¼0
ð3Þ
B.G. Sherlock, Y.P. Kakad / Journal of the Franklin Institute 339 (2002) 111–120
115
The second term on the right-hand side of Eq. (3), using the definitions of CðuÞ and SðuÞ and trigonometric identities, can be simplified as rffiffiffiffiffi N 1 X 2 2px ð2x þ 1Þup cos f ðxÞ cos a1 aðuÞ N N 2N x¼0 rffiffiffiffiffi
N 1 X 2 1 2px ð2x þ 1Þup ¼ a1 f ðxÞ cos aðuÞ N 2 N 2N x¼0 2px ð2x þ 1Þup þ þcos N 2N rffiffiffiffiffi
N 1 X 1 2 ð2x þ 1Þðu 2Þp þ 2p f ðxÞ cos ¼ a1 aðuÞ 2 N 2N x¼0 ð2x þ 1Þðu þ 2Þp 2p þcos 2N " rffiffiffiffiffi 1 X 1 p 2N ð2x þ 1Þðu 2Þp f ðxÞ cos ¼ a1 aðuÞ cos 2 N N x¼0 2N rffiffiffiffiffi N 1 p 2X ð2x þ 1Þðu 2Þp f ðxÞ sin sin N N x¼0 2N rffiffiffiffiffi N 1 p 2X ð2x þ 1Þðu þ 2Þp þ cos f ðxÞ cos N N x¼0 2N # rffiffiffiffiffi N 1 p 2X ð2x þ 1Þðu þ 2Þp þsin f ðxÞ sin N N x¼0 2N h 1 p Cðu 2Þ p Sðu 2Þ sin ¼ a1 aðuÞ cos 2 N aðu 2Þ N aðu 2Þ p Cðu þ 2Þ p Sðu þ 2Þ þ sin þ cos : N aðu þ 2Þ N aðu þ 2Þ Similarly, the third term on the right-hand side is simplified as rffiffiffiffiffi N 1 X 2 4px xð2x þ 1Þ up cos f ðxÞ cos aðuÞ a2 N N 2N x¼0 1 2p Cðu þ 4Þ 2p Sðu þ 4Þ þ sin ¼ a2 aðuÞ cos 2 N aðu þ 4Þ N aðu þ 4Þ 2p Cðu 4Þ 2p Sðu 4Þ þ cos sin : N aðu 4Þ N aðu 4Þ
ð4Þ
ð5Þ
116
B.G. Sherlock, Y.P. Kakad / Journal of the Franklin Institute 339 (2002) 111–120
Combining Eqs. (4) and (5) yields the final algorithm: 8 9 p Cðu þ 2Þ Cðu 2Þ > > > > > = < cos N aðu þ 2Þ þ aðu 2Þ > a1 Cw ðuÞ ¼ a0 CðuÞ þ aðuÞ > 2 p Sðu þ 2Þ Sðu 2Þ > > > > > þsin ; : N aðu þ 2Þ aðu 2Þ 8 9 2p Cðu þ 4Þ Cðu 4Þ > > > > > > þ cos < N aðu þ 4Þ aðu 4Þ = a2 ; u ¼ 0; 1; ::::; N 1: ð6Þ þ aðuÞ > 2 2p Sðu þ 4Þ Sðu 4Þ > > > > > þsin ; : N aðu þ 4Þ aðu 4Þ Analogously, we may obtain the following result for the DST: 8 9 p Sðu þ 2Þ Sðu 2Þ > > > > > = < cos N aðu þ 2Þ þ aðu 2Þ > a1 Sw ðuÞ ¼ a0 SðuÞ þ aðuÞ > 2 p Cðu 2Þ Cðu þ 2Þ > > > > > ; : þsin N aðu 2Þ aðu þ 2Þ 9 8 2p Sðu þ 4Þ Sðu 4Þ > > > > > > þ cos < N aðu þ 4Þ aðu 4Þ = a2 ; u ¼ 1; 2; ::::; N: þ aðuÞ > 2 2p Cðu 4Þ Cðu þ 4Þ > > > > > þsin ; : N aðu 4Þ aðu þ 4Þ
ð7Þ
Therefore, the DCT and DST of the non-windowed signal defined by Eqs. (1) and (2), respectively are used in applying the above equations to give the DCT and DST of the windowed signal. In applying Eqs. (6) and (7), terms such as Sðu þ 4Þ will result in arguments out of the range of u specified in the definitions of SðuÞ and CðuÞ: These may be resolved by applying the symmetries of the DCT and DST which are listed below for all integer u: CðuÞ ¼ CðuÞ ¼ Cð2N þ uÞ ¼ Cð2N uÞ ¼ Cð4N uÞ ¼ Cð4N þ uÞ; CðN uÞ ¼ CðN þ uÞ;
CðNÞ ¼ 0;
SðuÞ ¼ SðuÞ ¼ Sð2N uÞ ¼ Sð2N þ uÞ ¼ Sð4N uÞ ¼ Sð4N þ uÞ; SðN þ uÞ ¼ SðN uÞ;
Sð0Þ ¼ 0:
3.2. Implementation C language procedure implementing the update as described in Eqs. (6) and (7) is given in Figs. 1 and 2. This implementation is intended only to illustrate the algorithm and is not optimized for processing speed. A more efficient implementation would extend the DCT and DST once and for all before entering the main loop, rather than calculating the extended coefficients repeatedly for each value of u:
B.G. Sherlock, Y.P. Kakad / Journal of the Franklin Institute 339 (2002) 111–120
117
double extendDCT2(double C[], int u, int N){ // Given the DCT2 in array C[0..N-1], this function uses // the symmetries of the definition to extend the values // to all integer u, not just 0..N-1. int twoN,fourN, sign = 1; twoN=N+N; fourN= twoN+twoN; if(u<0)u = -u; u = u % fourN; if (u>twoN) u = fourN-u; if(u>N){ u = twoN-u; sign = -sign; } if (u==N){return(0.0);}else{return(sign*C[u]); } } double extendDST2(double S[], int u, int N){ // Given the DST2 in array S[1..N], this function uses // the symmetries of the definition to extend the values // to all integer u, not just 1..N. int twoN,fourN,sign = 1; twoN=N+N; fourN= twoN+twoN; if(u<0){ u = -u; sign = -sign;} u = u % fourN; if (u>=twoN){u -= twoN;sign = -sign;} if(u>N) u=twoN-u; if(u==0){return(0.0);}else{return(sign*S[u]);} } Fig. 1. C language procedures to extend the DCT2 and DST2 beyond the usual index ranges. These are used by the windowing procedure of Fig. 2.
4. Computational complexity The algorithms proposed have computational complexity of order N: The number of arithmetic operations required to window both the DCT and DST in the transform domain using the proposed algorithms is:
Hann :
16N þ 8 operations ð8N additions; 8N þ 8 multiplicationsÞ;
Hamming or Blackman :
30N þ 16 operations ð16N additions; 14N þ 8 multiplicationsÞ:
118
B.G. Sherlock, Y.P. Kakad / Journal of the Franklin Institute 339 (2002) 111–120
window_transform(double C[],double S[],double Cw[],double Sw[],int N, double a0, double a1, double a2){ // converts the DCT2 and DST2 of a length N signal into the // DCT2 and DST2 of the same signal with a window applied. // The window is w(x) = a0 + a1 cos 2 pi x/N + a2 cos 4 pi x/N // // C[0..N-1] DCT to be windowed (input) // S[1..N] DST to be windowed (input) // Cw[0..N-1] DCT windowed (output) // Sw[1..N] DST windowed (output) // N signal and transform length // a0,a1,a2 coefficients of window // // For Blackman window, a0=0.42, a1=-0.5, a2=0.08 // For Hamming window, a0=0.54, a1=-0.46, a2=0 // For Hann window, a0=0.5, a1=-0.5, a2=0 // int u; double double double double
Cum2,Cup2,Cum4,Cup4,Sum2,Sup2,Sum4,Sup4; cos_pivN,sin_pivN, root2, invroot2,cos_2pivN,sin_2pivN; Cterm1,Cterm2,Cterm3,Cterm4,Call4; Sterm1,Sterm2,Sterm3,Sterm4,Sall4;
cos_pivN = cos(M_PI/N); cos_2pivN = cos(2.0*M_PI/N); sin_pivN = sin(M_PI/N); sin_2pivN = sin(2.0*M_PI/N); root2=sqrt(2.0); invroot2=1.0/sqrt(2.0); for (u=0; u<=N;u++){ Cum2= extendDCT2(C,u-2,N); Cup2= extendDCT2(C,u+2,N); Sum2= extendDST2(S,u-2,N); Sup2= extendDST2(S,u+2,N); Cum4= extendDCT2(C,u-4,N); Cup4= extendDCT2(C,u+4,N); Sum4= extendDST2(S,u-4,N); Sup4= extendDST2(S,u+4,N); Cterm1 = cos_pivN*Cup2 + sin_pivN*Sup2; Cterm2 = cos_pivN*Cum2 - sin_pivN*Sum2; Cterm3 = cos_2pivN*Cup4 + sin_2pivN*Sup4; Cterm4 = cos_2pivN*Cum4 - sin_2pivN*Sum4; Sterm1 = cos_pivN*Sup2 - sin_pivN*Cup2; Sterm2 = cos_pivN*Sum2 + sin_pivN*Cum2; Sterm3 = cos_2pivN*Sup4 - sin_2pivN*Cup4; Sterm4 = cos_2pivN*Sum4 + sin_2pivN*Cum4; if(u== N-2){Cterm1 *= root2; Sterm1 *= root2;} if(u== 2){ Cterm2 *= root2; Sterm2 *= root2;} if(u== N-4){Cterm3 *= root2; Sterm3 *= root2;} if(u== 4){ Cterm4 *= root2; Sterm4 *= root2;} Call4 = (a1*(Cterm1+Cterm2) + a2*(Cterm3+Cterm4))/2.0; Sall4 = (a1*(Sterm1+Sterm2) + a2*(Sterm3+Sterm4))/2.0; if ((u==0)||(u==N)){Call4 *= invroot2; Sall4 *= invroot2;} Cw[u]= a0*C[u] + Call4; Sw[u]= a0*S[u] + Sall4; } } Fig. 2. C language implementation of the algorithm to perform windowing in the discrete cosine and sine transform domains, for the Blackman, Hamming, or Hann window. This procedure calls the subroutines of Fig. 1 to extend the DCT2 and DST2 beyond the usual index range.
B.G. Sherlock, Y.P. Kakad / Journal of the Franklin Institute 339 (2002) 111–120
119
The conventional approach requires N multiplications to apply the window function, and 2N log2 N operations for each transform if the fastest known DCT/ DST algorithms are used, where N is a power of 2. This yields a total of N þ 4N log2 N operations, for an overall computational complexity of order N log2 N: It is easily seen that our proposed approach is more efficient than the conventional approach when the transform length N is at least 4 for the Hann window, or 256 for the Hamming/Blackman windows. Here we have quoted the smallest powers of 2 for which our approach gives superior performance. For values of N that are not powers of 2, the computational advantage offered by our approach will be considerably greater, because of the additional computation required to calculate DCT/DST for such signal lengths. These algorithms would be of particular interest in applications where frequent updates are required in the presence of a sliding window, e.g. where new data points arrive at frequent intervals, and it is desired to maintain the DCT/DST of the N most recent data points.
5. Conclusion General, computationally efficient algorithms for applying Hann, Hamming, Blackman, and related windows to the DCT and DST directly in the transform domain data were analytically derived. These can be used in applications where the transforms must be updated in the presence of a window function, as new data is received in real time. Software implementations which demonstrate the correctness of developed algorithms in C language are also provided.
Acknowledgements This research was supported in part by the Office of Naval Research Grant Number N00014-99-1-0091. The authors would like to express their gratitude for the useful and detailed evaluation of the manuscript by the anonymous reviewers. The clarity and completeness of this paper has been substantially improved as a result of their suggestions.
References [1] N. Ahmed, T. Natarajan, K.R. Rao, Discrete cosine transform, IEEE Trans. Comput. C-23 (1974) 90–93. [2] H.C. Andrews, Multidimensional rotations in feature selection, IEEE Trans. Comput. C-20 (1971) 1045–1051. [3] W.K. Pratt, Generalized Wiener filtering computation techniques, IEEE Trans. Comput. C-21 (1972) 636–641.
120
B.G. Sherlock, Y.P. Kakad / Journal of the Franklin Institute 339 (2002) 111–120
[4] J. Pearl, Walsh processing of random signals, IEEE Trans. Electromagn. Components EMC-13 (1971) 137–141. [5] J. Pearl, H.C. Andrews, W.K. Pratt, Performance measures for transform data coding, IEEE Trans. Commun. Technol. COM-20 (1972) 411–415. [6] K.R. Rao, P. Yip, Discrete Cosine Transform: Algorithms, Advantages, Applications, Academic Press, Boston, 1990. [7] H. Kwakernaak, R. Sivan, Modern Signals and Systems, Prentice-Hall, Englewood Cliffs, NJ, 1991. [8] F.J. Harris, On the use of windows for harmonic analysis with the discrete fourier transform, Proc. IEEE 66 (1978). [9] P. Yip, K.R. Rao, On the shift property of DCT’s and DST’s, IEEE Trans. Acoust. Speech Signal Process. ASSP-35 (1987) 404–406. [10] B.G. Sherlock, D.M. Monro, Moving fast fourier transform, IEE Proc. F 139 (4) (1992) 279–282. [11] B.G. Sherlock, Windowed discrete fourier transform for shifting data, Signal Process. 74 (1999) 169– 177. [12] B.G. Sherlock, Y.P. Kakad, Windowing the discrete cosine transform in the transform domain, Fourth World Multiconference on Circuits, Systems, Communications and Computers (CSCC 2000), Athens, Greece, July 10–15, 2000, pp. 320–324. [13] H.B. Kekre, J.K. Solanki, Comparative performance of various trigonometric unitary transforms for transform image coding, Int. J. Electron. 44(3) (1978) 305–315. [14] Z. Wang, A fast algorithm for the discrete sine transform implemented by the fast cosine transform, IEEE Trans. Acoust. Speech Signal Process. ASSP-30 (1982). [15] R.C. Gonzalez, R.E. Woods, Digital Image Processing, Addison-Wesley, Reading, MA, 1992. [16] N.C. Geckinli, D. Yavuz, Some novel windows and a concise tutorial comparison of window families, IEEE Trans. Acoust. Speech Signal Process. ASSP-26 (1978) 501–507.