next up previous
Next: Multitapering Up: Autocorrelation Estimation Previous: Single Tapers

Non-parametric Estimation

Instead of presetting $ M$ to what is considered to be a reasonable value, we instead apply some constraints. The first assumption is that $ \rho_{xx}(\tau)>0$ for all $ \tau$. The second assumption is that the autocorrelation is monotonically decreasing. This means that low frequency components are favoured, which are widely accepted in the literature as being the most important to account for.

The autocorrelation is estimated using a standard unbiased estimator (equation 11) and then the best least squares fit that satisfies the constraint of monotonicity can be obtained using techniques from the literature of isotonic regression. The particular algorithm that we use is the Pool Adjacent Violators Algorithm (PAVA) (T. Robertson and Dykstra, 1988), which provides a unique, least squares fit under the constraint. Before using the algorithm, we set $ \hat{\rho}_{xx}(\tau)=r_{xx}(\tau)$ for $ \tau=1{,\ldots,}N/3$ and $ \hat{\rho}_{xx}(\tau)=0$ for $ \tau>N/3$. This is done partly because it reduces the amount of data the algorithm iterates over, and also because the raw autocorrelation estimate is very noisy for $ \tau \geq N/3$ (there is less data available to compute autocorrelations at high lags), and we do not expect significant autocorrelation at such high lags. Furthermore, the value of zero will propagate, eventually stopping at the lag which gives $ M$. For the purpose of the algorithm it is also necessary to define a weighting function $ w(\tau)=1$ for $ \tau=1{,\ldots,}N$. The algorithm then proceeds as follows:

  1. If $ \hat{\rho}_{xx}(\tau)$ is not isotonic there must exist a violator at $ k$ such that $ \hat{\rho}_{xx}(k-1)>\hat{\rho}_{xx}(k)$.
  2. Pool these two values, by replacing them both with their weighted average:

    $\displaystyle [\hat{\rho}_{xx}(k-1)w(k-1)+\hat{\rho}_{xx}(k)w(k)]/[w(k-1)+w(k)]$ (14)

  3. Replace $ w(k-1)$ and $ w(k)$ by $ w(k-1)+w(k)$
  4. Repeat until no more violators.

The algorithm was tested on artificial data which consisted of white noise of length $ N=200$, and had been low-pass filtered with a Gaussian of varying standard deviation. This highlighted a slight bias for white noise data, which was easily remedied by setting $ \hat{\rho}_{xx}(1)=0$ if $ \hat{\rho}_{xx}(1)<0.1$.


next up previous
Next: Multitapering Up: Autocorrelation Estimation Previous: Single Tapers
Mark Woolrich 2001-07-16