среда, 14 октября 2015 г.

Цифровая фильтрация на ПЛИС. Сглаживание и борьба с импульсными помехами

Когда возникают задачи сглаживания аналогового сигнала, к примеру, с датчиков температуры за сутки (хотя это и так достаточно инерционный параметр), или резкие колебания с датчиков давления около некоего среднего значения (рецессии) или что еще более жизненно – на ваши данные во времени накладывается шум или статический разряд, не поленитесь, загляните в любой математический справочник, к примеру в [1], и откройте для себя полиномы регрессии*, алгоритм скользящего среднего и прочие «вкусности». В случае с аналоговыми цепями все тривиально, фильтрация. А если сигнал уже оцифрован вместе с помехой? Можно-ли в том же микроконтроллере или ПЛИС («тупой» логике) собрать цифровой фильтр и проверить на работоспособность? Разумеется, да. А как схемотехнически реализовать подобный фильтр? Сегодня в продолжении нашего цикла [2, 3] мы рассмотрим данный вопрос в контексте сглаживания, скажем фильтрацию огибающей от импульсных помех в ПЛИС.

Краткий экскурс...

Следует отметить, что сглаживание осуществляется и для временной (ибо осуществляется над самими отсчетами сразу) и в неявном виде для частотной области. Данную операцию над сигналом можно осуществить несколькими способами: применением окон (фильтров) [4] с соответствующими коэффициентами (Блэкмана-Хэрриса, Хафмана и т.д.), использованием полиномов регрессии (экспоненциального, Гауссового, Лангранжа, показательного, гиперболического, параболического, линейного, логарифмического и т.п.), а также таких сглаживающих фильтрующих функций как «скользящее среднее» и медианная фильтрация. К примеру, запись для наиболее часто употребимого экспоненциального сглаживания имеет следующий вид (см. формулу 1).

V = Vp*(1-K)+V*K; (1)
где: V – результат функции, Vp – входное значение (аргумент), K – коэффициент сглаживания.

Основой этих полиномов является – вырождение точек (регрессия). Отсюда видно, что чем ближе к "1" коэффициент сглаживания, тем выше эффект сглаживания (см. рисунок 2).


Рис. 2. Тестовый проект визуализации использования полиномов регрессии.
Сигнал, шум с нормальным распределением и аддитивная смесь сигнала и шума
* Аппроксимация данных с учетом их статистических параметров (например, при обработке экспериментальных данных, полученных в результате измерений процессов или физических явлений) относится к задачам регрессии. Задачей регрессионного анализа является подбор математических формул, наилучшим образом описывающих экспериментальные данные.
Однако, использование полиномов регрессии в цифровой технике обусловлено рядом технических трудностей из-за большого количества проводимых арифметических операций, поэтому чаще ограничиваются более простыми в реализации скользящим усреднением и скользящими медианами. Последние по сути оптимальны для решения двух разных задач. В вопросах сглаживания и устранения многих типов шумов почти идеальны первые, но они совершенно не подходят для фильтрации импульсных помех. Почему? Для ответа на этот вопрос следует пояснить, что же такое это «скользящее среднее».

Простое скользящее среднее или арифметическое скользящее среднее (англ. Simple Moving Average) численно равно среднему арифметическому значений исходной функции за установленный период выборки. Для фильтра «скользящего среднего» (не взвешенного) доступен только один параметр – порядок или длина окна. От него зависит и частота среза и уровень пульсаций в полосе заграждения. Рассмотрим разные длины окон и их влияние на результат (см. таблицу 1).

Таблица 1. Варианты весовых окон «скользящего среднего» над входными данными


Как видно из таблицы 1, чем больше длина окна, тем более сглаженный результат будет на выходе. Однако эта идиллия продолжится до тех пор, пока значения входных отсчетов (сигнала) находятся в пределах, так называемого, нормального или рабочего диапазона. Что будет при появлении импульсного шума (помехи)? Предположим на входе фильтра «скользящее среднее» следующую последовательность: 2, 5, 9, 1, 0, 3, 9 и вдруг раз и появилось 857. Очевидно, что даже при большом интервале (окне) усреднения эту импульсную помеху в виде аномального значения 857 полностью не убрать. Вот тут в игру вступает медианная фильтрация. Но, пока не будем забегать вперед и обо всем по-порядку...

Используемое ПО и оборудование

Для организации работы нам понадобится следующее оборудование и программное обеспечение:
  1. Стендовая макетная плата c ПЛИС из второго цикла [3].
  2. 3-х вольтовый источник питания мощностью от пяти ватт.
  3. ПК с установленной инструментальной системой проектирования логических матриц Xilinx Foundation Series 3.1i/4.1i, ISE Webpack или Xilinx ISE Design Suite 14 (или выше) [5, 6].
  4. Программатор JTAG, к примеру Xilinx JTAG Download Parallel Cable или Xilinx Platform Cable USB, подключенный к JTAG- разъему программирования [7].
Создание проекта, сборка и прошивка ПЛИС

Запустив САПР Xilinx, создадим новый проект и выберем соответствующую модель матрицы. После чего  можем войти в схемотехнический редактор, нажав третью кнопку на вкладке «Design Entry» менеджера проектов Project Manager и приступить к моделированию (см. рисунок 3). Все счетчики и триггеры в проекте должны быть асинхронные.

Рис. 3. Создание проекта в среде Xilinx

По окончании моделирования проверяем и собираем проект по нажатию кнопки «Implementation» в окне Project Manager (см. рисунок 4).

Рис. 4. Сборка проекта (имплементация) в среде Xilinx

После сборки проекта среда создаст файл прошивки формата JED с сигнатурой матрицы XC95288XL7-TQ144, который и нужно будет прошить в CPLD. Процесс прошивки рассмотрен нами ранее подробно в предыдущих материалах и останавливаться на нем не будем.

Нижний уровень. Реализация алгоритма «скользящее среднее»

Каждое значение аргумента выходной функции «скользящего среднего» – есть не что иное, как сумма соседних (предыдущих значений), деленное на длину окна. Число точек W, учавствующих в расчете такого среднего, называют окном скользящего усреднения. Чем оно больше, тем больше данных учавствует в расчете среднего и тем более сглаженная кривая получается. Если обратиться к таблице 1 (см. выше), то общую формулу для данной весовой функции можно записать как выражение (2).
Yw = SUM (Yi) / W; (2)
где: Yw – результат скользящего среднего, SUM() – сумма значений всех отсчетов на длине окна W.

При малых значениях длины окна W сглаженная кривая будет практически повторять ход изменения данных, а при больших W эти кривые будут отражать лишь закономерность их медленных вариаций, что соответственно сказывается на поведении окна в ВЧ-области (об этом чуть далее).

Поскольку мы имеем дело с двоичной логикой, то длину окна целесообразно выбирать кратной степени двойки. Что это дает? А дает это простоту реализации операции деления: скажем, если необходимо поделить наш результат сложения на число 2N (N- вес разряда), то это будет равноценно сдвигу разрядов числа вправо на N- позиций. Этот сдвиг реализуется просто переподключением нужных весовых разрядов на схеме. Всего-то? Именно так, в случае с целочисленной арифметикой. Предположим к нам поступают данные постоянно (беспрерывно) к необходимо реализовать «скользящее» окно с накопительным буфером на два слова, т.е. к примеру, длина окна W = 2. Схемотехнически реализация будет следующей (см. рисунок 5).


Рис. 5. Реализация фильтра «скользящее среднее» с длиной окна W = 2

Описание принципа работы

На цепочке (каскадное соединение) двух восьмиразрядных асинхронных D- триггеров FD8 реализованы накопительный буфер и «скользящее окно». По фронту первого такта производится запись поступивших данных на нижний буфер (единократно и на верхний буфер за счет разрешения прохождения через нижний BUFT8, управляемый нулем). При этом на входах сумматора будут присутствовать одинаковые сигналы, которые после деления на два (обычным сдвигом на один разряд с учетом выхода переполнения сумматора) дадут тот же сигнал, что и на входе. Сделано это для адекватной работы «скользящего среднего» в первый момент (центрированный алгоритм расчета), однако по первому же спаду тактов сработает триггер-защелка FDС, которая переподключит верхнее накопительное плечо FD8 в рабочий режим и в дальнейшем, по приходу второго фронта данные с выхода нижнего накопительного плеча FD8 запишутся в верхнее плечо FD8. Пришедшее второе слово (данные) по второму такту запишется в первое плечо буфера. Таким образом, по фронту каждого такта данные (слова, байты) со входа INDATA будут сдвигаться и каждый раз на входах сумматора будут данные окна. По спаду тактов результат деления будет защелкиваться выходным буфером FD8 и поступать на выход OUTDATA. По входу CLR есть возможность сброса фильтра в начальное состояние по сигналам запуска, при необходимости. При желании, «скользящее среднее» легко наращивается до длины** окна кратной степени двойки (22, 23  и т.д).
** Следует отметить, что при длине окна более восьми, значительно возрастает время реакции звена, но с этим приходится смириться. Время реакции напрямую определяется длиной окна и скоростью поступления данных, сами же данные поступают синхронно (по фронту) с тактовыми импульсами.
Оценка оптимальности порядка фильтра «скользящее среднее»

Допустим задан сигнал вида sin(x) и некая случайная составляющая с нормальным распределением (см. выше рисунок 1), накладываемая на сигнал в той же временной области. Данный сигнал мы решили сгладить с помощью фильтра «простое скользящее среднее». Выбрали порядок фильтра N. Как оценить оптимальность выбранного порядка N? То есть, какая длина окна будет достаточна со строгой математической точки зрения? Очевидно, что следует определить частоту среза фильтра и задаться допустимым уровнем пульсаций в полосе пропускания.

Мы знаем, что от длины окна зависит и частота среза и уровень пульсаций в полосе заграждения. Данный фильтр можно рассматривать как фильтр с конечной импульсной характеристикой (КИХ) с одинаковыми коэффициентами, равными 1/N, где N – длина окна (в этом случае окно и есть импульсная характеристика фильтра). Частотная характеристика такого фильтра есть преобразование Фурье от единичного импульса (правильнее конечно использовать термины Z- преобразования), это будет функция вида sin(x)/x, в которой главный лепесток определяет полосу пропускания, а остальные боковые лепестки задают полосу заграждения. Разумеется, чем длиннее окно (больше порядок фильтра), тем уровень боковых лепестков будет ниже, но ширина главного лепестка будет равна 1/N, то есть с ростом порядка фильтра она будет уменьшаться, а значит частота среза будет меньше. Логично? В случае одиночного сигнала, оптимальным в нашем случае будет такой порядок фильтра, который даст максимально близкую к частоте вашего синуса частоту среза. Количественно это будет определяться формулой 3.
N = Fd/2/Fs; (3)
где: N – порядок фильтра, Fd – частота дискретизации, Fs – частота сигнала.

Правильнее будет сказать, что оптимальная длина (она же порядок фильтра) будет определяться следующим выражением (4).
Nопт = round(N); (4)
где: Nопт – оптимальная длина окна, round() – есть округление в нижнюю сторону, так как требуется, чтобы частота среза фильтра была выше частоты сигнала (нам то необходимо, чтобы сигнал заданной частоты прошел через фильтр).

Таким образом, при заданных исходных параметрах: частоте сигнала, частоте дискретизации и разных уровнях соотношения сигнал/шум (SNR), мы будем всегда иметь одну оптимальную величину N. Отметим, что фильтр всегда оптимален для определенных характеристик сигнала, характеристики изменились – фильтр уже неоптимален. Да, вот такой вот недостаток. Но что же вы хотели от простого фильтра с одним настраиваемым параметром? Его задача – лишь сглаживание данных.

Практическое приложение подобной оценки

Рассмотрим конкретный пример расчета подобной оценки, скажем для гитарного тюнера. Это такое устройство или программа для определения звучащей ноты, т.е. частоты нотного ряда, на основе спектроанализа [8...10]. Допустим имеем сигнал (набор данных, сэмплов) с частотой дискретизации 44100 Гц и разрешением 32 бита на сэмпл. Среднее из скольки сэмплов нужно взять, чтобы остались частоты ниже 1000 Гц? То есть стоит задача отфильтровать высокие ноты на гитаре.

Обратимся к формуле (4) и подставим исходные данные:
Nопт = round(Fd/2/Fs) = round(44100/2/1000) = round(22.05) = 22.
Как видим, порядок фильтра или длина окна должна быть равна 22 сэмпла. Получается, если нам необходимо убрать перед анализатором спектра все частоты выше 1000 Гц или что эквивалентно введению окна ФНЧ, то мы должны выдавать один сэмпл, сформированный как среднее из 22 входных сэмплов аудиосигнала. Вот так просто.

Нижний уровень. Одномерная медианная фильтрация

Медианой набора данных 'х1, х2, …, при нечетном является средний по значению элемент, получающийся при упорядочивании набора по возрастанию или убыванию. Для четных медиану обычно определяют как среднее арифметическое двух средних отсчетов упорядоченной последовательности. Медианный фильтр представляет собой оконный фильтр, последовательно скользящий по массиву сигнала, и возвращающий на каждом шаге один из элементов, попавших в окно фильтра. Таким образом, медианная фильтрация осуществляет замену значений отсчетов в центре окна медианным значением исходных отсчетов. На практике длина окна фильтра для упрощения обработки данных, как правило, устанавливается с нечетным числом отсчетов. Следовательно, можно выделить два класса медианных фильтров: фильтры с минимальной длиной окна равной трем и фильтры с длиной окна более трех (к примеру, 5, 7, 9, 11, 13, 15 и т.д.).

Поскольку данный метод устраняет локальные (единичные) выбросы внутри окна, то медианная фильтрация срезает импульсные помехи на раз-два, оставляя полезный сигнал (сохраняя перепады) почти без изменения (см. рисунок 6). И что важно: медианный фильтр является нелинейным КИХ- фильтром, который может полностью затирать частотную информацию в сигнале, импульсная реакция фильтра равна нулю, а отклик на ступенчатую функцию равен 1. Поэтому будьте внимательны при выборе данного вида фильтров применительно к вашим задачам.

Рис. 6. Скользящие медианы
(набор данных – окружности, окно W = 3 – сплошная кривая,  окно W = 7 – пунктирная кривая)

Для того, чтобы найти значение скользящей медианы в точке X, вычисляется медиана значений ряда во временном интервале [X - delta, X + delta]. Соответствующее значение называется (2delta + 1)-точечной скользящей медианой. Исходя из этого, реализуем наименьший возможный медианный фильтр с длиной окна равной трем. Алгоритм следующий:
  1. Имеем набор данных (8-ми битных) A, B, C.
  2. Производим первую сортировку по убыванию: if (A<B) and (A<C), то считаем A – наименьшим элементом и далее выбираем между B и С, наименьший элемент является медианой и выдаем его на выход.
  3. Если условие по первой сортировке не выполняется, то производим вторую сортировку по убыванию: if (B<A) and (B<C), то считаем B – наименьшим элементом и далее выбираем между A и С, наименьший элемент является медианой и выдаем его на выход.
  4. Если условие по первой и второй сортировке не выполняется, то производим третью сортировку по убыванию между A и B, наименьший элемент является медианой и выдаем его на выход.
В итоге наша схема (см. рисунок 5) преобразуется следующим образом (см. рисунок 7).

Рис. 7. Реализация фильтра «скользящая медиана» с длиной окна W = 3

Описание принципа работы

На цепочке (каскадное соединение) трех восьмиразрядных асинхронных D- триггеров FD8 реализованы накопительный буфер элементов A, B, C и «скользящее окно» длиной W = 3. По фронту первого такта производится запись поступивших данных на нижний (первый) буфер, по второму и последующим тактам производится сдвиг и запись последующих данных (см. таблицу 2).

Таблица 2. Процесс сдвига данных в накопительном буфере медианного фильтра c длиной окна W = 3


Сортировку осуществляем с помощью семи 8-ми битных компараторов чисел COMPM8, выдающих логическую "1" на выходе GT при выполнении условия по входу A>B и  логическую "1" на выходе LT при выполнении условия по входу B>A. Реализация условия (a<b) and (a<c) осуществляется двумя компараторами над набором данных A, B, C, результат которого по И-НЕ коммутирует низким уровнем буфер с тремя состояниями BUFT8 прохождение полученной медианы по итогам сравнения элементов B и C. Реализация условия (b<a) and (b<c) осуществляется двумя компараторами над набором данных A, B, C, результат которого по И-НЕ коммутирует низким уровнем буфер с тремя состояниями BUFT8 прохождение полученной медианы по итогам сравнения элементов A и C. При невыполнении одновременно двух условий (a<b) and (a<c) и (b<a) and (b<c), т.е. когда их логические уровни равны по И-НЕ коммутирует низким уровнем буфер с тремя состояниями BUFT8 прохождение полученной медианы по итогам сравнения элементов A и B. Полученные медианы сведены монтажным ИЛИ на вход выходного регистра, тактируемого коммутатором. Поскольку для работы фильтра требуется минимум 3 точки, то первичное накопление осуществляется с помощью узла коммутатора. В первый момент времени инверсный сигнал с выхода RS- триггера разрешает работу счетчика и активирует передачу третьего такта с выхода счетчика-делителя на три. По приходу спада третьего такта сигнал с выхода дешифратора весового разряда b11 = 3 защелкивает вход установки RS- триггера, после инверсии этот же сигнал выключает счетчик и активирует передачу тактов по спаду на выход коммутатора, минуя счетчик, для осуществления процесса записи результата медианного значения сдвигаемых данных в выходной буфер на регистре FD8. Таким образом, работа медианного фильтра начинается при поступлении третьего такта, что равнозначно получению третьей точки из выборки.
Лицензионное соглашение :)

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

При реализации медианного фильтра в системах ЦОС следует обратить внимание на то, что при появлении импульсной помехи на границах окна предпочтительно варьировать длину окна, к примеру вести расчет с двумя длинами окон одновременно, либо, к примеру, использовать алгоритм быстрой медианной фильтрации [11]. Также следует учитывать фактор задержки выходного сигнала фильтров как «скользящего среднего», так и «скользящих медиан», определяемый порядком фильтра.

Спецификации, тестовый проект визуализации использования полиномов регрессии (под TDL) и схемотехническая реализация рассмотренных модулей сглаживания и фильтрации в среде САПР Xilinx Foundation Series 4.1i (файл plis3_res.zip) вы можете загрузить с сайта нашего журнала http://radioliga.com (раздел «Программы»), а также с сайта разработчика [12]. Если тема представляет для вас интерес – пишите, задавайте вопросы.

Ресурсы
  1. И.Н. Бронштейн, К.А. Семендяев. Справочник по математике для инженеров и учащихся ВТУЗов. – М.: Наука, 1981, 720 с.
  2. Е. Бадло, С. Бадло. ПЛИС. Часть 1 или... Удаленный прошивальщик по сети TCP-IP. – Радиолюбитель, Минск, 2013, №12, с.34
  3. Е. Бадло, С. Бадло. ПЛИС. Часть 2 или... Реализация псевдо-UART. – Радиолюбитель, Минск, 2014, №2, с.45
  4. Использование оконных функций в задачах цифрового спектрального анализа. Примеры и рекомендации http://www.dsplib.ru/content/winex/winex.html, онлайн-версия от 01.10.2011
  5. Xilinx download http://www.xilinx.com/support/download.html#download.html?&_suid=138996081380407635947911615178
  6. ISE In-Depth Tutorial http://www.xilinx.com/support/documentation/sw_manuals/xilinx11/ise11tut.pdf
  7. С. Бадло. JTAG.Xilinx программатор. – Радиолюбитель, Минск, 2008, №7, с.38
  8. Е. Бадло, С. Бадло. Виртуальные приборы. Спектроанализатор своими руками. – Радиолюбитель, Минск, 2009, №3, с.32
  9. С. Бадло. Быстрое преобразование Фурье. Практика использования. Часть 2. – ПРОграммист, 2010, №1, с.33
  10. С. Бадло. Быстрое преобразование Фурье. Практика использования. Часть 1. – ПРОграммист, 2010, №8, с.46
  11. Б.В. Бардин. Быстрый алгоритм медианной фильтрации. Научное приборостроение, 2011, том 21,  №3, с.135-139 http://213.170.69.26/mag/2011/full3/Art16.pdf
  12. Ресурсы к проекту http://raxp.radioliga.com/cnt/s.php?p=plis3_res.zip

2 комментария:

  1. Алгоритм фильтрации с помощью скользящей средней очень грубый и дающий сильное запаздывание. Лучше всего с фильтрацией справляются алгоритмы Калмана (Kalman) и алгоритм Джурика (Jurik)

    ОтветитьУдалить
  2. Преимущества и недостатки скользящего среднего уже рассмотрены выше. Да, он дает запаздывание. Но замечу, что идеального фильтра не бывает, есть понятие оптимального фильтра под свою задачу, т.е. под частотно-временные параметры сигнала. К тому же, тот же Калман не такой уж и скоростной, напротив у него низкая пропускная способность. Взглянем:
    //prediction update
    //omit x = x
    p = p + q;

    //measurement update
    k = p / (p + r);
    x = x + k * (data - x);
    p = (1 - k) * p;

    Дохрена ведь операций.

    ОтветитьУдалить

В комментариях уважайте собеседника, внимательно читайте посты и не додумывайте. Просьбы и предложения из разряда: «можно ваш Skype/Viber/телефон», «напишите мне в vk/FB», а также другие им подобные — игнорируются. Выход новых версий ПО, внешняя ссылка, переставшая работать с течением времени и т.п. не является основанием для претензий. Желающие спокойно подискутировать и высказаться — Welcome. Желающие спонсировать блог — Donate. Нарушение этих простых правил ведет к бану и удалению комментариев без предупреждения.