Код фильтра и теста

Вначале представим то, что можно быстро скопировать и вставить. Структура и функция расчёта одной итерации фильтра:

typedef struct FirFrac15 {
  int16_t x, x_1, y, y_1;
} FirFrac15;

#define N 3 //задание полосы пропускания (целое число)
void FirFrac15Calc(FirFrac15 *Filter) {
  register int32_t Acc;
  register int16_t xAcc = Filter->x + Filter->x_1;
  Acc = (((int32_t)Filter->y_1 << (N + 1)) - ((int32_t)Filter->y_1 << 1) + (int32_t)xAcc) >> (N + 1);
  Filter->y = (int16_t)Acc;
  Filter->x_1 = Filter->x; // входной отсчёт и запаздывающий на один
  Filter->y_1 = Filter->y; // выходной отсчёт и запаздывающий на один
}

Скетч для тестирования:

#define Size 100 //размер тестового массива

int16_t Vhod[Size];
int16_t Vihid[Size];

typedef struct FirFrac15 {
  int16_t x, x_1, y, y_1;
} FirFrac15;

#define N 2 //задание полосы пропускания (целое число)
void FirFrac15Calc(FirFrac15 *Filter) {
  register int32_t Acc;
  register int16_t xAcc = Filter->x + Filter->x_1;
  Acc = (((int32_t)Filter->y_1 << (N + 1)) - ((int32_t)Filter->y_1 << 1) + (int32_t)xAcc) >> (N + 1);
  Filter->y = (int16_t)Acc;
  Filter->x_1 = Filter->x; // входной отсчёт и запаздывающий на один
  Filter->y_1 = Filter->y; // выходной отсчёт и запаздывающий на один
}
void FirFrac15Reset(FirFrac15 *Filter) {
  Filter->x = 0, Filter->x_1 = 0, Filter->y = 0, Filter->y_1 = 0;
}

#define Float2IQ15(f) (int16_t(float(f) * (1 << 15))) //переход от плавающей к фиксированной
#define IQ15ToFloat(i) (float(i) / float(1 << 15)) //переход от фиксированной к плавающей

void setup() {

  Serial.begin(9600);
  Serial.setTimeout(5);
  //Инициализация массива
  for (int i = 0; i < (sizeof(Vhod) / sizeof(int16_t)); i++) {
    //Vhod[i] = Float2IQ15(0.499); //проверка статического
    Vhod[i] = Float2IQ15(0.25) + random(-Float2IQ15(0.1),Float2IQ15(0.1)); // проверка сигналом
  }
}

uint8_t ctr = 0;
FirFrac15 Filter = { 0, 0, 0, 0 };  //фильтр с нулевыми нач условиями
void loop() {
  if (ctr < 50) { // 50 точек чтобы не было автообновления плоттера
    int16_t x = Vhod[ctr];
    Filter.x = x;
    FirFrac15Calc(&Filter);
    int16_t y = Filter.y;
    Serial.print(x);
    Serial.print(",");
    Serial.println(y);
    ctr++;
  } else {
    ctr = 0;  //повтор
    FirFrac15Reset(&Filter);
    delay(10000);
  }
}

Для просмотра выходных данных открыть плоттер в среде "Ардуино" создания скетчей.

Для чего необходимо

Иногда бываеттак, что необходимо что‑то сделать с входным сигналом АЦП, чтобы он не содержал всплесков, отфильтровать из него постоянную составляющую и многие другие преобразования. В основном на скорую руку используется простейший КИХ‑фильтр как среднее значение от суммы отсчётов. Недостатком является сложность получения большой постоянной времени. Предлагаемое решение основано на реализации БИХ‑фильтра первого порядка, при этом, вместо операции умножения, обычно реализуемой как MUL или в виде отдельной функции для контроллеров, не содержащих аппаратной реализации (или использующих внешний модуль умножения через регистры в памяти). Заместо этой «дорогостоящей» операции для простых контроллеров используется арифметический сдвиг вправо‑влево ASR/ASL, которые являются эффективными, также, предлагаемый метод может быть легко реализован на ПЛИС на уровне регистров. Сдвиг производится на заранее заданное количество разрядов, данный сдвиг фактически определяет частоту среза. Недостатком является фиксированное значение частоты среза фильтра. Однако, комбинируя данные фильтры можно реализовать некоторое приближение к желаемой характеристики.

Как это работает

Чтобы посмотреть расчёт фильтра, скопируйте и вставьте этот код в систему символьных вычислений Maxima (желательно по исполняемым группам после каждого комментария построчно) и нажмите Ctrl+Shift+R - пересчитать скрипт

	
	kill(all)$ /*очистить все выражения*/
/*Фильтр первого порядка в общем виде*/
	Filter:sum(b[i]*z^(-i),i,0,1)/(1+sum(a[i]*z^(-i),i,1,1));
factor(Filter); /*"упрощение"*/
/*подстановка для перехода в комплексно-частотную характеристику, fd - частота дискретизации*/
	ComEx:z=exp(%i*2*%pi*f/fd);
/*КЧХ*/
	Fc:Filter,ComEx;;
/*на нулевой частоте значение КЧХ равно*/
	Fdc:Fc,f=0;
/*на частоте Найквиста (половина частоты дискретизации) значение КЧХ равно*/
	FNf:Fc,f=fd/2;
/*Критерий ФНЧ - на нулевой частоте частотная характеристика равна 1, на Найквиста - ную*/
	Crt:[Fdc=1,FNf=0];
/*Находение коэффициентов, соответствующих данному критерию*/
	Snc:solve(Crt,[b[0],b[1]])[1];
/*Нахождение полюса. Коэффициент a[1] - фактически полоса пропускания, должен быть от -1 до 1 внутри единичной окружности*/
	solve(denom(factor(Filter)),z);
/*Принять значения коэффициента как степени двойки*/
	Cdv:[a[1]=2^(-N)-1];
/*Произвести подстановку коэффициента*/
	Sab:expand(append(Cdv,subst(Cdv,Snc)));
/*Рекуррентное выражение для фильтра*/
	RV:sum(x[i]*b[i],i,0,1)-sum(y[i]*a[i],i,1,1);
/*Неоптимизированная форма - прямая подстановка*/
	RVab:RV,Sab;
declare(N,integer);/*N - целое число*/
/*Варианты формы рекуррентного соотношения*/
	factor(RVab);expand(RVab);facsum(RVab,[x[0],x[1]]);factorfacsum(RVab,2^N);
/*Раскрыть скобки и подобрать подобные слагаемые*/
	C1:collectterms(expand(RVab),x[0],x[1],y[1]);
/*Выбор оптимальной формы - внутри скобок умножение, за скобками - деление*/
	collectterms(expand(factor(C1*2^(N+1))),x[0],x[1],y[1])/2^(N+1);factor(C1-%);
/*сформированные выражения*/;
  RVmx:Sx:(x[1]+x[0]);
  RVmy:2^(-N-1)*(y[1]*2^(N+1) - 2*y[1] + Sx);
  factor(RVab-RVmy);/*проверка*/
/*расчёт цифрового фильтра*/
	Npt:20$LY:makelist(makelist (0,i,1,Npt),j,1,4)$
	for m:1 thru 4 do  
	(
	  /**/
	  Cfn:[N=m], [yn,yn1,xn,xn1]:[0,0,0,0], /*Задать начальные условия*/
	  for n: 0 thru Npt -1 do  
	  (
	        tn:n,  xn:1*signum(n), /*реакция на ступеньку*/
	        x01:xn+xn1, /*xn - текущий отсчёт*/
	        yn:2^(-N-1)*(yn1*2^(N+1)-yn1*2+x01), /*формирование ответа - сравнить с коэффициентами*/
	        yn:float(subst(Cfn,yn)),
	        xn1:xn,
	        yn1:yn,
	        LY[m][n+1]:[tn,yn]
	  )
	)$
LY[1]; /*выходные точки*/
/*Переходной процесс */
	wxplot2d([[discrete, LY[1]],[discrete, LY[2]],[discrete, LY[3]],[discrete, LY[4]]],[t,0,20],[style, points, points,points,points],[legend,false]),wxplot_size=[480,320];;
assume(N>0,f>0,fd>0);
/*При заданных коэффициентах находится КЧХ, зависящая только от N*/
	Fcab:Fc,Sab;
/*Упрощение выражений для АЧХ и ФЧХ*/
	FcabA:factor(trigsimp(cabs(Fcab)));FcabP:factor(trigsimp(carg(Fcab)));
/*Построение семейства частотных характеристик*/
	Num:fd=10000;NumA:makelist(N=i,i,1,6);
/*АЧХ*/
	FAn:create_list (subst([i,Num],FcabA),i,NumA);
/*ФЧХ*/
	FPn:create_list (subst([i,Num],FcabP),i,NumA)$
    wxplot2d(FAn,[f,0,5000], [legend,false]),wxplot_size=[480,320];wxplot2d(FPn,[f,0,5000],[legend,false]),wxplot_size=[480,320];

Основные этапы:

• Функция фильтра в общем виде как коэффициенты при z^-1

• Подстановка комплексной экспоненты z=exp(%i2%pi*f/fd) %i - мнимая единица равная
  (-1)^(2^(-1)), %pi - длина замкнутой дуги единичной окружности

• Нахождение граничных точек на нулевой частоте (постоянный сигнал) и частоте Найквиста (точка бесконечной частоты для аналогичных непрерывных функций)

• Решение для коэффициентов относительно граничных точек

• Дискретизация коэффициента путём квантования его 2^N, при этом, полюса передаточной функции при изменении N лежат внутри единичной окружности (для первого порядка один действительный полюс)

• Оптимизация рекуррентного соотношения с выносом за скобки 2^M как умножения и
  2^-M как деления

Что получается в результате

При запуске скетча монитор порта должен вывести следующие графики

Семейство переходных процессов для различных N:

Семейство АЧХ для частоты дискретизации 10 кГц

Следует также отметить, что чем больше постоянная времени фильтра или его порядок, тем больше должна быть разрядность аккумулятора. Так, в DSP-процессорах имеется расширение аккумулятора, например с 32 до 40 и более бит, которые при операции REP MAC (повторить умножение с накоплением заданное количество раз) позволяют произвести суммирование небольших величин, расширяя динамический диапазон. Таким образом, если имеется ПИ-регулятор с постоянной времени часы и сутки и частотой квантования в килогерцы, то он должен иметь внутреннюю разрядность 32 бита несмотря, например, на 10-ти битный вход от АЦП. Пример эмуляции разрядности при заданной постоянной времени представлен на рисунке ниже.

• Эмуляциюможно осуществить с использованием функции floor — отсечение дробной части вещественного числа

• Общий вид функции эмуляции квантования по уровню как floor(x*2^N)/2^N, N — целое, имитирующее «разрядность процессора»

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