qm-dsp
1.8
|
00001 /* -*- c-basic-offset: 4 indent-tabs-mode: nil -*- vi:set ts=8 sts=4 sw=4: */ 00002 00003 /* 00004 QM DSP Library 00005 00006 Centre for Digital Music, Queen Mary, University of London. 00007 This file Copyright 2010 Chris Cannam. 00008 00009 This program is free software; you can redistribute it and/or 00010 modify it under the terms of the GNU General Public License as 00011 published by the Free Software Foundation; either version 2 of the 00012 License, or (at your option) any later version. See the file 00013 COPYING included with this distribution for more information. 00014 */ 00015 00016 #ifndef MEDIAN_FILTER_H 00017 #define MEDIAN_FILTER_H 00018 00019 #include <algorithm> 00020 #include <cassert> 00021 #include <cmath> 00022 #include <iostream> 00023 #include <vector> 00024 00025 template <typename T> 00026 class MedianFilter 00027 { 00028 public: 00029 MedianFilter(int size, float percentile = 50.f) : 00030 m_size(size), 00031 m_frame(new T[size]), 00032 m_sorted(new T[size]), 00033 m_sortend(m_sorted + size - 1) { 00034 setPercentile(percentile); 00035 reset(); 00036 } 00037 00038 ~MedianFilter() { 00039 delete[] m_frame; 00040 delete[] m_sorted; 00041 } 00042 00043 void setPercentile(float p) { 00044 m_index = int((m_size * p) / 100.f); 00045 if (m_index >= m_size) m_index = m_size-1; 00046 if (m_index < 0) m_index = 0; 00047 } 00048 00049 void push(T value) { 00050 if (value != value) { 00051 std::cerr << "WARNING: MedianFilter::push: attempt to push NaN" << std::endl; 00052 return; // nan 00053 } 00054 drop(m_frame[0]); 00055 const int sz1 = m_size-1; 00056 for (int i = 0; i < sz1; ++i) m_frame[i] = m_frame[i+1]; 00057 m_frame[m_size-1] = value; 00058 put(value); 00059 } 00060 00061 T get() const { 00062 return m_sorted[m_index]; 00063 } 00064 00065 int getSize() const { 00066 return m_size; 00067 } 00068 00069 T getAt(float percentile) { 00070 int ix = int((m_size * percentile) / 100.f); 00071 if (ix >= m_size) ix = m_size-1; 00072 if (ix < 0) ix = 0; 00073 return m_sorted[ix]; 00074 } 00075 00076 void reset() { 00077 for (int i = 0; i < m_size; ++i) m_frame[i] = 0; 00078 for (int i = 0; i < m_size; ++i) m_sorted[i] = 0; 00079 } 00080 00081 static std::vector<T> filter(int size, const std::vector<T> &in) { 00082 std::vector<T> out; 00083 MedianFilter<T> f(size); 00084 for (int i = 0; i < int(in.size()); ++i) { 00085 f.push(in[i]); 00086 T median = f.get(); 00087 if (i >= size/2) out.push_back(median); 00088 } 00089 while (out.size() < in.size()) { 00090 f.push(T()); 00091 out.push_back(f.get()); 00092 } 00093 return out; 00094 } 00095 00096 private: 00097 const int m_size; 00098 T *const m_frame; 00099 T *const m_sorted; 00100 T *const m_sortend; 00101 int m_index; 00102 00103 void put(T value) { 00104 // precondition: m_sorted contains m_size-1 values, packed at start 00105 // postcondition: m_sorted contains m_size values, one of which is value 00106 T *point = std::lower_bound(m_sorted, m_sortend, value); 00107 const int n = m_sortend - point; 00108 for (int i = n; i > 0; --i) point[i] = point[i-1]; 00109 *point = value; 00110 } 00111 00112 void drop(T value) { 00113 // precondition: m_sorted contains m_size values, one of which is value 00114 // postcondition: m_sorted contains m_size-1 values, packed at start 00115 T *point = std::lower_bound(m_sorted, m_sortend + 1, value); 00116 if (*point != value) { 00117 std::cerr << "WARNING: MedianFilter::drop: *point is " << *point 00118 << ", expected " << value << std::endl; 00119 } 00120 const int n = m_sortend - point; 00121 for (int i = 0; i < n; ++i) point[i] = point[i+1]; 00122 *m_sortend = T(0); 00123 } 00124 00125 MedianFilter(const MedianFilter &); // not provided 00126 MedianFilter &operator=(const MedianFilter &); // not provided 00127 }; 00128 00129 #endif 00130