IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 58, NO. 7, JULY 2010 3907 A Recursive Scheme for Computing Autocorrelation Functions of Decimated Complex Wavelet Subbands Bart Goossens, Jan Aelterman, Aleksandra Piˇ zurica, and Wilfried Philips Abstract—This correspondence deals with the problem of the exact computation of the autocorrelation function of a real or complex dis- crete wavelet subband of a signal, when the autocorrelation function or alternatively the power spectral density (PSD) of the signal in the time domain (or spatial domain) is either known or estimated using a separate technique. The solution to this problem allows us to couple time domain noise estimation techniques to wavelet domain denoising algorithms, which is crucial for the development of “blind” wavelet-based denoising techniques. Specifically, we investigate the Dual-Tree complex wavelet transform (DT-CWT), which has a good directional selectivity in 2-D and 3-D, is approximately shift-invariant and yields better denoising results than a discrete wavelet transform (DWT). The proposed scheme gives an analytical relationship between the PSD of the input signal/image and the PSD of each individual real/complex wavelet subband which is very useful for future developments. We also show that a more general technique, that relies on Monte Carlo simulations, requires a large number of input samples for a reliable estimate, while the proposed technique does not suffer from this problem. Index Terms—Autocorrelation functions, complex wavelets. I. INTRODUCTION Many noise estimation techniques, either designed for white or col- ored noise (e.g., [1]–[8]), estimate the noise power in the time (or spa- tial) domain, while this estimate is often used in a transformed domain, such as the wavelet domain [9]–[12]. An orthonormal linear transform preserves the noise variance of stationary white Gaussian noise. For nonorthogonal, overcomplete transforms or for colored noise, the sit- uation is more complicated. Suppose that the covariance matrix of a finite length signal is given by and the transform matrix is denoted as . The covariance matrix of the transformed signal becomes [13] (1) Even though this relationship is simple, its computation is not, because the dimensions of the matrices can become very large, e.g., in 2-D or 3-D, which makes it impractical in realistic situations due to high memory and computational requirements. Also, if is not known in advance and has to be estimated, a very large number of observations is required, which is also not feasible. Instead, it is common to assume weak-sense stationarity, which states that the correlation between two samples depends on their difference in position, but not on their ab- solute positions. This assumption significantly reduces the number of free parameters in to be estimated, however the memory and com- putational requirements are still remaining. For shift invariant transforms, such as the undecimated wavelet transform [13] or steerable pyramids [14], [15], the task of computing Manuscript received July 13, 2009; accepted March 04, 2010. Date of pub- lication April 05, 2010; date of current version June 16, 2010. The associate editor coordinating the review of this manuscript and approving it for publica- tion was Dr. Ivan W. Selesnick. The work of B. Goossens was supported by Special Research Fund (BOF) from Ghent University. The authors are with the Department of Telecommunications and Informa- tion Processing (TELIN-IPI-IBBT), Ghent University, B-9000 Gent, Belgium (e-mail: Bart.Goossens@telin.UGent.be; Jan Aelterman@telin.ugent.be; Alek- sandra.Pizurica@telin.UGent.be; Wilfried Philips@telin.UGent.be). Digital Object Identifier 10.1109/TSP.2010.2047392 autocorrelation functions of subbands can be easily done by filtering and transforming a Dirac impulse to the transform domain and by computing the covariance matrix (or autocorrelation function) in each subband [11]. However, for shift variant transforms (such as the DWT [16], the DT-CWT [17] and the recently proposed Marr wavelet trans- form [18]), the obtained covariance matrix would depend on the exact position of the Dirac impulse, hence, the result would be incorrect. A general approach would then be to use Monte Carlo simulations: first generate a noise signal with the known covariance matrix (or auto- correlation function), transform this noise signal to the wavelet domain and estimate the autocorrelation functions in each wavelet subband. This method is easy to implement, but has as major drawback that it suffers from estimation errors: often many input samples are needed to have a reliable estimate of the autocorrelation in each wavelet subband, as we will show in Section VI. Despite the fact that recent efficient wavelet-based denoising tech- niques (e.g., [11], [19]) entirely rely on knowledge of the correlation properties of the transform coefficients, the computation and study of the autocorrelations has only received limited attention compared to the many papers that appeared on the topic of denoising. Averkamp and Houdré [20] analyze stationary second-order random processes in the DWT domain. Fowler [21] studies the variance of additive white noise, in nontight undecimated DWT frames. Chaux et al. [22], investi- gate noise covariance properties for dual-tree wavelet decompositions. Their analysis starts from the (cross)correlation functions between dif- ferent continuous wavelet basis functions that have a closed form ex- pression in frequency domain. In this correspondence, we present an exact computation method for the autocorrelation functions in the real or complex wavelet domain, based on DSP theory, in a recursive manner similar to the Fast DWT scheme of Mallat based on iterated filter banks [13], with including dec- imations at every scale. This has the additional advantage that the com- putation is very fast and memory-friendly. Because our DSP approach directly uses the wavelet filter coefficients, we can also apply our tech- nique in case of wavelets that do not have a closed form expression in frequency domain. We develop the expressions for the Dual-Tree Com- plex Wavelet transform (DT-CWT) [17], [23]. We have concentrated on the DT-CWT because of its popularity in signal and image processing. This transform is two times redundant in 1-D (four times in 2-D) and the coefficient magnitudes are approximately shift invariant, which yields significantly better results than a DWT for many real applications, e.g., denoising and usually better results than a nondecimated wavelet trans- form due to a better orientation selectivity. As such the derived expres- sions are important in a number of denoising methods. The method can also be easily extended to other types of multiresolution represen- tations, based on similar derivations. This correspondence is organized as follows: Section II introduces general concepts that are used in the remainder of this correspondence. In Section III the computation method for the 1-D DT-CWT is de- scribed. Next, extensions to the 2-D and higher dimensional oriented DT-CWT are presented, respectively, in Sections IV and V. Results and a discussion are given in Section VI. Finally, Section VII concludes this correspondence. II. PRELIMINARIES In this section, we introduce some concepts that are used in the re- mainder of the correspondence. Let denote the -transform of a zero-mean signal , i.e., . The discrete autocorrelation function of is given by (2) 1053-587X/$26.00 © 2010 IEEE Authorized licensed use limited to: University of Gent. Downloaded on July 05,2010 at 16:43:31 UTC from IEEE Xplore. Restrictions apply.