Честно простой цифровой фильтр
Вы работаете с АЦП. Получаете результаты преобразования, один за одним. И замечаете, что эти результаты «скачут». А хотелось бы, чтобы стояли, как… Ну, короче, чтобы стояли! Есть много причин, почему отсчеты АЦП могут быть нестабильны. В своей заметке я не говорю об этих причинах. Я говорю о том, как успокоить показания, получая их 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 < EXP(-Ntau / N)
Задавшись двумя граничными значениями Ntau, зная требуемое время переходного процесса и период выборки, мы находим 2 значения Alfa, соответствующие «грубому» и «точному» решениям. Выбираем значение Alfa (ниже я покажу как) в диапазоне:
EXP(-6,908 / N) <= Alfa <= EXP(-2,996 / N) (2)
где <= означает «меньше или равно» А потом легко вычисляем второй коэффициент:
Beta = 1 — Alfa (3)
Вот и все. Коэффициенты фильтра вычислены. Поставив их в утилиту, описанную в заметке коллеги mc2, можете полюбоваться частотной характеристикой. И убедиться, что за требуемое вам время Тпп переходной процесс заканчивается. Если уж говорить о программе для расчета, то можно обойтись и без «моего» метода. Но тут уж вам решать. Расчет по формулам (1). (3) несложен, но гарантирует вам безытерационное решение для заданного переходного процесса.
ПРОГРАММНАЯ РЕАЛИЗАЦИЯ
Если мы работаем с МК, то очень хорошей практикой является использование целочисленной арифметики. Даже работа с 32-разрядными LONG переменными в 8-битном МК здесь выполняется намного быстрее, чем с FLOAT. Поэтому я рекомендую выбирать Alfa и Beta таким образом:
Alfa = Na / 2^k Beta = Nb / 2^k Na + Nb = 2^k Иначе говоря, вместо дробных Alfa и Beta нужно взять целые числа, такие, чтобы их сумма была равна целой степени двойки: 2, 4, 8, 16, 32, 64, 128, 256, 512 и т.д. Вот тогда процессору работы поменьше:
Y(n) = (Na*Y(n-1) + Nb*X(n)) >> k Микроконтроллер умножает 2 раза целое на целое, складывает, производит сдвиг вправо на k разрядов. Несколько микросекунд — и у вас новый отсчет выходной величины.
ПРИМЕР ДЛЯ ЦЕЛОЧИСЛЕННОЙ АРИФМЕТИКИ
Рассмотрим пример, в котором я задался N=100. Тогда условие (2) будет выглядеть так:
0,9333 <= Alfa <= 0,9705 Для целочисленной арифметики выберем однобайтную базу коэффициентов. Иначе говоря, пусть сумма Na и Nb будет равна 256. Тогда k=8 и получаем такие границы целочисленного эквивалента Alfa:
239 <= Alfa*256 <= 248 239 <= Na <= 248 Выбрав любое значение Na из этого диапазона, мы можем быть уверены, что за 100 выборок после подачи сигнала на вход выходной сигнал войдет в зону 0,1. 5 % от установившегося значения. Хотим именно 0,1% — берем меньшее значение Na из указанного диапазона, допускаем 5% — берем большее значение. И уже от выбранного значения Na расчитываем второй коэффициент:
Nb = 256 — Na Теперь формула для работы микроконтроллера такова:
На этом расчет фильтра с целочисленными коэффициентами закончен.
Коллега _pv подсказал в обсуждении хорошую запись вместо (4), полностью эквивалентую по результату:
Здесь имеем экономию на одно умножение (при произвольных соотношениях между коэффициентами) и формула приобретает вид приращения старого значения 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, то есть микроконтроллеру не нужно умножать:
Видим, что расчет (5) имеет на одну операцию умножения меньше, чем расчет (4).
РЕЗЮМЕ
Расчет простого цифрового фильтра нижних частот сводится к таким действиям:
Первое. Зададимся временем переходного процесса, периодом выборок и вычислим, за сколько выборок процесс должен закончиться (1):
N = Тпп / Т
Второе. Вычислим границы возможных значений коэффициента при Y (2):
EXP(-6,908 / N) <= Alfa <= EXP(-2,996 / N) Третье. Выбираем k (рекомендую 8, но можно и иное) и умножаем границы Alfa на 2^k. Выбираем любое из возможных значений в качестве первого коэффициента Na, а второй коэффициент находим как дополнение:
Nb = 2^k — Na
Результирующая формула обработки входной последовательности X(n):
Как показано в примере, иногда можно подобрать Na или Nb, равными 1, что упрощает вычисления.
БЛАГОДАРНОСТИ Хочу выразить искреннюю благодарность уважаемому коллеге angel5a за то, что он вернулся из детского садика (кто читал дискуссию, тот поймет) и очень помог мне с переделкой исходной заметки. Радикальной переделкой, хочу заметить. Также в доработке статьи помогли дискуссии с коллегами _pv и известным в наших кругах товарищем avreal, выступившим на ненашем форуме.