Банк цифровых фильтров
В операциях цифровой обработки сигналов особое внимание уделяется цифровой фильтрации, которая по объему вычислений в среднем занимает от 20 до 60%. Для того чтобы свести к минимуму возможные потери информации и повысить качество ее обработки, цифровые фильтры должны обеспечивать возможность быстрой работы с большими блоками данных. Одним из вариантов решения этой задачи является использование банка цифровых фильтров.
В узком смысле цифровой фильтр — это частотно-избирательная цепь, обеспечивающая селекцию цифровых сигналов по частоте [2]. После выполнения цифровой фильтрации мы, как правило, получаем интересующий нас сигнал, то есть сигнал, несущий нужную нам информацию в виде, удобном для последующей обработки. Соответственно, к параметрам цифровых фильтров в современных системах цифровой обработки сигналов начинают предъявляться повышенные требования. Частоты, на которых работают цифровые фильтры, нередко достигают нескольких сотен мегагерц и более. Соответственно, растет ширина полос цифровых фильтров. Это ведет к увеличению объема и скорости вычислений, а значит, и к резкому росту аппаратных затрат. Для того чтобы свести к минимуму возможные потери информации и повысить качество ее обработки, цифровые фильтры должны обеспечивать возможность быстрой работы с большими блоками данных. Одним из вариантов решения этой задачи является использование банка цифровых фильтров.
Структура банка цифровых фильтров
Банк цифровых фильтров предназначен для разбиения входного сигнала на несколько подканалов. В рассматриваемом случае банк цифровых фильтров — совокупность однотипных полосовых фильтров, перекрывающих весь исследуемый частотный диапазон.
Пусть исследуемая полоса:
где Fs — частота дискретизации входного комплексного сигнала.
Тогда центральная частота k-ого цифрового фильтра:
где K — число подканалов, равное числу фильтров. Выходные отсчеты k-ого канала (фильтра) определяются формулой [1]:
Все полосовые цифровые фильтры получены из исходного ФНЧ сдвигами его частотной характеристики (входного сигнала) (рис. 1). Такие сдвиги может обеспечить дискретное преобразование Фурье:
где K — количество отсчетов в выборке; k — номер гармоники, k = (0, K–1).
Повторяя преобразования (2) на каждом текущем отсчете, получим:
что соответствует формуле (1), когда h(i) = 1, i = (0, K–1). Теперь ДПФ (рис. 2) можно рассматривать как набор из K полосовых фильтров:
где K — k-номер цифрового фильтра (канала).
Такая частотная характеристика имеет три существенных недостатка: растекание в боковые лепестки, наложение соседних каналов и «эффект частокола», вызванный остроконечной формой главного лепестка. Попытки улучшения АЧХ стандартными окнами Ханнинга, Хемминга, Хана [2] и т. п. хотя и позволяют убрать боковые лепестки (растекание), но лишь за счет усиления эффекта наложения. Это объясняется тем, что во временной области все стандартные окна фактически сужают интервал анализа относительно исходного прямоугольного окна, что в частотной области приводит к обратному эффекту.
Вывод прост: для того чтобы частотные характеристики каналов не перекрывались, интервал, на котором происходит взвешивание сигнала, должен быть больше интервала ДПФ-анализа. Фактически, нужно сначала сформировать взвешивающим окном желаемую форму частотной характеристики, а потом проводить ДПФ.
Если снять ограничение на длину интервала взвешивания N = K и заменить его на более легкое — N = L×K, L = 2, 3, 4,…, то есть N больше, но кратно интервалу ДПФ-анализа, то подбором взвешивающего окна можно задать любую форму частотной характеристики цифрового фильтра. Это позволит обеспечить и отсутствие перекрытия соседних каналов, и максимально равномерную характеристику в полосе пропускания. Как показывают вычисления, для обеспечения перекрытия соседних каналов менее 5% при любом К длина окна N должна быть в 12–16 раз больше К.
Чтобы вернуться к выбранной длине интервала ДПФ-анализа, взвешенную последовательность длины N = L×K разбивают на L блоков по K отсчетов, после чего эти блоки накладывают друг на друга и поэлементно суммируют. Каждый r-й отсчет наложенной последовательности, полученной в момент времени t, zt (r) = zt (K–i), i = (0, K–1), определяется выражением:
где N = L×K, n — номер блока, n = (0, L–1).
Далее над полученными K отсчетами проводится ДПФ. Поэлементное сложение блоков длины K взвешенной последовательности допустимо, так как все используемые в ДПФ комплексные экспоненты укладываются в K отсчетах целое число периодов, поэтому каждый K-й отсчет умножается на одно и то же значение.
Отсчеты после ДПФ описываются выражением:
Если число каналов ДПФ K = 8 и форма АЧХ взвешивающего окна такая, как на рис. 1, то АЧХ банка цифрового фильтров будет такой, как показано на рис. 3. При этом номера каналов K–k согласно (5) соответствуют обозначениям на рис. 3.
Фактически взвешивающее окно — это импульсная характеристика КИХ цифрового фильтра желаемой формы для одного канала.
Оценим число операций типа сложение-умножение, приходящихся на один отсчет входного комплексного сигнала. Для КИХ цифрового фильтра оно равно 2K×L, для БПФ — 2K log2K. Тогда на вычисление одного выходного отсчета во всех каналах банка цифровых фильтров приходится 2(L+log2K) операций (множитель 2 учитывает, что на сложение или умножение комплексных отсчетов приходится 2 операции: обработка вещественной и мнимой частей). Поскольку при подобной фильтрации частота на выходе цифровых фильтров будет в К раз меньше частоты входного сигнала, следующее преобразование можно выполнять через К входных отсчетов. Фактически это эквивалентно децимации цифровых фильтрованного сигнала. На практике обычно имеет место перекрытие АЧХ соседних каналов. Перекрытие вызвано тем, что невозможно получить идеально прямоугольную форму АЧХ взвешивающего окна. Это означает, что частотная полоса в каждом канале будет несколько шире, чем Fs/K. Следовательно, после децимации в K раз выходной сигнал будет искажен (рис. 4) [1].
Поэтому для устранения нежелательных эффектов децимации будем проводить следующее преобразования не через К, а через К/2 входных отсчетов, таким образом создадим двукратный запас по частоте дискретизации выходного сигнала. Тогда общее число операций на один отсчет частоты дискретизации составит:
Например, при числе цифровых фильтров К = 1024 и L = 12 число операций на один отсчет равно 88, а при К = 32–68. Фактически, определяющим фактором вычислительной сложности является не ДПФ, а цифровой фильтр, от формы которого зависит L. Отметим, что при практической реализации операции над вещественными и мнимыми частями отсчетов проводятся параллельно, так же одновременно проводятся операции взвешивания следующей группы отсчетов и БПФ предыдущей при конвейерной организации алгоритма. Параллелизация и конвейеризация и позволяют значительно ускорить анализ.
Здесь следует отметить принципиальную разницу между традиционным БПФ-анализом и реализуемым банком цифровых фильтров. При обычном вычислении спектра временной интервал между выборками может выбираться произвольно и зависит только от скорости изменения спектра. Для банка цифровых фильтров базисные функции Фурье — это опорные частоты гетеродинов, транспонирующих спектр. Поэтому сдвиг во времени в исходном сигнале должен сопровождаться точно таким же сдвигом всех функций базиса. Этот сдвиг легко реализуется только для частных случаев, когда сдвиг равен К/2 или К/4.
Из (6) не следует, что число каналов ДПФ может быть сколь угодно велико. С ростом K все сложнее получить желаемую форму АЧХ. При достаточно больших K либо вообще невозможно получить взвешивающее окно, обеспечивающее желаемую форму АЧХ, либо его длина чрезмерно увеличивается за счет одновременного увеличения L и К. При практической реализации важно не только количество операций на один входной отсчет, но и длина взвешивающего окна N = L×K, так как в памяти нужно хранить как сами весовые коэффициенты окна, так и равное его длине количество комплексных входных отсчетов.
Уточняющее преобразование
Чтобы провести более детальный анализ сигнала, вместо увеличения количества каналов K можно использовать прием, который назовем «Уточняющее преобразование». Уточнение идеально подходит для случаев, когда требуется более подробно рассмотреть спектральные компоненты в одном из каналов, и нет необходимости детализировать остальную часть спектра. С вычислительной точки зрения наиболее эффективно разбиение К на две равные части, при котором К = m1×m2 = m2, где mi — число каналов ДПФ на каждом шаге. Однако это не жесткое условие, так как повторное преобразование производится уже над децимированными отсчетами и не сказывается на скорости обработки.
Уточняющее преобразование заключается в том, что сигнал выбранного канала может быть повторно подвергнут описанному выше БПФ-преобразованию с тем же или другим ФНЧ-фильтром. В этом случае анализ сигнала разбивается на два этапа: на первом этапе исследователю предъявляется амплитудный спектр (или спектральная плотность мощности) сигнала с достаточно грубым разбиением на каналы (например, на 64 канала). Выбранный по этому спектру канал вновь разбивается на составляющие. Их количество произвольно, но кратно двум. На этом этапе исследователь может наблюдать уже не только амплитудный спектр сигнала, но и любую выбранную частотную составляющую. Каждая из этих составляющих может быть выдана в двух формах: либо как предварительно гетеродинированная и отфильтрованная, либо просто в форме сигнала, пропущенного через узкополосный цифровой фильтр. Первая форма подразумевает децимацию исходного сигнала, вторая использует восстановление исходного сигнала по децимированным отсчетам.
Оценка скорости анализа сигнала
Рассмотрим реализацию банка цифровых фильтров на ПЛИС Spartan-3 фирмы Xilinx.
Для оценки реального быстродействия спектроанализатора разобьем все операции на две группы:
- операции взвешивания входных отсчетов (КИХ-фильтрация);
- операции БПФ-анализа.
Примем, как и ранее, длину выборки
где К = 2 — размерность БПФ-анализа (число спектральных коэффициентов), L = 10ч16 — коэффициент расширения интервала К, выбираемый исходя из требований к крутизне КИХ-фильтра.
Оценим временные затраты на операции первой группы.
Минимальную скорость обработки можно получить, если при взвешивании входных отсчетов КИХ-коэффициентами использовать всего два множительных устройства: одно для вещественной, другое — для мнимой части. В этом случае взвешивание всех N отсчетов займет K×L тактов. Поскольку повторение вычислений происходит со сдвигом (децимацией) на половину длины БПФ-анализа, то на один отсчет входного сигнала будет приходиться (K×L)/K/2 = 2L операций. Например, при L = 12 на один отсчет входного сигнала должно приходиться 24 такта тактовой частоты ПЛИС. Если эта частота равна 25 МГц, частота дискретизации сигнала может составлять 1 МГц.
Поднять предельную частоту анализа можно распараллеливанием операций умножения. SPARTAN-3 XC3S200 имеет 12 встроенных умножителей. Шесть из них заняты процедурой БПФ, остальные могут быть использованы для взвешивания входных отсчетов. Таким образом, на данной микросхеме увеличить быстродействие можно только в три раза. При этом частота дискретизации сигнала составит 2,5 МГц на тактовой частоте ПЛИС 25 МГц, и 6 МГц — на 70 МГц.
Существенного увеличения быстродействия можно добиться только при использовании ПЛИС нового поколения (Virtex 4), ориентированных на цифровую обработку сигналов, содержащих не менее 1 млн системных вентилей и более сотни умножителей. В стандартной библиотеке таких ПЛИС имеются оптимизированные алгоритмы БПФ, осуществляющие обработку за один такт.
При тактовой частоте порядка 400 МГц такая ПЛИС может обеспечить обработку сигнала с частотой дискретизации около 100 МГц.
Временные затраты БПФ-анализа
Алгоритмы БПФ, предлагаемые Xilinx, можно разделить на две группы. Первая группа позволяет за счет распараллеливания вычислений производить обработку таким образом, что за время загрузки очередных N отсчетов происходит выгрузка такого же числа отсчетов, рассчитанных по предыдущей выборке. Иными словами, БПФ вносит только задержку выдачи результата, но не влияет на скорость обработки.
Вторая группа алгоритмов не требует распараллеливания вычислений (а значит, и больших аппаратных затрат), использует порядка 6–10 множительных устройств, а операции загрузки, БПФ и выгрузки производятся последовательно. В этом случае количество тактов, требуемое для БПФ, составляет примерно (2ч3)К. Поскольку это число в три-четыре раза меньше времени умножения на КИХ-коэффициенты, БПФ-устройство даже при такой организации вычислений большую часть времени будет находиться в режиме ожидания.
Основной недостаток выбранного алгоритма БПФ заключается в сокращении разрядности выходных отсчетов. Стандартная процедура БПФ для исключения переполнений после каждой итерации сдвигает данные на один разряд, что приводит к потере log2K разрядов в выходном сигнале. При K = 1024 от исходной сетки 216 останется всего 6 разрядов. Такое масштабирование оправданно лишь при взвешивании исходного сигнала прямоугольным весовым окном. Когда используется взвешивание расширенным весовым окном, такое масштабирование становится избыточным. В разработанном проекте масштабирование проводилось лишь на нескольких промежуточных итерациях, благодаря чему динамический диапазон выходного сигнала составляет 212.
Библиотека стандартных алгоритмов предлагает и более оптимальное решение, когда на каждой итерации проводится коррекция коэффициентов в зависимости от значения старшего разряда результата. Для экономии ресурсов ПЛИС этот алгоритм не использовался.
Оценка требуемой логической емкости при реализации банка цифровых фильтров
В ПЛИС Spartan-3 -200 удалось реализовать максимальное число каналов банка цифровых фильтров, равное 128. Данные об использованных ресурсах ПЛИС приведены в таблице 1.
Попытка реализации в использованной ПЛИС банка цифровых фильтров на 256 каналов привела к загрузке ресурсов, показанной в таблице 2. Вместить проект в выбранную ПЛИС удалось только c исключением операции вычисления мнимой части выходных отсчетов.
Оценка параметров реализованного банка цифровых фильтров
Анализ искажения выходных сигналов во временной области
Оценка искажения мгновенных значений сигнала, связанная с нелинейностями преобразований, округлением, шумами и эффектом проникновения сигналов соседних каналов, может быть произведена по рис. 5, 6.
В каналах, в которых отсутствует сигнал в полосе пропускания, уровень выходного сигнала по СКО составляет 0,0008, что меньше –60 дБ.
В полосе пропускания синусоидальный гармонический сигнал имеет среднее значение –0,001, а СКО соответствует СКО гармонического сигнала с амплитудой 0,262 с погрешностью менее 1%.
Мгновенные значения выходного гармонического сигнала отличаются от тестового сигнала на величины, не превышающие 3%.
Анализ времени переходного процесса
На рис. 7 представлен график переходного процесса, полученный для нулевого канала. По времени группового запаздывания и характеру переходного процесса график в точности соответствует переходному процессу эталонного КИХ-фильтра.
Время реакции на скачкообразное изменение частоты гармонического сигнала показывают рис. 8, 9. Это время составляет не более 13 отсчетов.
Анализ скорости работы банка цифровых фильтров
Тестирование скорости работы банка цифровых фильтров проводилось при помощи программы, которая осуществляет передачу тестового сигнала в виде массива размером в 320 Мбайт, что составляет 83 886 080 комплексных отсчетов. Полученное время передачи составило около 91 с, что соответствует скорости обработки в 921 825 комплексных отсчетов в секунду.
Анализ АЧХ банка цифровых фильтров
Частотная характеристика первого канала представлена в линейном масштабе на рис. 10, а в логарифмическом масштабе — на рис. 12, красным цветом наложена эталонная характеристика, синтезированная в Matlab. На рис. 11 в линейном масштабе показано отличие синтезированной ЧХ от реальной (синяя кривая). Это отличие составляет менее 0,1 дБ. Как видно на рис. 12, в полосе прозрачности эти характеристики практически совпадают, а в полосе подавления проявляются эффекты округления, что видно в увеличенном масштабе на рис. 13. Наложение частотных характеристик смежных каналов иллюстрируют рис. 14 (линейный масштаб) и рис. 15 (логарифмический масштаб).
Динамический диапазон на выходе банка цифровых фильтров
При подаче на вход гармонического сигнала с максимальной амплитудой 32 767 выходная амплитуда составила 8690. Уменьшение динамического диапазона составило 12 дБ, что видно и по графикам частотных характеристик.
Литература
- www.rfel.com
- Солонина А. И., Улахович Д. А., Арбузов С. М., Соловьева Е. Б., Гук И. И. Основы цифровой обработки сигналов: Курс лекций. СПб.: БХВ-Петербург, 2003.
- Рабинер Л., Гоулд Б. Теория и применение цифровой обработки сигналов. М.: Мир, 1978.