Нечёткие измерения, или фильтр Калмана

[We Are]

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

И выглядит это примерно так:

 

А хочется получить ровный график. Для этого есть несколько способов:

Один из простейших - усреднение, а точнее метод "скользящего окна":
Создается массив длиной, скажем, 50 измерений, первые 50 просто заполняют его, а вот начиная со следующего - рассчитывается среднее значение этих 50 измерений, выводится, самое раннее значение из массива удаляется, остальные смещаются вниз, сверху добавляется новое измеренное - и так далее.

#define LEN 50
float arr[ LEN + 1 ];
int cnt = 0;

float F(float in){
  float x = in;
  if(cnt < LEN){
    arr[cnt] = in;
  }else{
    arr[LEN] = in;
    x = 0;
    for(int i=0; i<LEN; i++){
      x += arr[i];
      arr[i] = arr[i + 1];
    } 
    x /= LEN;
  }
  cnt++;
  return x;
}

 

Второй способ - фильтр Калмана. Расписывать подробно не буду, потому что там статистика, ковариационное исчисление, матрицы - всё это довольно сложно и много.

Итоговый расчет может быть примерно по вот такому алгоритму (скопировано как есть):

    class KalmanFilterSimple1D
    {
        public double X0 {get; private set;} // predicted state
        public double P0 { get; private set; } // predicted covariance

        public double F { get; private set; } // factor of real value to previous real value
        public double Q { get; private set; } // measurement noise
        public double H { get; private set; } // factor of measured value to real value
        public double R { get; private set; } // environment noise

        public double State { get; private set; }
        public double Covariance { get; private set; }

        public KalmanFilterSimple1D(double q, double r, double f = 1, double h = 1)
        {
            Q = q;
            R = r;
            F = f;
            H = h;
        }

        public void SetState(double state, double covariance)
        {
            State = state;
            Covariance = covariance;
        }

        public void Correct(double data)
        {
            //time update - prediction
            X0 = F*State;
            P0 = F*Covariance*F + Q;

            //measurement update - correction
            var K = H*P0/(H*P0*H + R);
            State = X0 + K*(data - H*X0);
            Covariance = (1 - K*H)*P0;            
        }
    }



    // Применение...

    var fuelData = GetData();
    var filtered = new List<double>();

    var kalman = new KalmanFilterSimple1D(f: 1, h: 1, q: 2, r: 15); // задаем F, H, Q и R
    kalman.SetState(fuelData[0], 0.1); // Задаем начальные значение State и Covariance
    foreach(var d in fuelData)
    {
        kalman.Correct(d); // Применяем алгоритм

        filtered.Add(kalman.State); // Сохраняем текущее состояние 
    }

 

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

Каждое измерение содержит в себе некоторое фактическое значение + случайную ошибку, значит новое измерение будет содержать новое фактическое значение, близкое к предыдущему + новую случайную ошибку, и значит новое значение должно быть больше похоже на предыдущее, чем на новое измеренное, что дает такой расчет:

float prev = 0;
float k = 0.01;

float F(float in){
  float x = in;
  if(cnt > 0){
    x = k * in + (1 - k) * prev;
  }
  cnt++;
  prev = x;
  return x;
}

А теперь - результат обработки входных данных всеми тремя способами:

Если плохо видно, здесь три графика: красный - скользящее окно, зеленый - фильтр Калмана, синий - упрощенный фильтр Калмана

Разница видна только на начальном этапе, потом ее практически незаметно.

Применимо к случаям, когда измеряемая величина меняется медленно по сравнению с частотой измерений.

[Реклама]

Реклама: ООО ТД Промэлектроника, ИНН: 6659197470, Тел: 8 (800) 1000-321

[colorad]

Как-то программно моделировал RC НЧ фильтр . Запоминается значение напряжения на С . Для первого включения можно взять напряжение на С равное входному .

Рассчитываться ток заряда С через разницу входного напряжения и напряжения на С . И корректируется напряжение на С .  

Делал для того,  чтоб напряжение измеренное на аккумуляторе не сильно дергалось .

[Реклама]

20% скидка на весь каталог электронных компонентов в ТМ Электроникс!

Акция "Лето ближе - цены ниже", успей сделать выгодные покупки!

Плюс весь апрель действует скидка 10% по промокоду APREL24 + 15% кэшбэк и бесплатная доставка!

Перейти на страницу акции

Реклама: ООО ТМ ЭЛЕКТРОНИКС, ИНН: 7806548420, info@tmelectronics.ru, +7(812)4094849

[donec]

В 11.02.2024 в 18:05, We Are сказал:

а точнее метод "скользящего окна":

cnt++ , в итоге будет превышение индекса массива.

Как вариант, но для чисел с плавающей запятой лучше не применять, ошибка накапливается

#define LEN 50
int arr[LEN] = {0};
int cnt = 0;
int x = 0;

int F(int in)
{
	x -= arr[cnt];
	arr[cnt] = in;
	x += arr[cnt++];
	if (cnt > (LEN -1))
		cnt = 0;
	return (x / LEN);
}

 

[Реклама]

Особенности хранения литиевых аккумуляторов и батареек

Потеря емкости аккумулятора напрямую зависит от условий хранения и эксплуатации. При неправильном хранении даже самый лучший литиевый источник тока с превосходными характеристиками может не оправдать ожиданий. Технология, основанная на рекомендациях таких известных производителей литиевых источников тока, как компании FANSO и EVE Energy, поможет организовать правильный процесс хранения батареек и аккумуляторов. Подробнее>>

Реклама: АО КОМПЭЛ, ИНН: 7713005406, ОГРН: 1027700032161

[We Are]

13 hours ago, donec said:

cnt++ , в итоге будет превышение индекса массива.

Не уловил логики: а зачем делать так чтобы было превышение индекса массива и ошибки накапливались? 

Изменено 13 февраля пользователем We Are

[Реклама]

Выбираем схему BMS для корректной работы литий-железофосфатных (LiFePO4) аккумуляторов

 Обязательным условием долгой и стабильной работы Li-FePO4-аккумуляторов, в том числе и производства EVE Energy, является применение специализированных BMS-микросхем. Литий-железофосфатные АКБ отличаются такими характеристиками, как высокая многократность циклов заряда-разряда, безопасность, возможность быстрой зарядки, устойчивость к буферному режиму работы и приемлемая стоимость. Но для этих АКБ, также как и для других, очень важен контроль процесса заряда и разряда, а специализированных микросхем для этого вида аккумуляторов не так много. Инженеры КОМПЭЛ подготовили список имеющихся микросхем и возможных решений от разных производителей. Подробнее>>

Реклама: АО КОМПЭЛ, ИНН: 7713005406, ОГРН: 1027700032161

[BARS_]

Главная проблема любого фильтра - быстродействие. Там либо хорошо фильтруем, либо быстро отрабатываем изменения сигнала.  Особенно заметно при резком изменении сигнала.

 

В 11.02.2024 в 18:05, We Are сказал:

напряжение дергается

Как правило, при грамотной разработке схемы, достаточно обычного усреднения.

 

В 11.02.2024 в 18:05, We Are сказал:

координаты прыгают

А их фильтрацией GNSS приемник занимается.

 

 

[We Are]

4 minutes ago, BARS_ said:

А их фильтрацией GNSS приемник занимается.

Все равно прыгают. +- до 10 метров, а то и больше

5 minutes ago, BARS_ said:

Как правило, при грамотной разработке схемы, достаточно обычного усреднения.

Ну вот например показометр напряжения питания ESP8266 - там вообще схема не трогается никак. Разброс показаний - неточность самого ADC чипа.

Разница только в том, как считать: делать штук 50 измерений и выдавать результат, или делать одно измерение и выдавать результат по нему и сохраненному предыдущему. Это время выполнения. А итог одинаковый.

7 minutes ago, BARS_ said:

Главная проблема любого фильтра - быстродействие.

Истинно так.

Поэтому область применения ограничена медленноизменяющимися величинами - температура (тепловая инерция не даст скакать туда-сюда), напряжение аккумулятора, трафик типа "сколько прошло людей в час" и подобное.

[BARS_]

1 минуту назад, We Are сказал:

Все равно прыгают. +- до 10 метров, а то и больше

В дешманских приемниках да, в более-менее профессиональных - нет. А если еще и RTK им скормить, то точность миллиметровая будет.

 

2 минуты назад, We Are сказал:

показометр напряжения питания ESP8266

Вот именно, что показометр. У меня ESP32 без каких-либо ухищрений измеряет нормально напряжение.

[donec]

43 минуты назад, We Are сказал:

а зачем делать так чтобы было превышение индекса массива

Это к вашему коду относится, ошибка там.

[We Are]

25 minutes ago, donec said:

Это к вашему коду относится, ошибка там.

Покажите? 

[donec]

В 11.02.2024 в 18:05, We Are сказал:

#define LEN 50
float arr[ LEN + 1 ];
int cnt = 0;

float F(float in){
  float x = in;
  if(cnt < LEN){
    arr[cnt] = in;
  }else{
    arr[LEN] = in;
    x = 0;
    for(int i=0; i<LEN; i++){
      x += arr[i];
      arr[i] = arr[i + 1];
    } 
    x /= LEN;
  }
  cnt++;
  return x;
}

Смотрите

[We Are]

3 minutes ago, donec said:

Смотрите

Ну и? Внимательно гляньте-то.

[donec]

Что, ну и..

Каждый вызов функции F , будет увеличивать cnt на единицу, в конце концов произойдет переполнение числа int и cnt станет нулем, и начнется все заново

[We Are]

Just now, donec said:

Каждый вызов функции F , будет увеличивать cnt на единицу, в конце концов произойдет переполнение числа int и cnt станет нулем, и начнется все заново

Ну разумеется. Это называется "ограничения области применения": этот код не может работать вечно, он ограничен пределом int.
Можно расширить ну скажем до unsigned long long - но и там будет свой предел, просто больше.

В данном конкретном случае cnt увеличивается каждый раз не просто так, он там дальше в счетчике используется - именно поэтому он там и увеличивается, потому что нужно количество.
Можно было бы тупо убрать его в скобки - но тогда счет остановился бы, и пришлось бы делать второй счетчик, внешний к F...

Это, кстати, к вопросу "взять чужой код не думая"

[donec]

3 минуты назад, We Are сказал:

Это, кстати, к вопросу "взять чужой код не думая"

Вот кто-то и возьмет, а потом будет голову ломать почему раз в сутки глюк происходит

[We Are]

5 minutes ago, donec said:

Вот кто-то и возьмет, а потом будет голову ломать почему раз в сутки глюк происходит

Классика же: "я использовал код из примера в ардуине, а через 49 дней она зависла! Ууу, китайское качество!"

Просто millis() перешел через 0 , а суть примера была в работе с чипом, а не в обучении работе с байтами

Изменено 13 февраля пользователем We Are

[Алексей Еремеев]

Есть еще такая штука как   "Экспоненциальное скользящее среднее"

Очень неприхотлива к памяти, не нужно хранить массив данных, только одно последнее значение

  float alpha = 0.7; // 0 < alpha <= 1.0
  value = (alpha * new_measurement) + ((1-alpha) * value);

где

new_measurement - новое значение измерения

alpha - коэффициент усреднения 

value - результат, он-же аккумулятор (static или global)

[IMXO]

2 минуты назад, Алексей Еремеев сказал:

Есть еще такая штука

очень хреновая штука.... особенно в целочисленной математике,  потому как (1-alpha) определяет порог нечувствительности фильтра.

[We Are]

20 minutes ago, Алексей Еремеев said:

value = (alpha * new_measurement) + ((1-alpha) * value);

x = k * in + (1 - k) * prev;

Собственно, то же самое, третий вариант. Сразу не увидел.

Изменено 13 февраля пользователем We Are

[Алексей Еремеев]

нормальная это штука, для своих применений вполне

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

#define FILTER_SHIFT 3

uint16_t          dist_filter;       // filter integrator

inline uint8_t update_filter(uint8_t value)
{
  dist_filter -= (dist_filter >> FILTER_SHIFT);  // Update the filter with the old input
  dist_filter += value;                          // Update filter with new input
  return (uint8_t)(dist_filter >> FILTER_SHIFT); // Scale output for unity gain.
}
 

не забываем "прогреть" фильтр для первого значения

inline uint8_t warmup_filter(uint8_t value)
{
  dist_filter = ((uint16_t)value << FILTER_SHIFT);
  return value;
}

Изменено 13 февраля пользователем Алексей Еремеев

[We Are]

17 minutes ago, IMXO said:

очень хреновая штука.... особенно в целочисленной математике,  потому как (1-alpha) определяет порог нечувствительности фильтра

Это вот как раз тот третий вариант, который "упрощенный фильтр Калмана".

По графикам - сходятся данные, т.е. результат тот же самый.

[Алексей Еремеев]

1 минуту назад, We Are сказал:

Это вот как раз тот третий вариант, который "упрощенный фильтр Калмана".

Прошу прощения, не все портянки кода смотрел. Ну, калманом здесь и не пахнет, конечно... а так да - то-же самое.

Поставьте к = 0.05, например, будет похоже на скользящее по 1/0.05 = 20 точкам, зато без буфера

[We Are]

12 minutes ago, Алексей Еремеев said:

Поставьте к = 0.05, например, будет похоже на скользящее по 1/0.05 = 20 точкам, зато без буфера

да, такое соотношение

12 minutes ago, Алексей Еремеев said:

Ну, калманом здесь и не пахнет, конечно...

ну, как называли ))

Но с тем кодом, который "по Калману", график тоже сходится.
Вообще подозреваю что многие вещи которые можно сделать просто - можно обосновать и просчитать через сложные формулы расчета - вопрос в глубины погружения в математику и том, насколько всеобьемлющим должно быть решение.

Например, простое усреднение не умеет работать с "выбросами" - но медиана уже умеет, а кто-то просто обрезает всё что не укладывается в +-10% - и это тоже работает

 

[Алексей Еремеев]

6 минут назад, We Are сказал:

а кто-то просто обрезает всё что не укладывается в +-10% - и это тоже работает

Посмотрите на это с другой стороны. Предположим, датчик света попал с тени на солнце, и старые данные типа 10-11-12-11-10 вдруг стали 120-123-119... отсечка по 10% уже будет отсекать вообще все новое без шансов, зато медианный выскочит после энной итерации. Все зависит от того что у вас есть и что надо получить. Если ваши данные ретроспективные то можно все усреднять до ровной полки, если это данные, приходящие в реальном времени, то "вот пришел этот пичок" это новый тренд или мусорный выброс? Насколько "ровный" вы хотите? Насколько близко к "реальности" ? Можно просто усреднить, поставить медианные пред-фильтры, и этого будет достаточно для задачи, или действительно углубится в матан+физику и разработать физ. модель для того-же фильтра Калмана.

[Алексей Еремеев]

А еще. Кроме физики и математики еще есть физиология. Вот чем хороши стрелочные приборы? Стрелка колеблется, дрожит но глазами мы все равно "усредняем" и получаем данные с нужной нам точностью, а последствия "цифровизации" таковы что на дисплее параметр скачет, а если на границе старших разрядов то как вам болтанка 999.9 - 1000.1 - на этож смотреть невозможно, не то что воспринимать, и усредняй-не усредняй а болтать будет, меньше болтанки больше задержка в актуальном значении. Хочется и побыстрее и поточнее, только никому это не надо, а надо чтоб "стрелка в зеленом секторе - все нормально"

[We Are]

У стрелок есть большой минус - их надо видеть.

А цифры - их давно в коридор значений загоняют и оценивают "в норме" - "не в норме". Еще и меры принимаются автоматически, с учетом тенденций и скорости их нарастания.

Графики просто для контроля адекватности намерянных цифр физической сущности явления.