ФИЛЬТРЫ ЗНАЧЕНИЙ ДЛЯ ARDUINO

Очень часто приходится работать с датчиками, и, соответственно, с получаемыми с них значениями. Не все датчики выдают стабильные значения, иногда значение “шумит” вокруг некоторой величины. Как избавиться от шума? Программный фильтр! Но какой? Приведённые ниже алгоритмы фильтров собраны в библиотеке GyverFilters для более удобного использования, подробнее о библиотеке читайте ниже.

Среднее арифметическое


Среднее арифметическое знакомо вам со школы: суммировать несколько измерений и делить сумму на их количество. Рассмотрим пример усреднения значения с аналогового пина:

/*
  Элементарная реализация среднего арифметического. Сложили NUM_READINGS измерений,
  затем разделили сумму на NUM_READINGS и всё!
  Является "частным случаем" предыдущего фильтра
  Время выполнения примерно равно: 10 значений 50 мкс, 50 значений 92 мкс, 100 значений 146 мкс
*/
#define NUM_READINGS 500
int average;
void setup() {
  Serial.begin(9600);
}
void loop() {
  long sum = 0;                                  // локальная переменная sum
  for (int i = 0; i < NUM_READINGS; i++) {      // согласно количеству усреднений
    sum += analogRead(0);                        // суммируем значения с любого датчика в переменную sum
  }
  average = sum / NUM_READINGS;                  // находим среднее арифметическое, разделив сумму на число измерений
  Serial.println(average);                       // для примера выводим в порт
}

Можно делать чуть по-другому: хранить массив нескольких предыдущих измерений, новые измерения записывать в старые ячейки, и делить на количество при добавлении нового значения!

/*
   Готовая функция для вычисления среднего арифметического
   Принимает новые значения, суммирует их в своём массиве
*/
#define NUM_AVER 10       // выборка (из скольки усредняем)
long average;             // перем. среднего
int valArray[NUM_AVER];   // массив
byte idx = 0;             // индекс
int middleArifm(int newVal) {       // принимает новое значение
  valArray[idx] = newVal;           // пишем каждый раз в новую ячейку
  if (++idx >= NUM_AVER) idx = 0;   // перезаписывая самое старое значение
  average = 0;                      // обнуляем среднее
  for (int i = 0; i < NUM_AVER; i++) {
    average += valArray[i];         // суммируем
  }
  average /= NUM_AVER;              // делим
  return average;                   // возвращаем
}

Вот так среднее арифметическое с размером буфера 10 справляется с равномерно растущим сигналом + случайные выбросы. Синий график – реальное значение, красный – фильтрованное.

Бегущее среднее


Бегущее среднее (англ. Running Average) – практически то же среднее арифметическое, но алгоритм гораздо проще – вместо массива у нас одно значение, которое обновляется из нового значения и коэффициента усиления, т.е. данный алгоритм имеет ещё и настройку “резкости” фильтрации. Классический алгоритм выглядит так: filteredValue = newValue * k + filteredValue * (1 – k), где filteredValue – фильтрованная величина, newValue – новая величина для фильтрации, а k – коэффициент фильтрации, от 0 до 1 (float). Данный алгоритм можно сократить до filtered_val += (val – filtered_val) * k. Важный момент! Фильтр должен работать по времени, т.е. при помощи задержек или таймера на millis().

/*
  Простейший фильтр: запаздывающий, бегущее среднее, "цифровой фильтр", фильтр низких частот - это всё про него любимого
  Имеет две настройки: постоянную времени FILTER_STEP (миллисекунды), и коэффициент "плавности" FILTER_COEF
  Данный фильтр абсолютно универсален, подходит для сглаживания любого потока данных
  При маленьком значении FILTER_COEF фильтрованное значение будет меняться очень медленно вслед за реальным
  Чем больше FILTER_STEP, тем меньше частота опроса фильтра
  Сгладит любую "гребёнку", шум, ложные срабатывания, резкие вспышки и прочее говно. Пользуюсь им постоянно
*/
#define FILTER_STEP 5
#define FILTER_COEF 0.05
int val;
float val_f;
unsigned long filter_timer;
void setup() {
  Serial.begin(9600);  
}
void loop() {
  
  if (millis() - filter_timer > FILTER_STEP) {
    filter_timer = millis();    // просто таймер
    // читаем значение (не обязательно с аналога, это может быть ЛЮБОЙ датчик)
    val = analogRead(0);
    // основной алгоритм фильтрации. Внимательно прокрутите его в голове, чтобы понять, как он работает
    val_f = val * FILTER_COEF + val_f * (1 - FILTER_COEF);
    // для примера выведем в порт
    Serial.println(val_f);
  }
}

У фильтра “бегущее среднее” не один настраиваемый параметр, как может показаться на первый взгляд. Помимо коэффициента фильтрации k очень большую роль играет время итерации, то есть период вызова фильтра. В реальном коде фильтр вызывается с определённым промежутком времени, чтобы фильтровать шумы. Я обычно настраиваю фильтр вручную по графику, который строится средствами Arduino IDE или программой Serial Port Plotter. Задаюсь периодом итерации и подгоняю k, пока не станет “хорошо”. Но есть и аналитический способ расчёта коэффициента фильтрации (или времени итерации).
При выборе значения коэффициента k необходимо отталкиваться от того, какие изменения сигнала нам интересны, а какие мы будем считать за шум. Сделать это можно с помощью следующего выражения:
t = dt * (1 / k – 1)
где t — период времени, который отделяет слишком быстрые изменения от требуемых; dt — время итерации (период вызова фильтра).
Например, если в нашем случае с потенциометром, k = 0,1, а время между двумя измерениями dt = 20 мс, то время t = (1-0.1) * 0,02 / 0.1 = 0,18 сек. То есть все изменения сигнала, которые длятся меньше 0,18 секунд будут подавляться. Во втором случае (при k = 0,3), мы получим t = 0,047 сек. Вникнув в эту связь, можно настроить фильтр, всего лишь глянув на график сырого значения!

Вот так бегущее среднее справляется с равномерно растущим сигналом + случайные выбросы. Синий график – реальное значение, красный – фильтрованное с коэффициентом 0.1, зелёное – коэффициент 0.5.

Медианный фильтр


Медианный фильтр – в отличие от среднего арифметического, находит “среднее среди средних”, а не среди всех, то есть отсеивает, именно отсеивает резкие выбросы!

/*
  Медианный фильтр — довольно простая и интересная штука. Берёт значения и выбирает из них среднее.
  Не усредняет, а именно ВЫБИРАЕТ, отбрасывает все сильно отличющиеся.
  Время выполнения близко к нулю мкс
  Простой пример, чем отличается медианный фильтр от среднего арифметического:
  Возьмём числа 3, 4, 50. Среднее арифметическое даст нам 19. Целью медианного фильтра является фильтрация
  резких скачков, и после фильтрации он даст нам 4, как среднее между 3 и 50, а 50 будет отброшено как скачок.
  В данном скетче реализована фильтрация по трём значениям. Если интересен вариант с фильтрацией более трёх значений,
  то добро пожаловать в исходную статью. Осторожно, жесть. http://tqfp.org/programming/mediannyy-filtr-na-sluzhbe-razrabotchika.html
*/
int val[3];
int val_filter;
byte index;
void setup() {
  Serial.begin(9600);
}
void loop() {
  if (++index > 2) index = 0; // переключаем индекс с 0 до 2 (0, 1, 2, 0, 1, 2…)
  val[index] = analogRead(0); // записываем значение с датчика в массив
  // фильтровать медианным фильтром из 3ёх ПОСЛЕДНИХ измерений
  val_filter = middle_of_3(val[0], val[1], val[2]);
  Serial.println(val_filter); // для примера выводим в порт
}
// медианный фильтр из 3ёх значений
int middle_of_3(int a, int b, int c) {
  int middle;
  if ((a <= b) && (a <= c)) {
    middle = (b <= c) ? b : c;
  }
  else {
    if ((b <= a) && (b <= c)) {
      middle = (a <= c) ? a : c;
    }
    else {
      middle = (a <= b) ? a : b;
    }
  }
  return middle;
}

На скриншоте представлен результат работы фильтра против выбросов: синий график – реальное значение, красный – фильтрованное. Фильтр просто отсекает выбросы. Чем выше размерность фильтра, тем лучше отсекаются частые выбросы.

Фильтр с размерностью больше трёх представляет собой гораздо более громоздкий алгоритм, почитать о нём можно здесь. Реализацию для Arduino смотрите в библиотеке GyverFilters.

Альфа-бета фильтр


Альфа-Бета фильтр является разновидностью фильтра Калмана для одномерного случая, подробнее читайте на Википедии. Реализация описанного алгоритма для Arduino выглядит вот так:

float dt;  // период дискретизации
float sigma_process;  // отклонение процесса
float sigma_noise;  // шум

// служебные переменные фильтра
float xk_1, vk_1, a, b;
float xk, vk, rk;
float xm;

void setup() {
  // расчёт коэффициентов A и B
  float lambda = (float)sigma_process * dt * dt / sigma_noise;
  float r = (4 + lambda - (float)sqrt(8 * lambda + lambda * lambda)) / 4;
  a = (float)1 - r * r;
  b = (float)2 * (2 - a) - 4 * (float)sqrt(1 - a);

  int value;  // value - фильтруемая величина
  // процесс фильтрации
  xm = value;
  xk = xk_1 + ((float) vk_1 * dt );
  vk = vk_1;
  rk = xm - xk;
  xk += (float)a * rk;
  vk += (float)( b * rk ) / dt;
  xk_1 = xk;
  vk_1 = vk;

  // xk_1 - результат фильтрации
}

void loop() {

}

Упрощённый фильтр Калмана


Данный алгоритм я нашёл на просторах Интернета, источник потерял. В фильтре настраивается разброс измерения (ожидаемый шум измерения), разброс оценки (подстраивается сам в процессе работы фильтра, можно поставить таким же как разброс измерения), скорость изменения значений (0.001-1, варьировать самому).

float errmeasure = 40; // разброс измерения
float errestimate = 40;  // разброс оценки
float q = 0.5;  // скорость изменения значений

float currentestimate = 0.0;
float lastestimate = 0.0;
float kalmangain = 0.0;

void setup() {

}

// функция фильтрации
float filter(int value) {
  kalmangain = errestimate / (errestimate + errmeasure);
  currentestimate = lastestimate + kalmangain * (value - lastestimate);
  errestimate =  (1.0 - kalmangain) * errestimate + fabs(lastestimate - currentestimate) * q;
  lastestimate = currentestimate;
  return currentestimate;
}

void loop() {

}

Вот так упрощённый Калман справляется с равномерно растущим сигналом + случайные выбросы. Синий график – реальное значение, красный – фильтрованное.

Линейная аппроксимация


Про теорию метода наименьших квадратов можно почитать тут. Что делает предложенная мной реализация: берёт два массива (например: один значения, второй – время) и рассчитывает параметры для уравнения прямой линии, которая проходит равноудалённо от всех значений массива. Уравнение прямой вида у = A*x + B. Реализация:

int32_t sumX, sumY, sumX2, sumXY;
float a, b, delta;

void setup() {
  int x_array[] = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10};
  int y_array[] = {1, 5, 2, 8, 3, 9, 10, 5, 15, 12};
  byte arrSize = 10;  // размер массива

  for (int i = 0; i < arrSize; i++) {    // для всех элементов массива
    sumX += x_array[i];
    sumY += (long)y_array[i];
    sumX2 += x_array[i] * x_array[i];
    sumXY += (long)y_array[i] * x_array[i];
  }
  a = (long)arrSize * sumXY;             // расчёт коэффициента наклона приямой
  a = a - (long)sumX * sumY;

  // получаем a, b и delta
  a = (float)a / (arrSize * sumX2 - sumX * sumX);
  b = (float)(sumY - (float)a * sumX) / arrSize;
  delta = a * arrSize;          // расчёт изменения

}

void loop() {
}

Результат работы из примера выше, сравнённый с расчётом линии тренда в Excel

Какой фильтр выбрать?


Выбор фильтра зависит от типа сигнала и ограничений по времени фильтрации. Среднее арифметическое – самый простой и быстрый, но и результат соответствующий, из настроек только ширина “окна” и период вызова. Для большинства ситуаций подходит бегущее среднее, он довольно быстрый и даёт хороший результат на правильной настройке. Медианный фильтр 3-го порядка тоже очень быстрый, но он может только отсеить выбросы, сам сигнал он не сгладит. Медиана большего порядка является довольно таки громоздким алгоритмом, но работает уже лучше. Очень хорошо работает медиана 3 порядка + бегущее среднее, получается сглаженный сигнал с отфильтрованными выбросами (сначала фильтровать медианой, потом бегущим). AB фильтр и фильтр Калмана – отличные фильтры, справляются с шумным сигналом не хуже связки медиана + бегущее среднее, но нуждаются в тонкой настройке, также они довольно громоздкие с точки зрения кода. Линейная аппроксимация – инструмент специального назначения, позволяющий буквально предсказывать поведение значения за период – ведь мы получаем уравнение прямой.

Для “побаловаться” предлагаю скетч, в котором генерируется синусоида со случайными шумами. В скетч уже встроено бегущее среднее для демонстрации:

void setup() {
  Serial.begin(9600);
}

float fil = 0;
int degree = 0;

void loop() {
  degree += 5;
  int valNoise = (float)sin(radians(degree)) * 100 + (random(5) == 0 ? 1 : 0) * (random(0, 2) ? -1 : 1) * 20;

  fil += (valNoise - fil) * 0.3;
  Serial.print(valNoise);
  Serial.print(' ');
  Serial.println(fil);
  delay(20);
}

БИБЛИОТЕКА GYVERFILTER

GyverFilters v1.5

Страница на GitHub
свежий релиз

Библиотека с некоторыми удобными фильтрами для Arduino:

  • GFilterRA – компактная альтернатива фильтра экспоненциальное бегущее среднее (Running Average)
  • GMedian3 – быстрый медианный фильтр 3-го порядка (отсекает выбросы)
  • GMedian – медианный фильтр N-го порядка. Порядок настраивается в GyverFilters.h – MEDIAN_FILTER_SIZE
  • GABfilter – альфа-бета фильтр (разновидность Калмана для одномерного случая)
  • GKalman – упрощённый Калман для одномерного случая (на мой взгляд лучший из фильтров)
  • GLinear – линейная аппроксимация методом наименьших квадратов для двух массивов

Поддерживаемые платформы: все Arduino

ДОКУМЕНТАЦИЯ


// ============== Бегущее среднее ==============
GFilterRA();                                // инициализация фильтра
GFilterRA(float coef, uint16_t interval);   // расширенная инициализация фильтра (коэффициент, шаг фильтрации)
void setCoef(float coef);                   // настройка коэффициента фильтрации (0.00 - 1.00). Чем меньше, тем плавнее
void setStep(uint16_t interval);            // установка шага фильтрации (мс). Чем меньше, тем резче фильтр
float filteredTime(int16_t value);          // возвращает фильтрованное значение с опорой на встроенный таймер	
float filtered(int16_t value);              // возвращает фильтрованное значение
float filteredTime(float value);            // возвращает фильтрованное значение с опорой на встроенный таймер	
float filtered(float value);                // возвращает фильтрованное значение
// ============== Медиана из трёх ==============
GMedian3();                     // инициализация фильтра
uint16_t filtered(uint16_t);    // возвращает фильтрованное значение
// ============== Медиана из MEDIAN_FILTER_SIZE (настраивается в GyverFilters.h) ==============
GMedian();                      // инициализация фильтра
uint16_t filtered(uint16_t);    // возвращает фильтрованное значение
// ============== Альфа-Бета фильтр ==============
GABfilter(float delta, float sigma_process, float sigma_noise);
// период дискретизации (измерений), process variation, noise variation
void setParameters(float delta, float sigma_process, float sigma_noise);
// период дискретизации (измерений), process variation, noise variation
float filtered(float value);      // возвращает фильтрованное значение
// ============== Упрощённый Калман ==============
GKalman(float mea_e, float est_e, float q);
// разброс измерения, разброс оценки, скорость изменения значений
GKalman(float mea_e, float q);
// разброс измерения, скорость изменения значений (разброс измерения принимается равным разбросу оценки)
void setParameters(float mea_e, float est_e, float q);
// разброс измерения, разброс оценки, скорость изменения значений
void setParameters(float mea_e, float q);
// разброс измерения, скорость изменения значений (разброс измерения принимается равным разбросу оценки)
float filtered(float value);     // возвращает фильтрованное значение
// ============== Линейная аппроксимация ==============
void compute(int *x_array, int *y_array, int arrSize);		// аппроксимировать
float getA();		// получить коэффициент А
float getB();		// получить коэффициент В
float getDelta();	// получить аппроксимированное изменение

ПРИМЕРЫ


#include "GyverFilters.h"
GFilterRA analog0;    // фильтр назовём analog0
void setup() {
Serial.begin(9600);
// установка коэффициента фильтрации (0.0... 1.0). Чем меньше, тем плавнее фильтр
analog0.setCoef(0.01);
// установка шага фильтрации (мс). Чем меньше, тем резче фильтр
analog0.setStep(10);
}
void loop() {
Serial.println(analog0.filteredTime(analogRead(0)));
}
/*
Пример использования быстрого медианного фильтра 3 порядка
*/
#include "GyverFilters.h"
GMedian3 testFilter;
void setup() {
Serial.begin(9600);
}
void loop() {    
int value = analogRead(0);
// добавляем шум "выбросы"
value += random(2) * random(2) * random(-1, 2) * random(50, 250);
Serial.print("$"); 
Serial.print(value);
Serial.print(" ");  
value = testFilter.filtered(value);
Serial.print(value);
Serial.println(";"); 
delay(80);
}
/*
Пример использования медианного фильтра.
Порядок фильтра настраивается в GyverHacks.h - MEDIAN_FILTER_SIZE
*/
#include "GyverFilters.h"
GMedian testFilter;
void setup() {
Serial.begin(9600);
}
void loop() {
delay(80);
int value = analogRead(0);
// добавляем шум "выбросы"
value += random(2) * random(2) * random(-1, 2) * random(50, 250);
Serial.print("$"); 
Serial.print(value);
Serial.print(" ");  
value = testFilter.filtered(value);
Serial.print(value);
Serial.println(";"); 
}
/*
Пример альфа-бета фильтра
*/
#include "GyverFilters.h"
// параметры: период дискретизации (измерений), process variation, noise variation
GABfilter testFilter(0.08, 40, 1);
void setup() {
Serial.begin(9600);
}
void loop() {
delay(80);
int value = analogRead(0);
value += random(2) * random(-1, 2) * random(10, 70);
Serial.print("$");
Serial.print(value);
Serial.print(" ");
value = testFilter.filtered((int)value);
Serial.print(value);
Serial.println(";");
}
/*
Пример простого одномерного фильтра
*/
#include "GyverFilters.h"
// параметры: разброс измерения, разброс оценки, скорость изменения значений
// разброс измерения: шум измерений
// разброс оценки: подстраивается сам, можно поставить таким же как разброс измерения
// скорость изменения значений: 0.001-1, варьировать самому
GKalman testFilter(40, 40, 0.5);
// также может быть объявлен как (разброс измерения, скорость изменения значений)
// GKalman testFilter(40, 0.5);
void setup() {
Serial.begin(9600);
}
void loop() {
delay(80);
int value = analogRead(0);
value += random(2) * random(-1, 2) * random(10, 70);
Serial.print("$");
Serial.print(value);
Serial.print(" ");
value = testFilter.filtered((int)value);
Serial.print(value);
Serial.println(";");
}
/*
Сравнение калмана и бегущего среднего
*/
#include "GyverFilters.h"
// параметры: разброс измерения, разброс оценки, скорость изменения значений
// разброс измерения: шум измерений
// разброс оценки: подстраивается сам, можно поставить таким же как разброс измерения
// скорость изменения значений: 0.001-1, варьировать самому
GKalman kalman(90, 90, 0.5);
GFilterRA average(0.5, 80);
void setup() {
Serial.begin(9600);
}
void loop() {  
int value = analogRead(0);
value += random(2) * random(-1, 2) * random(50, 100);
Serial.print("$");
Serial.print(value);
Serial.print(" ");  
Serial.print((int)kalman.filtered(value));
Serial.print(" ");
Serial.print((int)average.filtered(value));
Serial.println(";");
delay(80);
}
/*
Пример линейной аппроксимации методом наименьших квадратов
Два массива: по оси Х и по оси У
Линейная аппроксимация повозоляет получить уравнение прямой,
равноудалённой от точек на плоскости ХУ. Удобно для расчёта
роста изменяющейся шумящей величины. Уравнение вида у = A*x + B
В папке с данным примером есть скриншот из excel,
иллюстрирующий работу аппроксимации с такими же исходными
*/
// два массива с данными (одинаковой размероности и размера)
int x_array[] = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10};
int y_array[] = {1, 5, 2, 8, 3, 9, 10, 5, 15, 12};
#include 
GLinear test;
void setup() {
Serial.begin(9600);
// передаём массивы и размер одного из них
test.compute((int*)x_array, (int*)y_array, sizeof(x_array));
// Уравнение вида у = A*x + B
Serial.println(test.getA());  // получить коэффициент А
Serial.println(test.getB());  // получить коэффициент В
Serial.println(test.getDelta());  // получить изменение (аппроксимированное)
}
void loop() {
}
/*
Пример линейной аппроксимации методом наименьших квадратов
Два массива: по оси Х и по оси У
Наполнение массивов осуществляется динамически: сдвигом и записью в крайнюю ячейку,
то есть аппроксимация по последним ARRAY_SIZE изменениям!!
*/
#define ARRAY_SIZE 10   // размер пространства для аппроксимации
// два массива с данными (одинаковой размероности и размера)
int x_array[ARRAY_SIZE] = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10};  // ось x от 1 до 10, допустим СЕКУНД
int y_array[ARRAY_SIZE];    // значения по оси У будем брать с датчика
#include 
GLinear test;
void setup() {
Serial.begin(9600);
}
void loop() {
for (byte i = 0; i < ARRAY_SIZE - 1; i++) {    // счётчик от 0 до ARRAY_SIZE
y_array[i] = y_array[i + 1];    // сдвинуть массив давлений КРОМЕ ПОСЛЕДНЕЙ ЯЧЕЙКИ на шаг назад
}
// последний элемент массива теперь - новое значение (просто с аналог. датчика)
y_array[ARRAY_SIZE - 1] = analogRead(0);
// передаём массивы и размер одного из них
test.compute((int*)x_array, (int*)y_array, sizeof(x_array));
// по нашим исходным данным это будет производная, т.е. "изменение единиц в секунду"
Serial.println(test.getDelta());  // получить изменение (аппроксимированное)
delay(1000);  // секундная задержка
}
2019-09-15T15:26:02+03:00