Jump to content

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


We Are
 Share

Recommended Posts

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

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

967388041_Screenshotfrom2024-02-1115-44-55.png.707924afb66a180f341d3be7b837bca6.png

 

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

Один из простейших - усреднение, а точнее метод "скользящего окна":
Создается массив длиной, скажем, 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;
}

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

1663432508_Screenshotfrom2024-02-1115-36-05.png.af70b484c26f7e4f3dc7a1f178de7e8b.png

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

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

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

Link to comment
Share on other sites

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

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

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

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

Link to comment
Share on other sites

В 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);
}

 

Link to comment
Share on other sites

Новый аккумулятор EVE серии PLM для GSM-трекеров, работающих в жёстких условиях (до -40°С)

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

Эти аккумуляторы поддерживают заряд при температуре от -40/-20°С (сниженным значением тока), безопасны (не воспламеняются и не взрываются) при механическом повреждении (протыкание и сдавливание), устойчивы к вибрации. Они могут применяться как для автотранспорта (трекеры, маячки, сигнализация), так и для промышленных устройств мониторинга, IoT-устройств. Подробнее параметры и результаты тестов новой серии PLM по ссылке.

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

13 hours ago, donec said:

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

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

Edited by We Are
Link to comment
Share on other sites

Как измерить внутреннее сопротивление литиевого аккумулятора. Измеряем правильно

Внутреннее сопротивление – один из наиболее значимых параметров, от которого напрямую зависят другие характеристики литиевого аккумулятора. От этого параметра зависят такие характеристики источника питания, как напряжение на выходе под нагрузкой, максимальный выходной ток и коэффициент полезного действия (КПД). Рассмотрим как измерить действительное значение внутреннего сопротивления АКБ, используя определенные методики. Подробнее>>

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

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

 

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

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

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

 

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

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

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

 

 

Link to comment
Share on other sites

4 minutes ago, BARS_ said:

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

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

5 minutes ago, BARS_ said:

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

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

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

7 minutes ago, BARS_ said:

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

Истинно так.

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

Link to comment
Share on other sites

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

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

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

 

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

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

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

Link to comment
Share on other sites

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

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

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

Link to comment
Share on other sites

В 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;
}

Смотрите

Link to comment
Share on other sites

Что, ну и..

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

Link to comment
Share on other sites

Just now, donec said:

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

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

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

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

Link to comment
Share on other sites

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

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

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

Link to comment
Share on other sites

5 minutes ago, donec said:

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

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

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

Edited by We Are
Link to comment
Share on other sites

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

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

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

где

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

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

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

Link to comment
Share on other sites

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

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

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

Link to comment
Share on other sites

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

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

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

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

Edited by We Are
Link to comment
Share on other sites

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

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

#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;
}

Edited by Алексей Еремеев
Link to comment
Share on other sites

17 minutes ago, IMXO said:

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

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

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

1068142689_Screenshotfrom2024-02-1115-36-05.png.4e44b7dbefa4f6b505f049e19a3ebaf5.png

Link to comment
Share on other sites

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

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

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

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

Link to comment
Share on other sites

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

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

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

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

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

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

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

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

 

Link to comment
Share on other sites

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

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

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

Link to comment
Share on other sites

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

Link to comment
Share on other sites

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

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

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

Link to comment
Share on other sites

Join the conversation

You can post now and register later. If you have an account, sign in now to post with your account.
Note: Your post will require moderator approval before it will be visible.

Guest
Unfortunately, your content contains terms that we do not allow. Please edit your content to remove the highlighted words below.
Reply to this topic...

×   Pasted as rich text.   Restore formatting

  Only 75 emoji are allowed.

×   Your link has been automatically embedded.   Display as a link instead

×   Your previous content has been restored.   Clear editor

×   You cannot paste images directly. Upload or insert images from URL.

Loading...
 Share

  • Recently Browsing   0 members

    • No registered users viewing this page.
  • Сообщения

    • А если " WiFi реле " придет не рабочее и на него плю́нешь.    Не совсем удобно, нужно заходить в приложение для управления реле, манипуляций будет больше.
    • Чявой? За схемотехнику благодарю) 
    • Микроконтроллеры и микросхемы стандартной логики тоже? Возможно для нагрева воды в системе отопления, а требования поддержания заданной температуры. 
    • Так если на контакты поступает напряжение, то можно на эти контакты поставить реле, а контакты реле уже использовать для замыкания, там как раз и не будет напряжения.
    • Про две группы написали, а про третью. Похоже вы сами не совсем понимаете режим работы вашего котла. Поставьте три отдельных регулятора и настройте их на каждую температуру.  А температура у вас падает с увеличением расхода воды. Только это будет не дешевле чем покупать например ПЛК-107 (Овен). А в ПЛК можно как раз использовать несколько реле и контролировать по расходу или по давлению. Если это производство или коммерция, то самоделки запрещено использовать без сертификации. Откуда такая экономия. https://termoshkaf.com/product/uko/ https://phantom-stab.ru/catalog/termoregulyatory/elektricheskie/termoregulyator-cifrovoy-TK-125
    • Может и не дурак, но сделал довольно топорно. Или выхода иного не было: По теме: двухканальные термостаты достаточно доступны по ценам. Можно добавить ещё один одноканальный, либо хорошенько подумать, сформулировать внятное ТЗ и решить - а нужно ли три канала? Может, и двух хватит. Допустим, "холодный" котёл включается на полную мощность и быстро набирает температуру 50 градусов. При достижении её часть ТЭНов отключается совсем, а оставшаяся работает на "подогрев" до 60, снижая нагрузку на электросеть. При критическом падении температуры (большой расход), снова включаем на полную мощность и т.д.  
×
×
  • Create New...