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 }