Простейший цифровой ФНЧ, для конечного массива точек Опции

[Ivan Kuznetzov]

Есть массив, в котором записаны значения с АЦП. Напряжения формируют некую "огибающую"

Платформа: STM32, сигнал - массив значений с АЦП в вольтах.

1)Подскажите фильтр (алгоритм на Си) чтобы сгладить(усреднить) значения в массиве?
2) Как определить характерные места у огибающей (резкое увеличение значения, резкий спад)



Сообщение отредактировал Ivan Kuznetzov - Jan 20 2011, 11:26

Эскизы прикрепленных изображений

 

[Demeny]

Есть массив, в котором записаны значения с АЦП. Напряжения формируют некую "огибающую"

1)Подскажите фильтр чтобы сгладить(усреднить) значения в массиве?

Для каждого элемента массива посчитайте среднее значение от окружающих его элементов.
Это и есть простейший цифровой ФНЧ. Чем шире окрестность, в которой считается среднее - тем сильнее "сглаживание", но и искажение формы исходного сигнала также сильнее.
Если не вдаваться в теорию - вот как-то так.

[Ivan Kuznetzov]

Demeny, спасибо! функция уже заметно красивей стала!!!

Код

for (x=1;x<ARRAYSIZE;x++) Pulsearray[x] = (Pulsearray[x-1] + Pulsearray[x+1])/2;

Сообщение отредактировал Ivan Kuznetzov - Jan 20 2011, 12:02

Эскизы прикрепленных изображений

 

[AHTOXA]

Код

for (x=1;x<ARRAYSIZE;x++) Pulsearray[x] = (Pulsearray[x-1] + Pulsearray[x+1])/2;

Наверное лучше

Код

Pulsearray[x] = (Pulsearray[x-1] + Pulsearray[x] + Pulsearray[x+1])/3;

[Ivan Kuznetzov]

AHTOXA, спасибо, действительно красивее!

Вот еще не пойму, как на синем графике найти указанные стрелками точки? дельту измерять какую-нибудь?

Эскизы прикрепленных изображений

 

[firstvald]

Есть метод скользящего среднего, считается так:

(текущий обработанный результат)=( (предыдущий обработанный результат)*(N-1)/N )+(текущий код АЦП)/N

Под кодом АЦП может выступать текущий отсчет необработанной вашей реализации. Собственно точки тогда имеют железную привязку

N=1....и до бесконечности

на практике достаточно N до 32, ну и естественно с ростом N теряются подробности сигнала.


Так характерные места - величины первой производной. Смотрите постоянно разницу предыдущего отсчета (возможно даже с некоторой глубиной, т е не последний, а предпоследний) и текущего. Ну стравнивайте полученную дельту с порогом или за знаком следите. Возможно понадобится и саму производную усреднять, если у вас по производной какие то критичные вычисления делаются.

[AHTOXA]

AHTOXA, спасибо, действительно красивее!

Да, но так всё равно неправильно Получается, что для первой точки мы взяли начальные данные, а для остальных - одна из точек (которая x-1) - уже отфильтрованная. Надо либо занычивать её в переменной, либо фильтровать в другой массив.
А лучше сделать скользящее среднее, как написал firstvald. Там можно регулировать степень сглаживания.

[SasaVitebsk]

Код

/************************************************************************
*                                                                        *
*                   Библиотека для цифровой фильтрации.                    *
*                            Версия:     1.01.                                *
*                                                                        *
* Файл: fir.c        Дата создания: 08.07.2010г.                            *
*              Последние изменения: 08.07.2010г.                            *
* Сапего Александр Леонидович. (sapegoal@mail.ru)                        *
*                                                                        *
************************************************************************/

#include "stdint.h"

// Фильтр 2-го порядка исходя из формулы y0^ = a0*(x0+x2)+a1*x1 + b1*y0 + b2*y1 + c1*y0
// где a и b - коэффициенты, c - целая часть y0
// После выполнения производится сдвиг y.
// y1 -> y2, y0 -> y1, y0^ -> y0;
//
// В процедуру передаётся 2 параметра
// 1 - указатель на значения, в последовательности: x0,x1,x2,y0,y1,y2
// 2 - указатель на значения коэффициентов в последовательности: a0,a1,b1,b2,c1
// По финишу результат делится на 32768

struct    data_s
{
  int16_t        x[3],y[3];                                // Данные фильтра
};

struct    coef_s
{
  int16_t        a0, a1, b1, b2, c1;                        // Коэффициенты фильтра
};

void fir2_16(uint8_t * data, uint8_t * coef)
{
  struct data_s    *d_s;
  struct coef_s    *c_s;
  int32_t        acc;
  
  d_s = (struct data_s *) data;
  c_s = (struct coef_s *) coef;
  acc = (int32_t)(d_s->x[0] + d_s->x[2]) * (int32_t)c_s->a0;
  acc += (int32_t)d_s->x[1] * (int32_t)c_s->a1;
  acc += (int32_t)d_s->y[0] * (int32_t)c_s->b1;
  acc += (int32_t)d_s->y[1] * (int32_t)c_s->b2;
  acc >>= 15;
  acc += d_s->y[0] * c_s->c1;
  d_s->y[2] = d_s->y[1];
  d_s->y[1] = d_s->y[0];
  d_s->y[0] = acc;
}

Код

// Фильтр Баттерворта 4 порядка 1000 -> 100; 200  = -30 дб
//
// 1 звено y0 = 0,061885*(x0+x2) + 0,123770*x1 + 1,048600*y1 - 0,296140*y2
// Коэффициенты a0,a1,b1,b2,c1 = 2028, 4056, 1593, -9704, 1
// 2 звено y0 = 0,077956*(x0+x2) + 0,155913*x1 + 1,320910*y1 - 0,632739*y2
// Коэффициенты a0,a1,b1,b2,c1 = 2554, 5109, 10516, -20734, 1
//
int16_t            fltr_c1[5] = {2028, 4056, 1593, -9704, 1},
                fltr_c2[5] = {2554, 5109, 10516, -20734, 1};

Код

    fir2_16((uint8_t *)&lin[cnt_line].x[0],(uint8_t *)fltr_c1);// Первое звено фильтра
    fir2_16((uint8_t *)&lin[cnt_line].y[0],(uint8_t *)fltr_c2);// Второе звено фильтра

Так например. ))

[Ivan Kuznetzov]

SasaVitebsk, поясните пожалуйста как пользоваться? Есть массив, 200 точек. Входные параметры фильтра - это data?

[demiurg_spb]

Пробегало давненько:

Код

u16_t Yavg;

for(;;)
{
   Yavg -= Yavg/256;
   Yavg += ADCH;
}

Запуск и готовность АЦП за Вами.
Утверждается, что в старшем байте Yavg получим фильтр НЧ с частотой среза = частоте запуска АЦП/256
с погрешностью обработки <0.5 LSB, по АФЧХ эквивалентный RC цепочке с такой же частотой среза.

Почему? Хемминг уже рассказал, у меня лучше не получится.

[777777]

Платформа: STM32, сигнал - массив значений с АЦП в вольтах.

1)Подскажите фильтр (алгоритм на Си) чтобы сгладить(усреднить) значения в массиве?

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

[_pv]

1)Подскажите фильтр (алгоритм на Си) чтобы сгладить(усреднить) значения в массиве?
2) Как определить характерные места у огибающей (резкое увеличение значения, резкий спад)

1) demiurg_spb и firstvald уже подсказали простейший IIR фильтр:

Код

y[0] = x[0];
for (int i = 1; i < len; i++) y[i] = y[i-1] + (x[i] - y[i-1]) * K

K = 0.0 .. 1.0, временная постоянная фильтра. Ну и конечно можно (и нужно) заменить на целочисленное умножение и битовый сдвиг вправо.

Хотя если надо потом точно определять положение "характерных мест", то fir фильтр (скользящее среднее - частный случай) лучше из-за линейной фазовой характеристики. Кстати, то что предложил SasaVitebsk - fir фильтром не является.

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

разводить плату надо было так, чтобы программисту не пришлось фильтровать результаты АЦП

таким сигнал вполне может быть изначально, и совсем не обязательно вызван шумами на плате.

[demiurg_spb]

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

Ага и всем-всем заказчикам отрубать руки если они датчик (TC или RTD например) подключат проводом без экрана и не дай Бог вблизи силовых кабелей или или или.......... :-)
В любом нормальном промышленном измерительном оборудовании имеется возможность регулировки степени фильтрации сигнала так-сказать по месту применения. Так-что схемотехника+трассировка - отдельно, грамотный софт - отдельно и перекладывать с больной головы на здоровую не стоит.

[GetSmart]

В Саратове живут суровые профэссианалы

[777777]

В Саратове живут суровые профэссианалы

Увы - не все. Почему-то считается немерянной крутизной развести плату минимального размера, сэкономив на конденсаторах по питанию, землях и размещении элементов - главное упихать их поплотнее, а то, что АЦП шумит на 4 разряда - это фигня, программисты отфильтруют.

[Ivan Kuznetzov]

Внесу ясность в суть дела. Речь идет о цифровом измерителе артериального давления.
И делается это так - момент появления тонов в манжете - это систолическое давление, момент когда тоны вдруг резко спадают и начиают идти примерно одинакого уровня - диастолическое давление. Проблема в том что интенсивность тонов у всех людей разная и "порог" тут не катит по-видимому...

[SasaVitebsk]

SasaVitebsk, поясните пожалуйста как пользоваться? Есть массив, 200 точек. Входные параметры фильтра - это data?

Там же всё видно. ))
Давайте я упрощу.

1. Объявляем

Код

// Каналы АЦП
struct    data_ADC
{
  int16_t        x[3],y[3],z[3];                            // Данные фильтра
};

struct data_ADC lin[8]; // Объявляем массив данных.

2. Объявляем массив констант коэффициентов фильтра так, как это в предыдущем моём посте показано.

3. Читаем АЦП

Код

    lin[cnt_line].x[2] = lin[cnt_line].x[1];            // Коэффициенты x переписать
    lin[cnt_line].x[1] = lin[cnt_line].x[0];            // Коэффициенты x переписать
    lin[cnt_line].x[0] = ADC;                            // Прочитать АЦП
    fir2_16((uint8_t *)&lin[cnt_line].x[0],(uint8_t *)fltr_c1);// Первое звено фильтра
    fir2_16((uint8_t *)&lin[cnt_line].y[0],(uint8_t *)fltr_c2);// Второе звено фильтра

4. Результат получаем в lin[cnt_line].z[0]


Если у вас фильтр 6 порядка, то вам надо добавить в структуру data_ADC ещё один массив, и вызвать ещё одно звено фильтра с соответствующими коэффициентами.

PS: Есть и IIR фильтр если хотите. ))

[Ivan Kuznetzov]

SasaVitebsk, хотим! чем отличается от этого?

Сообщение отредактировал Ivan Kuznetzov - Jan 25 2011, 19:32

[SasaVitebsk]

Кстати, то что предложил SasaVitebsk - fir фильтром не является.

)) Во-первых Вы мне льстите. Это предложил не я, а Баттерворт. ))
Во-вторых, а чем по вашему мнению является уравнение вида "y0 = 0,061885*(x0+x2) + 0,123770*x1 + 1,048600*y1 - 0,296140*y2"?
(где x1 = x(n-1) и так далее)
Я не математик, если честно, и не теоретик, но некоторые источники типа фильтрсолюшн, QED и другие, просто вводят нас в заблуждение, а Вы нам сейчас откроете всю правду. ))
Ждём.
===
2 777777 по поводу

разводить плату надо было так, чтобы программисту не пришлось фильтровать результаты АЦП

. А причём здесь АЦП? Иногда входной сигнал требует фильтрации. Или обработанный. То есть даже если этот сигнал идеально снят АЦП.

[firstvald]

Внесу ясность в суть дела. Речь идет о цифровом измерителе артериального давления.
И делается это так - момент появления тонов в манжете - это систолическое давление, момент когда тоны вдруг резко спадают и начиают идти примерно одинакого уровня - диастолическое давление. Проблема в том что интенсивность тонов у всех людей разная и "порог" тут не катит по-видимому...

Тут все сурово. Шумов много. Один из основных - шуршание манжеты ( и с увеличением возраста манжеты растет). Причем спектрально он практически в полосе тонов короткова. Немножко спасают специальные микрофоны с мелой чувствительностью к низкочастотным шумам. Я думаю вам нужно как фнч, так и фвч. А решение о начале и конце тонов принимать только после того как всю раелизацию записали - т е пороги принятия решения плавающие.

[_pv]

)) Во-первых Вы мне льстите. Это предложил не я, а Баттерворт. ))
Во-вторых, а чем по вашему мнению является уравнение вида "y0 = 0,061885*(x0+x2) + 0,123770*x1 + 1,048600*y1 - 0,296140*y2"?
(где x1 = x(n-1) и так далее)
Я не математик, если честно, и не теоретик, но некоторые источники типа фильтрсолюшн, QED и другие, просто вводят нас в заблуждение, а Вы нам сейчас откроете всю правду. ))
Ждём.

я всегда считал что, соответсвенно названию, FIR фильтры (с конечной импульсной характеристикой) это Y = A1*X1 + A2 * X2 + ..., одна из особенностей которых - линейная фазовая характеристика, т.е. просто задержка, а это может быть важно для автора исходного сообщения, он хочет определять после фильтрации когда именно спад на сигнале произошол.
а IIR - с обратной связью и, соответственно, с бесконечной импульсной характеристикой это Y = A1*X1 + A2*X2 + ... + B1*Y1 + B2*Y2 + ...
и, соответственно, с нелинейной фазовой характеристикой
не думаю что qed и фильтрсолюшен могут их перепутать, посмотрите внимательнее.

[firstvald]

Навсякий случай: FIR все же постабильнее, хотя и длиннее, IIR может так зазвенеть, что не задушишь не убъешь. Если честно, топику надо простейшими методами добиться удовлетворительной работы, иначе утопнет в этих фильтрах.

[SasaVitebsk]

не думаю что qed и фильтрсолюшен могут их перепутать, посмотрите внимательнее.

Я вот даже не буду с Вами спорить. Просто посмотрите сами любую прогу из перечисленных. Лучше фильтрсолюшен - там прямо формулу можно посмотреть. И задайте любой фильтр 2 порядка, к примеру. Фильтр ведь можно рассматривать как последовательное соединение звеньев.

Навсякий случай: FIR все же постабильнее, хотя и длиннее, IIR может так зазвенеть, что не задушишь не убъешь. Если честно, топику надо простейшими методами добиться удовлетворительной работы, иначе утопнет в этих фильтрах.

Если не на обум фильтр выбирать, то ничего там звенеть не будет. И не надо пугать человека. Пусть человек возьмёт и выведет результирующий график с моего примера. Хотя, этот фильтр не под него делался. Характеристики его я привёл. А если захочет, то может и свой фильтр построить с любыми характеристиками. Работы на подбор фильтра - 15 минут.

[ig_z]

Кстати, то что предложил SasaVitebsk - fir фильтром не является.

Во-вторых, а чем по вашему мнению является уравнение вида "y0 = 0,061885*(x0+x2) + 0,123770*x1 + 1,048600*y1 - 0,296140*y2"?
(где x1 = x(n-1) и так далее)
Я не математик, если честно, и не теоретик, но некоторые источники типа фильтрсолюшн, QED и другие, просто вводят нас в заблуждение, а Вы нам сейчас откроете всю правду. ))
Ждём.

Так все таки, каким типом является этот фильтр? Если по вашему мнению утверждение "fir фильтром не является" есть ложное, то это все же fir фильтр?
Я тоже не математик, и не теоретик, как и многие здесь. Поэтому вы заинтриговали.
Ждем. И если не затруднит, на пальцах, без ссылок на проги. Это ведь вопрос принципа, а не реализации

[ViKo]

Так все таки, каким типом является этот фильтр? И если не затруднит, на пальцах, без ссылок на проги. Это ведь вопрос принципа, а не реализации

Из приведенных формул не совсем понятно, какие индексы подразумеваются для отсчетов y. Но если в выражении используются не только входные данные (x), но и выходные (y), рассчитанные чуть раньше, то это будет фильтр с бесконечной импульсной характеристикой (БИХ, по-английски IIR). Если бы использовались только входные данные, это был бы фильтр с конечной импульсной характеристикой (КИХ, по-английски FIR).

[SasaVitebsk]

Извиняюсь. Действительно IIR.

[Ivan Kuznetzov]

Еще вопрос по теме тонометра. Я сейчас использую датчик давление-напряжение, пульсации давления вылавливаю так: убираю постоянную составляющую и усиливаю пульсации операционником.
Есть вариант использовать датчик давление-частота (1 МГц). Как улавливать малейшие колебания частоты (тоны короткова) при частотном методе?

[=AK=]

1)Подскажите фильтр (алгоритм на Си) чтобы сгладить(усреднить) значения в массиве?

Можно сделать доморощенное грубое подобие фильтра Калмана. Его достоинство - в простоте реализации: не нужно запоминать массив данных, фильтр Калмана фильтрует данные "на лету", по мере их поступления.

Пусть Vin - результат на выходе АЦП, а Vout - отфильтрованный текущий результат, n - некое число, при помощи которого оценивается "достоверность" текущего результата. Вычисление ведется по простой формуле:

Vout = Vout - (Vout/n) + (Vin/n), где n по ходу фильтрации увеличивается от 1 до некоторого разумного (т.е. не очень большого) значения, определяемого ожидаемыми свойствами сигнала.

Пример:

1. Самое первое измерение принимаем как данность, поскольку сравнивать и усреднять не с чем, n=1:
Vout = Vin

2. У второго измерения "достоверность" такая же, какая была у первого, так что принимаем n=2:
Vout = Vout - (Vout/2) + (Vin/2)

3. Для третьего измерения принимаем n=3:
Vout = Vout - (Vout/3) + (Vin/3)

...

16. Ограничимся, например, n=16, иначе все вообще усреднится нафиг. Это и все последующие измерения считаем так:
Vout = Vout - (Vout/16) + (Vin/16)


Если говорить упрощенно, то у настоящего фильтра Калмана n меняется динамически в зависимости от сигнала. Когда сигнал меняется быстро, вес новых измерений увеличивается (т.е. n уменьшается), когда же изменения сигнала невелики, то вес каждого нового измерения уменьшается (n увеличивается).

[firstvald]

Не , там напряжение ли или частота - все едино - переводим в давление. Если частота, можно период мерять и тогда получим реализацию с хорошим временным разрешением. Так тоны короткова ловить - в этом самый смысл. Не зря покупные тонометры все врут по-своему , алгоритмы разные. Ну сделайте несколько огроменных массивов. В одном чистая реализация, в другом отфильтрованная , в третьем первая производная, в четвертом отфильтрованная производная. Если на них посмотреть - какой алгоритмик и придумаете. Я уже сказал - порог принятия решения должен быть переменным. Из поступающей реализации надо вычитать медленную устредненную - тогда получаем собственно изменения без постоянной составляющей. Но тут все лепять по своему, вряд ли кто будет так запросто делиться. Хотя лет 30 назад выходила книжка по кардиомониторам - мужики там все алгоритмы расписали, включая фильтрацию. Где теперь те мужики.

[ViKo]

http://search.analog.com/search/default.as...CG&local=en

[Guest_orthodox_*]

Из приведенных формул не совсем понятно, какие индексы подразумеваются для отсчетов y. Но если в выражении используются не только входные данные (x), но и выходные (y), рассчитанные чуть раньше, то это будет фильтр с бесконечной импульсной характеристикой (БИХ, по-английски IIR). Если бы использовались только входные данные, это был бы фильтр с конечной импульсной характеристикой (КИХ, по-английски FIR).

Несколько FIR звеньев, включенных последовательно , по этому определению попали бы в IIR?

[ViKo]

Несколько FIR звеньев, включенных последовательно , по этому определению попали бы в IIR?

С чего это вдруг? Каждое звено использует только то, что приходит к ним на вход. И не использует то, что сами рассчитали.
А вот если охватить все эти звенья обратной связью, подать выходные данные на вход - получим IIR.
Годится?