Алгоритм скользящего среднего (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, желтая клетка).
    На следующем шаге процедура повторяется.

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

    Similar posts

    Ads
    AdBlock has stolen the banner, but banners are not teeth — they will be back

    More

    Comments 15

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

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

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

      Трюк с вычитанием хорош только если вы уверены что не произойдёт никаких катастрофических округлений.
        –1
        признаю, по скорости работы код далеко не оптимален и оптимизировать есть куда. Совместными усилиями можем поправить ))
        0
        А зачем Вы окно изменяете?
        Еще можно бы добавить как считать — задним окном, центром или наперед.
        И, мне кажется, не стоит начало и конец отдельно считать, просто значения которые меньше окна не выводить вообще.
          +1
          по условиям задачи входной массив должен быть равен выходному, поэтому не выводить данные не получится
          0
          Подходит только для квазистатических сигналов. Представьте сигнал у которого амплитуда меняется в разы или у которого собственная частота очень близка к частоте его дискретизации.
            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/ приведены графики, частично похожие на приведеные здесь. Можно ли представленный тут алгоритм, или может существует другой, повернуть так, чтобы он правильно идентефицировал тип графика? Типа острые пики, пологие пики, «забор»?
                  • UFO just landed and posted this here

                  Only users with full accounts can post comments. Log in, please.