Экспоненциальные фильтры скользящего среднего Печать
Добавил(а) microsin   

Экспоненциальный фильтр скользящего среднего (Exponential Moving Average Filter, EMA) это дискретный фильтр низкой частоты (low-pass, ФНЧ), с бесконечной импульной характеристикой (БИХ, или по-английски Infinite-Impulse Response, IIR). Он придает больший вес последним данным, дисконтируя старые данные экспоненциальным образом, и ведет себя аналогично дискретному низкочастотному RC-фильтру первого порядка.

Примечание: это перевод статьи [1].

В отличие от фильтра простого скользящего среднего (Simple Moving Average, SMA), большинство EMA-фильтров не оконные, и у них следующее значение зависит от всех предыдущих значений. Таким образом большинство фильтров EMA это форма IIR-фильтра, а SMA это фильтр с конечной импульсной характеристикой (КИХ, или по-английски Finite Impulse Response, FIR). Есть исключения, и вы действительно можете построить оконный экспоненциальный фильтр скользящего среднего, в котором коэффициенты взвешиваются экспоненциально.

EMA-фильтр благодаря своей простоте очень хорош для использования вместе с микроконтроллерами.

[Формула EMA]

Дифференциальное уравнение для exponential moving average фильтра:

y[i] = α ⋅ x[i] + (1 − α) ⋅ y[i − 1]

Здесь: 

y это выходное значение фильтра ([i] обозначает номер выборки выходного сигнала)
x это входной сигнал
α это константа, которая задает частоту среза фильтра (значение между 0 и 1)

Обратите внимание, что вычисления в формуле не требуют хранения нескольких предыдущих значений x, нужно запоминать только предыдущее значение y, что значительно упрощает вычисления и экономит память (это особенно важно для микроконтроллеров). Требуется только одно сложение, одно вычитание и два умножения.

Константа α определяет агрессивность фильтра. Она может меняться от 0 до 1 (включительно). По мере приближения α к нулю фильтр становится все агрессивнее. При α = 0 вход никак не влияет на выход (если фильтр запущен в такой конфигурации, на выходе всегда будет 0). При α стремящемся к 1 фильтр все меньше и меньше влияет на входные данные, пропуская их на выход, и при α = 1 фильтрация вообще не работает (данные проходят со входа на выход без изменения).

Фильтр называется экспоненциальным, потому что взвешенное влияние на выход предыдущих входных значений уменьшается экспоненциально, чем дальше был ввод во времени. Это можно увидеть в уравнении разности, если мы подставим в предыдущие входные данные:

y[i] = α ⋅ x[i] + (1-α) ⋅ y[i-1]
y[i] = α ⋅ x[i] + (1-α) ⋅ (α ⋅ x[i-1] + (1-α) ⋅ y[i-2])
y[i] = α ⋅ x[i] + (1-α) ⋅ (α ⋅ x[i-1] + (1-α) ⋅ (α ⋅ x[i-2] + (1-α) ⋅ y[i-3]))
...

EMA eq02                                (2)

Следующий код реализует базовый EMA-фильтр на C++, подходящий для микроконтроллеров и других встраиваемых устройств:

#include < stdio.h>
#include < cstdint>

class EmaFilter
{
public:
    EmaFilter(double alpha) :
        m_alpha(alpha), m_lastOutput(0.0) {}

    double Run(double input)
    {
        m_lastOutput = m_alpha * input + (1 - m_alpha) * m_lastOutput;
        return m_lastOutput;
    }
private:
    double m_alpha;
    double m_lastOutput;
};

int main() {
    EmaFilter emaFilter(0.5);

    for(uint32_t i = 0; i < 20; i++)
    {
        double input = 1.0;
        double output = emaFilter.Run(input);
        printf("%f\n", output);
    }
}

При настройке alpha=0.5, если жестко подать в итерациях цикла на вход фильтра значение 1.0 (что соответствует скачку входного сигнала от 0.0 до 1.0 в момент ремени t=0), то этот код напечатает следующее:

0.500000
0.750000
0.875000
0.937500
0.968750
0.984375
0.992188
0.996094
0.998047
0.999023
0.999512
0.999756
0.999878
0.999939
0.999969
0.999985
0.999992
0.999996
0.999998
0.999999

Примечание: вы можете поэкспериментировать с фильтром с помощью платформы https://replit.com/, см. [2].

[Частотная характеристика]

Частотную характеристику EMA-фильтра можно найти с помощью Z-преобразования. Начнем с выражения EMA-фильтра в домене времени:

y[i] = α ⋅ x[i] + (1-α) ⋅ y[i-1]                           (3)

Выполним над ней Z-преобразование:

Y(z) = αX(z) + (1-α)z-1Y(z)                               (4)

Найдем передаточную функцию H(z):

        Y(z)           α
H(z) = ------ = ---------------                           (5)
        X(z)     1 - (1-α)z-1

Эта передаточная функция может быть использована для создания гистограмм амплитуды и фазовой характеристики фильтра EMA. На графике ниже показана реакция фильтра EMA с α=0.25. Частота по оси x - нормализованная частота, в единицах Гц/выборка, благодаря чему график подойдет для любой частоты дискретизации.

EMA bode plot fig01

Рис. 1. Графики характеристик АЧХ и ФЧХ фильтра EMA для α=0.25.

Частоту среза (по уровню -3dB) EMA-фильтра можно вычислить по формуле:

       fs               α2
fc = ------cos-1[1 - --------]                          (6)
      2*pi             2(1-α)

Здесь fs это частота дискретизации, Гц.

Важное следствие из этого - ответ фильтра (АЧХ, ФЧХ) зависит как от α, так и от частоты дискретизации. Помните, что если вы измените частоту дискретизации, то понадобится также поменять α, чтобы получить такую же характеристику фильтра.

Имейте в виду, что не все допустимые значения α дадут допустимую частоту среза. Это когда отклик фильтра настолько плавный, что он не падает до -3db или ниже, пока не достигнет частоты Найквиста. Математически это происходит, когда:

       α2
|1 - --------| > 1                                       (7)
      2(1-α)

и, таким образом, обратная функция косинуса не определена. Выражение выше становится больше 1, когда α больше или равно:

α=0.82843                                                (8)

​[Импульсная характеристика]

Дискретная единичная функция выборки определяется как:

δ[n] = 1 при n=0
     = 0 при n!=0                                        (9)

Если мы будем использовать её в качестве входных данных для формулы (2), и учитывая, что большинство входных значений будет равно 0, то формулу (2) можно упростить до:

y[i] = α(1−α)n                                          (10)

​Исходя из этого, мы можем построить график того, как будет выглядеть ответ для импульса на входе. Как вы можете видеть на следующем графике, вывод начинается с y[0]=α, и затем спадает до 0. Сем больше α, тем быстрее падение.

EMA impulse response fig02

Рис. 2. Ответ на входной импульс EMA-фильтра при разных значениях α.

[Постоянная времени]

В отношении частотной характеристики и импульсной характеристики можно определить постоянную времени для фильтра EMA. Постоянная времени - это время, необходимое для достижения выхода 1 - 1/e (63.2%) от значения установившегося режима. Рис. 3 показывает количество итераций, требуемых для достижения заданной доли устойчивого состояния для различных α. Это предполагает тупенчатое изменение уровня на входе от 0 до 1 в момент t=0.

EMA time constant fig03

Рис. 3. Количество итераций, требуемых для достижения заданного дробного значения (т. е. постоянные времени) для разных значений α.

[Прайминг фильтра]

При инициализации фильтра важно начальное значение y. Если установить его в 0, то фильтру потребуется некоторое время, чтобы "разогреться" и выйти на правильный выход. Если установить его на первое входное значение, то фильтр получит "head start". Это может быть важным или неважным для вашего приложения. Часто используют двоичный флаг наподобие firstTime, который устанавливается в true, когда фильтр впервые создается, и затем устанавливается в false после обработки первого ввода. На первом проходе фильтр просто устанавливает выхо y в значение входа x. Лучше это называть заполнением (праймингом) фильтра, чтобы не перепутать с термином "инициализация" фильтра (т. е. с его созданием в конструкторе).

Ниже показан тот же пример класса C++, переработанный для прайминга:

class EmaFilterWithPriming
{
public:
    EmaFilterWithPriming(double alpha) :
        m_alpha(alpha), m_lastOutput(0.0), m_firstRun(true) {}
    double Run(double input)
    {
        if (m_firstRun)
        {
            m_firstRun = false;
            m_lastOutput = input; // Bypass filter and prime with input
        }
        else
        {
            m_lastOutput = m_alpha * input + (1 - m_alpha) * m_lastOutput;
        }
        return m_lastOutput;
    }
private:
    double m_alpha;
    double m_lastOutput;
    bool m_firstRun;
};

[Ссылки]

1. Exponential Moving Average (EMA) Filters site:blog.mbedded.ninja.
2. https://replit.com/@gbmhunter/cpp-ema-filter
3. Простой цифровой фильтр EMA.