annotate base/MovingMedian.h @ 1886:f803d3c33f76 tip

Switch off copious debug in soft synth driving
author Chris Cannam
date Fri, 14 Aug 2020 10:44:44 +0100 (2020-08-14)
parents f04038819c26
children
rev   line source
Chris@1573 1 /* -*- c-basic-offset: 4 indent-tabs-mode: nil -*- vi:set ts=8 sts=4 sw=4: */
Chris@1573 2
Chris@1573 3 /*
Chris@1573 4 Sonic Visualiser
Chris@1573 5 An audio file viewer and annotation editor.
Chris@1573 6 Centre for Digital Music, Queen Mary, University of London.
Chris@1573 7 This file copyright 2007-2015 Particular Programs Ltd,
Chris@1573 8 copyright 2018 Queen Mary University of London.
Chris@1573 9
Chris@1573 10 This program is free software; you can redistribute it and/or
Chris@1573 11 modify it under the terms of the GNU General Public License as
Chris@1573 12 published by the Free Software Foundation; either version 2 of the
Chris@1573 13 License, or (at your option) any later version. See the file
Chris@1573 14 COPYING included with this distribution for more information.
Chris@1573 15 */
Chris@1573 16
Chris@1573 17 #ifndef SV_MOVING_MEDIAN_H
Chris@1573 18 #define SV_MOVING_MEDIAN_H
Chris@1573 19
Chris@1573 20 #include <bqvec/Allocators.h>
Chris@1573 21 #include <bqvec/VectorOps.h>
Chris@1573 22
Chris@1573 23 #include <algorithm>
Chris@1573 24 #include <iostream>
Chris@1573 25 #include <stdexcept>
Chris@1573 26
Chris@1573 27 /**
Chris@1573 28 * Obtain the median (or other percentile) of a moving window across a
Chris@1573 29 * time series. Construct the MovingMedian object, then push() each
Chris@1573 30 * new value in the time series and get() the median of the most
Chris@1573 31 * recent window. The size of the window, and the percentile
Chris@1573 32 * calculated, can both be changed after construction.
Chris@1573 33 *
Chris@1573 34 * Note that for even-sized windows, the "median" is taken to be the
Chris@1573 35 * value at the start of the second half when sorted, e.g. for size 4,
Chris@1573 36 * the element at index 2 (zero-based) in the sorted window.
Chris@1573 37 *
Chris@1573 38 * Not thread-safe.
Chris@1573 39 */
Chris@1573 40 template <typename T>
Chris@1573 41 class MovingMedian
Chris@1573 42 {
Chris@1573 43 public:
Chris@1573 44 MovingMedian(int size, double percentile = 50.f) :
Chris@1573 45 m_size(size),
Chris@1573 46 m_percentile(percentile) {
Chris@1573 47 if (size < 1) throw std::logic_error("size must be >= 1");
Chris@1573 48 m_frame = breakfastquay::allocate_and_zero<T>(size);
Chris@1573 49 m_sorted = breakfastquay::allocate_and_zero<T>(size);
Chris@1573 50 calculateIndex();
Chris@1573 51 }
Chris@1573 52
Chris@1573 53 ~MovingMedian() {
Chris@1573 54 breakfastquay::deallocate(m_frame);
Chris@1573 55 breakfastquay::deallocate(m_sorted);
Chris@1573 56 }
Chris@1573 57
Chris@1573 58 MovingMedian(const MovingMedian &) =delete;
Chris@1573 59 MovingMedian &operator=(const MovingMedian &) =delete;
Chris@1573 60
Chris@1573 61 void setPercentile(double p) {
Chris@1573 62 m_percentile = p;
Chris@1573 63 calculateIndex();
Chris@1573 64 }
Chris@1573 65
Chris@1573 66 void push(T value) {
Chris@1573 67 if (value != value) {
Chris@1573 68 std::cerr << "WARNING: MovingMedian: NaN encountered" << std::endl;
Chris@1573 69 value = T();
Chris@1573 70 }
Chris@1573 71 drop(m_frame[0]);
Chris@1573 72 breakfastquay::v_move(m_frame, m_frame+1, m_size-1);
Chris@1573 73 m_frame[m_size-1] = value;
Chris@1573 74 put(value);
Chris@1573 75 }
Chris@1573 76
Chris@1573 77 T get() const {
Chris@1573 78 return m_sorted[m_index];
Chris@1573 79 }
Chris@1573 80
Chris@1573 81 int size() const {
Chris@1573 82 return m_size;
Chris@1573 83 }
Chris@1573 84
Chris@1573 85 void reset() {
Chris@1573 86 breakfastquay::v_zero(m_frame, m_size);
Chris@1573 87 breakfastquay::v_zero(m_sorted, m_size);
Chris@1573 88 }
Chris@1573 89
Chris@1573 90 void resize(int target) {
Chris@1573 91 if (target == m_size) return;
Chris@1573 92 int diff = std::abs(target - m_size);
Chris@1573 93 if (target > m_size) { // grow
Chris@1573 94 // we don't want to change the median, so fill spaces with it
Chris@1573 95 T fillValue = get();
Chris@1573 96 m_frame = breakfastquay::reallocate(m_frame, m_size, target);
Chris@1573 97 breakfastquay::v_move(m_frame + diff, m_frame, m_size);
Chris@1573 98 breakfastquay::v_set(m_frame, fillValue, diff);
Chris@1573 99 m_sorted = breakfastquay::reallocate(m_sorted, m_size, target);
Chris@1573 100 for (int sz = m_size + 1; sz <= target; ++sz) {
Chris@1573 101 put(m_sorted, sz, fillValue);
Chris@1573 102 }
Chris@1573 103 } else { // shrink
Chris@1573 104 for (int i = 0; i < diff; ++i) {
Chris@1573 105 drop(m_sorted, m_size - i, m_frame[i]);
Chris@1573 106 }
Chris@1573 107 m_sorted = breakfastquay::reallocate(m_sorted, m_size, target);
Chris@1573 108 breakfastquay::v_move(m_frame, m_frame + diff, target);
Chris@1573 109 m_frame = breakfastquay::reallocate(m_frame, m_size, target);
Chris@1573 110 }
Chris@1573 111
Chris@1573 112 m_size = target;
Chris@1573 113 calculateIndex();
Chris@1573 114 }
Chris@1573 115
Chris@1573 116 void checkIntegrity() const {
Chris@1573 117 check();
Chris@1573 118 }
Chris@1573 119
Chris@1573 120 private:
Chris@1573 121 int m_size;
Chris@1573 122 double m_percentile;
Chris@1573 123 int m_index;
Chris@1573 124 T *m_frame;
Chris@1573 125 T *m_sorted;
Chris@1573 126
Chris@1573 127 void calculateIndex() {
Chris@1573 128 m_index = int((m_size * m_percentile) / 100.f);
Chris@1573 129 if (m_index >= m_size) m_index = m_size-1;
Chris@1573 130 if (m_index < 0) m_index = 0;
Chris@1573 131 }
Chris@1573 132
Chris@1573 133 void put(T value) {
Chris@1573 134 put(m_sorted, m_size, value);
Chris@1573 135 }
Chris@1573 136
Chris@1573 137 static void put(T *const sorted, int size, T value) {
Chris@1573 138
Chris@1573 139 // precondition: sorted points to size-1 sorted values,
Chris@1573 140 // followed by an unused slot (i.e. only the first size-1
Chris@1573 141 // values of sorted are actually sorted)
Chris@1573 142 //
Chris@1573 143 // postcondition: sorted points to size sorted values
Chris@1573 144
Chris@1573 145 T *ptr = std::lower_bound(sorted, sorted + size - 1, value);
Chris@1573 146 breakfastquay::v_move(ptr + 1, ptr, int(sorted + size - ptr) - 1);
Chris@1573 147 *ptr = value;
Chris@1573 148 }
Chris@1573 149
Chris@1573 150 void drop(T value) {
Chris@1573 151 drop(m_sorted, m_size, value);
Chris@1573 152 }
Chris@1573 153
Chris@1573 154 static void drop(T *const sorted, int size, T value) {
Chris@1573 155
Chris@1573 156 // precondition: sorted points to size sorted values, one of
Chris@1573 157 // which is value
Chris@1573 158 //
Chris@1573 159 // postcondition: sorted points to size-1 sorted values,
Chris@1573 160 // followed by a slot that has been reset to default value
Chris@1573 161 // (i.e. only the first size-1 values of sorted are actually
Chris@1573 162 // sorted)
Chris@1573 163
Chris@1573 164 T *ptr = std::lower_bound(sorted, sorted + size, value);
Chris@1573 165 if (*ptr != value) {
Chris@1573 166 throw std::logic_error
Chris@1573 167 ("MovingMedian::drop: value being dropped is not in array");
Chris@1573 168 }
Chris@1573 169 breakfastquay::v_move(ptr, ptr + 1, int(sorted + size - ptr) - 1);
Chris@1573 170 sorted[size-1] = T();
Chris@1573 171 }
Chris@1573 172
Chris@1573 173 void check() const {
Chris@1573 174 bool good = true;
Chris@1573 175 for (int i = 1; i < m_size; ++i) {
Chris@1573 176 if (m_sorted[i] < m_sorted[i-1]) {
Chris@1573 177 std::cerr << "ERROR: MovingMedian::checkIntegrity: "
Chris@1573 178 << "mis-ordered elements in sorted array starting "
Chris@1573 179 << "at index " << i << std::endl;
Chris@1573 180 good = false;
Chris@1573 181 break;
Chris@1573 182 }
Chris@1573 183 }
Chris@1573 184 for (int i = 0; i < m_size; ++i) {
Chris@1573 185 bool found = false;
Chris@1573 186 for (int j = 0; j < m_size; ++j) {
Chris@1573 187 if (m_sorted[j] == m_frame[i]) {
Chris@1573 188 found = true;
Chris@1573 189 break;
Chris@1573 190 }
Chris@1573 191 }
Chris@1573 192 if (!found) {
Chris@1573 193 std::cerr << "ERROR: MovingMedian::checkIntegrity: "
Chris@1573 194 << "element in frame at index " << i
Chris@1573 195 << " not found in sorted array" << std::endl;
Chris@1573 196 good = false;
Chris@1573 197 break;
Chris@1573 198 }
Chris@1573 199 }
Chris@1573 200 for (int i = 0; i < m_size; ++i) {
Chris@1573 201 bool found = false;
Chris@1573 202 for (int j = 0; j < m_size; ++j) {
Chris@1573 203 if (m_sorted[i] == m_frame[j]) {
Chris@1573 204 found = true;
Chris@1573 205 break;
Chris@1573 206 }
Chris@1573 207 }
Chris@1573 208 if (!found) {
Chris@1573 209 std::cerr << "ERROR: MovingMedian::checkIntegrity: "
Chris@1573 210 << "element in sorted array at index " << i
Chris@1573 211 << " not found in source frame" << std::endl;
Chris@1573 212 good = false;
Chris@1573 213 break;
Chris@1573 214 }
Chris@1573 215 }
Chris@1573 216 if (!good) {
Chris@1573 217 std::cerr << "Frame contains:" << std::endl;
Chris@1573 218 std::cerr << "[ ";
Chris@1573 219 for (int j = 0; j < m_size; ++j) {
Chris@1573 220 std::cerr << m_frame[j] << " ";
Chris@1573 221 }
Chris@1573 222 std::cerr << "]" << std::endl;
Chris@1573 223 std::cerr << "Sorted array contains:" << std::endl;
Chris@1573 224 std::cerr << "[ ";
Chris@1573 225 for (int j = 0; j < m_size; ++j) {
Chris@1573 226 std::cerr << m_sorted[j] << " ";
Chris@1573 227 }
Chris@1573 228 std::cerr << "]" << std::endl;
Chris@1573 229 throw std::logic_error("MovingMedian failed integrity check");
Chris@1573 230 }
Chris@1573 231 }
Chris@1573 232 };
Chris@1573 233
Chris@1573 234 #endif
Chris@1573 235