qm-dsp  1.8
MedianFilter.h
Go to the documentation of this file.
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