c@389: /* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
c@390: /*
c@390:     QM DSP Library
c@390: 
c@390:     Centre for Digital Music, Queen Mary, University of London.
c@390:     This file Copyright 2010 Chris Cannam.
c@390: 
c@390:     This program is free software; you can redistribute it and/or
c@390:     modify it under the terms of the GNU General Public License as
c@390:     published by the Free Software Foundation; either version 2 of the
c@390:     License, or (at your option) any later version.  See the file
c@390:     COPYING included with this distribution for more information.
c@390: */
c@390: 
cannam@489: #ifndef QM_DSP_MEDIAN_FILTER_H
cannam@489: #define QM_DSP_MEDIAN_FILTER_H
c@389: 
c@389: #include <algorithm>
c@389: #include <cassert>
c@389: #include <cmath>
c@389: #include <iostream>
c@392: #include <vector>
c@389: 
c@389: template <typename T>
c@389: class MedianFilter
c@389: {
c@389: public:
c@389:     MedianFilter(int size, float percentile = 50.f) :
cannam@477:         m_size(size),
cannam@477:         m_frame(new T[size]),
cannam@477:         m_sorted(new T[size]),
cannam@477:         m_sortend(m_sorted + size - 1) {
c@389:         setPercentile(percentile);
cannam@477:         reset();
c@389:     }
c@389: 
c@389:     ~MedianFilter() { 
cannam@477:         delete[] m_frame;
cannam@477:         delete[] m_sorted;
c@389:     }
c@389: 
c@389:     void setPercentile(float p) {
c@389:         m_index = int((m_size * p) / 100.f);
c@389:         if (m_index >= m_size) m_index = m_size-1;
c@389:         if (m_index < 0) m_index = 0;
c@389:     }
c@389: 
c@389:     void push(T value) {
c@394:         if (value != value) {
c@407:             std::cerr << "WARNING: MedianFilter::push: attempt to push NaN, pushing zero instead" << std::endl;
c@407:             // we do need to push something, to maintain the filter length
c@407:             value = T();
c@394:         }
cannam@477:         drop(m_frame[0]);
cannam@477:         const int sz1 = m_size-1;
cannam@477:         for (int i = 0; i < sz1; ++i) m_frame[i] = m_frame[i+1];
cannam@477:         m_frame[m_size-1] = value;
cannam@477:         put(value);
c@389:     }
c@389: 
c@389:     T get() const {
cannam@477:         return m_sorted[m_index];
c@389:     }
c@389: 
c@393:     int getSize() const {
c@393:         return m_size; 
c@393:     }
c@393: 
c@389:     T getAt(float percentile) {
cannam@477:         int ix = int((m_size * percentile) / 100.f);
c@389:         if (ix >= m_size) ix = m_size-1;
c@389:         if (ix < 0) ix = 0;
cannam@477:         return m_sorted[ix];
c@389:     }
c@389: 
c@389:     void reset() {
cannam@477:         for (int i = 0; i < m_size; ++i) m_frame[i] = 0;
cannam@477:         for (int i = 0; i < m_size; ++i) m_sorted[i] = 0;
c@389:     }
c@389: 
c@392:     static std::vector<T> filter(int size, const std::vector<T> &in) {
c@392:         std::vector<T> out;
c@392:         MedianFilter<T> f(size);
c@392:         for (int i = 0; i < int(in.size()); ++i) {
c@392:             f.push(in[i]);
c@392:             T median = f.get();
c@392:             if (i >= size/2) out.push_back(median);
c@392:         }
c@392:         while (out.size() < in.size()) {
c@392:             f.push(T());
c@392:             out.push_back(f.get());
c@392:         }
c@392:         return out;
c@392:     }
c@392: 
c@389: private:
c@389:     const int m_size;
c@389:     T *const m_frame;
c@389:     T *const m_sorted;
c@389:     T *const m_sortend;
c@389:     int m_index;
c@389: 
c@389:     void put(T value) {
cannam@477:         // precondition: m_sorted contains m_size-1 values, packed at start
cannam@477:         // postcondition: m_sorted contains m_size values, one of which is value
cannam@477:         T *point = std::lower_bound(m_sorted, m_sortend, value);
cannam@477:         const int n = m_sortend - point;
cannam@477:         for (int i = n; i > 0; --i) point[i] = point[i-1];
cannam@477:         *point = value;
c@389:     }
c@389: 
c@389:     void drop(T value) {
cannam@477:         // precondition: m_sorted contains m_size values, one of which is value
cannam@477:         // postcondition: m_sorted contains m_size-1 values, packed at start
cannam@477:         T *point = std::lower_bound(m_sorted, m_sortend + 1, value);
cannam@477:         if (*point != value) {
cannam@477:             std::cerr << "WARNING: MedianFilter::drop: *point is " << *point
cannam@477:                       << ", expected " << value << std::endl;
cannam@477:         }
cannam@477:         const int n = m_sortend - point;
cannam@477:         for (int i = 0; i < n; ++i) point[i] = point[i+1];
cannam@477:         *m_sortend = T(0);
c@389:     }
c@390: 
c@390:     MedianFilter(const MedianFilter &); // not provided
c@390:     MedianFilter &operator=(const MedianFilter &); // not provided
c@389: };
c@389: 
c@389: #endif
c@389: