Тема нерекурсивные частотные цифровые фильтры недостаточно овладеть премудростью, нужно уметь пользоваться ею

ЦИФРОВАЯОБРАБОТКАСИГНАЛОВ

Digital signals processing. Digital nonrecursive frequency filters.

Тема 7. НЕРЕКУРСИВНЫЕ ЧАСТОТНЫЕ ЦИФРОВЫЕ ФИЛЬТРЫ

Недостаточно овладеть премудростью, нужно уметь пользоваться ею.

Марк Туллий Цицерон. О высшем благе и высшем зле.

Римский сенатор и философ, 1 в.д.н.э.

Мало пользы от теории бокса, пока сам не научишься махать кулаками.

Евгений Буцко. Идеология белых воротничков.

Радиоинженер, геофизик Уральской школы, ХХ в.

Содержание

1. Общие сведения. Типы фильтров. Методика расчетов нерекурсивных цифровых фильтров. Фильтры с линейной фазовой характеристикой.

2. Идеальные частотные фильтры. Импульсная реакция фильтров.

3. Конечные приближения идеальных фильтров. Ограничение окна операторов фильтров. Применение весовых функций для нейтрализации явления Гиббса. Основные весовые функции. Весовая функция Кайзера.

4. Гладкие частотные цифровые фильтры. Принцип синтеза фильтров.

5. Дифференцирующие цифровые фильтры. Передаточная функция. Точность дифференцирования. Применение весовых функций. Фильтры с линейной групповой задержкой.

6. Альтернативные методы расчета НЦФ. Оптимизационные методы. Метод частотной выборки.

Введение

Нерекурсивные фильтры реализуют алгоритм свертки двух функций: yk = hn ③ xk-n, где xk – массив входных данных фильтра, hn – оператор (ядро, импульсный отклик) фильтра, k и n – нумерация числовых значений массива данных и числовых значений коэффициентов фильтра, k = 0, 1, 2, … ,K; n = 0, 1, 2, … ,N; K ≥ N. Значения выходных отсчетов свертки yk для любого аргумента k определяются текущим и "прошлыми" (до k-N) значениями входных отсчетов. Такой фильтр называется нерекурсивным цифровым фильтром (НЦФ). Интервал [0-N] оператора получил название "окна" фильтра. Окно фильтра составляет N+1 отсчет, фильтр является односторонним каузальным, т.е. причинно обусловленным текущими и "прошлыми" значениями входного сигнала, и выходной сигнал не опережает входного. В общем случае, каузальный фильтр меняет в спектре сигнала состав гармоник, их амплитуды и фазы.

Каузальный фильтр может быть реализован физически в реальном масштабе времени. Начало фильтрации возможно только при задании определенных начальных условий – N значений отсчетов для точек x(k-n) при k<n. Как правило, в качестве начальных условий принимаются нулевые значения, тренд сигнала или значения отсчета х(0), т.е. продление отсчета x(0) назад по аргументу.

При обработке данных на ЭВМ ограничение по каузальности снимается. В программном распоряжении фильтра могут находиться как "прошлые", так и "будущие" (k+n, до k+N') значения входной последовательности отсчетов относительно текущей точки вычислений k, при этом для завершения свертки (аналогично началу) требуется N' точек конечных условий при (k+n)>K. При N' = N и h(-n) = h(n) фильтр называется двусторонним симметричным фильтром. Симметричные фильтры, в отличие от односторонних, не изменяют фазы обрабатываемого сигнала.

7.1. Общие сведения.

Основное свойство любого фильтра – его частотная (frequency response) и фазовая характеристики. Они показывают, какое влияние фильтр оказывает на амплитуду и фазу различных гармоник обрабатываемого сигнала.

К наиболее известным типам нерекурсивных цифровых фильтров (НЦФ) относятся частотные фильтры, алгоритм которых для симметричных НЦФ, не изменяющих фазу сигналов, имеет вид:

yk =hn sk-n.

Типы фильтров. В зависимости от вида частотной характеристики выделяют три основных группы частотных фильтров: ФНЧ - фильтры низких частот (low-pass filters) - пропускание низких и подавление высоких частот во входном сигнале, ФВЧ - фильтры высоких частот (high-pass filters) - пропускание высоких и подавление низких частот, и ПФ - полосовые фильтры, которые пропускают (band-pass filters) или подавляют (band-reject filters) сигнал в определенной частотной полосе. Среди последних в отдельную группу иногда выделяют РФ - режекторные фильтры, понимая под ними фильтры с подавлением определенной гармоники во входном сигнале, и СФ – селекторные фильтры, обратные РФ.

Если речь идет о подавлении определенной полосы частот во входном сигнале, то такие фильтры называют заградительными. Ни теоретического, ни практического интереса к методам их расчета обычно не проявляется, так как их частотная характеристика обычно задается инверсией характеристики полосового фильтра (1-Hп()) и каких-либо дополнительных особенностей при своем проектировании не имеет.

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

Рис. 7.1.1. Типы основных частотных фильтров.

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

Рис. 7.1.2.

Так, например, фильтр высоких частот g(n) может быть получен инверсией фильтра низких частот h(n) - вычислением разности между исходным сигналом и результатом его фильтрации низкочастотным НЦФ:

y(k) = s(k) –h(n) s(k-n).

Отсюда, условие инверсии симметричного низкочастотного фильтра в высокочастотный:

g(0) = 1-h(0), g(n) = -h(n) при n0.

Пример обращения и спектры фильтров приведены на рис. 7.1.2 (в правой части главных диапазонов).

Рис. 7.1.3.

Применяется также способ получения фильтров высоких частот из низкочастотных фильтров путем реверса частоты в передаточной функции низкочастотного фильтра, т.е. заменой переменной  на переменную ' = (при t = 1). Для симметричных фильтров, содержащих в передаточной функции только косинусные члены аргумента , в результате такой операции будем иметь:

cos n(-) = cos n cos n = (-1)n cos n.

Последнее означает смену знака всех нечетных гармоник передаточной характеристики фильтра и, соответственно, всех нечетных членов фильтра:

g(n) = h(n) при n = ±1, ±3, …

Пример частотного реверса приведен на рис. 7.1.3. Физическую сущность такой операции инверсии спектра легко понять на постоянной составляющей сигнала. При изменении на противоположный знака каждого второго отсчета постоянной величины это постоянной значение превращается в "пилу", частота которой равна частоте Найквиста главного частотного диапазона (отсчеты по амплитудным значениям этой частоты), равно как и наоборот, отсчеты гармоники сигнала на частоте Найквиста (знакочередующиеся в силу сдвига по интервалам дискретизации на ) превращаются в постоянную составляющую.

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

vk = hн(n) ③ s(k-n), yk = hв(n) ③ vk = hн(n) ③ hв(n) ③ s(k-n).

Так как операция свертки коммутативна, то вместо отдельных массивов коэффициентов ФНЧ и ФВЧ их сверткой может быть определен непосредственно массив коэффициентов полосового фильтра: hn = hн(n) ③ hв(n).

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

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

Методика расчетов НЦФ. Обычно при фильтрации сигналов задается требуемая частотная характеристика фильтра. Задачей является построить фильтр, отвечающий заданным требованиям и провести фильтрацию. Зачастую бывает невозможно построить в точности заданный фильтр, и выполняется фильтр, близкий по характеристикам к заданному.

Существует много способов построения фильтров с заданной частотной характеристикой. Наиболее простой из них – проектирование фильтров с линейной фазой с помощью весовых окон. Этот способ является универсальным и позволяет получить фильтр с любой заданной частотной характеристикой. Отметим, однако, что с помощью других, математически более строгих и совершенных методов, иногда удается построить фильтр меньшей длины, удовлетворяющий тем же требованиям к частотной характеристике.

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

1. Задание идеальной амплитудно-частотной характеристики передаточной функции фильтра. Термин идеальности понимается здесь в том смысле, что на характеристике указываются полосы пропускания и подавления частот с коэффициентами передачи 1 и 0 соответственно без переходных зон.

2. Расчет функции импульсного отклика идеального фильтра (обратное преобразование Фурье частотной характеристики фильтра). При наличии скачков функций на границах пропускания/подавления импульсный отклик содержит бесконечно большое количество членов.

3. Ограничение функции отклика до определенного количества членов, при этом на передаточной характеристике фильтра возникает явление Гиббса – осцилляции частотной характеристики с центрами на скачках.

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

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

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

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

(7.1.1)

Оно выполняется, если импульсная характеристика фильтра имеет положительную симметрию:

h(n) = h(N-n-1), n = 0, 1, 2, …, (N-1)/2, N – нечетное (тип 1);

n = 0, 1, 2, …, (N/2)-1, N – четное (тип 2).

При этом фазовая характеристика будет определяться длиной фильтра:

(N-1)/2.

Частотная характеристика фильтра:

H() = |H()| exp(j()), (7.1.2)

где модуль |H()| задается аналогично АЧХ симметричных фильтров. Следует также учитывать, что частотную характеристику типа 2 нельзя использовать для проектирования фильтров верхних частот, т.к. она всегда равна нулю на частоте Найквиста.

Собственно методика расчета каузальных фильтров, за исключением использования (7.1.2) для задания частотной характеристики, не отличается от методики расчета симметричных фильтров, включая необходимость использования весовых функций для нейтрализации явления Гиббса. Это позволяет применять чисто практический метод расчетов – вычислить и отработать сначала симметричный фильтр на N-точек (тип 1), а затем превратить его в каузальный сдвигом вправо на (N-1)/2 точек в область только положительных значений n ≥ 0.

7.2. Идеальные частотные фильтры.

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

Импульсная реакция фильтра (коэффициенты оператора) находится обратным преобразованием Фурье заданной передаточной функции H(). В общем случае:

h(nt) = (1/2)H() exp(jnt) d

Для получения вещественной функции импульсного отклика фильтра действительная часть передаточной функции должна быть четной, а мнимая - нечетной. Цифровые фильтры задаются в главном частотном диапазоне, границы которого (частота Найквиста N) определяются интервалом дискретизации данных (N = /t), подлежащих фильтрации, и соответственно определяют интервал дискретизации оператора фильтра (t = /N). Для фильтров с нулевым фазовым сдвигом мнимая часть передаточной функции должна быть равна нулю, при этом оператор фильтра определяется косинусным преобразованием Фурье:

h(nt)= (1/)H() cos(n/N) dn = 0, 1, 2,... (7.2.1)

Для идеального полосового фильтра H()=1 в полосе частот от н до в, и интеграл (7.2.1) вычисляется в этих пределах. Идеальные фильтры низких и высоких частот, как частные случаи идеальных ПФ, интегрируются в диапазоне от 0 до в для низкочастотного и от н до N для высокочастотного фильтра.

При интервале дискретизации данных t, условно принимаемым за 1, главный частотный диапазон передаточных функций ограничивается значением частоты Найквиста от - до . Если на практике интервал дискретизации данных в физических единицах отличается от 1, то это сказывается только на изменении масштаба частотной шкалы передаточных функций.

Пример 1. t = 0.1 сек. fN = 1/2t = 5 Гц. N =/t = 10 .

Пример 2. x = 10 метров. fN = 0.05 м-1. N= 0.1 .

Во всех дальнейших выражениях значение t, если это специально не оговорено, будем принимать равным 1.

При H()=A=1 в полосе пропускания (н, в), и H()=0 за ее пределами, для идеальных симметричных полосовых НЦФ из (7.2.1) с границами интегрирования, соответственно, от н до в в общем виде получаем:

h(n) = (А/) [в sinc(nв) - н sinc(nн)], (7.2.2)

ho = (в - н)/, h(n) = (sin nв - sin nн)/(n).

где sinc(n) = sin(n)/(n) - функция интегрального синуса (функция отсчетов), бесконечная по координате .

При инверсии частотной характеристики в заградительный фильтр:

ho = (1-(н - в))/, h(n) = (sin nн - sin nв)/(n).

Рис. 7.2.1. Входные сигналы. Рис. 7.2.2. Спектр сигнала и границы фильтра.

На рис. 7.2.1 приведен пример сигнала однотональной балансной амплитудной модуляции (чистого сигнала – вверху, с наложенными шумами внизу, мощность шумов равна мощности сигнала). Если информация заключена в частоте и амплитуде модулирующего сигнала, то полосовой фильтр выделения сигнала из шумов, спектр которого для одной модулирующей частоты приведен на рис. 7.2.2, в идеальном случае должен иметь плоскую частотную характеристику в границах возможных вариаций модулирующей частоты (от н до в).

Размер оператора фильтра определяется приблизительно из следующих соображений. Чем больше размер оператора, тем круче будет переходная зона и меньше ее размер, т.е. тем ближе будет фактически реализованная передаточная функция фильтра к идеальной. Обычно сначала стоит попробовать построить фильтр достаточно большого размера, оценить его соответствие заданной частотной характеристике и в дальнейшем попытаться уменьшить. Значение N для симметричных НЦФ должно быть нечетным числом.

Рис. 7.2.3. Оператор фильтра.

На рис. 7.2.3 приведен оператор полосового фильтра, вычисленный по (7.2.2) для приведенных выше условий, с ограничением по числу коэффициентов оператора до N=100. Как видно из рисунка, оператор затухает достаточно медленно и явно усечен, что должно сказаться на форме частотной характеристики фильтра. Все дальнейшие вычисления будут проводиться на продолжении данного примера.

7.3. Конечные приближения идеальных фильтров /24/.

Оператор идеального частотного НЦФ, как это следует из выражения (7.2.2), представляет собой бесконечную затухающую числовую последовательность, реализующую заданную передаточную функцию:

H() =h(n) cos n. (7.3.1)

Ограничение окна операторов фильтров. На практике бесконечный ряд (7.3.1) всегда приходится ограничивать определенным числом членов его конечного приближения

H'() =h(n) cos n,

при этом передаточная функция осложняется явлением Гиббса, и появляется переходная зона между полосами пропускания и подавления сигнала (рис. 7.3.1, пунктирная кривая при N=100). Явление Гиббса формирует первые выбросы передаточной функции на расстоянии /(2(N+1)) от скачков (разрывов первого рода). Если ширину переходной зоны p в первом приближении принять по расстоянию между первыми выбросами по обе стороны от скачка функции H(), то ее значение будет ориентировочно равно /(N+1) =p

Рис. 7.3.1. Передаточные функции полосового фильтра.

Применение весовых функций. Если уровень пульсаций передаточной функции, определяемый явлением Гиббса, не удовлетворяет поставленным задачам фильтрации данных, рекомендуется использование сглаживающих весовых функций. С учетом того, что при применении весовых функций происходит расширение переходных зон примерно в два раза, значение ширины переходной зоны будет равным p = 2/N. Отсюда можно определить минимальное число членов усеченного ряда по заданному размеру переходной зоны:

N = 2/p. (7.3.2)

Для примера на рис. 7.3.1 значение N принято равным 200, при этом крутизна переходной зоны увеличилась (тонкая кривая H'(), N=200), создавая запас на последующее сглаживание весовой функцией.

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

hn = pn h(n). (7.3.3)

Рис. 7.3.2. Полосовая фильтрация (вверху – входной сигнал, внизу – выходной).

Подстановкой коэффициентов (7.3.3) в (7.3.1) рекомендуется произвести построение полученной передаточной характеристики фильтра и непосредственно по ней оценить пригодность фильтра для поставленных задач. Это наглядно видно на рис. 7.3.1, где для нашего примера была применена весовая функция Гаусса. Передаточная функция Hp() имеет практически такую же крутизну, как и функция H'() при N=100 и практически плоскую вершину в интервале спектра сигнала. Качество работы фильтра для сигнала, приведенного на рис. 7.2.1, можно видеть на рис. 7.3.2.

При необходимости более точной оценки полученной передаточной функции можно рекомендовать увеличение ее частотного разрешения в 2-4 раза перед выполнением преобразования Фурье, что можно выполнить путем увеличения размеров оператора hn дополнением нулями.

Основные весовые функции. Ниже в таблицах приведены формулы и основные спектральные характеристики наиболее распространенных весовых окон. Носители весовых функций, в принципе, являются неограниченными и при использовании в качестве весовых окон действуют только в пределах окна и обнуляются за его пределами. Для упрощения записи формулы приводятся в аналитической форме с временным окном 2, симметричным относительно нуля (0). При переходе к дискретной форме окно 2 заменяется окном 2N+1, а значения t – дискретами t = nt. Большинство весовых функций на границах окна (n = N) принимают нулевые или близкие к нулевым значения. Последнее исключается, если принять 2= (2N+3)t, при этом близкие к нулю значения перемещаются за границы окна.

Основные весовые функции.

Временное окно	Весовая функция	Фурье-образ
Естественное (П)	П(t) = 1, |t|П(t) t	П() = 2 sinc[]
Бартлетта ()	b(t) = 1-|t|/	B() =  sinc2(/2).
Хеннинга, Ганна	p(t) = 0.5[1+cos(t/)]	0.5П()+0.25П(+/)+0.25П(-/)
Хемминга	p(t) = 0.54+0.46 cos(t/)	0.54П()+0.23П(+/)+0.23П(-/)
Карре (2-е окно)	p(t) = b(t) sinc(t/)	·B()*П(), П() = 1 при ||</
Лапласа-Гаусса	p(t) = exp[-2(t/)2/2]	[(/) exp(-22/(22))] ③ П()
Кайзера-Бесселя	p(t) =, Jo[x] =[(x/2)k/k!]2	Вычисляется преобразованием Фурье. Jo[x] - модифицированная функция Бесселя нулевого порядка

Характеристики спектров весовых функций.

Параметры	Ед. изм.	П- Окно	Барт- летт	Лан-цош	Хен- нинг	Хемминг	Кар- ре	Лаплас	Кайзер
Амплитуда: Главный пик 1-й выброс(-) 2-й выброс(+) Ширина Гл. пика Положения: 1-й нуль 1-й выброс 2-й нуль 2-й выброс	 %Гл.п. - “ - / / / / /	2 0.217 0.128 0.60 0.50 0.72 1.00 1.22	1 - 0.047 0.89 1.00 - - 1.44	1.18 0.048 0.020 0.87 0.82 1.00 1.29 1.50	1 0.027 0.0084 1.00 1.00 1.19 1.50 1.72	1.08 0.0062 0.0016 0.91 1.00 1.09 1.30 1.41	0.77 - - 1.12 - - - -	0.83 0.0016 0.0014 1.12 1.74 1.91 2.10 2.34	0.82 .00045 .00028 1.15 1.52 1.59 1.74 1.88

Весовая функция Кайзера. Наибольшее распространение при расчетах частотных НЦФ получила весовая функция Кайзера:

p(n) = .

Это объясняется тем, что параметры функции Кайзера могут устанавливаться непосредственно по техническим требованиям к передаточным функциям проектируемых фильтров – допустимой ширине переходной зоны p и значению коэффициента шума фильтра (максимальным значениям осцилляций передаточной функции в единицах коэффициента передачи в полосе пропускания).

Кайзером установлено, что для заданного значения  произведение количества членов оператора НЦФ на ширину переходной зоны является величиной постоянной. Оно получило название D-фактора:

D = N·p/.

С другой стороны, установлены следующие эмпирические соотношения между D-фактором и параметром  функции Кайзера:

D = (А-7.95)/14.36 при А>21.

= 0.9222 при А<21.

 = 0.1102(A-8.7) при А>50.

= 0 при А<21.

= 0.5842(A-21)0.4 + 0.07886(A-21), 21<А<50.

где: А = -20 log  - затухание в децибелах.

Приведенные выражения позволяют по заданному значению коэффициента шума  определить параметр  функции Кайзера, а через D-фактор число членов фильтра:

N = D/p.

При проектировании полосовых фильтров проверка передаточной функции полученного оператора НЦФ исходному заданию по значению коэффициента шума является обязательной. Это объясняется тем, что поскольку полоса пропускания полосового фильтра ограничена двумя скачками, на передаточной характеристике возникают два центра осцилляций, при этом наложение осцилляций может как уменьшить, так и увеличить амплитуду суммарных осцилляций. Если за счет наложения произойдет увеличение амплитуды осцилляций, то расчет НЦФ следует повторить с уменьшением исходного значения .

Пример расчета полосового фильтра.

Произвести расчет ПФ при следующих исходных параметрах: н = 0.3, в = 0.6, p = 0.1,  = 0.02.

1. А = -20 log . А = 34.

2. N =  (A-7.95)/(14.36 p). N = 18.

3.  = 0.5842(A-21)0.4 +0.07886(A-21).= 2.62.

4. hо = (в-н)/. hо = 0.3

5. h(n) = (sin nв-sin nн)/(n). h(n)= 0.04521, -0.24490, -0.09515, ... , 0.02721.

6. pn= Jo{} / Jo{}. pn = 1.00, 0.997, 0.9882, .......

7. Оператор фильтра: hn = pn h(n), n = 0, 1, 2,..., N. h-n = hn. hn = 0.3000, 0.04508, -0.2420, ....

8. Проверка по формуле: H() =hn cos n, 0  .

Для оценки формы передаточной функции количество точек спектра в интервале 0- достаточно задать равным 2N, т.е. с шагом /36.

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

Ошибки отклонения () частотной характеристики относительно заданной при проектировании кроме разрядности коэффициентов В в битах зависит также от размеров N оператора фильтра и в первом приближении может оцениваться по формулам:





ln N)/3

Выражение (7.3.4) наиболее пессимистично и предполагает наихудшие ситуации вычислений. Два других выражения носят более реальный характер по статистическим данным.

Курсовая работа 8-07. Разработка программы расчета универсального частотного НЦФ (низкочастотный, высокочастотный, полосовой) и его использования при обработке цифровых сигналов.

7.4. Гладкие частотные фильтры /24/.

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

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

Допустим, мы имеем симметричный НЦФ с передаточной функцией:

H() = hо+2hn cos n. (7.4.1)

Как известно, cos n равен полиному по cos  степени n, при этом выражение (7.4.1) можно записать в виде:

H() =gn (cos )n =gn xn, (7.4.2)

где переменная х=cos  изменяется от 1 до -1 (поскольку  изменяется от 0 до ). Преобразование переменной представляет собой нелинейное растяжение оси абсцисс с поворотом на 180o (по переменной х передаточные функции ФНЧ похожи на ФВЧ, и наоборот) с выражением функции через степенной полином. Последнее примечательно тем, что синтез гладких функций на базе степенных полиномов затруднений не представляет.

Так, например, для конструирования ФНЧ в качестве исходной может быть принята степенная функция вида:

g(x)= (1+x)z (1-x)r, (7.4.3)

где z и r - параметры.

Рис. 7.4.1. Примеры синтеза гладких фильтров.

Функция (7.4.3) имеет нули порядка z и r в точках соответственно х = -1 и х = 1 (рис. 7.4.1), при этом значения параметров z и r характеризуют степень касания функцией оси абсцисс. Чем больше порядок, тем медленнее функция отходит ("отрывается") от оси абсцисс.

Если выражение функции (7.4.3) проинтегрировать в пределах от -1 до х и нормировать на значение интеграла от -1 до 1 , то будет получена гладкая передаточная характеристика низкочастотного фильтра. На рисунке 7.4.1 приведены передаточные функции для двух пар параметров z и r, вычисленные по формуле:

H(x)=g(x)dx /g(x)dx. (7.4.4)

Рис. 7.4.2. Схема возврата к ряду Фурье.

Функция H(x) имеет перегиб в точке (z-r)/(z+r) и переходную зону, крутизна которой тем больше, чем больше значения z и r. Подстановкой x=cos  осуществляется возврат к частотной переменной  с сохранением монотонности функции.

В заключение, для определения коэффициентов фильтра hn требуется осуществить обратное преобразование от степенной формы (7.4.2) к ряду Фурье (7.4.1). Выполнение данной операции достаточно просто производится рекурсивным способом, показанным на рис. 7.4.2. Подробное обоснование рекурсии приведено в /24/.

Пример расчета гладкого фильтра.

Произвести расчет ФНЧ с гладкой частотной характеристикой с перегибом характеристики в точке /3. За исходную функцию принять функцию (7.4.3).

1. x= cos(/3)= 0.5= (z-r)/(z+r). Принято: z=3, r=1.

Исходный многочлен: g(x) = (1-x)(1+x)3 = 1+2x-2x3-x4.

2. H(x)=g(x)dx = C+x+x2-0.5 x4-0.2 x5. При х= -1, H(-1)= 0, откуда С=0.3. При х=1, H(1)=1.6.

Отсюда: H(x)= (3+10x+10x2-5x4-2x5)/16. gn = {3/16, 10/16, 10/16, 0, -5/16, -2/16}.

3. Применяя рекурсивное преобразование, получаем: hn= {(98, 70, 20, -5, -5, -1)/256}.

Для расчетов гладких фильтров высоких частот в выражении (7.4.4) достаточно поменять местами пределы интегрирования. Гладкие полосовые фильтры получаются комбинацией ФНЧ и ФВЧ с перекрытием частот пропускания.

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

7.5. Дифференцирующие цифровые фильтры.

Передаточная функция. Из выражения для производной d(exp(jt))/dt = j exp(jt)

следует, что при расчете фильтра производной массива данных необходимо аппроксимировать рядом Фурье передаточную функцию вида H() = j. Поскольку коэффициенты такого фильтра будут обладать нечетной симметрией (h-n = -hn) и выполняется равенство

hn [exp(jn)-exp(-jn)] = 2j hn sin n,

то передаточная характеристика фильтра имеет вид:

H() = 2j(h1 sin + h2 sin 2+ ... + hN sin N),

т.е. является мнимой нечетной, a сам фильтр является линейной комбинацией разностей симметрично расположенных относительно sk значений функции. Уравнение фильтрации:

yn =hn(sk+n - sk-n).

Если дифференцированию подлежит низкочастотный сигнал, а высокие частоты в массиве данных представлены помехами, то для аппроксимации в пределах главного частотного диапазона задается (без индекса мнимости) передаточная функция фильтра вида:

H() = в, H() = 0, в< N.

Оператор дифференцирующего фильтра:

h(n) = (2/)H() sin(n/N) dn = 0,1,2,... (7.5.1)

Принимая, как обычно, N =  (t = 1) и решая (7.5.1) при H() = , получаем:

hn = (2/)[sin(nв)/n2 - в cos(nв)/n], (7.5.2)

hо = 0, h-n = -hn.

Частотная характеристика:

Im(H()) =hn sin n = 2hn sin n (7.5.3)

Точность дифференцирования. На рис. 7.5.1 приведен пример расчета коэффициентов дифференцирующего фильтра на интервал частот {0-0.5} при t=1 (в = ). Операторы дифференцирующих фильтров, как правило, затухают очень медленно и, соответственно, достаточно точная реализация функции (7.5.3) весьма затруднительна.

Рис. 7.5.1. Коэффициенты оператора фильтра.

Ряд (7.5.3) усекается до N членов, и с помощью весовых функций производится нейтрализация явления Гиббса. Явление Гиббса для дифференцирующих фильтров имеет весьма существенное значение, и может приводить к большим погрешностям при обработке информации, если не произвести его нейтрализацию. Примеры ограничения оператора, приведенного на рис. 7.5.1, и соответствующие передаточные функции H'() усеченных операторов показаны на рис. 7.5.2.

Для оценки возможных погрешностей дифференцирования усеченными операторами произведем расчет фильтра при в = . По формулам (7.5.2) определяем:

h0-10 = 0, 0.3183, 0.25, -0.0354, -0.125, 0.0127, 0.0833, -0.0065, -0.0625, 0.0039, 0.05.

Рис. 7.5.2. Частотные функции фильтров.

Произведем проверку работы фильтра на простом массиве данных sn = n, производная которого постоянна и равна 1. Для массива с постоянной производной фильтр может быть проверен в любой точке массива, в том числе и в точке n=0, для которой имеем:

у =hn so-n = 2n hn,

при этом получаем: у=0.5512 при N=5, у=1.53 при N=10.

Рис. 7.5.3. Погрешность дифференцирования.

Такое существенное расхождение с действительным значением производной объясняется тем, что при =0 тангенс угла наклона реальных передаточных функций фильтра, как это видно на рисунке 7.5.2, весьма существенно отличается от тангенса угла наклона аппроксимируемой функции H() = . На рис. 7.5.3 приведены частотные графики относительной погрешности дифференцирования  = Hн'()/Hн() с вычислением значений на нулевой частоте по пределам функций при N → ∞На рис. 7.5.4 приведен пример операции дифференцирования гармоники s с частотой  оператором с N=10 в сопоставлении с точным дифференцированием ds/dk.

Рис. 7.5.4. Пример операции дифференцирования.

Применение весовых функций. Применим для нейтрализации явления Гиббса весовую функцию Хемминга. Результат нейтрализации для фильтра с N=10 приведен на рис. 7.5.5. Повторим проверочный расчет дифференцирования на массиве sn = n и получим результат у=1.041, т.е. погрешность дифференцирования уменьшается порядок.

Рис. 7.5.5. Дифференцирование с применением весовой функции.

Аналогично производится расчет и полосовых дифференцирующих фильтров с соответствующим изменением пределов интегрирования в (7.5.1) от н до в. При этом получаем:

hn = (нcos nн-вcos nв)/(n) + (sin nв-sin nн)/(n2).

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

(7.5.4)

где иконстанты.

Оно выполняется, если импульсная характеристика фильтра имеет положительную симметрию:

h(n) = -h(N-n-1), n = 0, 1, 2, …, (N-1)/2, N – нечетное (тип 1);

n = 0, 1, 2, …, (N/2)-1, N – четное (тип 2).

При этом фазовая характеристика будет определяться длиной фильтра:

(N-1)/2,  = /2.

Частотная характеристика фильтра:

H() = |H()| exp(j()), (7.5.4)

где модуль |H()| задается нечетным. Оба типа фильтров вводят в выходной сигнал сдвиг фазы на 90о. Кроме того, частотная характеристика фильтра типа 1 всегда равно нулю на частоте Найквиста, что определяется знакопеременностью левой и правой части главного диапазона спектра с учетом периодизации спектра дискретных функций.

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

7.6. АЛЬТЕРНАТИВНЫЕ МЕТОДЫ РАСЧЕТА нцф [43].

Метод прямого расчета НЦФ по частотной характеристики понятен и прост для применения. Недостаток метода – отсутствие гибкости. Он не позволяет проектировать фильтры с разной степенью неравномерности частотной характеристики в полосах пропускания и подавления, а степень неравномерности не зависит от количества членов фильтра и не может изменяться. Максимальные осцилляции частотной характеристики всегда наблюдаются в области полосовых границ и уменьшаются при удалении от них, но при близких границах могут наблюдаться явления интерференции осцилляций. Более гибкими в проектировании являются альтернативные методы: оптимизационные,

Оптимизационные методы позволяют проектировать экономные по размерам операторы фильтров с оптимальными (по Чебышеву) осцилляциями частотных характеристик. Они основаны на понятии полос равных колебаний.

Рис. 7.6.1. Оптимальный фильтр низких частот

Частотная характеристика оптимального фильтра низких частот приведена на рис. 7.6.1. В полосе пропускания реальная характеристика фильтра осциллирует с постоянными амплитудными колебаниями между значениями 1-p и 1+p. В полосе подавления осцилляции постоянной амплитуды находятся в интервале 0-s. Разность между идеальной и практической характеристиками представляет собой функцию ошибок E(f). Оптимальный метод позволяет определить коэффициенты фильтра h(n), для которых значение максимальной взвешенной ошибки минимизируется

min[max(E(f))]

в полосе пропускания и в полосе подавления, при этом характеристика фильтра будет иметь равные колебания в пределах полос пропускания и подавления, а количество экстремумов колебаний у фильтров с линейной фазовой характеристикой обычно прямо связано с количеством коэффициентов фильтра (N+1)/2.

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

Метод частотной выборки представляет собой вариант метода расчета фильтра по частотной характеристике без применения весовых функций и может применяться для расчетов как частотно-избирательных фильтров, так и фильтров с произвольной частотной характеристикой.

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

Рис. 7.6.2. Задание параметров НЦФ.

Допустим, нам требуется достаточно простой симметричный низкочастотный фильтр с шириной переходной зоны порядка 0.2 главного частотного диапазона (при k=1 для фильтра, fN = 0.5 Гц для спектра и ширина переходной зоны 0.2 х 0.5 = 0.1 Гц). Минимальный размер фильтра при идеальной характеристике для обеспечения такого перехода 2N+1 = 2(1+1/0.1) = 11 точек. С учетом расширения переходной зоны при уменьшении осцилляций на границе зон примем для начала N=8. Частотная характеристика проектируемого фильтра (правая половина) приведена на рис. 7.6.2 с границей раздела зон между 3 и 4 отсчетами спектра. Расчет оператора фильтра проводим обратным преобразованием Фурье, а по полученным отсчетам оператора вычисляем фактическую частотную характеристику этого оператора с уменьшением шага по частоте в 4-6 раз, что позволяет выявить осцилляции и определить погрешность фильтра (по максимумам осцилляций).

Рис. 7.6.3. Подбор отсчетов переходной зоны НЦФ.

На рис. 7.6.3. показан результат подбора частотных значений характеристики фильтра в районе переходной зоны (2 точки), что позволяет более чем в 30 раз снизить осцилляции частотной характеристики.

Рис. 7.6.4. НЦФ с точкой подбора на границе.

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

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

Метод частотных выборок допускает также рекурсивную реализацию фильтров /43/.

литература

24. Хемминг Р.В. Цифровые фильтры. – М.: Недра, 1987. – 221 с.

43. Айфичер Э., Джервис Б. Цифровая обработка сигналов. Практический подход. / М., "Вильямс", 2004, 992 с.

Главный сайт автора ~ Лекции по ЦОС ~ Практикум

О замеченных опечатках, ошибках и предложениях по дополнению: davpro@yandex.ru.

Copyright ©2008 Davydov А.V.