В процессе цифровой обработки сигналов нередко возникает задача фильтрации сигнала в очень узком частотном диапазоне. Примером такой задачи может служить ситуация, когда сигнал содержит несколько составляющих на близких частотах, и требуется выделить одну из них. Для этого необходимо спроектировать цифровой фильтр с очень узкими относительными полосами пропускания и перехода. В зависимости от расположения составляющих сигнала, это может быть фильтр нижних, верхних частот, полосовой или заграждающий фильтр. Прямое проектирование таких устройств часто приводит к фильтрам очень высоких порядков, практическая реализация которых либо нерациональна, либо невозможна. Рассмотрим следующий пример. Пусть имеется дискретный сигнал
u(n) = sin(2πf1nT) + sin(2πf2nT) + e(n) ,
где f1 = 90 Гц и f2 = 105 Гц - частоты составляющих; T = 1/fs = 1/8000 с - период дискретизации; fs = 8000 Гц - частота дискретизации; e(n) - гауссовский шум с нулевым средним и дисперсией σ2 = 0,1. Требуется выделить синусоидальную составляющую с частотой f1. Для этого необходимо построить фильтр нижних частот с граничными значениями частот полосы пропускания fpb и полосы задерживания fsb, удовлетворяющими условию f1 < fpb< fsb < f2 . Поскольку частота Найквиста fN = fs/2 = 4000 Гц, для нормализованных значений fpb = fpb/fN и fsb = fsb/fN это условие будет выглядеть так :
Следовательно, переходная полоса Δ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 раз.
После нажатия на кнопку Apply появится предупреждение, что оценочное значение порядка проектируемого фильтра 5022 слишком велико, и дальнейшая процедура расчёта может дать непредсказуемые результаты. Однако даже если такой фильтр построен, для его применения потребуется очень большой вычислительный ресурс. Так как коэффициенты фильтра обладают симметрией и общее их количество равно 5023, процессор должен выполнять
(5022/2 + 1) MACs/sample х 8000 samples/sec = 20096000 MACs/sec
(Multiply-And-aCcumulate operations per second - число операций умножения-накопления в секунду, показатель быстродействия процессора). Для хранения коэффициентов потребуется 2512 ячеек памяти. Альтернативный путь решения рассматриваемой задачи заключается в использовании идей многочастотной дискретизации - последовательном прохождении сигнала через полифазные фильтры. Как было показано, каждый из них выполняет две функции: фильтрацию в заданной полосе и изменение частоты дискретизации выходного сигнала. Прежде всего, для восстановления отфильтрованного сигнала [0, 90] Гц (сигнал содержит только составляющую на частоте 90 Гц), его достаточно дискретизировать с частотой, большей или равной fs1 = 2·90 = 180 Гц. Для выделения составляющей с частотой fs выполним фильтрацию сигнала, используя два полифазных фильтра-дециматора и два полифазных фильтра-интерполятора. Пусть первый фильтр имеет полосу пропускания
[0, 1/20–0,025]·FN = [0, 100] Гц
и полосу задерживания
[1/20, 1]·FN = [200, 4000] Гц.
АЧХ в полосе пропускания этого фильтра должна находиться в пределах [0,995, 1,005], а в полосе задерживания - не превышать величину 10-4. Исходя из этих требований, коэффициент децимации М положим равным 20. При этом частота дискретизации на выходе первого фильтра дециматора будет равна 8000/20 = 400 Гц. Расчёт коэффициентов первого фильтра можно выполнить опять с использованием GUI sptool. Вводимые данные для расчёта:
8000 - в поле Sampling Frequency;
100 и 200 - в поля Fp и Fs, соответственно;
20*log10(1,005)–20*log10(0,995) или 0,08686 - в поле Rp;
80 - в поле Rs, что соответствует ослаблению в 10000 раз.
Результаты расчёта дают приблизительную оценку порядка фильтра n1 = 270, однако полученные значения Actual Rp = 0,1436 и Actual Rs = 75,58 (они отображаются в правой части окна Filter Designer) указывают на необходимость увеличить порядок фильтра и повторить вычисления. В нашем примере мы выбрали значение n1 = 289, что соответствует Actual Rp = 0,08281 и Actual Rs = 80,16. Второй фильтр дециматор характеризуется следующими параметрами:
400 в поле Sampling Frequency;
95 и 100 в полях Fp и Fs, соответственно;
0,08686 в поле Rp;
80 в поле Rs, что соответствует ослаблению в 10000 раз.
Результаты вычислений дают предварительную оценку порядка фильтра n2 = 270 и фактические значения Rp и Rs для этого порядка, которые опять-таки не укладываются в заданные требования (пульсации недопустимо большие). Выбор n2 = = 275 обеспечивает выполнение требований. Чтобы коэффициенты спроектированных фильтров стали доступны для дальнейшего использования, их надо экспортировать в рабочее пространство MATLAB. Для этого необходимо открыть меню File главной панели sptool и выбрать раздел Export… Если экспортированные структуры имеют имена filt1 и filt2, извлечь соответствующие коэффициенты можно с помощью команд
b1=filt1.tf.num; b2=filt2.tf.num;
Применение двух фильтров дециматоров, реализуемых в виде полифазных структур, требует выполнения
и хранения 145 + 138 = 283 коэффициентов. Для восстановления исходной частоты дискретизации 8000 Гц следует к выходу второго фильтра последовательно подключить два фильтра-интерполятора, причём коэффициент интерполяции первого фильтра равен 2, а второго - 20. С учётом того, что эти фильтры выполняют столько же операций в секунду, сколько и предшествующие им фильтры- дециматоры, общее количество операций умножения с накоплением в секунду будет равно 85600·2 = 171200, а число коэффициентов фильтров - 283.
