analysis/app/kza.go

// Copyright 2015 The Go Authors. All rights reserved.
// Use of this source code is governed by a BSD-style
// license that can be found in the LICENSE file.

package app

import "math"

// TODO: This all assumes that data is sampled at a regular interval
// and there are no missing values. It could be generalized to accept
// missing values (perhaps represented by NaN), or generalized much
// further by accepting (t, x) pairs and a vector of times at which to
// evaluate the filter (and an arbitrary window size). I would have to
// figure out how that affects the difference array in KZA.

// TODO: These can generate a lot of garbage. Perhaps the caller
// should pass in the target slice? Or these should just overwrite the
// input array and leave it to the caller to copy if necessary?

// MovingAverage performs a moving average (MA) filter of xs with
// window size m. m must be a positive odd integer.
//
// Note that this is filter is often described in terms of the half
// length of the window (m-1)/2.
func MovingAverage(xs []float64, m int) []float64 {
	if m <= 0 || m%2 != 1 {
		panic("m must be a positive, odd integer")
	}
	ys := make([]float64, len(xs))
	sum, n := 0.0, 0
	for l, i, r := -m, -(m-1)/2, 0; i < len(ys); l, i, r = l+1, i+1, r+1 {
		if l >= 0 {
			sum -= xs[l]
			n--
		}
		if r < len(xs) {
			sum += xs[r]
			n++
		}
		if i >= 0 {
			ys[i] = sum / float64(n)
		}
	}
	return ys
}

// KolmogorovZurbenko performs a Kolmogorov-Zurbenko (KZ) filter of xs
// with window size m and k iterations. m must be a positive odd
// integer. k must be positive.
func KolmogorovZurbenko(xs []float64, m, k int) []float64 {
	// k is typically small, and MA is quite efficient, so just do
	// the iterated moving average rather than bothering to
	// compute the binomial coefficient kernel.
	for i := 0; i < k; i++ {
		// TODO: Generate less garbage.
		xs = MovingAverage(xs, m)
	}
	return xs
}

// AdaptiveKolmogorovZurbenko performs an adaptive Kolmogorov-Zurbenko
// (KZA) filter of xs using an initial window size m and k iterations.
// m must be a positive odd integer. k must be positive.
//
// See Zurbenko, et al. 1996: Detecting discontinuities in time series
// of upper air data: Demonstration of an adaptive filter technique.
// Journal of Climate, 9, 3548–3560.
func AdaptiveKolmogorovZurbenko(xs []float64, m, k int) []float64 {
	// Perform initial KZ filter.
	z := KolmogorovZurbenko(xs, m, k)

	// Compute differenced values.
	q := (m - 1) / 2
	d := make([]float64, len(z)+1)
	maxD := 0.0
	for i := q; i < len(z)-q; i++ {
		d[i] = math.Abs(z[i+q] - z[i-q])
		if d[i] > maxD {
			maxD = d[i]
		}
	}

	if maxD == 0 {
		// xs is constant, so no amount of filtering will do
		// anything. Avoid dividing 0/0 below.
		return xs
	}

	// Compute adaptive filter.
	ys := make([]float64, len(xs))
	for t := range ys {
		dPrime := d[t+1] - d[t]
		f := 1 - d[t]/maxD

		qt := q
		if dPrime <= 0 {
			// Zurbenko doesn't specify what to do with
			// the fractional part of qt and qh, so we
			// interpret this as summing all points of xs
			// between qt and qh.
			qt = int(math.Ceil(float64(q) * f))
		}
		if t-qt < 0 {
			qt = t
		}

		qh := q
		if dPrime >= 0 {
			qh = int(math.Floor(float64(q) * f))
		}
		if t+qh >= len(xs) {
			qh = len(xs) - t - 1
		}

		sum := 0.0
		for i := t - qt; i <= t+qh; i++ {
			sum += xs[i]
		}
		// Zurbenko divides by qh+qt, but this undercounts the
		// number of terms in the sum by 1.
		ys[t] = sum / float64(qh+qt+1)
	}

	return ys
}