fix #22 improve getAverage(nMedian) (#23)

- add **getMedianAverage(nMedians)**  removing bias - #22
  - thanks to Peter Kowald
- add example **RunningMedian_getMedianAverage.ino**
- extended performance test
- update readme.md

CHANGELOG.md

@@ -6,6 +6,14 @@ The format is based on [Keep a Changelog](http://keepachangelog.com/)
 and this project adheres to [Semantic Versioning](http://semver.org/).
 
 
+## [0.3.8] - 2023-07-11
+- add **getMedianAverage(nMedians)**  removing bias - #22
+  - thanks to Peter Kowald
+- add example **RunningMedian_getMedianAverage.ino**
+- extended performance test
+- update readme.md
+
+
 ## [0.3.7] - 2022-10-28
 - Add RP2040 support to build-CI.
 - Add CHANGELOG.md
@@ -14,7 +22,6 @@ and this project adheres to [Semantic Versioning](http://semver.org/).
 - add setSearchMode() for selecting fastest median search mode.  
   EXPERIMENTAL, select between LINEAR or BINARY search.  
 
-
 ## [0.3.6} - 2022-06-06
 - bump version for platformio
 

README.md

@@ -43,7 +43,7 @@ is large. For most applications a value much lower e.g. 19 is working well, and
 is performance wise O(100x) faster in sorting than 255 elements.
 
 
-### Note: Configurable Options
+#### Note: Configurable Options
 
 There are several options that can be configured via defines at compile time, those being:
 - **RUNNING_MEDIAN_USE_MALLOC**: bool
@@ -57,10 +57,22 @@ There are several options that can be configured via defines at compile time, th
   - Static: The buffer stores at most this many items.
 
 
+#### Related
+
+- https://github.com/RobTillaart/Correlation
+- https://github.com/RobTillaart/GST - Golden standard test metrics
+- https://github.com/RobTillaart/Histogram
+- https://github.com/RobTillaart/RunningAngle
+- https://github.com/RobTillaart/RunningAverage
+- https://github.com/RobTillaart/RunningMedian
+- https://github.com/RobTillaart/statHelpers - combinations & permutations
+- https://github.com/RobTillaart/Statistic
+
+
 ## Interface
 
 
-### Constructor
+#### Constructor
 
 - **RunningMedian(const uint8_t size)** Constructor, dynamically allocates memory.
 - **~RunningMedian()** Destructor.
@@ -69,7 +81,7 @@ There are several options that can be configured via defines at compile time, th
 - **bool isFull()** returns true if the internal buffer is 100% filled.
 
 
-### Base functions
+#### Base functions
 
 - **clear()** resets internal buffer and variables, effectively empty the buffer.
 - **add(const float value)** adds a new value to internal buffer, 
@@ -78,21 +90,47 @@ optionally replacing the oldest element if the buffer is full.
 - **float getAverage()** returns average of **all** the values in the internal buffer.
 - **float getAverage(uint8_t nMedian)** returns average of **the middle n** values. 
 This effectively removes noise from the outliers in the samples.
+The function is improved in 0.3.8 to correct a bias, see #22.
+- **float getMedianAverage(uint8_t nMedian)** almost same as above, 
+except it compensates for alignment bias, see #22. 
+This is done by adjusting the nMedian parameter (-1 or +1) if needed.
 - **float getHighest()** get the largest values in the buffer.
 - **float getLowest()** get the smallest value in the buffer.
 - **float getQuantile(const float quantile)** returns the Quantile value from the buffer. 
 This value is often interpolated.
 
 
-### Less used functions
+#### getMedianAverage(nMedian)
+
+**getAverage(nMedian)** and **getMedianAverage(uint8_t nMedian)** differ.
+When nMedian is odd and count is even or vice versa, the middle N are not 
+perfectly in the middle. 
+By auto-adjusting nMedian (-1 +1) this balance is restored.
+
+Assume an internal size of 7 elements \[0..6] then 
+- **getAverage(4)** will average element 1, 2, 3, 4
+- **getMedianAverage(4)** will adjust nMedian and average element 2, 3, 4.
+
+The example **RunningMedian_getMedianAverage.ino** shows the difference.
+
+The implementation of **getMedianAverage(uint8_t nMedian)** is experimental 
+and might change in the future. 
+Idea is taking top and bottom elements only for 50% if needed, however that 
+implies at least 2 extra float multiplications.
+
+It is possible that the name **getMedianAverage(uint8_t nMedian)**
+will change in the future to be more descriptive.
+
+
+#### Less used functions
 
 - **float getElement(const uint8_t n)** returns the n'th element from the values in time order.
 - **float getSortedElement(const uint8_t n)** returns the n'th element from the values in size order (sorted ascending).
 - **float predict(const uint8_t n)** predict the maximum change of median after n additions, 
 n must be smaller than **getSize()/2**.
 
 
-### SearchMode optimization
+#### SearchMode optimization
 
 Since 0.3.7 the internal sort has been optimized.
 It is now possible to select between LINEAR (=0) and BINARY (=1) insertion sort.
@@ -125,3 +163,4 @@ See examples.
   - get the median without (full) sorting. QuickSelect()
 - move all code to .cpp file
 
+

RunningMedian.cpp

@@ -1,7 +1,7 @@
 //
 //    FILE: RunningMedian.cpp
 //  AUTHOR: Rob Tillaart
-// VERSION: 0.3.7
+// VERSION: 0.3.8
 // PURPOSE: RunningMedian library for Arduino
 //
 //  HISTORY: see changelog.md
@@ -125,6 +125,47 @@ float RunningMedian::getAverage(uint8_t nMedians)
 }
 
 
+//  nMedians is the spread, or the middle N
+//  this version compensated for bias #22
+
+float RunningMedian::getMedianAverage(uint8_t nMedians)
+{
+  //  handle special cases.
+  if ((_count == 0) || (nMedians == 0)) return NAN;
+  if (_count == 1) return _values[0];
+  if (_count == 2) return (_values[0] + _values[1]) * 0.5;
+
+  //  nMedians can not be larger than current nr of elements.
+  if (_count <= nMedians) return getAverage();
+
+  //  _count is at least 3 from here
+
+  //  Eliminate the bias when the nMedians would fall slightly
+  //  to the left or right of the centre.
+  //  If count and nMedians are not both odd or both even reduce
+  //  the spread by 1 to make them the same.
+  //  If nMedians becomes 0 correct this. to 2.
+  if ((_count & 0x01) != (nMedians & 0x01))
+  {
+     --nMedians;
+     //  nmedians can not become 0
+     if (nMedians == 0) nMedians = 2;
+  }
+
+  uint8_t start = (_count - nMedians) / 2;
+  uint8_t stop = start + nMedians;
+
+  if (_sorted == false) sort();
+
+  float sum = 0;
+  for (uint8_t i = start; i < stop; i++)
+  {
+    sum += _values[_sortIdx[i]];
+  }
+  return sum / nMedians;
+}
+
+
 float RunningMedian::getElement(const uint8_t n)
 {
   if ((_count == 0) || (n >= _count)) return NAN;

RunningMedian.h

@@ -3,7 +3,7 @@
 //    FILE: RunningMedian.h
 //  AUTHOR: Rob Tillaart
 // PURPOSE: RunningMedian library for Arduino
-// VERSION: 0.3.7
+// VERSION: 0.3.8
 //     URL: https://github.com/RobTillaart/RunningMedian
 //     URL: http://arduino.cc/playground/Main/RunningMedian
 // HISTORY: See RunningMedian.cpp
@@ -12,7 +12,7 @@
 
 #include "Arduino.h"
 
-#define RUNNING_MEDIAN_VERSION        (F("0.3.7"))
+#define RUNNING_MEDIAN_VERSION        (F("0.3.8"))
 
 
 //  fall back to fixed storage for dynamic version => remove true
@@ -61,6 +61,9 @@ class RunningMedian
   float   getAverage();
   //  returns average of the middle nMedian values, removes noise from outliers
   float   getAverage(uint8_t nMedian);
+  //  returns average of the middle nMedian values, removes noise from outliers
+  //  Bias compensated see #22.
+  float   getMedianAverage(uint8_t nMedian);
 
   float   getHighest() { return getSortedElement(_count - 1); };
   float   getLowest()  { return getSortedElement(0); };

examples/RunningMedian_getMedianAverage/RunningMedian_getMedianAverage.ino

@@ -0,0 +1,66 @@
+//
+//    FILE: RunningMedian_getMedianAverage.ino
+//  AUTHOR: Rob Tillaart
+// PURPOSE: test sketch to show the difference between
+//          getAverage(nMedian) and getMedianAverage(nMedian).
+//     URL: https://github.com/RobTillaart/RunningMedian
+
+
+#include <RunningMedian.h>
+
+RunningMedian samples = RunningMedian(11);
+
+
+void test_compare()
+{
+  Serial.println(__FUNCTION__);
+  for (int x = 0; x < 9; x++)
+  {
+    for (int y = 0; y <= x; y++)
+    {
+      float a = samples.getAverage(y);
+      float b = samples.getMedianAverage(y);
+      Serial.print(x);
+      Serial.print('\t');
+      Serial.print(y);
+      Serial.print('\t');
+      Serial.print(a, 4);
+      Serial.print('\t');
+      Serial.print(b, 4);
+      Serial.print('\t');
+      Serial.println(a - b, 4);
+    }
+    Serial.println();
+
+    samples.add(x);
+  }
+}
+
+
+void setup()
+{
+  Serial.begin(115200);
+  Serial.print("Running Median Version: ");
+  Serial.println(RUNNING_MEDIAN_VERSION);
+
+  test_compare();
+
+  samples.clear();
+  samples.add(1);
+  samples.add(2);
+
+  delay(1000);
+  uint32_t start = micros();
+  float f = samples.getMedianAverage(1);
+  uint32_t stop = micros();
+  Serial.println(stop - start);
+  Serial.println(f);
+}
+
+
+void loop()
+{
+}
+
+
+//  -- END OF FILE --

examples/RunningMedian_performance/RunningMedian_performance.ino

@@ -45,10 +45,16 @@ void setup()
 #endif
 
   samples.setSearchMode(0);
-  test1();
+  test_median();
 
   samples.setSearchMode(1);
-  test1();
+  test_median();
+
+  samples.setSearchMode(0);
+  test_average();
+
+  samples.setSearchMode(1);
+  test_average();
 
   Serial.println("\ndone..\n");
 }
@@ -57,14 +63,15 @@ void loop()
 {
 }
 
-void test1()
+void test_median()
 {
   uint32_t start = 0;
   uint32_t stop  = 0;
   uint32_t total = 0;
 
   samples.clear();
   Serial.println();
+  Serial.println(__FUNCTION__);
   Serial.print(F("Allocated size = "));
   Serial.println(samples.getSize());
   Serial.print(F("nr of elements = "));
@@ -92,7 +99,7 @@ void test1()
   stop = micros();
   Serial.print(F(" 1 x get: "));
   Serial.print(stop - start);
-  Serial.println(F("       == sorting"));
+  Serial.println(F("\t\t== sorting"));
   Serial.print(  "  median: ");
   Serial.println(result);
   delay(100);
@@ -103,11 +110,65 @@ void test1()
   stop = micros();
   Serial.print(F(" 1 x get: "));
   Serial.print(stop - start);
-  Serial.println(F("       == no sorting"));
+  Serial.println(F("\t\t== no sorting"));
   Serial.print(  "  median: ");
   Serial.println(result);
   delay(100);
+}
+
+
+void test_average()
+{
+  uint32_t start = 0;
+  uint32_t stop  = 0;
+  uint32_t total = 0;
+
+  samples.clear();
+  Serial.println();
+  Serial.println(__FUNCTION__);
+  Serial.print(F("Allocated size = "));
+  Serial.println(samples.getSize());
+  Serial.print(F("nr of elements = "));
+  Serial.println(samples.getCount());
+  Serial.print(F("    searchMode = "));
+  Serial.println(samples.getSearchMode());
+  Serial.println();
+  delay(50);
+
+  for (uint8_t i = 0; i < sourceSize; i++)
+  {
+    start = micros();
+    samples.add(sourceData[i]);
+    total += (micros() - start);
+  }
+  Serial.print(F("50 x add: "));
+  Serial.println(total);
+  Serial.print(F("     avg: "));
+  Serial.println(total / 50.0);
+  delay(100);
+
+  // time to access the data
+  start = micros();
+  float result = samples.getAverage(20);
+  stop = micros();
+  Serial.print(F(" 1 x get: "));
+  Serial.print(stop - start);
+  Serial.println(F("\t\t== sorting"));
+  Serial.print(  " average: ");
+  Serial.println(result);
+  delay(100);
 
+  // time to access the data
+  start = micros();
+  result = samples.getAverage(20);
+  stop = micros();
+  Serial.print(F(" 1 x get: "));
+  Serial.print(stop - start);
+  Serial.println(F("\t\t== no sorting"));
+  Serial.print(  " average: ");
+  Serial.println(result);
+  delay(100);
 }
 
-// -- END OF FILE --
+
+//  -- END OF FILE --