Построение цифрового фильтра с конечной импульсной характеристикой

Вступление издалека

Недавно передо мной встала достаточно интересная задача, с которой я раньше никогда не сталкивался — борьба с шумом. Мы принимали сигнал с датчиков на аналогово-цифровой преобразователь (АЦП)
А так как данная тема для меня была (хотя и сейчас есть кое-где) темным лесом, я пошел мучить вопросами гугл, мне показалось освещена эта тема не очень подробно и доступно, поэтому решил написать статью с примером разработки и готовым исходником.

Ближе к делу

Цифровые фильтры могут быть двух видов – с конечной и с бесконечной импульсной характеристикой (КИХ и БИХ). Для решения моей задачи подходит КИХ-фильтр, поэтому про него и расскажу.

Для начала посмотрим как же он работает:



Здесь показан пример фильтра нижних частот, как видно на рисунке, этот фильтр пропускает нижние частоты, а все остальные старается отсечь (подавление), или хотя бы ослабить (переход). Отклонения в полосе пропускания и полосе подавления выбираются в зависимости от принимаемого сигнала, но при использовании различных весовых функций, на них могут накладываться определенные ограничения. Например, если используется весовая функция Хэмминга, то эти отклонения будут равны между собой.
Ширина полосы перехода ∆F зависит от длины фильтра и от весовой функции (для функции Блэкмена ∆F=5,5|N).

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



Она реализуется через цикл, но постойте, а где же взять нужные коэффициенты? Вот тут-то как раз и зарыта собака (и не одна).

Параметры фильтра

Естественно для разных фильтров нужны разные коэффициенты, и для этого нужно определиться с параметрами фильтра, это обычно сначала делается теоретически (с умным видом прикидываем какая у нашего сигнала частота, потом частоты, которые надо отсеивать), а потом изучаем АЧХ реальных измерений (и осознаем, как сильно мы ошибались).
По этим АЧХ мы определяемся с идеальной частотной характеристикой (какие частоты проходят свободно, какие мы убираем и как сильно), теперь нам нужна идеальная импульсная характеристика её можно посчитать как Фурье-образ от идеальной частотной:



где H_D(w) – идеальная характеристика.

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




где fc и wc – частота среза.

Итак, осталось уже немного идеал идеалом, а мы имеем дело с практикой, и нам нужна «реальная» импульсная характеристика. Для её расчета нам понадобится весовая функция w(n), их есть несколько разновидностей, в зависимости от требований к фильтру (Хэмминга, Хеннинга, Блэкмена, Кайзера, о них не говорю, ибо статья и так большая), в нашем случае я использую функцию Блэкмена:



где N – длина фильтра, т.е. количество коэффициентов.

Теперь надо перемножить идеальную импульсную характеристику и весовую функцию:



Финишная прямая

Теперь мы готовы рассчитать выходные значения, по формуле фильтра, она самая первая в этой статье, ну вот и всё, в завершение привожу исходный код фильтра:
void Filter (const double in[], double out[], int sizeIn)
{
const int N = 20; //Длина фильтра
long double Fd = 2000; //Частота дискретизации входных данных
long double Fs = 20; //Частота полосы пропускания
long double Fx = 50; //Частота полосы затухания

long double H [N] = {0}; //Импульсная характеристика фильтра
long double H_id [N] = {0}; //Идеальная импульсная характеристика
long double W [N] = {0}; //Весовая функция

//Расчет импульсной характеристики фильтра
double Fc = (Fs + Fx) / (2 * Fd);

for (int i=0;i<N;i++)
{
if (i==0) H_id[i] = 2*M_PI*Fc;
else H_id[i] = sinl(2*M_PI*Fc*i )/(M_PI*i);
// весовая функция Блекмена
W [i] = 0.42 - 0.5 * cosl((2*M_PI*i) /( N-1)) + 0.08 * cosl((4*M_PI*i) /( N-1));
H [i] = H_id[i] * W[i];
}

//Нормировка импульсной характеристики
double SUM=0;
for (int i=0; i<N; i++) SUM +=H[i];
for (int i=0; i<N; i++) H[i]/=SUM; //сумма коэффициентов равна 1 

//Фильтрация входных данных
for (int i=0; i<sizeIn; i++)
{
out[i]=0.;
for (int j=0; j<N-1; j++)// та самая формула фильтра
if(i-j>=0)
out[i]+= H[j]*in[i-j];
}
}


При подготовке статьи использовались:
Основные характеристики и параметры фильтров. analogiu.ru/6/6-5-2.html
Айфичер Э. Джервис Б. Цифровая обработка сигналов. Практический подход. 2-е издание
Поделиться публикацией
Комментарии 31
    0
    Это какой фильтр? ФНЧ, ФВЧ или полосовой?
      0
      Извините, увидел.
      +1
      Вы бы добавили, как фильтром убрать шум, для полноты картины.
        0
        И добавьте пожалуста что х — это свертка а не умножение.
          0
          а вот х это всё-таки умножение, самое простое и всем известное, или я чего-то не понимаю?
            +1
            Ой, да, извините. В прошлом году был курс по обработке сигналов — вылетело из головы.
              0
              Всё в порядке, ошибиться может каждый.
          0
          это я дополню, но боюсь только завтра) спасибо за совет
          +1
          А у вас конечный кадр данных? То есть вы кладете все в огромный массив, потом фильтрует? Я это к тому, что если прием данных происходит непрерывно и вы его искуственно нарезаете на части, то в местах стыков у вас будут проблемы…
            0
            Да кадр конечный, согласен, если сбор данных непрерывный этот вариант будет не очень приемлем.
              0
              А как тогда осуществлять фильтрацию, если идет тот же непрерывный поток данных с АЦП, а надо ослабить 50 Гц помеху? Сейчас делаю обычным усреднением, но подозреваю, что это крайне не эффективно.
                0
                Такого я ещё не делал, но могу посоветовать рекурсивную реализацию цифрового фильтра (есть в книге Айфичера (см. подвал статьи)), ещё видел реализацию непрерывной фильтрации через обратную связь (измеряется шум, далее из получаемого сигнала вычитается шум, и изменение шума подается снова на вход, ссылки под рукой нет, но кажется довольно быстро это можно найти в гугле), и есть шум довольно маленький по сравнению с сигналом — можно попробовать метод «скользящего среднего».
                  +1
                  При программной реализации в случае покадровой обработки самое простое — подавать данные в функцию фильтрации описанную выше «в нахлест». То есть к основному кадру данных пристегивать кусочек данных впереди и кусочек позади. Размер этого «кусочка» взять равным половине порядка фильтра. Естественно, что эффективно это будет только для случаев, когда размер кадра значительно больше размера фильтра, иначе большие накладные расходы.
                  Другой способ — поштучная обработка пришедших данных. То есть в рабочем массиве храните только число отсчетов равное порядку фильтра. На каждом новом сэмпле новый отсчет в массив добавляете, самый старый выкидываете и запускаете процедуру однократно (для одного сэмпла). Так чаще делается в DSP. Ну и подобная же архитектура используется в аппаратной реализации фильтра (ASIC/FPGA).
                    0
                    По первому способу — после обработки эти «приклеенные кусочки» нужно выкинуть, т.к. они содержат некорректные результаты, ну а сами «ядра» кадров можно склеивать и отдавать дальше.
                    +1
                    Мне нравится подход с построением модели в пространстве состояний (пример построения). Получается по сути система дискретных линейных уравнений, которую легко использовать в режиме реального времени. К тому же строится эта модель может на основе элементарных динамических звеньев. К примеру звено вида 1/Ts-1, где Т — постоянная интегрирования (1/Т — частота среза), а s — оператор Лапласса, является простейшим звеном задержки (простейший ФНЧ). Перемножая несколько элементарных звеньев можно управлять частотами среза, нелинейностями в и во вне полос пропускания и т.д.
                +1
                Для завершенности статьи как-то не хватает достоинств и недостатков FIR-фильтров, а также описания, как на основе ФНЧ построить остальные типы фильтров. В принципе могу и сам написать в комментарии, если у вас нет времени на доработку или кому-то не терпится. Там не много.
                  0
                  Я когда начал тоже подумал добавить ещё и это, но потом решил что статья будет великовата. Если вас не затруднит — добавьте, буду признателен, время на доработку раньше чем завтра не появится уж точно.
                  +2
                  Достоинства:
                  — простота теоретического анализа;
                  — простота практической реализации;
                  — наглядная связь коэффициентов фильтра с отсчетами его импульсной характеристики;
                  — устойчивость фильтра;
                  — при условии симметричности ИХ фильтра, фазовая характеристика линейная, что позволяет уменьшить искажения фронтов
                  Недостатки:
                  — для обеспечения хорошей АЧХ необходим высокий порядок фильтра (доходит до тысяч)

                  Расчет остальных типов фильтров базируется на теореме суммирования преобразований Фурье. За основу берется всепропускающий фильтр, который пропускает без ослабления все частоты. Такой фильтр имеет лишь один коэффициент, который не равен нулю (a0=1). В итоге расчет ФВЧ сводится к разнице:

                  ak, ФВЧ = ak, ВФ — ak, ФНЧ

                  image
                  Получаем ФВЧ с такой же частотой среза, какая была у рассчитанного ранее ФНЧ.
                  Подход к расчету полосового фильтра аналогичный, только берется разница двух ФНЧ:

                  ak, ПФ = ak, ФНЧ1 — ak, ФНЧ2

                  И последнее — режекторный фильтр:

                  ak, РФ = ak, ВФ — ak, ПФ
                    +1
                    > — для обеспечения хорошей АЧХ необходим высокий порядок фильтра (доходит до тысяч)
                    Почему?.. Есть полиномы (многочлены) Баттерворта, дающие гладкие АЧХ даже для порядков знаменателя передаточной функции в пределах пятого.
                      0
                      Здесь имеется ввиду то, что при высоких порядках АЧХ будет приближаться к идеальной — прямоугольной функции, т.е. строго убирать все шумы и пропускать то что надо, в таком фильтре полоса перехода равна 0 (ну это идеал конечно же).
                        0
                        > в таком фильтре полоса перехода равна 0
                        Имеется в виду крутизна фронта/спада, как я понял. В общем случае да, крутизна повышается с увеличением порядка, но на мой взгляд порядок модели 1000+ это уже очень много. Возможно просто потому, что я не сталкивался с задачами где это необходимо.
                        0
                        Посмотрел графики полиномов Баттерворта, впечатляет, как же здорово комментарии на Хабре расширяют кругозор, спасибо.
                          +1
                          Хорошей крутизной отличаются полиномы Чебышева, если мне не изменяет память.
                            +1
                            А также Кайзера. Его методика мне почему-то больше нравится.
                              +2
                              Можно поподробнее? Ссылку?
                                +2
                                  0
                                  Спасибо! Пробежал пока поверхностно. Пока кажется сложнее того что я использую, но там синтез ведется из критерия линейности смещения фаз. Мне пока в этом не было необходимости (гарантировать отсутствие смещения фаз в некоторой полосе). При использовании простых моделей (асимптотически устойчивых заведомо) итак это гарантирует. На досуге изучу потщательнее.
                      +1
                      Алгоритм Ремеза
                      MatLab, octave: функция (g)remez.
                      И, да: Рабинер и Голд «Теория и практика цифровой обработки сигналов». С третьей главы и далее о КИХ фильтрах и весовых функциях подробно и с картинками.
                        0
                        Спасибо огромное за ссылку на алгоритм Ремеза, не мог нигде найти достаточно подробное описание по нему (по крайней мере такое, чтобы я его ещё и понял). Насчет картинок — я вообще думаю подробно написать про весовые функции отдельно, а то когда видишь много-много вариантов — глаза разбегаются.
                          +2
                          Да не за что, мне частенько приходится реализовывать цифровые фильтры на железе (на FPGA или на DSP TI), работа такая. Доедаю не первую собаку, связанную с нелинейностью, эффектом Гиббса и прочим. Обычно для реализации хватает алгоритма из octave функция remez (там, кажется, частный случай Ремеза — Parks-McClellan algorithm и исходники есть). Но бывают и фейлы: например, если попробовать посчитать ФНЧ скажем 1024 порядка с частотой среза 0.5w, то remez из octave выдаёт ошибочные результаты (впрочем причины вроде понятны). А вот MatLab'овский gremez не чудит, по-крайней мере не встречались такие условия.

                          Ещё интересный момент связан с конечной точностью вычислений. В частности, если фильтр реализуется в целочисленной арифметике (по очевидным причинам это быстрее и менее затратно по ресурсам на FPGA к примеру), выбор разрядности конечных и промежуточных результатов, метода округления — тоже целая наука.
                        0
                        Мне тоже скоро надо будет разбираться в нелинейностях и эффекте Гиббса, потому как первый этап — написать хоть какой-нибудь работающий цифровой фильтр (а пару-тройку недель назад я о их существовании и не знал) вроде бы пройден успешно.
                        Ещё раз спасибо, пойду-ка я изучать исходники функции remez.

                        Только полноправные пользователи могут оставлять комментарии. Войдите, пожалуйста.

                        Самое читаемое