Банк цифровых фильтров Текст научной статьи по специальности «Электротехника, электронная техника, информационные технологии»

Леонид АЗАРЕНКОВ Иван КАНАТОВ, к. т. н. Дмитрий КАПЛУН

mitya_kapl@front.ru

Банк цифровых фильтров

В операциях цифровой обработки сигналов особое внимание уделяется цифровой фильтрации, которая по объему вычислений в среднем занимает от 20 до 60%. Для того чтобы свести к минимуму возможные потери информации и повысить качество ее обработки, цифровые фильтры должны обеспечивать возможность быстрой работы с большими блоками данных. Одним из вариантов решения этой задачи является использование банка цифровых фильтров.

В узком смысле цифровой фильтр — это частотно-избирательная цепь, обеспечивающая селекцию цифровых сигналов по частоте [2]. После выполнения цифровой фильтрации мы, как правило, получаем интересующий нас сигнал, то есть сигнал, несущий нужную нам информацию в виде, удобном для последующей обработки. Соответственно, к параметрам цифровых фильтров в современных системах цифровой обработки сигналов начинают предъявляться повышенные требования. Частоты, на которых работают фильтры, нередко достигают нескольких сотен мегагерц и более. Соответственно, растет ширина полос фильтров. Это ведет к увеличению объема и скорости вычислений, а значит, и к резкому росту аппаратных затрат. Для того чтобы свести к минимуму возможные потери информации и повысить качество ее обработки, цифровые фильтры должны обеспечивать возможность быстрой работы с большими блоками данных. Одним из вариантов решения этой задачи является использование банка цифровых фильтров.

Структура банка фильтров

Банк цифровых фильтров предназначен для разбиения входного сигнала на несколько подканалов. В рассматриваемом случае банк фильтров — совокупность однотипных полосовых фильтров, перекрывающих весь исследуемый частотный диапазон.

Пусть исследуемая полоса:

/е [-Л/2, Л/2М0, Л],

где Рв — частота дискретизации входного комплексного сигнала.

Тогда центральная частота к-ого фильтра:

А = (кхРв)/К = кх/р

где К — число подканалов, равное числу фильтров. Выходные отсчеты к-ого канала (фильтра) определяются формулой [1]:

лг-1

Ук(і) = £ ВДх(*-і)ехр(;2я/*х*) =

1=0

N-1

=^іІі(і)х(і-ї)ехр(]2п/1кхҐ). (1)

(=0

Все полосовые фильтры получены из исходного ФНЧ сдвигами его частотной характеристики (входного сигнала) (рис. 1). Такие сдвиги может обеспечить дискретное преобразование Фурье:

Х(к) = ^ х(і) ехр(- /), (2)

і=0 Л

где К — количество отсчетов в выборке; к — номер гармоники, к = (0, К-1).

Повторяя преобразования (2) на каждом текущем отсчете, получим:

Х(к, і) = ^ х(і-і) ехр(- ^ і), (3)

(=0 К

что соответствует формуле (1), когда к(і) = 1, і = (0, К-1). Теперь ДПФ (рис. 2) можно рассматривать как набор из К полосовых фильтров:

Ук-кМ = х(к, г),

где К — к-номер фильтра (канала).

Такая частотная характеристика имеет три существенных недостатка: растекание в боко-

Номер фильтра, к

Рис. 2. Эффекты наложения и растекания в ДПФ

А(0

[ X X X л,

1 2М И1 1М

Рис. 1. АЧХ банка фильтров

гл гл гл гл гл гл гл г ■\ п п п п г АПППП Г 1ППППППП"

/ о’/1\/ 2\/ 3\/ 4\/ 5\/ 6\/ 7\/ < А /Д/з\j2\j 1\/< /« )У1У2УзУ4У5УбУД;р

' 0 і 0 і 0

-Ре -?$/2 Ре/2 р5

Рис. 3. Эквивалентные распределения каналов

вые лепестки, наложение соседних каналов и «эффект частокола», вызванный остроконечной формой главного лепестка. Попытки улучшения АЧХ стандартными окнами Хан-нинга, Хемминга, Хана [2] и т. п. хотя и позволяют убрать боковые лепестки (растекание), но лишь за счет усиления эффекта наложения. Это объясняется тем, что во временной области все стандартные окна фактически сужают интервал анализа относительно исходного прямоугольного окна, что в частотной области приводит к обратному эффекту.

Вывод прост: для того чтобы частотные характеристики каналов не перекрывались, интервал, на котором происходит взвешивание сигнала, должен быть больше интервала ДПФ-анализа. Фактически, нужно сначала сформировать взвешивающим окном желаемую форму частотной характеристики, а потом проводить ДПФ.

Если снять ограничение на длину интервала взвешивания N = К и заменить его на более легкое — N = 1хК, Ь = 2, 3, 4,..., то есть N больше, но кратно интервалу ДПФ-анали-за, то подбором взвешивающего окна можно задать любую форму частотной характеристики фильтра. Это позволит обеспечить и отсутствие перекрытия соседних каналов, и максимально равномерную характеристику в полосе пропускания. Как показывают вычисления, для обеспечения перекрытия соседних каналов менее 5% при любом К длина окна N должна быть в 12-16 раз больше К.

Чтобы вернуться к выбранной длине интервала ДПФ-анализа, взвешенную последовательность длины N = 1хК разбивают на I блоков по К отсчетов, после чего эти блоки накладывают друг на друга и поэлементно суммируют. Каждый г-й отсчет наложенной последовательности, полученной в момент времени % (г) = (К-1), I = (0, К-1), определяется выражением:

Ф) = 2^(К-1) =

1-1

= ^1х^-пК-г)Н(пК+г); (4)

п=0

где N = 1хК, п — номер блока, п = (0,1-1).

Далее над полученными К отсчетами проводится ДПФ. Поэлементное сложение блоков длины К взвешенной последовательности допустимо, так как все используемые в ДПФ комплексные экспоненты уклады-

ваются в К отсчетах целое число периодов, поэтому каждый К-й отсчет умножается на одно и то же значение.

Отсчеты после ДПФ описываются выражением:

к-1

у,(К-к) = £ г{0ехр(- ]2п/кх{);

1=0

Ук = Уо- (5)

Если число каналов ДПФ К =8 и форма АЧХ взвешивающего окна такая, как на рис. 1, то АЧХ банка фильтров будет такой, как показано на рис. 3. При этом номера каналов К-к согласно (5) соответствуют обозначениям на рис. 3.

Фактически взвешивающее окно — это импульсная характеристика КИХ фильтра желаемой формы для одного канала.

Оценим число операций типа сложение-умножение, приходящихся на один отсчет входного комплексного сигнала. Для КИХ фильтра оно равно 2Кх1, для БПФ — 2Кlog2K. Тогда на вычисление одного выходного отсчета во всех каналах банка фильтров приходится 2(1+^2К) операций (множитель 2 учитывает, что на сложение или умножение комплексных отсчетов приходится 2 операции: обработка вещественной и мнимой частей). Поскольку при подобной фильтрации частота на выходе фильтров будет в К раз меньше частоты входного сигнала, следующее преобразование можно выполнять через К входных отсчетов. Фактически это эквивалентно децимации фильтрованного сигнала. На практике обычно имеет место перекрытие АЧХ соседних каналов. Перекрытие вызвано тем, что невозможно получить идеально прямоугольную форму АЧХ взвешивающего окна.

Рис. 4. Иллюстрация эффекта наложения при децимации: а) дискретный спектр исходного сигнала, б) дискретный спектр сигнала после децимации в 2 раза

Это означает, что частотная полоса в каждом канале будет несколько шире, чем Рх/К. Следовательно, после децимации в К раз выходной сигнал будет искажен (рис. 4) [1].

Поэтому для устранения нежелательных эффектов децимации будем проводить следующее преобразования не через К, а через К/2 входных отсчетов, таким образом создадим двукратный запас по частоте дискретизации выходного сигнала. Тогда общее число операций на один отсчет частоты дискретизации составит:

е2К*1+2К1оё2К) л/т .

^------- 62 ; = 4(1+ logj.fi). (6)

Например, при числе фильтров К = 1024 и I = 12 число операций на один отсчет равно 88, а при К = 32-68. Фактически, определяющим фактором вычислительной сложности является не ДПФ, а цифровой фильтр, от формы которого зависит I. Отметим, что при практической реализации операции над вещественными и мнимыми частями отсчетов проводятся параллельно, так же одновременно проводятся операции взвешивания следующей группы отсчетов и БПФ предыдущей при конвейерной организации алгоритма. Параллелизация и конвейеризация и позволяют значительно ускорить анализ.

Здесь следует отметить принципиальную разницу между традиционным БПФ-ана-лизом и реализуемым банком фильтров. При обычном вычислении спектра временной интервал между выборками может выбираться произвольно и зависит только от скорости изменения спектра. Для банка фильтров базисные функции Фурье — это опорные частоты гетеродинов, транспонирующих спектр. Поэтому сдвиг во времени в исходном сигнале должен сопровождаться точно таким же сдвигом всех функций базиса. Этот сдвиг легко реализуется только для частных случаев, когда сдвиг равен К/2 или К/4.

Из (6) не следует, что число каналов ДПФ может быть сколь угодно велико. С ростом К все сложнее получить желаемую форму АЧХ. При достаточно больших К либо вообще невозможно получить взвешивающее окно, обеспечивающее желаемую форму АЧХ, либо его длина чрезмерно увеличивается за счет одновременного увеличения I и К. При практической реализации важно не только количество операций на один входной отсчет,

но и длина взвешивающего окна N = 1хК, так как в памяти нужно хранить как сами весовые коэффициенты окна, так и равное его длине количество комплексных входных отсчетов.

Уточняющее преобразование

Чтобы провести более детальный анализ сигнала, вместо увеличения количества каналов К можно использовать прием, который назовем «Уточняющее преобразование». Уточнение идеально подходит для случаев, когда требуется более подробно рассмотреть спектральные компоненты в одном из каналов, и нет необходимости детализировать остальную часть спектра. С вычислительной точки зрения наиболее эффективно разбиение К на две равные части, при котором К = т1 х т2 = т2, где т{ — число каналов ДПФ на каждом шаге. Однако это не жесткое условие, так как повторное преобразование производится уже над децимированными отсчетами и не сказывается на скорости обработки.

Уточняющее преобразование заключается в том, что сигнал выбранного канала может быть повторно подвергнут описанному выше БПФ-преобразованию с тем же или другим ФНЧ-фильтром. В этом случае анализ сигнала разбивается на два этапа: на первом этапе исследователю предъявляется амплитудный спектр (или спектральная плотность мощности) сигнала с достаточно грубым разбиением на каналы (например, на 64 канала). Выбранный по этому спектру канал вновь разбивается на составляющие. Их количество произвольно, но кратно двум. На этом этапе исследователь может наблюдать уже не только амплитудный спектр сигнала, но и любую выбранную частотную составляющую. Каждая из этих составляющих может быть выдана в двух формах: либо как предварительно гетеродинированная и отфильтрованная, либо просто в форме сигнала, пропущенного через узкополосный фильтр. Первая форма подразумевает децимацию исходного сигнала, вторая использует восстановление исходного сигнала по децимированным отсчетам.

Оценка скорости анализа сигнала

Рассмотрим реализацию банка фильтров на ПЛИС 8раЛап-3 фирмы Хіііпх.

Для оценки реального быстродействия спектроанализатора разобьем все операции на две группы:

• операции взвешивания входных отсчетов

(КИХ-фильтрация);

• операции БПФ-анализа.

Примем, как и ранее, длину выборки

N = ЬхК,

где К =2 — размерность БПФ-анализа (число спектральных коэффициентов), I = 10^16 —

коэффициент расширения интервала К, выбираемый исходя из требований к крутизне КИХ-фильтра.

Оценим временные затраты на операции первой группы.

Минимальную скорость обработки можно получить, если при взвешивании входных отсчетов КИХ-коэффициентами использовать всего два множительных устройства: одно для вещественной, другое — для мнимой части. В этом случае взвешивание всех N отсчетов займет Кх1 тактов. Поскольку повторение вычислений происходит со сдвигом (децимацией) на половину длины БПФ-анализа, то на один отсчет входного сигнала будет приходиться (Кх1)/К/2 = 21 операций. Например, при I = 12 на один отсчет входного сигнала должно приходиться

24 такта тактовой частоты ПЛИС. Если эта частота равна 25 МГц, частота дискретизации сигнала может составлять 1 МГц.

Поднять предельную частоту анализа можно распараллеливанием операций умножения. 8РАИТА^3 ХС38200 имеет 12 встроенных умножителей. Шесть из них заняты процедурой БПФ, остальные могут быть использованы для взвешивания входных отсчетов. Таким образом, на данной микросхеме увеличить быстродействие можно только в три раза. При этом частота дискретизации сигнала составит 2,5 МГц на тактовой частоте ПЛИС

25 МГц, и 6 МГц — на 70 МГц.

Существенного увеличения быстродействия можно добиться только при использовании ПЛИС нового поколения (Уй1ех 4), ориентированных на цифровую обработку сигналов, содержащих не менее 1 млн системных вентилей и более сотни умножителей. В стандартной библиотеке таких ПЛИС имеются оптимизированные алгоритмы БПФ, осуществляющие обработку за один такт.

При тактовой частоте порядка 400 МГц такая ПЛИС может обеспечить обработку сигнала с частотой дискретизации около 100 МГц.

Временные затраты БПФ-анализа

Алгоритмы БПФ, предлагаемые ХШих, можно разделить на две группы. Первая группа позволяет за счет распараллеливания вычислений производить обработку таким образом, что за время загрузки очередных N отсчетов происходит выгрузка такого же числа отсчетов, рассчитанных по предыдущей выборке. Иными словами, БПФ вносит только задержку выдачи результата, но не влияет на скорость обработки.

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

Таблица 1. Распределение ресурсов ПЛИС для банка фильтров на 128 каналов

Логическое использование Занятых Доступных Занятость, %

Общее количество регистров 1932 3840 50

Количество регистров, использованных в качестве flip flop регистров 1927

Количество регистров, использованных в качестве защелок 5

Количество регистров для 4 входных LUT-блоков 1686 3840 43

Логическое распределение

Количество занятых вентилей 1588 1920 82

Количество вентилей, включающих только связную логику 1588 1588 100

Количество вентилей включающих несвязную логику 0 1588 0

Общее количество для 4 входных LUT-блоков 2231 3840 58

Количество блоков, использованных как логические 1686

Количество блоков, использованных как трассировочные 204

Количество блоков, использованных как сдвигающие регистры 341

Количество граничных блоков ввода/вывода 31 173 17

Количество блоков памяти ОЗУ 11 12 91

Количество умножителей 18x18 8 12 66

Количество блоков глобальной синхронизации 2 8 25

Таблица 2. Распределение ресурсов ПЛИС для банка фильтров на 256 каналов

Логическое использование Занятых Доступных Занятость, %

Количество регистров, использованных в качестве flip flop регистров 1882 3840 49

Количество регистров для 4 входных LUT-блоков 2013 3840 52

Логическое распределение

Количество занятых вентилей 1880 1920 97

Количество вентилей, включающих только связную логику 1880 1880 100

Количество вентилей включающих несвязную логику 0 1880 0

Общее количество для 4 входных LUT-блоков 3248 3840 84

Количество блоков, использованных как логические 2013

Количество блоков, использованных как трассировочные 130

Количество блоков, использованных как сдвигающие регистры 1105

Количество граничных блоков ввода/вывода 31 173 17

Количество блоков памяти ОЗУ 12 12 100

Количество умножителей 18x18: 8 12 66

Количество блоков глобальной синхронизации 2 8 25

Рис. 5. Сигнал на выходе банка фильтров

при подаче тестового гармонического сигнала с частотой 5 Гц

-

■ ^ , .

0 5000 10000 15000

- -

0 5000 10000 15000

■

. . ^ ■

0 5000 10000 15000

■ “

. _

5000 сИаппеїз і гот 1 іо 4 10000

Рис. 6. Сигнал на выходе банка фильтров при подаче тестового сигнала с линейно нарастающей частотой от 0 Гц до 0,2xFs (при Fs = 1 МГц)

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

Основной недостаток выбранного алгоритма БПФ заключается в сокращении разрядности выходных отсчетов. Стандартная процедура БПФ для исключения переполнений после каждой итерации сдвигает данные на один разряд, что приводит к потере ^2К разрядов в выходном сигнале. При К = 1024 от исходной сетки 216 останется всего 6 разрядов. Такое масштабирование оправданно лишь при взвешивании исходного сигнала прямоугольным весовым окном. Когда используется взвешивание расширенным весовым окном, такое масштабирование становится избыточным. В разработанном проекте масштабирование проводилось лишь на нескольких промежуточных итерациях, благодаря чему динамический диапазон выходного сигнала составляет 212.

Библиотека стандартных алгоритмов предлагает и более оптимальное решение, когда на каждой итерации проводится коррекция коэффициентов в зависимости от значения стар-

шего разряда результата. Для экономии ресурсов ПЛИС этот алгоритм не использовался.

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

В ПЛИС ЗраГап-З -200 удалось реализовать максимальное число каналов банка фильтров, равное 128. Данные об использованных ресурсах ПЛИС приведены в таблице 1.

Попытка реализации в использованной ПЛИС банка фильтров на 256 каналов привела к загрузке ресурсов, показанной в таблице 2. Вместить проект в выбранную ПЛИС удалось только с исключением операции вычисления мнимой части выходных отсчетов.

Оценка параметров реализованного банка фильтров

Анализ искажения выходных сигналов во временной области

Оценка искажения мгновенных значений сигнала, связанная с нелинейностями преобразований, округлением, шумами и эффектом проникновения сигналов соседних каналов, может быть произведена по рис. 5, 6.

В каналах, в которых отсутствует сигнал в полосе пропускания, уровень выходного сигнала по СКО составляет 0,0008, что меньше -60 дБ.

В полосе пропускания синусоидальный гармонический сигнал имеет среднее значение -0,001, а СКО соответствует СКО гармонического сигнала с амплитудой 0,262 с погрешностью менее 1%.

Мгновенные значения выходного гармонического сигнала отличаются от тестового сигнала на величины, не превышающие 3%.

Анализ времени переходного процесса На рис. 7 представлен график переходного процесса, полученный для нулевого канала. По времени группового запаздывания и характеру переходного процесса график в точности соответствует переходному процессу эталонного КИХ-фильтра.

Время реакции на скачкообразное изменение частоты гармонического сигнала показывают рис. 8, 9. Это время составляет не более 13 отсчетов.

Анализ скорости работы банка фильтров Тестирование скорости работы банка фильтров проводилось при помощи программы, которая осуществляет передачу тестового сигнала в виде массива размером в 320 Мбайт, что составляет 83 886 080 комплексных отсчетов. Полученное время передачи составило около 91 с, что соответствует скорости обработки в 921 825 комплексных отсчетов в секунду.

Анализ АЧХ банка фильтров

Частотная характеристика первого канала представлена в линейном масштабе на рис. 10, а в логарифмическом масштабе — на рис. 12, красным цветом наложена эталонная характеристика, синтезированная в Ма11аЬ. На рис. 11 в линейном масштабе показано отличие синтезированной ЧХ от реальной (синяя кривая). Это отличие составляет менее 0,1 дБ. Как видно на рис. 12, в полосе прозрачности эти характеристики практически совпадают, а в полосе подавления проявляются эффекты округления, что видно в увеличенном масштабе на рис. 13.

Рис. 7. Переходный процесс в нулевом канале банка фильтров

°о \ /\ /і 1 ■

-0,2 V \/ . .

500 1000 1500 2000 2500

0,2 ' . А 'А ■

-0,2 ■ 1/\А .

500 1000 1500 2000 2500

0,2 ■

-0,2 ■ . \/ V/. .

500 1000 1500 2000 2500

0,2 —,А А,- ■

-0,2 _ У.УГ .

500 1000 1500 channels from 1 to 4 2000 2500

Рис. 8. Сигнал на выходе банка фильтров (каналы 1—4)

при подаче тестового гармонического сигнала скачкообразно изменяющего частоту каждые 2 с с шагом Fs/K начиная с частоты 1 Гц (при Fs = 10000, К = 128)

0,2 " J\ /\ 7L_ ■

-0,2 _ _

500 1000 1500 2000 2500

0,2 -

-0,2 ■ . . . іЛ/\і . -

500 1000 1500 2000 2500

0,2 -

-0,2 ■ . . . \J.\J -

500 1000 1500 2000 2500

0,2 -

-0,2 ■ . . . . V\J -

500 1000 1500 channels from 5 to 8 2000 2500

Рис. 9. Сигнал на выходе банка фильтров (каналы 5—8) при подаче

тестового гармонического сигнала, скачкообразно изменяющего частоту каждые 2 с

с шагом Fs/K начиная с частоты 1 Гц (при Fs = 10000, К = 128)

Рис. 10. Сравнение практически полученной АЧХ первого канала банка фильтров и теоретически рассчитанной АЧХ (с учетом сдвига и масштабирования по оси у)

Рис. 11. Сравнение практически полученной АЧХ первого канала банка фильтров и теоретически рассчитанной АЧХ (с учетом сдвига и масштабирования по оси у), увеличенный фрагмент

Рис. 12. Сравнение практически полученной АЧХ первого канала банка фильтров и теоретически рассчитанной АЧХ (с учетом сдвига и масштабирования по оси у), полулогарифмический масштаб

Рис. 13. Сравнение практически полученной АЧХ первого канала банка фильтров и теоретически рассчитанной АЧХ (с учетом сдвига и масштабирования по оси у), полулогарифмический масштаб, увеличенный фрагмент

Рис. 14. АЧХ банка фильтров — 3 соседних канала (1,2и3)

Overlapped MFR for channels 0,1 and 2

Рис. 15. АЧХ банка фильтров в полулогарифмическом масштабе — 3 соседних канала (1, 2 и3)

Рис. 16. Сигнал на выходе банка фильтров (каналы 33—36). Тестовый сигнал — синусоиды с частотами 25050, 25 250, 24750 и 50 100 Гц (при Fs = 100 кГц)

Рис. 18. Сигнал на выходе банка фильтров (каналы 33—36).

Зашумленный тестовый сигнал — синусоиды с частотами 25 050, 25 250, 24750 и 50 100 Гц (при Fs = 100 кГц и соотношении «сигнал-шум» 0,622)

Наложение частотных характеристик смежных каналов иллюстрируют рис. 14 (линейный масштаб) ирис. 15 (логарифмический масштаб).

Динамический диапазон на выходе банка фильтров

При подаче на вход гармонического сигнала с максимальной амплитудой 32 767 выходная амплитуда составила 8690. Уменьшение динамического диапазона составило 12 дБ, что видно и по графикам частотных характеристик. ■

Литература

1. www.rfel.com

2. Солонина А. И., Улахович Д. А., Арбузов С. М., Соловьева Е. Б., Гук И. И. Основы цифровой обработки сигналов: Курс лекций. СПб.: БХВ-Петер-бург, 2003.

3. Рабинер Л., Гоулд Б. Теория и применение цифровой обработки сигналов. М.: Мир, 1978.