Во многих цифровых устройствах для преобразования аналоговых сигналов используется АЦП. Часто аналоговые сигналы содержат нежелательный высокочастотный шум.
Чтобы "очистить" сигнал от этих шумов применяются аналоговые RC фильтры низких частот, которые устанавливаются после источника сигнала. Такой подход не всегда идеален и практичен. Например, для больших постоянных времени требуются большие значения R и C.
В качестве альтернативы, можно "очистить" зашумленный сигнал с помощью цифрового эквивалента аналогового RC фильтра нижних частот.
По сути, программа этого цифрового фильтра состоит всего из двух строчек на Си.:
Dacc = Dacc + Din — Dout
Dout = Dacc/K
где Dout — выходное значение фильтра, Din — входное значение фильтра, K — постоянный коэффициент, который рассчитывается по формуле:
K = T x SPS
где T — это постоянная времени фильтра, SPS — частота дискретизации АЦП.
Dacc и Dout должны сохранять свои значения, после выполнения алгоритма. Если алгоритм реализовать в виде функции, то эти переменные можно просто сделать статическими.
Для 8-ми разрядных входных данных алгоритм цифрового фильтра в Си коде может выглядеть так:
Для наглядности рассмотрим реальный пример использования этого алгоритма в микроконтроллере AVR atmega16. Допустим мы хотим сымитировать RC фильтр с постоянной времени 0.001 c. Схема представлена на рисунке ниже.
R1 = 10 кОм
C1 = 0.1 мкФ
Trc = R1*C1 = 10000 Ом * (0.1/1000000) Ф = 0.001 с
F = 1/(2*Pi*R1*C1) = 1/(6.28 * Trc) =
Тактовая частота модуля АЦП в микроконтроллерах AVR зависит от его тактовой частоты и внутреннего предделителя. Допустим, наш микроконтроллер тактируется от внутреннего генератора с частотой 8 МГц, а предделитель в модуле АЦП установлен равным 64. Тогда тактовая частота модуля АЦП будет равна:
Fadc = Fcpu/Pre = 8000000/64 = 125000 Гц
Из этой частоты можно рассчитать частоту дискретизации АЦП при работе в режиме непрерывного преобразования. Она равна отношению тактовой частоты АЦП к количеству тактов, которые требуются для выполнения одного преобразования. По даташиту можно узнать, что одно преобразование выполняется за 13 тактов (если это не первое преобразование).
Fs = Fadc/n = 125000/13 =
Итак, частота дискретизации равна 9600 Гц, а постоянная времени 0.001 с. Коэффициент фильтра будет равен:
K = SPS x T = 9600 * 0.001 = 9.6 =
Теперь все данные известны и можно написать тестовую программу для проверки алгоритма.
Я сделал эту программу очень простой. АЦП работает в режиме непрерывного преобразования. В прерывании 8-ми разрядный результат преобразования обрабатывается алгоритмом и записывается в порт C. К порту C подключен схема R-2R ЦАПа на основе резисторов, чтобы можно было сравнить полученный сигнал с сигналом от аналогового RC фильтра. Тактовая частота микроконтроллера atmega16 — 8МГц, коэффициент предделителя АЦП — 64. Тестовая схема показана на рисунке ниже.
Код программы
Результат моделирования
Результат моделирования программы в Proteus`e показан на рисунке ниже. Красный сигнал — это входной меандр частотой 200 Гц, синий — сигнал на выходе RC фильтра с постоянной времени 0.001 с, а желтый — сигнал обработанный микроконтроллером. Он имеет ступенчатую форму, так как после ЦАПа не подвергался фильтрации.
Как видно из рисунка форма сигнала микроконтроллера достаточно точно повторяет сигнал с выхода RC фильтра.
Для наибольшего быстродействия коэффициент К лучше выбирать кратным степени 2 (например 2, 4, 8..), тогда компилятор будет заменять операцию деления сдвигами. В противном случае при высокой частоте дискретизации, микроконтроллер может не успевать рассчитывать следующее выходное значение фильтра. Я столкнулся с этим при моделировании.
Также необходимо учесть тот момент, что при больших значениях коэффициента К, переменная Dacc должна иметь достаточную разрядность, чтобы не наступало ее переполнение.
Вступление издалека
Недавно передо мной встала достаточно интересная задача, с которой я раньше никогда не сталкивался — борьба с шумом. Мы принимали сигнал с датчиков на аналогово-цифровой преобразователь (АЦП)
А так как данная тема для меня была (хотя и сейчас есть кое-где) темным лесом, я пошел мучить вопросами гугл, мне показалось освещена эта тема не очень подробно и доступно, поэтому решил написать статью с примером разработки и готовым исходником.
Ближе к делу
Цифровые фильтры могут быть двух видов – с конечной и с бесконечной импульсной характеристикой (КИХ и БИХ). Для решения моей задачи подходит КИХ-фильтр, поэтому про него и расскажу.
Для начала посмотрим как же он работает:
Здесь показан пример фильтра нижних частот, как видно на рисунке, этот фильтр пропускает нижние частоты, а все остальные старается отсечь (подавление), или хотя бы ослабить (переход). Отклонения в полосе пропускания и полосе подавления выбираются в зависимости от принимаемого сигнала, но при использовании различных весовых функций, на них могут накладываться определенные ограничения. Например, если используется весовая функция Хэмминга, то эти отклонения будут равны между собой.
Ширина полосы перехода ∆F зависит от длины фильтра и от весовой функции (для функции Блэкмена ∆F=5,5|N).
Работает фильтр довольно просто: фильтр получает значения, с помощью коэффициентов преобразует их и выдаёт выходную последовательность, тогда с формулой самого фильтра всё понятно:
Она реализуется через цикл, но постойте, а где же взять нужные коэффициенты? Вот тут-то как раз и зарыта собака (и не одна).
Параметры фильтра
Естественно для разных фильтров нужны разные коэффициенты, и для этого нужно определиться с параметрами фильтра, это обычно сначала делается теоретически (с умным видом прикидываем какая у нашего сигнала частота, потом частоты, которые надо отсеивать), а потом изучаем АЧХ реальных измерений (и осознаем, как сильно мы ошибались).
По этим АЧХ мы определяемся с идеальной частотной характеристикой (какие частоты проходят свободно, какие мы убираем и как сильно), теперь нам нужна идеальная импульсная характеристика её можно посчитать как Фурье-образ от идеальной частотной:
где H_D(w) – идеальная характеристика.
Но можно пойти и по более простому пути – есть уже заранее вычисленные идеальные импульсные характеристики, например для фильтра нижних частот формула выглядит следующим образом:
где fc и wc – частота среза.
Итак, осталось уже немного идеал идеалом, а мы имеем дело с практикой, и нам нужна «реальная» импульсная характеристика. Для её расчета нам понадобится весовая функция w(n), их есть несколько разновидностей, в зависимости от требований к фильтру (Хэмминга, Хеннинга, Блэкмена, Кайзера, о них не говорю, ибо статья и так большая), в нашем случае я использую функцию Блэкмена:
где N – длина фильтра, т.е. количество коэффициентов.
Теперь надо перемножить идеальную импульсную характеристику и весовую функцию:
Финишная прямая
Теперь мы готовы рассчитать выходные значения, по формуле фильтра, она самая первая в этой статье, ну вот и всё, в завершение привожу исходный код фильтра:
Вы работаете с АЦП. Получаете результаты преобразования, один за одним. И замечаете, что эти результаты «скачут». А хотелось бы, чтобы стояли, как… Ну, короче, чтобы стояли!
Есть много причин, почему отсчеты АЦП могут быть нестабильны. В своей заметке я не говорю об этих причинах. Я говорю о том, как успокоить показания, получая их AS IS. И как сделать это максимально просто. При этом, возможно, не имея ни малейшего понятия о науке под названием «цифровая обработка сигналов».
Эта заметка написана в качестве полной замены предыдущей заметки. Ту лучше не читать 🙂
ПОСТАНОВКА ЗАДАЧИ
Имеется последовательность кодов, дискретно во времени представляющих физический сигнал. Мы будем говорить о последовательности кодов с АЦП. Физический сигнал и полученная последовательность кодов имеют шумы, выбросы, «болтанку» — назвать можно как угодно. Наша задача — сгладить входную последовательность, то есть, выдать выходную последовательность такую, чтобы влияние шумов было уменьшено.
При этом мы стремимся выполнить задачу максимально простыми программными средствами обычного микроконтроллера.
Более того, у нас поставлено еще одно условие. Допустим, что нам неудобно накапливать данные, а потом обрабатывать их и выдавать результат один раз на N полученных кодов. То ли буфер для данных негде организовать, то ли темп выдачи результатов должен совпадать с темпом получения кодов от АЦП, то ли еще что. Но условие поставлено.
И тут нам на помощь приходит цифровой рекурсивный фильтр, самой простой реализацией которого является фильтр первого порядка. Вот его и будем делать.
ТЕРМИНОЛОГИЯ
Мы хотим реализовать блок программной обработки кодов АЦП, на который поступают входные отсчеты, обрабатываются и преобразуются в выходные отсчеты. Вот этот блок обработки и будем называть цифровым фильтром: 
В каждый из моментов времени t1-2Т , t1-T, t1, t1+T и так далее на входе появляется очередной отсчет последовательности Х, а цифровой фильтр выдает новый отсчет последовательности Y. Если на каждый входной отсчет фильтр выдает и выходной отсчет, то это и есть классический цифровой фильтр.
Иногда поступают не так: набирают некоторое количество входных отсчетов и затем обрабатывают их. Таким образом, на несколько входных получают один выходной отсчет. Это, как правило, уже задача получения того или иного интегрального значения (например, среднего). Она настолько распространена при фильтрации шумов, что часто воспринимается как «обычная». И поэтому классический цифровой фильтр, выдающий отсчеты с каждым входным, называют "скользящим усреднением", как нечто не совсем обычное для усреднения. Что ж, значит мы рассматриваем скользящее усреднение. Но понимаем, что это обычный цифровой фильтр 🙂
Важной особенностью цифрового фильтра является то, что входные отсчеты приходят с постоянным темпом, в нашем случае — с периодом Т. Обратная величина называется частотой дискретизации Fд и играет огромную роль во всех параметрах цифрового фильтра. Самое простое, что необходимо себе четко представлять: частотная характеристика цифрового фильтра симметрична относительно Fд/2 и полностью повторяется за пределами Fд. Из чего следует, что фильтр нижних частот, который мы проектируем сейчас, не способен подавлять сигналы с частотой Fд, 2Fд, 3Fд, и т.д. Этого понимания пока достаточно для нашей задачи.
Математически цифровой фильтр 1-го порядка описывается различными способами. Мы будем использовать такое представление:
Y(n) = Alfa*Y(n-1) + Beta*X(n)
То есть, очередной отсчет Y(n) получаем путем взвешенного сложения предыдущего выходного отсчета Y(n-1) и нового входного кода X(n). При этом обычно коэффициент «усиления» фильтра желательно иметь равным 1. Для этого нужно, чтобы выполнялось
Alfa + Beta = 1
В рамках данной заметки задача расчета цифрового рекурсивного фильтра 1-го порядка состоит в нахождении коэффициентов Alfa иBeta с учетом удобства их использования в микроконтроллере (МК) для цифровой фильтрации отсчетов.
РАСЧЕТ ФИЛЬТРА
Цифровой фильтр (равно как и не цифровой) имеет как бы 2 лица: частотную характеристику и переходную (временнУю) характеристику. Задачу расчета можно ставить как получение требуемой частотной либо требуемой переходной характеристики. Я люблю исходить из заданного поведения во временной области. И объясню почему.
Дело в том, что частотная характеристика фильтра первого порядка… как бы это выразиться… простенькая. Возьмем 2 фильтра с частотой дискретизации 200 Гц. Первый с частотой среза по уровню -3 дБ, равной 5 Гц, второй — с частотой 1 Гц:
Как видим, в большей части частотной характеристики подавление шумов составляет всего 15. 30 дБ, а у второго — 20..40 дБ. Вспомним, я отмечал выше, что за пределами Fд/2 (у нас это 100 Гц) характеристика симметрична и в районе 200 Гц подавления шумов снова нет. Если нам нужно «взрослое зубастое подавление», то необходимо строить более серьезные фильтры.
Но все же, второй фильтр лучше давит помехи! А что, если частоту среза еще понизить?
И тут оказывается, что в погоне за парочкой децибел дополнительного подавления мы совсем забыли о временнОм «лице» фильтра. А давайте глянем, как реагируют эти два фильтра на ступенчатое входное воздействие (сигнала не было, а потом он вдруг стал равен 1 и таким и остался).
Что мы видим? Если говорить о времени переходного процесса, то первый фильтр (который похуже фильтрует) отрабатывает входной скачек примерно за 160 мс, а второй, наш передовик с подавлением 20. 40 дБ — почти за 800 мс. Вот так-то. Лучше фильтрация — хуже переходной процесс. Поэтому и нужно выбрать некий оптимум.
Вот, понимая это, я и предлагаю: исходить из требований по быстродействию. Задавшись временем реакции на ступенчатое входное воздействие, мы получим параметры фильтра (вот те самые Alfa иBeta), а параметры фильтрации примем уж какие получатся. Парочка децибел туды-сюды уже мало что изменят, а быстродействие будет известно.
Вот почему в качестве исходного требования я выбираю обеспечение заданного времени переходного процесса при подаче на вход ступенчатого воздействия. Рассмотрим чуть пристальнее, как же фильтр отрабатывает такой сигнал. Если дать скачек от нуля до некоторого значения, обозначенного 100%, то на выходе фильтра увидим:
Как видим, на 8-м отсчете фильтр уже отработал 90% входного воздействия, на 10-м — 95%, и так далее. Обычно принято говорить о «недоработанном», т.е. о тех 10 или 5 процентах, на которые фильтр еще «врет» к какому-то отсчету. Говорят еще, что погрешность установления выходного сигнала составляет столько-то процентов. Далее я привожу формулы для погрешности установления от 5 до 0,1%. В первом случае переходной процесс закончился «на глазок», а точности 0,1% обычно достаточно для того, чтобы считать процесс полностью законченным.
Разница между этими 5%-ым и 0,1%-ым фильтрами не столь уж велика. Я предполагаю, что все фильтры, которые будут разработаны по описываемой методике, находятся в континиуме между этими двумя крайними точками. В качестве характеристики фильтра по степени «законченности» переходного процесса введем такой параметр: Ntau — и вычислим его крайние значения:
LN(1/5%) = 2,996
LN(1/0,1%) = 6,908
Ntau = 2,996. 6,908
Смысл Ntau примерно таков: чем более жесткие требования мы предъявляем к завершенности переходного процесса, тем больше это число. Так что нижней границе значений Ntau соответствует «грубый» процесс, когда мы спешим считать отработку законченной, а верхней границе — «точный» переходной процесс.
Итак, разработчик задался временем переходного процесса Тпп, за которое фильтр отработает скачек на 95. 99,9%. Что еще нужно знать? Время выборки — тот период, с которым на вход фильтра поступают выборки и с таким же темпом с фильтра уходят результаты. Он у нас обозначен Т.
Ясно, что за время переходного процесса будет обработано много отсчетов. Сколько?
N = Тпп / Т (1)
И наша задача — выбрать значения Alfa иBeta так, чтобы уложиться в Тпп.
Оказывается, все очень просто. Должно выполняться условие:
Alfa > k
Микроконтроллер умножает 2 раза целое на целое, складывает, производит сдвиг вправо на k разрядов. Несколько микросекунд — и у вас новый отсчет выходной величины.
ПРИМЕР ДЛЯ ЦЕЛОЧИСЛЕННОЙ АРИФМЕТИКИ
Рассмотрим пример, в котором я задался N=100. Тогда условие (2) будет выглядеть так:
0,9333 > 8 (4)
На этом расчет фильтра с целочисленными коэффициентами закончен.
Коллега _pv подсказал в обсуждении хорошую запись вместо (4), полностью эквивалентую по результату:
Y(n) = Y(n-1) + [Nb*(X(n) — Y(n-1)) >> 8]
Здесь имеем экономию на одно умножение (при произвольных соотношениях между коэффициентами) и формула приобретает вид приращения старого значения Y(n) на разность между входом и выходом, умноженную на коэффициент Beta. Чем меньше этот коэффициент, тем медленнее следит выход за входом — более сильная фильтрация.
Приведу сишный код, соответствующий формуле от _pv:
Здесь учтено, что при делении (сдвиге вправо) в целочисленной арифметике теряются разряды. Поэтому промежуточный результат сохраняется в переменной z в масштабе суммы (до деления), а точнее, в масштабе выходной величины, деленной на Beta. Кстати, здесь и ответ на вопрос о необходимом диапазоне представления z: по максимальному значению выходной величины, деленному на Beta. Например, однобайтные переменные и коэффициент Nb, равный 1 (при последующем сдвиге на 8 бит, как в показанной функции) требуют двухбайтного z.
Еще один фокус — и я откланяюсь 🙂
Если не стремиться к конкретной цифре погрешности установления (0,1% или 5% или что-то еще), то можно попытаться выбросить одно умножение, выбрав меньшее из чисел Na и Nb равным 1.
В нашем примере меньшее из чисел именно Nb и легко посчитать, что оно может быть в диапазоне 8. 17. Вычислим следующие границы:
256 / Nb = 32. 16
Здесь числа равны степени двойки случайно. Важно то, что из диапазона допустимых значений 256 / Nb мы выбираем одно из чисел, равное именно степени двойки. Например, в нашем случае, берем 16.
Тогда 16 — это сумма «удобных» целочисленных Na и Nb, а меньшее из них равно 1, то есть микроконтроллеру не нужно умножать:
Y(n) = (15*Y(n-1) + X(n)) >> 4 (5)
Видим, что расчет (5) имеет на одну операцию умножения меньше, чем расчет (4).
РЕЗЮМЕ
Расчет простого цифрового фильтра нижних частот сводится к таким действиям:
Первое. Зададимся временем переходного процесса, периодом выборок и вычислим, за сколько выборок процесс должен закончиться (1):
N = Тпп / Т
Второе. Вычислим границы возможных значений коэффициента при Y (2):
EXP(-6,908 / N) > k
Как показано в примере, иногда можно подобрать Na или Nb, равными 1, что упрощает вычисления.
БЛАГОДАРНОСТИ
Хочу выразить искреннюю благодарность уважаемому коллеге angel5a за то, что он вернулся из детского садика (кто читал дискуссию, тот поймет) и очень помог мне с переделкой исходной заметки. Радикальной переделкой, хочу заметить.
Также в доработке статьи помогли дискуссии с коллегами _pv и известным в наших кругах товарищем avreal, выступившим на ненашем форуме.