Алгоритм скользящего среднего (Simple Moving Average)

    Возникла задача реализовать на С++ алгоритм скользящего среднего, который находит широкое применение в обработке сигналов и статистике.
    image
    За основу была взята функция smooth в MATLAB.
    Данную функцию можно использовать для фильтрации сигналов. В качестве входных параметров определяются массив данных и окно усреднения.
    Кому интересно, прошу под кат


    Итак, есть несколько реализаций данного алгоритма. Рассмотрим самый простой из них:
    предположим, размер окна усреднения будет равен 5, тогда на каждом шаге усреднения берется текущее значение, к нему прибавляются 4 предыдущих и результат делится на 5. Очевидная проблема здесь в инициализации алгоритма, сначала нужно накопить определенное количество данных, не меньшее, чем окно усреднения.

    В MATLAB алгоритм фильтрации с помощью скользящего среднего реализован в функции smooth
    Пример использования smooth(input,window),
    где input — массив входящих данных
    window — окно усреднения.
    Изменяя параметр window можно получить большую или меньшую степень сглаживания данных:
    image

    Исходник, реализующий данный пример представлен ниже:
    clear all;
    %% Параметры
    t=1;% длительность сигнала
    Fd=512;% Частота дискретизации (Гц)
    A1=1;% Амплитуда синусоиды
    F1=10;% Частота синусоиды (Гц)
    SNR=0.3;% Соотношение сигнал/шум
    T=0:1/Fd:t;% Массив отсчетов
    Noise=A1*SNR*randn(1,length(T));%Сгенерированный шум
    Signal=A1*sind((F1*360).*T);%Сгенерированный сигнал
    SN=Signal+Noise;
    figure(1);
    subplot(4,1,1);
    plot(SN);
    subplot(4,1,2);
    plot(smooth(SN,3));
    subplot(4,1,3);
    plot(smooth(SN,8));
    subplot(4,1,4);
    plot(smooth(SN,20));
    


    Для компенсации задержки в обработке сигнала, MATLAB использует динамически изменяемый размер окна, суть метода иллюстрирует следующий пример:
    image

    Собственную реализацию данного алгоритма я сначала написал в MATLAB, а затем перенес на С++

    MATLAB edition:

    %my_smooth
    
    %в случае, если размер окна четный, увеличиваем его на 1 для симметрии;
    window = 5;
    if(mod(window,2)==0)
        window=window+1;
    end
    hw=(window-1)/2; %размах окна влево и вправо от текущей позиции
       
    n=length(Signal);
    result=zeros(n,1);
    result(1)=SN(1);  %первый элемент берем из исходного массива SN как есть
    
    for i=2:n              %организовываем цикл по числу элементов
        init_sum = 0;
        if(i<=hw)        %если индекс меньше половины окна, мы находимся в начале массива, 
                               %нужно брать окно меньшего размера
            k1=1;         %в качестве начала окна берем первый элемент
            k2=2*i-1;    %конец окна 
            z=k2;          %текущий размер окна
        elseif (i+hw>n) %если индекс+половина окна больше n - мы приближаемся к концу массива и размер окна
                                %также нужно уменьшать
            k1=i-n+i;     %начало окна 
            k2=n;          %конец окна - последний элемент массива
            z=k2-k1;     %размер окна
        else                 %если первые два условия не выполняются, мы в середине массива
            k1=i-hw;
            k2=i+hw;
            z=window;
        end
        
        for j=k1:k2                          %организуем цикл от начала до конца окна
           init_sum=init_sum+SN(j); %складываем все элементы
        end
        result(i)=init_sum/(z);        %и делим на текущий размер окна
    end
    


    Результат работы данной программы абсолютно совпадает с матлабовским smooth при нечетных размерах окна и несколько отличается при четном его значении (четные окна матлаб считает чуть иначе)
    Исходный m-файл можно взять тут

    С++ Edition

    void smooth(double *input, double *output, int n, int window)
    {
       int i,j,z,k1,k2,hw;
       double tmp;
       if(fmod(window,2)==0) window++;
       hw=(window-1)/2;
       output[0]=input[0];
    
       for (i=1;i<n;i++){
           tmp=0;
           if(i<hw){
               k1=0;
               k2=2*i;
               z=k2+1;
           }
           else if((i+hw)>(n-1)){
               k1=i-n+i+1;
               k2=n-1;
               z=k2-k1+1;
           }
           else{
               k1=i-hw;
               k2=i+hw;
               z=window;
           }
    
           for (j=k1;j<=k2;j++){
               tmp=tmp+input[j];
           }
           output[i]=tmp/z;
       }
    


    спасибо за внимание, конструктивная критика приветствуется
    p.s.:
    Алгоритм можно оптимизировать по скорости работы изменив подсчет суммы:
    image
    Видно, что для подсчета суммы элементов на 4-м шаге нужно из суммы на третьем шаге вычесть 1-й элемент массива (2, отмечен красным) и прибавить 6-й элемент (8, желтая клетка).
    На следующем шаге процедура повторяется.

    Данный подход будет эффективным при большом размере окна усреднения
    Поделиться публикацией
    AdBlock похитил этот баннер, но баннеры не зубы — отрастут

    Подробнее
    Реклама
    Комментарии 15
    • 0
      <source=c>if(fmod(window,2)==0)
      fmod для целых чисел да ещё и сравнение на строгое равенство чисел с плавающей запятой…

      > C++ edition
      и где тут C++? Это си, не более.

      На каждый элемент результирующего массива лишние ветвления: особые случаи для начала и конца нужны только в начале и в конце, а основной объём данных должен обрабатываться без этого.

      Трюк с вычитанием хорош только если вы уверены что не произойдёт никаких катастрофических округлений.
      • –1
        признаю, по скорости работы код далеко не оптимален и оптимизировать есть куда. Совместными усилиями можем поправить ))
      • 0
        А зачем Вы окно изменяете?
        Еще можно бы добавить как считать — задним окном, центром или наперед.
        И, мне кажется, не стоит начало и конец отдельно считать, просто значения которые меньше окна не выводить вообще.
        • +1
          по условиям задачи входной массив должен быть равен выходному, поэтому не выводить данные не получится
        • 0
          Подходит только для квазистатических сигналов. Представьте сигнал у которого амплитуда меняется в разы или у которого собственная частота очень близка к частоте его дискретизации.
          • 0
            Цифровая обработка такого сигнала дрлжна вестись с оглядкоц на теорему Котельникова. Тоесть оцифровав таким образом сигнал вы должны понимать что то что вы получили в цифре может не совсем соответствовать оригиналу.
          • 0
            Для сигнала с выборосами логично будет использовать сглаживание скользящей медианой. А вообще, естественно, эти «скользящие» алгоритмы применяются лишь для обработки сигналов в реальном времени, либо для обработки сигналов, у которых не прослеживается четкой структуры.
          • 0
            Мне кажется, на практике для сглаживания временного ряда чаще используют реккурентную формулу

            output[i] = alpha*input[i] + (1-alpha)*output[i-1]
            0<=alpha<=1

            Она очень хорошо выравнивает, вычислительно проста, имеет всего 1 параметр, и этот параметр имеет четкий физический смысл в терминах теории автоматического управления.
            • +1
              Это другой метод, не Simple Moving Average, а Exponential Moving Average.
            • +4
              Кто-то упомянул C++? Вот мой вариант:
              #include <iostream>
              #include <iterator>
              #include <algorithm>
              #include <queue>
              #include <vector>
              
              class MovingAverage {
                  std::queue<double> window; // окно
                  size_t period; // максимальный размер окна
                  double sum; // сумма элементов в окне
              public:
                  MovingAverage(size_t period) : period(period) {}
                  double operator()(double num) {
                      sum += num;
                      window.push(num);
                      if (window.size() > period) {
                          sum -= window.front();
                          window.pop();
                      }
                      return sum / window.size();
                  }
              };
              
              int main() {
                  using namespace std;
                  double indata[] = {1, 2, 3, 2, 4, 5, 4, 5, 4, 5, 6, 2, 5, 6, 6, 7};
                  size_t size = sizeof(indata) / sizeof(double);
                  vector<double> out(size);
              
                  // применение функтора MovingAverage к исходному массиву
                  transform(indata, indata + size, out.begin(), MovingAverage(5));
              
                  // вывод результирующего массива
                  copy(out.begin(), out.end(), ostream_iterator<double>(cout, "\n"));
              }


              Кроме того, такой функтор удобно использовать в классической ситуации, когда данные приходят в реальном времени и по новому значению нужно достроить график скользящей.
              • 0
                Не соответствует стандарту. std::transform не обязан в функтор кормить элементы по порядку.

                Effects: Assigns through every iterator i in the range [result,result + (last1 — first1)) a new corresponding value equal to op(*(first1 + (i — result)) or binary_op(*(first1 + (i — result)), *(first2 + (i — result))).

                • +1
                  Что интересно, эта цитата из стандарта C++11 не накладывает ограничения на порядок обхода. По крайней мере не так явно, как в старой версии стандарта:
                  25.2.3 Transform [lib.alg.transform]
                  Requires: op and binary_op shall not have any side effects.

                  Мне не очень понятно, зачем нужно было убирать. Вот оставшиеся ограничения:
                  Requires: op and binary_op shall not invalidate iterators or subranges, or modify elements in the ranges [first1,last1], [first2,first2 + (last1 — first1)], and [result,result + (last1 — first1)].

                  Parallel STL, впрочем, пока никто не отменял, поэтому спорить не буду. Пусть будет так:
                  ...
                      std::vector<double> &out;
                  
                  public:
                      MovingAverage(std::vector<double> &out, size_t period) : out(out), period(period) {}
                  
                      void operator()(double num) {
                          sum += num;
                          window.push(num);
                          if (window.size() > period) {
                              sum -= window.front();
                              window.pop();
                          }
                          out.push_back(sum / window.size());
                      }


                      vector<double> out;
                      out.reserve(size);
                      for_each(indata, indata + size, MovingAverage(out, 5));
                  
              • 0
                Я начал писать статьи в разделе habrahabr.ru/blogs/bioinformatics/, в последней статье понял, что многие алгоритмы описаны на хабре, но с ходу я не знаю подходящий он или нет.

                Хочу просить помощи. В этой статье habrahabr.ru/blogs/bioinformatics/137453/ приведены графики, частично похожие на приведеные здесь. Можно ли представленный тут алгоритм, или может существует другой, повернуть так, чтобы он правильно идентефицировал тип графика? Типа острые пики, пологие пики, «забор»?
                • НЛО прилетело и опубликовало эту надпись здесь

                Только полноправные пользователи могут оставлять комментарии. Войдите, пожалуйста.