Многоскоростные фильтры и банки фильтров (фильтры и банки фильтров с многочастотной дискретизацией) определили самостоятельное направление в теории и практике цифровой обработки сигналов (ЦОС). В этой связи достаточно упомянуть квадратурно-зеркальные фильтры — новый класс многополосных фильтров, банки многоскоростных фильтров для реализации вейвлет-преобразований, полифазные фильтры. Многоскоростная фильтрация нашла широкое практическое применение в задачах компрессии речи, звука и изображений, построения эффективных систем фильтрации, очистки сигнала от помех, оптимизации вычислительных ресурсов при реализации алгоритмов ЦОС.

Для решения задач анализа (моделирования) и синтеза (проектирования) систем с многочастотной дискретизацией пакет MATLAB предоставляет широкие возможности. В рамках проекта "MATLAB для DSP" рассмотрим частную, но практически важную задачу проектирования и моделирования узкополосного фильтра. Для этих целей, как это и оговорено в проекте, используются два наиболее дружественных для пользователя инструмента MATLAB - Simulink и GUI.
Для удобства читателя в первом разделе статьи приведены краткие теоретические сведения по существу затрагиваемых вопросов. При этом предполагается, что читатель знаком с предметной областью. В последующих разделах на примере решения конкретной задачи - синтез и моделирование узкополосного ФНЧ - излагаются правила и практические рекомендации по использованию инструментов MATLAB.

Основные соотношения теории многочастотной дискретизации.

Основы теории многоскоростной фильтрации (многоскоростной обработки сигналов) и описание приложений можно найти в [1-4].

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

Таким образом, частота дискретизации f_s сигнала x(n) связана с частотой дискретизации f'_s соотношением Mf'_s = f'_s , или соотношением периодов дискретизации T' = MT. Частоты дискретизации сигналов x(n) (частота f_s ) и y_l(n) (частота f_s ) при интерполяции связаны соотношениями f'_s = Lf_s или T = LT'.

Пример децимации сигнала x(n) при M = 2: отсчёты сигнала до (а) и после (б) децимации.

Рис 3.

Пример интерполяции сигнала x(n) при L = 2: отсчёты сигнала до (а) и после (б) интерполяции

Рис 4.

При децимации и интерполяции сигнала происходит деформация спектров.
Для децимации

где z = ejw , а

и, следовательно,

где нормированная частота w = wT', а T' - период дискретизации после децимации. Он будет в M раз больше, так что в шкале ненормированных частот получим

Учитывая, что T' = MT , можно записать окончательное соотношение между спектрами децимированного сигнала Y_D(e^jw) и исходного X(e^jw):

Для интерполяции Y_I(z) = X(z^L) или Y_I(e^jw) = X(e^jwL) или Y_I(e^jw) = X(e^jwL).
Таким образом, спектр децимированного сигнала является взвешенной суммой исходного спектра X(e^jw) и его (M–1) сдвинутых по частоте копий (отражений) с шагом 2Пw_s/M. Спектр же интерполированного сигнала является спектром исходного сигнала с изменённым периодом по частоте. Период увеличивается в L раз.
Учитывая упомянутые выше свойства спектра, необходимо перед децимацией ставить фильтр децимации, чтобы исключить наложения спектра, а для интерполяции - фильтр интерполяции для устранения отражений, то есть тех дополнительных компонент спектра, которые попали в рабочую полосу [0,F_N] за счёт увеличения периода спектра.
Замечательные тождества. При построении систем с многочастотной дискретизацией очень полезны преобразования, изображённые на рисунках. Они полезны во многих случаях при реализации фильтров, что будет продемонстрировано в следующем разделе.

Рис 5.

Рис 6.

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

может быть представлено суммой

H(z) = h(0) + h(2)z-2 + ...
+ h(1)z-1 + h(3)z-3 + ... или
H(z) = h(0) + h(2)z-2 + ...
+ z-1 (H(1) + h(3)z-2 + ...) = E0(z2) + E1(z2) .

Смысл многочленов E₀ и E₁ понятен из контекста.Если E₀ и E₁ рассматривать как передаточные функции КИХ-фильтров, то нетрудно заметить, что базовым элементом задержки таких фильтров является z^-2, обеспечивающий задержку на два такта. Следовательно, фильтры E₀ и E₁ могут работать на частоте дискретизации, в два раза меньшей исходной. Если использовать разложения по степеням z^-3 или z^-4 и так далее, то можно получить блоки фильтров, работающие на ещё более низких частотах дискретизации. Итоговая схема фильтра, когда H(z) = E0(z2) + z-1E1(z2) ,показана на рисунке 7.

Рассмотренное разбиение называется полифазным, а реализующие его схемы - полифазными фильтрами.В качестве примера, иллюстрирующего построение полифазного фильтра и использование замечательных тождеств, покажем, как можно эффективно реализовать КИХ-фильтр дециматор.

Рис 8.

В схеме рис. 8 для вычисления каждого отсчёта необходимо выполнить N+1 умножений и N сложений. В то же время, за счёт децимации (M = 2) половину результатов отсчётов мы отбрасываем и, следовательно, используем ресурсы вычислителя неэффективно.Построим фильтр на основе полифазного разбиения (полифазный фильтр, рис. 7) и воспользуемся замечательным тождеством (рис. 5). Последовательность преобразований при этом показана на рис. 9.

Рис. 9

В результате фильтры E₀(z²) и E₁(z²), работающие на частоте f_s, заменяются фильтрами E₀(z) и E₁(z), работающими на частоте f_s/2. Таким образом мы получим двукратную эконом ию в скорости вычислений. При коэффициенте децимации M применение полифазных фильтров позволяет получить экономию в M раз. Аналогичные результаты получаются и при полифазном построении фильтров интерполяторов [1,2]. Таким образом, полифазные реализации в совокупности с рациональными преобразованиями схем фильтров позволяют строить эффективные вычислители в задачах ЦОС.

Расчёт узкополосного низкочастотного фильтра.

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

u(n) = sin(2πf₁nT) + sin(2πf₂nT) + e(n) ,

где f₁ = 90 Гц и f₂ = 105 Гц - частоты составляющих; T  = 1/f_s = 1/8000 с - период дискретизации; f_s = 8000 Гц - частота дискретизации; e(n) - гауссовский шум с нулевым средним и дисперсией σ² = 0,1. Требуется выделить синусоидальную составляющую с частотой f₁. Для этого необходимо построить фильтр нижних частот с граничными значениями частот полосы пропускания f_pb и полосы задерживания fsb, удовлетворяющими условию f₁ < f_pb< f_sb < f₂ . Поскольку частота Найквиста f_N = f_s/2 = 4000 Гц, для нормализованных значений f_pb = f_pb/f_N и f_sb = f_sb/f_N это условие будет выглядеть так :

Следовательно, переходная полоса Δf амплитудно-частотной характеристики фильтра (АЧХ) не должна превышать величину Δf < 0,02625 – 0,0225 = 0,00375. Эти требования схематически показаны на рис. 10.
Для решения сформулированной задачи попробуем спроектировать нерекурсивный фильтр с полосой пропускания от 0 до 95 Гц и полосой задерживания от 100 до 4000 Гц, то есть до частоты Найквиста. Кроме того, потребуем, чтобы АЧХ в полосе пропускания находилась в пределах [0,99, 1,01], а в полосе задерживания не превышала значения 0,0001. Для этого нам необходимо воспользоваться функциями пакета MATLAB, оценивающими по заданным требованиям порядок фильтра и рассчитывающими его коэффициенты. Эти действия удобно выполнять с помощью графической программы (GUI) sptool, входящей в библиотеку signal пакета MATLAB и описанную в [5]. Загрузим эту программу с помощью инструкции sptool и на появившейся главной панели нажмём кнопку New Design, находящуюся под окном Filters. После этого на экране появится панель Filter Designer.
Чтобы рассчитать фильтр, надо ввести данные в поля Sampling Frequency: Fp (верхнее граничное значение полосы пропускания), Rp (максимально допустимая величина пульсаций в полосе пропускания), Fs (нижнее граничное значение полосы задерживания) и Rs (максимально допустимая величина пульсаций в полосе задерживания). Как видно, соответствующие значения пульсаций, вводимые в поля Rp и Rs, должны быть заданы в децибелах (вводятся абсолютные величины). Введём данные, необходимые для расчёта нашего фильтра:

• 8000 - в поле Sampling Frequency;
• 95 и 100 - в поля Fp и Fs, соответственно;
• 20*log10(1,01)–20*log10(0,99) - в поле Rp;
• 80 - в поле Rs, что соответствует ослаблению в 10000 раз.