28-01-2013, 03:54 PM
Wavelet Thresholding Techniques for Power Spectrum Estimation
Wavelet Thresholding.pdf (Size: 61.63 KB / Downloads: 37)
ABSTRACT
Estimation of the power spectrum S( f ) of a stationary random process can be viewed as a
nonparametric statistical estimation problem. We introduce a nonparametric approach based on a
wavelet representation for the logarithm of the unknown S( f ). This approach offers the ability to
capture statistically significant components of lnS( f ) at different resolution levels and guarantees
nonnegativity of the spectrum estimator. The spectrum estimation problem is set up as a problem
of inference on the wavelet coefficients of a signal corrupted by additive non-Gaussian noise. We
propose a wavelet thresholding technique to solve this problem under specified noise/resolution
tradeoffs and show that the wavelet coefficients of the additive noise may be treated as independent
random variables. The thresholds are computed using a saddle-point approximation to the
distribution of the noise coefficients.
INTRODUCTION
We consider the problem of estimating the power spectrum {S( f ) , - 1/2 < f £ 1/2} of a
discrete-time, wide-sense stationary, real, Gaussian random process
{x(n) , n = 0 , ±1 , ±2 , . . . } from a finite set of observations. Consistent estimates may be
obtained by suitable processing of the empirical spectrum estimates (periodogram). This may be
done using sieves [1][2], spline smoothing [3-5], or kernel smoothing [6, §5.4][7][8]. These
methods all require the choice of a certain resolution parameter, often called bandwidth in the
spectrum estimation literature [7][8]. Various techniques for optimal bandwidth selection have
been studied in [2][4-8]. They produce estimates that have a good overall bias vs variance tradeoff.
However, the bias vs variance tradeoff is generally not optimal locally. Sharp peaks and
troughs in the spectrum require a narrow spectral bandwidth, otherwise they are oversmoothed
and the bias is large; in contrast, smooth components of the spectrum require a wide bandwidth in
order to achieve a significant noise reduction. These concerns have long been recognized [8,
§7.2], so various ideas have been introduced to remedy this situation. Jenkins’ technique of window
closing produces a family of estimates using different bandwidths [8]. Unfortunately, this
technique does not provide a mechanism for combining the information obtained at different resolutions.
Capon’s method (originally developed for wavenumber analysis with arrays [9]) produces
estimates of the spectrum power in different frequency bands using adaptive filtering techniques,
but is not designed to produce an estimate of the spectral density itself.
Fixed-Bandwidth Smoothing
The drawbacks of applying single-resolution (fixed-bandwidth) smoothing techniques to the
scaled log-periodogram are illustrated in Fig. 2. In Fig. 2a and b, we show the scaled logperiodogram
for Example 1, convolved with Hamming windows with bandwidths 128 ´ ( 0. 5/ N)
and 32 ´ ( 0. 5/ N), respectively. For the wide bandwidth, the peak at the pole frequency is
smeared out. For the narrow bandwidth, the peak is captured correctly, but the spectral estimates
far from the peak are too wiggly. In the following section, we show how wavelet techniques can
be used to achieve a local adaptation of the bandwidth to the data.
Discrete Wavelet Transform
We assume that N is a multiple of 2J and consider an orthonormal, discrete wavelet
transform (DWT) for signals on the interval [ 0 , N - 1 ]. The transform is implemented using an
iterated bank of subsampled lowpass and highpass quadrature-mirror, finite-impulse-response
filters {g(l)} and {h(l)}, respectively. We denote by y(x) the (continuous) wavelet associated
with the filter bank [29, §5.6]. In order to keep the exposition of basic concepts clear, the theoretical
results which motivate our approach are presented for the simple case of periodic wavelets.
This is equivalent to assuming that the signals are periodic beyond the boundaries at l = 0 and
l = N - 1. A more accurate treatment of the boundaries is proposed in Section 3.5.
Inference on Wavelet Coefficients
If ln S is smooth enough, the wavelet coefficients {aj l} decay rapidly at fine scales [32].
This behavior is to be contrasted with that of the wavelet coefficients for the noise {ej l}, which
have scale-independent variance se
2. This property is illustrated in Fig. 5, where the DWT
coefficients for the log-spectrum and for the scaled log-periodogram in Examples 1 and 2 are
displayed at various scales. The DWT is computed using the RD6 filters. It can be seen that the
noise corrupts wavelet coefficients at all scales equally, while the signal contributions appear
either at coarse scales or in the vicinity of abrupt changes in the spectrum.
Boundary Conditions
The easiest solution to the boundary problem is to construct a periodic extension of the signal.
Unfortunately, this introduces a discontinuity at the boundaries. After thresholding of the
DWT coefficients, the spectral estimates exhibit typical artifacts due to the Gibbs phenomenon.
A possible remedy is to use special wavelets near the boundaries [29, §10.7][34], a solution
analogous to the use of special boundary kernels in linear smoothing [33, §4.4]. Another solution,
which has the advantage of using standard wavelet software, consists in applying the
periodic DWT to a full period of the scaled log-periodogram (2N points). The spectral estimates
are obtained by thresholding the 2N DWT coefficients, applying the inverse DWT, and displaying
the N positive-frequency samples. This approach is equivalent to introducing extra coefficients to
handle the boundaries. It should be noticed that the assumption of asymptotic mutual independence
of the wavelet coefficients breaks down near the boundaries, so the thresholding technique
is not expected to perform as well as in the interior of the interval.
Saddle-Point Approximation
A good approximation to Pj (. ) is needed at fine and moderate scales. The Fourier technique
for computing Pj (x) performs poorly for large x due to numerical instabilities (see Appendix).
Cramer’s large deviations upper bound [35, §3] is not tight, in the sense that the relative
error on Pj (x) is unbounded. A more suitable technique is the saddle-point approximation [22-
24]. The saddle-point approximation is the first term in an asymptotic series expansion for Pj (. ),
as described in the Appendix. It is displayed on a log scale in Fig. 6 at different scales, for a
variety of wavelets. The size of the next-order correction term is typically no greater than 30% at
the finest scale and decreases at coarser scales.
CONCLUSION
We have shown how wavelet techniques may be applied to perform a nonlinear smoothing
of the log-periodogram and produce an automatic adjustment of the bandwidth to the data. We
have proposed a wavelet thresholding scheme for which the wavelet coefficients of the noise may
be treated as independent random variables. The presence of log-spectrum components is decided
from binary hypothesis tests. The significance level of the test determines the desired
noise/resolution tradeoff, and the thresholds are determined from the distribution of the noise
coefficients. Central-limit-theorem results apply at coarse scales; however, the non-Gaussian
character of the distribution should be taken into account at fine and moderate scales. The
saddle-point approximation has been shown to be an effective tool for computing thresholds that
attain a specified probability of false alarm.