Всем привет! Написать эту статью меня побудило выступление на семинарах по цифровой обработке сигналов, где слушатели всегда заостряли интерес к методике вычисления корректирующих FIR-фильтров, несмотря на то, что эту тему я затрагивал поверхностно и по большей части рассказывал об этом в ознакомительных чертах. Если публика желает получить тайные знания, то почему бы ими не поделиться. В этой статье я постараюсь в доступной форме изложить алгоритм расчета корректирующих КИХ фильтров, который необходим для выравнивания АЧХ в полосе пропускания после звеньев CIC фильтров в задачах децимации и интерполяции сигналов. В частности, рассмотрим проектирование фильтров на современных ПЛИС Xilinx. Как обычно, в конце статьи будет ссылка на полезные скрипты для расчета различных фильтров и получение файла коэффициентов фильтра-корректора.

Предполагается, что читатель знаком с основами цифровой обработки сигналов и имеет представление о CIC и FIR фильтрах. Приступим.

Введение

Как известно, CIC фильтры часто используются в задачах децимации и интерполяции сигналов, поскольку обладают одним очень важным свойством — простота реализации, что делает этот класс фильтров дешевыми в исполнении. Особенность CIC-фильтров — отсутствие операций умножения для вычисления отклика на выходе фильтра. CIC фильтры применяются повсеместно, где требуется работа на нескольких или разных скоростях, то есть в тех задачах, где требуется децимация и интерполяция потока данных. Напомню, что под децимацией понимается процесс понижения частоты дискретизации сигнала (уменьшение скорости потока данных), а интерполяция — это процесс повышения частоты дискретизации сигнала (увеличение скорости потока данных). CIC фильтр — неотъемлемая часть узлов DDC (Digital Down Converter) и DUC (Digital Up Converter).


АЧХ CIC фильтра эквивалентна частотной характеристике фильтра нижних частот (ФНЧ). На форму АЧХ и положение нулей характеристики влияют некоторые параметры фильтра: порядок- N, коэффициент децимации/интерполяции — R, величина задержки в дифференцирующем звене — M. И именно во влиянии этих параметров на ЧХ кроется главный недостаток CIC фильтров. При больших порядках фильтра N и значениях коэффициента децимации/интерполяции R существенно ухудшается результирующий спектр в полосе пропускания, а именно искажается форма главного лепестка АЧХ. Для многих задач многочастотной и широкополосной обработки сигналов необходимо иметь максимально равномерный и прямоугольный спектр в полосе пропускания. Но если взглянуть на график АЧХ CIC фильтра, то мы видим резко убывающую кривую в полосе пропускания. Такая неравномерность недопустима и приводит к потере энергии полезного сигнала и искажении его формы. В связи с этим актуальна задача по выравниванию АЧХ CIC фильтров в полосе пропускания при минимизации помех в полосе ослабления.

Постановка задачи

Пусть есть CIC фильтр, для которого заданы порядок N и коэффициент изменения частоты дискретизации R. Величина задержки, если это не фильтр скользящего среднего, на практике M = 1 или 2 и определяется регистрами в узлах умножителей DSP48 ПЛИС (независимо от вендора — Altera или Xilinx). Необходимо рассчитать корректирующий фильтр для выравнивания АЧХ после прохождения сигнала через звено CIC фильтра.

Решение

Корректирующий фильтр легче всего сделать на базе FIR фильтра. Для расчета этого фильтра необходимо заранее определить его параметры:

• порядок фильтра или длина импульсной характеристики — NFIR,
• эффективная полоса, в пределах которой требуется выравнивание — Fr,
• разрядность коэффициентов фильтра — Bc,
• наличие оконной фильтрации — WIN.

Порядок — NFIR. Этот параметр определяет «качество» частотной характеристики фильтра. Разработчик ЦОС на ПЛИС постоянно сталкивается с задачей выбора оптимального порядка фильтра, ведь чем больше порядок — тем лучше с точки зрения частотных свойств, но тем больше ресурсов кристалла ПЛИС необходимо потратить, чтобы этот фильтр реализовать. Оптимальные значения порядка фильтра-корректора находятся в диапазоне от 32 до 128. В некоторых задачах значение порядка фильтра может быть выше и достигать цифры 256, но на практике это не имеет смысла, и я таких длинных фильтров-корректоров не встречал. Корректирующий фильтр должен быть небольшим и не затрагивать вычислительные ресурсы, которые нужны для реализации более сложных вещей. Порядок фильтра для данной реализации лучше задавать четным числом и кратным степени двойки, если речь идёт о применении фильтра на базе ПЛИС, но это требование необязательное.

Под эффективной полосой пропускания Fr я понимаю значение относительной нормализованной частоты среза после децимации или интерполяции. Для задач децимации вычисляется по формуле Fr = Fo/R, где Fo — параметр нормированной частоты (например, от 0.1 до 0.5), а R — коэффициент децимации. Значение Fo определяет полосу пропускания фильтра. Следует отметить, что при Fr = 0.5, порядок фильтра должен быть нечетным, чтобы ноль АЧХ фильтра-корректора не попал в полосу полезного сигнала. Для других значений Fr порядок фильтра может быть любым. Это связано с особенностями вычисления ИХ встроенной функции FIR2 в САПР Matlab.

Разрядность коэффициентов для описываемого алгоритма определяет неравномерность АЧХ в полосе пропускания и полосе заграждения. На практике разрядность коэффициентов задается от 16 и до 27 бит, что связано с особенностями узла DSP48 современных FPGA. В практических целях достаточно задать разрядность коэффициентов равной 16-18 бит и это в должной мере обеспечит требуемую равномерность АЧХ. Для фильтров высоких порядков (например, N = 256) разрядности в 18 бит недостаточно и начинают проявляться эффекты квантования, в частности в полосе заграждения. Поэтому для фильтров больших порядков необходимо повысить разрядность коэффициентов, благо структура узлов обработки DSP48 в ПЛИС это позволяет (см. рис. узла DSP48).



Оконная фильтрация позволяет сгладить пульсации результирующей АЧХ в полосе пропускания после применения фильтра-корректора. Оконная функция накладывается на импульсную характеристику рассчитываемого FIR фильтра (математически это свёртка спектров). Исходя из личного опыта, отмечу, что функция Кайзера — это наилучший вариант оконной функции. С помощью одного параметра BETA можно менять частотную характеристику фильтра под требуемую задачу. Функция Кайзера рассчитывается через модифицированные функции Бесселя нулевого порядка I0, но во многих стандартных пакетах является встроенной (Matlab, GNU Octave, MathCAD). Физически параметр BETA определяет долю энергии, сосредоточенную в пределах главного лепестка спектра. Чем больше параметр BETA, тем больше сосредоточено энергии и тем шире главный лепесток, но меньше уровень боковых лепестков. В практических целях параметр BETA = 3-11.

Подробную методику расчета FIR фильтра я приводить не буду, об этом можно почитать в моей предыдущей статье. Сложного ничего нет — нужно любым удобным способом получить коэффициенты FIR фильтра. Сделать это можно с помощью MathCAD, GNU Octave, Matlab (утилита FDATool), приложения ScopeFIR, программы LabView, либо написать свою методику расчета на базе известных алгоритмов.

Алгоритм

Опишем алгоритм расчета корректирующего КИХ фильтра. Ниже будут использованы вставки кода на языке скриптов Matlab для лучшего понимания процесса реализации.

1 шаг: задаем исходные параметры CIC-фильтра —

R = 8; % Decimation factor
N = 4; % Number of stages
M = 1; % Differential delay (only 1)

2 шаг: задаем параметры FIR-фильтра (разрядность коэффициентов, порядок, окно) и нормированную частоту среза АЧХ —

NFIR  = 128; % Filter order, must be odd when Fo = 0.5 !
Bc    = 16;  % Coef. Bit-width
Fo    = 0.3; % Normalized Cutoff: 0.2 < Fo < 0.5;
BETA  = 8;   % BETA parameter for Kaiser window (if IS_WIND = 'Y') 

3 шаг: Выбираем шаг «дискретизации» векторов для расчета характеристик. В соответствии с этим шагом создаем требуемые массивы чисел. Рассчитываем CIC-фильтр по известной формуле и переводим результат в децибелы.

HCIC = (R^-N*abs(1*M*sin(pi*M*R*ff) ./ sin(pi*ff)).^N);
HCICdb = 20 * log10(abs(HCIC));

4 шаг: В соответствии с параметром эффективной полосы пропускания Fo, делим вектор отсчетов частоты на две части: fp — вектор частот в полосе пропускания, fs — вектор частот в полосе ослабления. Производим расчет идеальной характеристики компенсирующего фильтра с частотой среза Fo/R по следующей формуле:


где f – отсчеты нормированной частоты в диапазоне полосы пропускания от 0 до Fo/R. Остальные отсчеты равны нулю.

% Calculate ideal response
Mp = ones(1, length(fp)); % Pass band response; Mp(1) = 1
Mp(2:end) = abs(M * R * sin(pi*fp(2:end)/R) ./ sin(pi*M*fp(2:end))).^(N);
Mf = [Mp zeros(1, length(fs))];

Вид полученной характеристики представлен на рисунке ниже:


5 шаг: По идеальной частотной характеристике производим расчет импульсной характеристики компенсирующего FIR фильтра с заданными параметрами (весовое окно, если есть, разрядность коэффициентов и порядок фильтра).

Вид импульсной характеристики фильтра-корректора (расчет производится с помощью встроенной функции FIR2, а оконная функция считается через функцию FIR1):


Исходный код для расчета ИХ корректирующего фильтра:

% Calculate FIR
hFIR = fir2(NFIR-1, f, Mf); % Filter length NFIR
hFIR = hFIR / max(hFIR); % Double coefficients
hCOE = round(hFIR*(2^(Bc-1)-1)); % Fixed point coefficients

% Windowed FIR (Kaiser with BETA) 
if (IS_WIND == 'Y')
  WIND = kaiser(NFIR, BETA); % KAISER WINDOW IS USED!
  hWIND = fir1(NFIR-1, Fo/R, 'low', WIND);
  hNEW = hCOE .* hWIND;% conv2(hCOE,Hwind);
  hCOE = hNEW;
end

На этом расчет фильтра заканчивается. Результатом является файл коэффициентов фильтра в удобном для пользователя виде. Как видно, ничего сложного нет и весь процесс разбивается на 4-5 этапов:

• Определение входных параметров CIC и FIR фильтров,
• Формирование идеальной частотной характеристики,
• Расчет коэффициентов FIR фильтра (функция FIR2),
• Построение АЧХ компенсирующего фильтра для проверки,
• Расчет неравномерности в полосе пропускания.

Следует отметить, что для схем с понижением частоты дискретизации (задача децимации) компенсирующие фильтры обычно стоят после CIC фильтров. Для задач повышения частоты дискретизации (интерполяция) компенсирующие фильтры стоят до CIC фильтров. Такая расстановка фильтров наилучшим образом подходит с точки зрения объема занимаемых ресурсов кристалла, поскольку корректирующий FIR фильтр в этом случае работает на пониженной частоте (после децимации или до интерполяции).

На следующем рисунке приведены АЧХ CIC и FIR фильтров и результирующая АЧХ после компенсации при следующих параметрах:

• Коэффициент децимации R = 11;
• Порядок CIC фильтра N = 4;
• Задержка в цепи дифф. звена M = 1;
• Порядок FIR-фильтра NFIR = 64;
• Разрядность коэффициентов Bc = 16;
• Нормированная частота среза Fo = 0.4;
• Оконная функция — не используется

Как видно, АЧХ CIC фильтра скомпенсирована, а результирующая АЧХ визуально обладает достаточной прямоугольностью (кривая красного цвета).

Неравномерность АЧХ в полосе пропускания

К сожалению, на этом выравнивание АЧХ не заканчивается. Увеличим график в пределах полосы пропускания. На следующем графике можно заметить небольшие биения в полосе.


Почему так? Неравномерность АЧХ в полосе пропускания тем выше, чем больше порядок CIC-фильтра N. Коэффициент децимации/интерполяции R не влияет на неравномерность. Также на неравномерность влияет порядок компенсирующего КИХ-фильтра. Чем меньше порядок фильтра, тем больше неравномерность. Но чем больше порядок фильтра, тем больше частота «биений» в полосе пропускания. Именно по этой причине в практических целях не применяются компенсирующие фильтры высоких порядков NFIR > 128! Эффективная полоса тоже вносит незначительный вклад: чем меньше выравниваемая полоса, тем меньше неравномерность.

Для борьбы с биениями применяются оконные функции, о которых я говорил ранее. Посмотрим, как выглядит график в укрупненном масштабе без оконной фильтрации и с применением функции Кайзера (BETA = 8). Остальные параметры остаются неизменными.

• Синий — результат без оконной фильтрации,
• Красный — применение оконной функции Кайзера.

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

Второй способ улучшения свойств частотной характеристики — разбиение обработки на несколько стадий. Например, децимация проводится не в один этап, а с помощью нескольких звеньев понижения частоты дискретизации: R = R1 * R2 *… Rn.

С целью экономии ресурсов ПЛИС в некоторых задачах возможно совмещение корректирующего и формирующего FIR-фильтров. Для этого нужно перемножить их импульсные характеристики или выполнить свертку спектров.

Результат

Для удобства расчета фильтра-корректора написан m-скрипт, позволяющий получить результат в наглядной форме. Скрипт умеет выводить данные коэффициентов в одном из нескольких популярных форматов, но при желании его можно дописать и модифицировать под свои цели.

• COE – формат Xilinx для загрузки коэффициентов в Core Generator,
• H – файл-header для подключения в проекты на С/C++ (Code Composer Studio),

Для запуска скрипта необходим САПР Matlab или GNU Octave (я отлаживался в последнем). Далее полученный файл коэффициентов загружается в формате *.COE, в Xilinx Core Generator, где реализуется компенсирующий FIR-фильтр.



Схематичный вид КИХ фильтра в среде Xilinx Vivado (вкладка post-synthesis):

Альтернативный метод

Ещё один способ расчета корректирующих фильтров предложил мой коллега по работе — Антонов А.Е. Он написал интерактивное приложение, которое пошагово в ручном режиме подбирает требуемую неравномерность АЧХ фильтра и выгружает коэффициенты для дальнейшей работы. Предлагаемая программа — c_koeff.exe — предназначена для корректировки коэффициентов FIR фильтра, стоящего за системой из двух CIC-фильтров, для компенсации неравномерности АЧХ всей системы в полосе пропускания, вносимой CIC-фильтрами. Предполагается, что коэффициенты FIR фильтра получены в среде FDATool пакета MATLAB и заданы в виде C-header файла. Программа написана на C++ Builder. Для отображения рассчитанных коэффициентов скомпенсированного FIR фильтра и его АЧХ используется программа ISVI 6 и более версии, разработанная в ЗАО «ИнСис» и распространяемая бесплатно. Дабы не нарушать правила ресурса, рекламировать достижения фирмы и особенности прикладного софта я не буду.

Результат работы программы изображен на следующем рисунке. Путём перебора параметра Np, получаем спектральную функцию компенсирующего фильтра. В приложении ISVI также можно наблюдать импульсную характеристику фильтра.


Описание приложения находится в word-файле среди исходников на github (ссылка ниже).

Исходный код

• m-скрипты для Matlab / GNU Octave (расчёт FIR фильтров и корректирующих фильтров)
• Приложение c_koeff для интерактивного расчета фильтра-корректора

Литература

• Цифровая обработка сигналов на ПЛИС — 1
• Цифровая обработка сигналов на ПЛИС — 2
• Altera DS (Understanding CIC Compensation Filters)
• Использование CIC фильтров в задачах децимации и интерполяции
• Э.С. Айфичер, Барри У. Дервис., Цифровая обработка сигналов. Практический подход
• Рабинер Л., Голд Б., Теория и применение цифровой обработки сигналов
• An economical class of digital filters for decimation and interpolation

Спасибо за внимание!