0,0
рейтинг
11 декабря 2011 в 04:31

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

Данный подход будет эффективным при большом размере окна усреднения
Коваленко Александр @alk0v
карма
79,0
рейтинг 0,0
Пользователь
Реклама помогает поддерживать и развивать наши сервисы

Подробнее
Реклама

Самое читаемое Разработка

Комментарии (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/ приведены графики, частично похожие на приведеные здесь. Можно ли представленный тут алгоритм, или может существует другой, повернуть так, чтобы он правильно идентефицировал тип графика? Типа острые пики, пологие пики, «забор»?
    • 0
      skewness считай и kurtosis

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