From 2de77865c1a0008b38075cc8f410614cc1f9338f Mon Sep 17 00:00:00 2001
From: Ahmet Inan <inan@aicodix.de>
Date: Thu, 9 Apr 2020 10:58:02 +0200
Subject: [PATCH] added SMA3 and SMA4

Pairwise reduction accelerator idea lifted from:
https://github.com/xdsopl/sma/blob/master/pairwise.hh
---
 README.md |  5 +++++
 sma.hh    | 46 ++++++++++++++++++++++++++++++++++++++++++++++
 2 files changed, 51 insertions(+)

diff --git a/README.md b/README.md
index 5980ba4..86bcfde 100644
--- a/README.md
+++ b/README.md
@@ -99,6 +99,11 @@ Implemented are the following [trigger functions](https://en.wikipedia.org/wiki/
 
 The [simple moving average](https://en.wikipedia.org/wiki/Moving_average#Simple_moving_average) gives us the mean of the last N data points.
 
+* SMA1 computes the sum of its internal history buffer at each new input from scratch is and therefore very slow.
+* SMA2 updates its internal sum using only the new input and the oldest value in the history buffer. It is therefore the fastest of all and works perfect with integers but suffers from drift when used with floats on sequences having a high dynamic range.
+* SMA3 is based on SMA2 but uses the [Kahan summation algorithm](https://en.wikipedia.org/wiki/Kahan_summation_algorithm) to reduce drift significantly.
+* SMA4 uses a tree and only update nodes that depend on the new input value and is slower than SMA3 but it has no drift.
+
 ### [calculus.hh](calculus.hh)
 
 Some [calculus](https://en.wikipedia.org/wiki/Calculus) functions:
diff --git a/sma.hh b/sma.hh
index c84cf41..ca284f9 100644
--- a/sma.hh
+++ b/sma.hh
@@ -6,6 +6,8 @@ Copyright 2019 Ahmet Inan <inan@aicodix.de>
 
 #pragma once
 
+#include "kahan.hh"
+
 namespace DSP {
 
 template <typename TYPE, typename VALUE, int NUM>
@@ -61,5 +63,49 @@ public:
 	}
 };
 
+template <typename TYPE, typename VALUE, int NUM>
+class SMA3
+{
+	TYPE hist_inp[NUM];
+	Kahan<TYPE> hist_sum;
+	int hist_pos;
+public:
+	SMA3() : hist_sum(0), hist_pos(0)
+	{
+		for (int i = 0; i < NUM; ++i)
+			hist_inp[i] = 0;
+	}
+	TYPE operator () (TYPE input)
+	{
+		hist_sum(-hist_inp[hist_pos]);
+		hist_inp[hist_pos] = input;
+		if (++hist_pos >= NUM)
+			hist_pos = 0;
+		return hist_sum(input) / VALUE(NUM);
+	}
+};
+
+template <typename TYPE, typename VALUE, int NUM>
+class SMA4
+{
+	TYPE tree[2 * NUM];
+	int leaf;
+public:
+	SMA4() : leaf(NUM)
+	{
+		for (int i = 0; i < 2 * NUM; ++i)
+			tree[i] = 0;
+	}
+	TYPE operator () (TYPE input)
+	{
+		tree[leaf] = input;
+		for (int child = leaf, parent = leaf / 2; parent; child = parent, parent /= 2)
+			tree[parent] = tree[child] + tree[child^1];
+		if (++leaf >= 2 * NUM)
+			leaf = NUM;
+		return tree[1] / VALUE(NUM);
+	}
+};
+
 }