blob: 9ad6fa62dba7e20959e08bd2536defe16f72a508 [file] [log] [blame]
// 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
}