Сигналы и системы

Лекция 9 Сигналы и системы 6 – й семестр 2017/18 – й учебный год доц. Щетинин Ю.И. Синтез фильтров с конечной импульсной характеристикой (КИХ – фильтров) Новосибирский государственный технический университет Факультет автоматики и вычислительной техники Кафедра систем сбора и обработки данных Общие сведения КИХ-фильтры или нерекурсивные фильтры – это цифровые фильтры с уравнением вида M y[n]  b[n]  x[n]   bk x[n  k ]  h[0]x[n]  h[1]x[n  1]  ..  h[ M ]x[n  M ] , (1) k 0 М - порядок, bk=h[k]- коэффициенты фильтра, N=M+1 – число коэффициентов (длина) фильтра. Уравнение фильтра (1)– это дискретная свертка. Такие фильтры называют также трансверсальными (transversal filters) или линиями задержки с отводами (tapped delay line). Пример структурной схемы фильтра 4-го порядка с 5 коэффициентами h[0] z -1 z -1 h[1] h[2] z -1 z -1 h[4] h[3] y[n] Для реализации фильтра необходимы элементы (операции) задержки, умножения и сложения. Применяя Z-преобразование к свертке, получим уравнение КИХ - фильтра в Z - плоскости Y ( z )  b0 X ( z )  b1z 1 X ( z )  ...  bM z  M X ( z ) . H  z  Передаточная функция такого фильтра M Y ( z)   bk z  k  b0  b1z 1  ...  bM z  M X ( z) k 0 (2) не имеет полюсов, а имеет только нули. Из уравнений(1) и (2) очевидно, что значения импульсной характеристики и коэффициенты передаточной функции фильтра совпадают. M Частотная характеристика КИХ – фильтра H (e j )  H ( )   bk e j k . k 0 Достоинства КИХ – фильтров: всегда устойчивы, достаточно легко получить линейную ФЧХ, эти фильтры более устойчивы к погрешностям квантования коэффициентов и округления операций умножения и сложения, чем БИХ фильтры. Задачей синтеза (проектирования) КИХ - фильтра является получение его коэффициентов h[k ]  bk . 1 Параметры спецификации АЧХ ФНЧ: fp (fpass) – граничная частота ПП, fs(fstop) – граничная частота ПЗ, δP(дБ) – допустимое отклонение АЧХ от 1 (0дБ) в ПП, δS (дБ) - минимальное затухание в ПЗ, Затухание (пульсации) в полосе пропускания в дБ определяются по выражению Затухание (ослабление) в полосе задерживания в дБ . 2δ1Спецификация 1+δP 1 1-δP δS ФНЧ ПП ПЗ f Основные этапы разработки КИХ – фильтров и БИХ – фильтров fP fS практически одинаковы: 1. Выбирается метод синтеза (проектирования) фильтра. 2. Определяется порядок и передаточная функция фильтра. 3. Выбирается структура фильтра, реализующая требуемую передаточную функцию. 4. Анализируется влияние конечной разрядности процессора на характеристики фильтра. 5. Реализация (implementation) фильтра: разрабатывается программный код или специальное аппаратное устройство в виде цифрового процессора или специальной СБИС. Для синтеза КИХ – фильтров разработано достаточно много различных методов. Они реализованы, в частности, в соответствующих функциях и процедурах пакетов – расширений Signal Processing Toolbox и Filter Design and Analysis Toolbox MATLAB. Наиболее часто используются следующие методы синтеза:  методы оптимизации частотных характеристик, основанные на процедурах минимизации нормы разности между желаемой и аппроксимирующей характеристиками фильтра.  оконный метод, основанный на разложении частотной характеристики целевого фильтра в ряд Фурье с последующим взвешиванием импульсной характеристики фильтра с помощью специальных взвешивающих оконных функций,  метод частотных выборок, В данном курсе (лекции) рассматриваются только первый и второй методы синтеза. Метод частотных выборок опускается из – за отсутствия времени. Типы КИХ – фильтров с линейной фазовой характеристикой Линейная фазовая характеристика необходима, например, для качественных аудиосистем, систем передачи данных, систем обработки изображений, обработки биомедицинских сигналов и др. задач. Для этого достаточно, как будет доказано ниже, чтобы импульсная характеристика фильтра имела четную или нечетную симметрию. Имеется 4 типа импульсных характеристик КИХ – фильтров с линейной ФЧХ. Здесь N – длина (число коэффициентов) фильтра. Тип 2 Тип 1 Тип 3 Тип 4 Таким образом, для фильтров с четной симметрией необходимо: hk=h N-k, k=0,1,…,N Для фильтров с нечетной симметрией (антисимметрия): hk= -h N-k, k=0,1,…,N . 2 Докажем, что ЦФ будет иметь линейную фазовую характеристику, если его импульсная характеристика симметрична относительно средней точки, т.е. (см. рис. выше). Доказательство. Возьмем каузальную импульсную характеристику типа 1 , Сдвинем h[n] на n1 = N/2 отсчетов влево. Получим симметричную относительно нуля последовательность Найдем ДВПФ h1[n], т.е. её частотную характеристику (3) Перейдем теперь от h1[n] к h[n], т.е. сдвинем h1[n] вправо на отсчетов. По свойству временного сдвига ДВПФ , поэтому , т.е. АЧХ при сдвиге не изменяется, изменяется только ФЧХ фильтра. Из полученного выражения следует, что частотная характеристика при такой импульсной характеристике имеет линейную фазу . ω -ωN/2 Производная ФЧХ по частоте носит название групповой задержки. Из вышеприведенного выражения следует, что групповая задержка для фильтров с симметричной ФЧХ равна N/2. Фазовые искажения сигнала из-за нелинейности ФЧХ фильтра Покажем, что линейная ФЧХ не изменяет форму сигнала, а нелинейная – искажает. В лекции № 7 было показано, что выходной сигнал дискретного (цифрового) фильтра при входном сигнале в виде гармоники с частотой ω0 , начальной фазой θ и амплитудой А имеет вид y (nT )  A H (e jω0T ) cos(ω0 nT  θ  φ(ω0 )) . j T В этом выражении H (e 0 ) - значение АЧХ на частоте ω0,  (0 ) - значение ФЧХ фильтра на той же частоте, при этом Т – интервал дискретизации (отсчетов) сигналов. Если ФЧХ фильтра – линейная, т.е. имеет вид  ()  k , то выход фильтра y[n]  A H cos((n  k )) , т.е. все частотные компоненты сигнала пройдут на выход с одинаковой задержкой. При H  1 выходной сигнал по сравнению с входным просто задержится на k отсчетов. Этот случай иллюстрируется ниже на первом и втором рис. для сигнала из двух гармоник с нормализованными цифровыми частотами 0.06π рад и 0.18π рад: Вход x=cos(0.06*pi*n)-0.5*cos(0.18*pi*n) Выход при задержке k= 10: y1=cos(0.06*pi*(n-10))-0.5*cos(0.18*pi*(n-10)). Видно, что в данном случае при линейной ФЧХ выходной сигнал представляет собой просто задержанный на 10 отсчетов входной сигнал. На третьем рис. показан результат, когда ФЧХ фильтра нелинейная и обе гармоники входного сигнала независимо от частоты задерживаются по фазе на одинаковый угол π/2 y2=cos(0.06*pi*n- pi/2)-0.5*cos(0.18*pi*n-pi/2); В этом случае сдвиг каждой гармоники входа не пропорционален частоте и выходной сигнал не просто имеет задержку, а претерпевает фазовые искажения (phase distortion). В высококачественной аудиообработке, в задачах видеообработки сигналов, обработке биомедицинских сигналов и многих других приложениях фазовые искажения не допустимы. Поэтому должны применяться фильтры с линейной ФЧХ, в частности, КИХ – фильтры. 3 Входной сигнал 2 -2 10 20 30 40 50 60 50 60 50 60 Линейный сдвиг фазы на 10 отсчетов 2 -2 10 20 30 40 Нелинейный фазовый сдвиг 2 -2 10 20 30 n 40 Оптимальные чебышевские КИХ – фильтры Среди множества методов проектирования КИХ – фильтров имеется класс т.н. оптимальных методов синтеза. В них специальными математическими методами минимизируется отклонение действительной частотной характеристики фильтра H(ω) от желаемой (desired) АЧХ D(ω). В качестве меры отклонения действительной АЧХ от желаемой используется т.н. р – норма ошибки, определяемая выражением 1/𝑃 𝜔 𝐿𝑝 = [∫0 𝑑 𝑊(𝜔)|𝐷 (𝜔) − 𝐻(𝜔)|𝑃 𝑑𝜔] , (4) в котором D(ω) – желаемая АЧХ, H(ω)- АЧХ реального фильтра, W(ω) - неотрицательная весовая функция частоты, позволяющая нужным образом увеличивать или, наоборот, уменьшать на тех или иных частотах (диапазонах) влияние разностей D ( )  H ( ) на значение нормы ошибки Lp. P Т.о., весовая функция W(ω) позволяет придавать разную значимость различным участкам разности в зависимости от частоты ω. Чаще всего используются два вида норм: р = 2 (квадратичная норма ошибки) и р = ∞ (максимальное абсолютное значение ошибки). При р = ∞ норма ошибки соответствует максимальному абсолютному отклонению действительной характеристики фильтра от желаемой. Процесс проектирования при применении такой нормы заключается в минимизации выражения ошибки Е(ω) (5) E()  W ()[D(e j )  H (e j )] на множестве значений коэффициентов КИХ – фильтра. Математически это можно записать в виде min[max|E(ω)|], (6) т.е. отыскивается минимальное значение максимума значения ошибки. Эта процедура в математике известна как задача чебышевской минимаксной аппроксимации. Получаемые в результате её решения фильтры называют чебышевскими оптимальными фильтрами или минимаксными фильтрами. Алгоритм решения применительно к синтезу КИХ – фильтров разработан учеными Парксом (T.W. Parks) и МакКлеланом (J. H. McClellan). Этот алгоритмв свою очередь использует для проектирования фильтра т.н. обменный алгоритм Ремеза (Remez) и теорию чебышевской аппроксимации. Получаемые фильтры оптимальны в смысле минимума максимальной ошибки между желаемой и действительной (actual) АЧХ фильтра. В результате спроектированные фильтры имеют равноволновое (equiripple) поведение АЧХ в частотной области и поэтому их называют также равноволновыми фильтрами (equiripple filters). 4 Сам алгоритм достаточно сложен, требует определенной математической подготовки и в этом курсе не рассматривается. Ниже на рис. поясняются особенности характеристик фильтров, получаемых методом Паркса – МакКлелана. На верхнем графике показывается, что АЧХ такого фильтра имеет равные колебания (пульсации) в полосе пропускания и полосе задерживания. Но сама амплитуда (уровень) этих пульсаций в ПП (δр) и ПЗ (δs) может быть разной. 1+δP |H(ejω)| реальная АЧХ 1 1-δP |Hd(ejω)| идеальная АЧХ ω δS На среднем рис. изображены ошибки между идеальной и действительной АЧХ δP δS -δP -δS Когда ошибки в ПП и ПЗ – разные, то сама оптимизационная процедура определения минимакса резко усложняется, поэтому для реализации процедуры необходимо уравнять (сбалансировать) амплитуды ошибок в ПП и ПЗ. Результат показан на нижнем графике ω Балансирование влияния ошибок можно выполнить с помощью весовых коэффициентов ошибок в полосе пропускания WP и полосе задерживания WS. Определение значений WP и WS происходит в процедуре проектирования. Алгоритм Паркса - МакКлелана проектирования оптимального КИХ – фильтра содержит следующие этапы: 1. Задание требований к АЧХ: задаются граничные частоты ПП, ПЗ, допустимый уровень пульсаций в ПП, ослабление (затухание) в ПЗ, порядок фильтра и частота дискретизации. Граничные угловые частоты нормируются на частоту Найквиста Fs/2. 2. Абсолютные значения отклонений (пульсаций) в ПП δP и затухания (ослабления) в ПЗ δS по исходным значениям δP и δS , выражаемым в дБ, вычисляются по выражениям:   P  10  P dB 20   S  10  1,  S dB 20 . Вычисленные значения определяют числитель и знаменатель в следующем выражении W P numerator   S . При этом значения Ws и Wp имеют смысл весовых множителей  S deno min ator WP ошибок для ПЗ и ПП фильтра. 4. Для вычисления коэффициентов фильтра используется специальный алгоритм Ремеза (Remez algorithm), входящий в алгоритм Паркса - МакКлелана. 5. Если получившийся фильтр не удовлетворяет спецификации, то порядок фильтра увеличивается, и шаги 1-4 повторяются. В данном курсе сама процедура алгоритма Ремеза не рассматривается из-за её громоздкого характера. Для её реализации в MATLAB имеется встроенная функция remez(). Рассмотрим пример проектирования оптимального ФНЧ со следующей спецификацией: Порядок фильтра = 50. Полоса пропускания фильтра (ПП) 0 – 800 Гц. Полоса задерживания 1000 – 4000 Гц. Пульсации в ПП должны быть максимум 1дБ. Затухание в ПЗ минимум 40 дБ. Частота отсчетов 8000 Гц. 3. 5 Решение. Выполним нормализацию частот на частоту Найквиста 4000 Гц. Нормализованная граничная частота ПП: 800/4000=0,2. Желаемая АЧХ в ПП д.б. равна 1. Граничные нормированные частоты ПЗ: 1000/4000=0,25 4000/4000 =1. Желаемая АЧХ в ПЗ = 0. Определяем по п.2 алгоритма абсолютные значения отклонений в ПП и затухания в ПЗ: . Определяем веса в ПП и ПЗ (п.3) . Получаем Ws = 12 и Wp = 1. Для вычисления коэффициентов фильтра h[k] применяем функцию remez() пакета MATLAB. Соответствующая программа MATLAB и её результаты приведены ниже Fs=8000; % частота отсчетов f=[0 0.2 0.25 1]; % нормированные частоты m=[1 1 0 0]; % Значения идеальной (желаемой) АЧХ w=[1 12]; % Весовые коэффициенты ошибок h=remez(50,f,m,w); % Вычисление коэффициентов фильтра freqz(h,1,256,Fs); % построение АЧХ figure, impz(h,1)% построение ИХ фильтра Примечание. Функция remez в настоящее время считается устаревшей, хотя и сохранена в MATLAB. Вместо remez сейчас применяется функция firpm. Синтез оптимальных чебышевских фильтров в среде MATLAB В настоящее время в MATLAB имеются 3 функции для синтеза КИХ – фильтров на основе алгоритма Паркса – МакКлелана: firpm, cfirpm и firgr. Ниже рассмотрена наиболее простая из них функция b = firpm(n, fn, a0, w). Эта функция в современных версиях MATLAB заменила функцию remez, которая пока сохранена, но считается устаревшей. Оценка порядка фильтра n и остальных входных параметров для функции firpm может быть выполнена с помощью специальной функции [n, fn, a0, w] = firpmord(f,a,dev,Ft). Эта функция определяет порядок фильтра приближенно, с точностью до 2-3 единиц. Поэтому после проектирования фильтра необходимо проверять соответствие АЧХ спроектированного фильтра исходной спецификации и при необходимости уточнить порядок и перепроектировать фильтр. Функция firpmord предназначена для определения порядка многополосного нерекурсивного фильтра, т.е. фильтра, у которого может быть (но необязательно) несколько полос пропускания и полос задерживания. При этом параметры f и a задают кусочно – постоянную АЧХ идеального (желаемого) фильтра. Вектор частот f должен содержать возрастающую последовательность значений частот, определяющих границы полос фильтра. Вектор а определяет постоянные желаемые значения АЧХ в соответствующих полосах частот. Первая полоса начинается с нулевой частоты, последняя частота из вектора f равна частоте Найквиста. Длина вектора f равна length(f)=2*length(a) - 2. Пример соответствия между векторами f и a показан ниже на рис. a f = [0 .3 1.0 a = [0 1 .4 .6 .7 .9] 0.5 0.5] 0.5 Нормализованная частота f 0. 0.1 0.2 0.3 0.4 0.5 0.6 6 0.7 0.8 0.9 1.0 Параметр dev в функции – это вектор максимально допустимых отклонений АЧХ фильтра от желаемой АЧХ на частотах, специфицированных в векторе f. Они связаны с соответствующими Rp 20  S  10  Rs 20 значениями Rp и Rs в децибелах выражениями:  P  10  1 и . Параметр (вектор) dev должен иметь такую же длину, что и параметр а. Параметр Ft в функции - частота дискретизации в Гц и не является обязательным. Указанные параметры передаются в основную функцию синтеза фильтра firpm() при её вызове. Функция синтеза фильтра firpm() имеет 9 вариантов синтаксиса. В одном из наиболее простых вариантов b = firpm(n,f,a)функция принимает основные выходные параметры функции порядка firpmord и возвращает вектор – строку b, содержащую (n+1)коэффициентов КИХ – фильтра порядка n , АЧХ которого соответствует значениям векторов f, a, dev. В функции синтеза фильтра [b, err, opt] = firpm(n, fn, a0, w) выходные параметры функции синтеза фильтра: b – вектор коэффициентов передаточной функции, он же является вектором значений импульсной характеристики, err - максимальное отклонение АЧХ от идеальной частотной характеристики на рассматриваемом интервале частот. Параметр opt - это массив дополнительных записей с 9 полями, содержащими дополнительные значения, полезные в уточненных процедурах синтеза фильтров. Выходные параметры err и opt - не обязательные, их можно опускать, т.е. использовать функцию в виде b = firpm(n, fn, a0, w) . Но для уточненных процедур синтеза, более точных, они нужны. Аналогичное назначение имеет группа функций B = cfirpm(N,F,A,W). В отличие от функции b = firpm(n, fn, a0, w) эта функция предназначена для синтеза не только действительных, но и комплексных фильтров, а также фильтров с линейной и нелинейной ФЧХ. Пример. Рассчитаем оптимальный чебышевский полосовой фильтр с граничными частотами полосы пропускания 2400 и 3200 Гц, полосы задерживания 2000 и 3600 Гц, допустимыми пульсациями (отклонениями) в полосе пропускания 1 дБ и минимальным затуханием в полосе задерживания 50 дБ. % Процедура проектирования % фильтра с помощью firpm Ft=10000; % частота отсчетов Ap=1; % Отклонение в ПП в дБ As=50;% Затухание в ПЗ в дБ a=[0 1 0]; % вектор значений желаемой АЧХ % вектор граничных значений частот ПП и ПЗ f=[2000, 2400, 3200, 3600]; % Абсолютные отклонения в ПП и ПЗ Dp=10^(Ap/20)-1; Ds=10^(-As/20); dev=[DsDpDs]; % вектор допустимых ошибок АЧХ % Определение порядка фильтра [n, f0, a0, w]=firpmord(f,a,dev,Ft); % Вычисление коэффициентов фильтра b=firpm(n+4, f0, a0, w); fvtool(b,1) % загрузка фильтра в среду fvtool АЧХ фильтра Импульсная характеристика фильтра 7 Результаты проектирования с помощью процедуры визуализации фильтров fvtoool выведены в виде набора графиков АЧХ и значений импульсной характеристики (коэффициентов фильтра). Синтез фильтра имеет итерационный характер. Вначале рассчитывается передаточная функция при порядке фильтра 36, найденном по функции оценки порядка firpmord (). После построения АЧХ спроектированного фильтра видим, что АЧХ не удовлетворяет исходным требованиям по необходимому ослаблению в ПЗ. Поэтому в этом примере дважды порядок увеличивается на 2. При порядке n =40 фильтр удовлетворяет спецификации и является окончательным порядком. Коэффициенты фильтра (в формате вывода short): b= Columns 1 through 7 0.0035 0.0029 -0.0124 0.0004 0.0230 -0.0105 -0.0270 Columns 8 through 14 0.0223 0.0172 -0.0198 -0.0008 -0.0101 0.0008 0.0604 Columns 15 through 21 -0.0387 -0.1006 0.1118 0.0966 -0.1860 -0.0400 0.2175 Columns 22 through 28 -0.0400 -0.1860 0.0966 0.1118 -0.1006 -0.0387 0.0604 Columns 29 through 35 0.0008 -0.0101 -0.0008 -0.0198 0.0172 0.0223 -0.0270 Columns 36 through 41 -0.0105 0.0230 0.0004 -0.0124 0.0029 0.0035 Главным преимуществом чебышевских фильтров по сравнению с оконными фильтрами является существенно более низкий порядок фильтра при сопоставимой АЧХ. По сравнению с рассматриваемыми далее оконными фильтрами порядок такого фильтра обычно ниже примерно в два раза. Проектирование КИХ – фильтров методом разложения желаемой частотной характеристики в ряд Фурье (метод оконного взвешивания) Это второй популярный метод проектирования КИХ – фильтров. Другое название – оконный метод или метод оконного взвешивания. Он используется в лаб. работе № 8. j Пусть H d (e ) - желаемая (идеальная) частотная характеристика КИХ–фильтра с нормализо- ванной частотой среза ωС H d (e j 1, 0    C ) , (7) 0, C     Она имеет значение 1 в ПП, и 0 в ПЗ. 1 -π -ωС ωС π ω Напомним, что нормализованная или нормированная частота – это угловая частота Ω =2πF в рад/с, деленная на частоту отсчетов FT, т.е.    / FT  2 F / FT рад. Hd(ejω) Нормированная частота ω -2π -ωС 0 ωС 2π Для ЦФ такая частотная характеристика – периодическая функция ω с периодом 2π для Т = 1 (см. лекцию № 7, для Т ≠ 1 период равен 2π/T). Поэтому она может быть разложена по переменной ω  H d (e j )   h[n]e j n в комплексный ряд Фурье вида (8) n  1  j j n hd [n]  с коэффициентами (9)  H d e e d ,    n   . 2  Обратное ДВПФ (ОДВПФ) от частотной характеристики – это импульсная характеристика фильтра h[n]. Таким образом, коэффициенты ряда Фурье желаемой частотной характеристики дают импульсную характеристику дискретного фильтра.   8 Вычисляя обратное ДВПФ частотной характеристики идеального фильтра импульсную характеристику: , получим его (10) Вид характеристики для нормированной частоты среза -M M Импульсная характеристика симметрична относительно n=0, поэтому фильтр имеет нулевую фазу (нулевую ФЧХ). Но полученный фильтр не является КИХ - фильтром, так как его ИХ длится от n= до n = . Для получения КИХ - фильтра ограничим импульсную характеристику значениями , т.е. просто усечем её M отсчетами слева и справа от 0, как показано выше. . В примере на рис. M= 20. Получаем некаузальную ИХ. Сдвинем импульсную характеристику вправо на отсчетов (точек), чтобы она стала каузальной (физически осуществимой), т.е. начиналась с n = 0. Выражение усеченной и задержанной импульсной характеристики В результате сдвига импульсной характеристики АЧХ фильтра не изменяется, а фазовая характеристика становится равной . Усечение импульсной характеристики можно рассматривать как умножение желаемой характеристики hd[n] на прямоугольное окно wR , где В частотной области произведению соответствует свертка Свертка ω ω В результате свертки идеальной частотной характеристики 𝐻𝑑 (𝑒 𝑗𝜔 ) и частотной характеристики прямоугольного окна W(e jω) (вспомните геометрический смысл свертки!) получается реальная 9 колебания в ПП и ПЗ. Они объясняются явлечастотная характеристика H(e jω). У неё появились нием Гиббса. Фильтр уже не является идеальным. Вид реальной частотной характеристики фильтра с прямоугольным окном для нормированной частоты среза ωС=2 при длине окна N = 8 и N=32. Переходная полоса N=32 Явление Гиббса N=32 N=8 У реальной АЧХ фильтра появляется переходная полоса и пульсации в ПП и ПЗ. Они нежелательны и возникают из-за явления Гиббса. С увеличением размера окна N частотная характеристика фильтра становится все более близкой к идеальной прямоугольной. Но… из-за явления Гиббса пульсации (резкие выбросы) порядка 9% вблизи переходов от полосы пропускания к полосе задерживания будут сохраняться при любом N. Т.е. увеличение порядка фильтра не решает проблему. Для уменьшения величины этих пульсаций на практике вместо прямоугольного окна используют другие более гладкие функции окон, стремящиеся к нулю более плавно, чем прямоугольное окно. Примером такого окна является окно Хэмминга. Здесь для коррекции АЧХ использовано окно Хэмминга. При этом уменьшаются по величине нежелательные пульсации АЧХ в полосе пропускания и полосе задерживания. Характеристика становится плавной без заметных пульсаций в ПП. Но…расширяется переходная полоса фильтра, как «плата» за уменьшение пульсаций. Таблица идеальных импульсных характеристик для стандартных КИХ - фильтров Ниже в табл.1 приведены выражения для идеальных импульсных характеристик цифровых КИХ – фильтров: ФНЧ, ФВЧ, ПФ и РФ. Получение этих выражений здесь опущено. Тип фильтра ФНЧ ФВЧ ПФ РФ Идеальная импульсная характеристика фильтра h[n] (некаузальная ИХ)  C , n  0  h[n]    sin(C n) ,  M  n  M , n  0  n   C , n  0   h[n]    sin(C n) ,  M  n  M , n  0 n   H L , n  0   h[n]    sin(H n)  sin(L n) ,  M  n  M , n  0 n  n   H L , n  0   h[n]    sin(H n)  sin(L n) ,  M  n  M , n  0 n n  Каузальная импульсная характеристика получается сдвигом любой из этих идеальных характеристик h[n] вправо на М отсчетов. Передаточная функция идеального фильтра: H ( z )  b0  b1z 1  b2 z 2  ...  b2 M z 2 M , bn  h[n  M ], n  0,1,..., 2M 10 Оконные функции При проектировании КИХ – фильтров методом оконного взвешивания для улучшения АЧХ идеальная (желаемая) импульсная характеристика фильтра hd[n] умножается на выбранную оконную функцию w[n]: h [n]  hd [n]  w[n] . Название функции w[n] происходит от английского слова window - окно. Желательны окна с возможно более узким главным лепестком частотного спектра и как можно более низкими боковыми лепестками спектра окна. Предложено и известно несколько десятков подобных окон. Примеры наиболее распространенных окон 1. Прямоугольное окно – это исходное окно, оно использовалось выше wR [n]  1,  M  n  M 2. (11) Окно Бартлетта (треугольное) wBart [n]  1  n / M,  M  n  M (12) 3. Окно Хэмминга  n  wham [n]  0,54  0, 46 cos  ,  M  n  M. M  4. 5. (13) Окно Хэннинга  n  whan [n]  0,5  0,5cos  ,  M  n  M. M  (14) Окно Блэкмана  n wB {n]  0, 42  0,5cos  M   2n   0, 08cos    M  ,  M  n  M.  (15) При проектировании КИХ – фильтров часто используется также окно Кайзера из-за его дополнительных свойств. Оно имеет 2 параметра, выбор значений которых обеспечивает более гибкое управление уровнем пульсаций и шириной переходной полосы фильтра. Это окно будет рассмотрено позже. Для формирования оконных функций в MATLAB имеются специальные функции: BARTLETT, BARTHANNWIN, BLACKMAN, BLACKMANHARRIS, BOHMANWIN, CHEBWIN, GAUSSWIN, HAMMING, HANN, KAISER, NUTTALLWIN, RECTWIN, TRIANG, TUKEYWIN. На практике применяют и другие оконные функции. Графики четырех временных окон и их спектров в последовательности: 1. прямоугольное окно, 2. окно Бартлетта (треугольное), 3. окно Хэмминга, 4. окно Блэкмана. Графики получены для длины окна N = 2М+1= 31 командой wvtool(rectwin(31), bartlett(31), hamming(31), blackman(31)). Главный лепесток Боковые лепестки Основными спектральными (частотными) характеристиками окон являются ширина главного лепестка, уровень боковых лепестков и их относительный уровень относительно главного лепестка. Чем шире главный лепесток окна, тем шире переходная полоса фильтра. При увеличении длины окна (порядка фильтра) ширина главного лепестка окна уменьшается, но это не влияет на уровни 11 боковых лепестков. Параметры основных окон, используемые при проектировании фильтров, приведены в таблице2, при этом порядок фильтра M = N-1. Таблица 2. Тип окна Прямоугольное Бартлетта Хэннинга Хэмминга Блэкмана Ширина главного лепестка 4π / M 8π / M 8π / M 8π / M 12π / M Относительный уровень боковых лепестков относительно главного лепестка в дБ -13 -25 -31 -41 -57 Ослабление (пульсации) в полосе задерживания -20lg δS, дБ 21 25 44 53 74 Пульсации в полосе пропусканияδ P, дБ 0,7416 0,0546 0,0194 0,0017 Эквивалентное окно Кайзера, значение β 1,33 3,86 4,86 7,04 Окно Блэкмана обеспечивает самый низкий уровень боковых лепестков по сравнению с другими, приведенными в табл.2, но самый широкий главный лепесток спектра. Это означает, что фильтр с данным окном будет иметь при одинаковом порядке самую широкую переходную полосу, но наибольшее затухание в полосе задерживания. Окна Хэмминга и Хэннинга (другое название – окно Ханна) имеют примерно одинаковую ширину главного лепестка, но окно Хэмминга обеспечивает большее затухание в полосе задерживания, чем окно Хэннинга. Применение окон позволяет уменьшить влияние эффекта Гиббса, заключающегося в появлении нежелательных пульсаций АЧХ фильтра в полосах пропускания и задерживания фильтра. За счет оконного взвешивания достигается снижение пульсаций (колебаний) в полосах пропускания и задерживания. Но … это влечет за собой нежелательные последствия, приводящие к расширению переходной полосы фильтра. Алгоритм проектирования (синтеза) КИХ – фильтра методом взвешивающих окон Некоторые пояснения к алгоритму проектирования:  Проектирование (расчет) такого фильтра заключается в выборе типа окна, определении порядка фильтра, и вычислении его коэффициентов (импульсной характеристики).  В оконных КИХ – фильтрах частота среза fc определяется как среднее арифметическое fpass и fstop: т.е. fc = (fpass + fstop)/2.  Минимальный порядок фильтра М определяется из выражения связи ширины главного лепестка Δωm соответствующего окна с длиной (размером) окна. Например, для прямоугольного окна M  4 / m . Формула для порядка окон Хемминга и Ханна M  8 / m , для окна Блэкмана M  12 / m , здесь Δωm- ширина главного лепестка окна (переходной полосы фильтра).  Пульсации в полосе пропускания в дБ определяются по выражению  P dB  20  lg(1   p ) (см. раздел лекции, относящийся к спецификации фильтра). Ослабление (затухание) в полосе задерживания в дБ  S dB  20  lg( S ) . Алгоритм проектирования КИХ – фильтра с использованием оконного взвешивания: 1. Определить «идеальную» или желаемую частотную характеристику фильтра H D ( ) . 2. С помощью обратного ДВПФ от H D ( ) найти импульсную характеристику hD [n] желаемого фильтра для целевого фильтра. 3. Выбрать тип оконной функции, которая удовлетворяет требованиям к поведению фильтра в полосе пропускания и/ или задерживания. 4. Найти значения коэффициентов h[n] реального оконного фильтра h[n]  hD [n] w[n] , где w[n] - значения весовых коэффициентов выбранного окна. 12 Пример расчета простого оконного фильтра Рассмотрим простейший КИХ - фильтр нижних частот всего с тремя коэффициентами (N=3) с использованием окна Хэмминга. Это демонстрация особенностей расчета оконных фильтров. Пусть этот фильтр имеет частоту среза 500 Гц и частоту дискретизации 5000 Гц. 1. Вычислим нормированную частоту среза фильтра рад. В этом примере N=2M+1=3, М=1. Определим коэффициенты идеального ФНЧ по выражению Отсюда коэффициенты фильтра . Исходя из симметрии импульсной характеристики h[-1]=h[1]=0,1871. Найдем значения отсчетов окна Хэмминга по выражению : . Умножая значения импульсной характеристики на соответствующие значения окна, получим Таким образом, передаточная функция фильтра . Следовательно, уравнение фильтра . Графики АЧХ и ФЧХ фильтра с помощью функции freqz(). Обратите внимание: фазовая характеристика – линейная. Импульсные характеристики Проектирование КИХ – фильтров методом взвешивающих окон в Signal Processing Toolbox MATLAB Для проектирования КИХ – фильтров с линейной фазой методом оконного взвешивания в SPT MATLAB предназначена группа функций fir1() с несколькими вариантами синтаксиса. Здесь используется синтаксис b=fir1(n, Wn, ‘ftype’, window). Функция возвращает вектор b, содержащий (n+1) коэффициентов фильтра с линейной фазой. Параметр n – порядок фильтра. Нормированная на частоту Найквиста частота среза Wn задается числом (для ФНЧ и ФВЧ) или вектором (для ПФ и РФ) с элементами от 0 до 1, где 1 соответствует частоте Найквиста. Для полосового фильтра нормированная Wn задается в виде вектора с двумя элементами Wn=[w1, w2], где w1, w2 - нормированные граничные частоты полосы пропускания. Для режекторного фильтра w1 и w2 – это нормированные граничные частоты полосы задерживания. Для ФВЧ символьный параметр ftype=’high’, для РФ этот параметр ftype=’stop’. Для ФНЧ и ПФ параметр ftype по умолчанию опускается, но может и вводиться. Параметр window определяет тип окна длиной (n+1) отсчетов. Если он не задан, то по умолчанию используется окно Хэмминга. Тип окна выбирается по минимальному значению допустимых пульсаций АЧХ в полосе пропускания δP и задерживания δS, т.е. исходя из значения   min( P ,  S ) . Функций для определения необходимого порядка оконных фильтров, кроме окна Кайзера, в MATLAB нет. Поэтому необходимый порядок оконных фильтров определяется «вручную» по вышеприведенным в таблице 2 выражениям. 13 Пример. Рассчитаем КИХ – фильтр нижних частот с граничной частотой полосы пропускания fP =2000 Гц, граничной частотой полосы задерживания fS = 2400 Гц, максимальным ослаблением в полосе пропускания – AP = 1 дБ, минимальным затуханием в полосе задерживания - AS = 70 дБ и с частотой дискретизации FT = 10 кГц, следовательно, частота Найквиста 5 кГц. Определим необходимые параметры для функции fir1(). Для обеспечения необходимого ослабления в полосе задерживания, равного 70 дБ, выбираем окно Блэкмана, поскольку для этого типа окна пульсации в ПЗ составляют 74 дБ. Это следует из вышеприведенной таблицы параметров оконных функций. Частота среза оконного фильтра – это полусумма ( f P  f S ) / 2  2200 . Нормированная на частоту Найквиста она равна ( f P  f S ) / 2 / 5000  0,44 . Нормированная ширина переходной полосы   2 ( f S  f P ) / FT  0,08 рад. Порядок фильтра с окном Блэкмана определяем по формуле M  12 /   47,74 , т.е. получаем М=48, N=M+1=49. Скрипт – файл MATLAB: b=fir1(48, 0.44, blackman(49)); % Коэффициенты фильтра [h,f]=freqz(b,1,256,10000); % расчет частотной характеристики plot(f, 20*log10(abs(h))), grid % графикАЧХ title('ФНЧ КИХ - фильтр с окном Блэкмана') xlabel('Частота, Гц') figure(2), impz(b,1) % график импульсной характеристики Вид АЧХ спроектированного фильтра. Импульсная характеристика Резюме. Фильтр удовлетворяет исходным требованиям. Достоинства оконного метода: простота и наглядность. Недостатки: невозможно точно задать граничные частоты полосы пропускания и полосы задерживания. Следует понимать, что для большинства применяемых оконных функций максимальная амплитуда пульсаций АЧХ фильтра не зависит от длины (порядка) фильтра. Пульсации зависят от выбора типа окна. 14 Проектирование фильтра с окном Кайзера При проектировании КИХ – фильтров оконным методом желаемая (идеальная) импульсная характеристика фильтра hd[n] умножается на функцию окна . У большинства окон единственным параметром является размер (длина) окна N. В этом смысле одним из наилучших считается окно Кайзера (Kaizer window). Оно имеет специальный параметр β для выбора соотношения между шириной главного лепестка и площадью (уровнем) боковых лепестков. Функция окна Кайзера Вид окна для М= 37 (17) где - функция Бесселя первого рода нулевого порядка. У окна Кайзера имеется два управляемых параметра: длина окна n и параметр формы β. Параметр формы β управляет минимальным ослаблением (затуханием) в полосе задерживания фильтра δS. Этот параметр обычно, но не всегда, выбирается в диапазоне . Увеличение β приводит к увеличению затухания в полосе задерживания, и расширяет полосу пропускания. Приближенная формула для расчета порядка фильтра (𝐴𝑠 −8) 𝑁 = 2,285 , 𝐴𝑆 = −20 lg(𝛿𝑆 ) , 𝛥𝜔 = 𝜔𝑆 − 𝜔𝑃 (18) β = 0 - прямоугольное окно, β = 10 - гауссово окно В MATLAB базовой процедурой для проектирования КИХ – фильтров является рассмотренная ранее fir1(). При использовании окна Кайзера для вычисления порядка фильтра в SPT MATLAB имеется группа функций kaiserord(). Основной синтаксис функции [n, Wn, beta, ftype] =kaiserord(f, a, dev, Ft). Эта функция вычисляет ориентировочный порядок n фильтра, частоту среза Wn и значение коэффициента затухания beta. Параметр ftype – тип фильтра. Параметр Ft – это частота отсчетов. Вектор f при этом содержит граничные частоты полос фильтра, вектор a – желаемые значения АЧХ полос, вектор dev определяет максимально допустимые отклонения АЧХ синтезированного фильтра от желаемой АЧХ. Размер вектора dev совпадает с размером вектора a. После расчета порядка фильтра используется функция fir1() с полученными в kaiserord() параметрами. В результате определяются коэффициенты фильтра h, удовлетворяющие условиям, заданным для функции kaiserord(). При необходимости порядок фильтра уточняется. Код MATLAB для проектирования фильтра: Fs=10000; % частота дискретизации fcuts=[700 1000 2000 2300]; % граничные частоты mags=[0 1 0]; % желаемые значения АЧХ devs=[0.01 0.05 0.01]; % допустимые отклонения в ПЗ и ПП % проектирование фильтра [n,Wn,beta,ftype]=kaiserord(fcuts,mags,devs,Fs); h=fir1(n,Wn,ftype,kaiser(n+1,beta),'noscale'); % Вывод АЧХ и ИХ [H,f]=freqz(h,1,1024,Fs); subplot(1,2,1), plot(f,abs(H)), grid title('АЧХ фильтра') xlabel('Частота, Гц') subplot(1,2,2), impz(h), grid title('Импульсная характеристика') xlabel('Номер отсчета') ylabel('Амплитуда') 15 В качестве примера рассчитаем полосовой КИХ – фильтр с граничными частотами ПП 1000 и 2000 Гц и граничными частотами ПЗ 700 и 2300 Гц, частотой отсчетов 10000 Гц. Допустимое отклонение АЧХ в ПП 0,05 (20*lg(1.05) = 0.4238 дБ ≈ 0,4 дБ), затухание в ПЗ 0,01 (40 дБ). Порядок спроектированного фильтра n = 75. Структуры реализации КИХ – фильтров После проектирования (синтеза) фильтра необходимо выбрать структуру его реализации. 1. Основной структурой является прямая форма реализации КИХ – фильтра. Она реализуется по основному уравнению КИХ – фильтра, т.е. уравнению дискретной свертки . Здесь M– порядок фильтра. Пример структуры для M = 4 x[n] h[0] Прямую форму называют также трансверсальным фильтром (transversal filter) и линией задержки с отводами (tapped delay line). КИХ – фильтр порядка N при прямой реализации имеет N+1 коэффициентов и N+1 умножителей. Эти фильтры более устойчивы к ошибкам квантования и округления, чем БИХ - фильтры. Поэтому прямая форма в большинстве случаев достаточна для практических целей. 2. Второй формой структуры КИХ – фильтра порядка N является каскадная, состоящая из каскада звеньев второго порядка и, возможно, дополнительного звена первого порядка Для разложения H(z) на квадратичные звенья (секции) можно применить функцию tf2sos () MATLAB. Структура отдельного квадратичного (биквадратного) звена каскадной схемы. Пример. >> b=fir1(7,0.35) % Расчет фильтра 7 – го порядка b = -0.0051 0.0134 0.1474 0.3444 0.3444 0.1474 0.0134 -0.0051 Преобразование в каскадную структуру >>sos=tf2sos(b,1) sos = -0.0051 0.0393 0 1.0000 1-й каскад 1.0000 0.8695 -0.1305 1.0000 0 2-й каскад 1.0000 2.9794 2.4732 1.0000 3-й каскад 1.0000 1.2047 0.4043 1.0000 4-й каскад Каскадная форма реализации КИХ – фильтров также часто используется на практике. 3. Третьей формой является структура для фильтров с линейной фазой, учитывающая симметрию импульсной характеристики КИХ - фильтра с линейной ФЧХ. Пример такого фильтра 6 – го порядка Эта структура требует 4-х умножений вместо 6. Для больших N возможна экономия почти в 2 раза. 4. Наконец, для ускорения вычислений при больших N может использоваться быстрая свертка или свертка в частотной области (4-я форма), основанная на применении БПФ. Схема БПФ вычисление свертки в частотной области требует операций умножений, при этом L = 2 m≥ M+N -1 , где M,N – размеры векторов x[n], h[n]. h[n] ]] (см. лабораторную работу № 2). 16