Нерекурсивные частотные фильтры

Тема 4: нерекурсивные частотные фильтры.

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

 Цицерон. О высшем благе и высшем зле.

 (Римский сенатор и философ)

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

Евгений Буцко. Назидания идиотам.

 (Ленинградский инженер и геофизик)

Содержание: 4.1. Общие сведения. Типы фильтров. Методика расчетов. 4.2. Идеальные частотные фильтры. Импульсная реакция фильтров. 4.3. Конечные приближения идеальных фильтров. Применение весовых функций. Весовая функция Кайзера. 4.4. Дифференцирующие цифровые фильтры.   4.4. Гладкие частотные фильтры.  Литература.

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

Рекомендуемые материалы

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

y_k  =h_n  s_k-n .

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

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

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

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

                                   h_в(0) = 1- h_н(0),      h_в(n) = -h_н(n)   при n0.                       (4.1.1)

Применяется также способ получения фильтров высоких частот из низкочастотных фильтров путем реверса частоты в передаточной функции низкочастотного фильтра, т.е. заменой переменной w на переменную w' = p-w (при Dt = 1). Для симметричных фильтров, содержащих в передаточной функции только косинусные члены аргумента w, в результате такой операции будем иметь: 

cos n(p-w) = cos np cos nw = (-1)ⁿ cos nw.

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

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

v_k = h_н(n) _* s(k-n),      y_k = h_в(n) _* v_k = h_н(n) _* h_в(n) _* s(k-n).

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

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

Методика расчетов НЦФ в самом общем виде включает:

1. Задание идеальной частотной передаточной функции фильтра.

2. Расчет функции отклика идеального фильтра (обратное преобразование Фурье передаточной функции фильтра).

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

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

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

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

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

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

h(nDt) = (1/2p)H(w) exp(jwnDt) dw.

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

                            h(nDt)=  (1/p)H(w) cos(npw/w_N) dw,  n = 0,1,2,...          (4.2.1)

Для идеального полосового фильтра H(w)=1 в полосе частот от w_н до w_в и интеграл (4.2.1) вычисляется в этих пределах. Идеальные фильтры низких и высоких частот можно считать частными случаями идеальных полосовых фильтров с полосой пропускания от 0 до w_в для низкочастотного и от w_н до w_N для высокочастотного фильтра.

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

  Пример 1.  Dt = 0.1 сек.  f_N = 1/2Dt = 5 Гц.  w_N = p/Dt = 10 p.

  Пример 2.  Dx = 10 метров.  f_N = 0.05 м^-1.  w_N= 0.1 p.

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

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

                                 h(n) = (А/p) [w_в sinc(nw_в) - w_н sinc(nw_н)],                          (4.2.2)

h_o = (w_в - w_н)/p,    h(n) = (sin nw_в - sin nw_н)/(np).

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

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

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

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

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

                                                           H(w) =h(n) cos nw.                                       (4.3.1)

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

                                                           H'(w) =h(n) cos nw,

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

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

                                                                       N = 2p/D_p.                                                    (4.3.2)

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

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

                                                        h_n = h(n)·p_n.                                                (4.3.3)

Подстановкой коэффициентов (4.3.3) в (4.3.1) рекомендуется произвести построение полученной передаточной характеристики фильтра и непосредственно по ней оценить пригодность фильтра для поставленных задач. Это наглядно видно на рис. 4.3.1, где для нашего примера была применена весовая функция Гаусса. Передаточная функция H_p(w) имеет практически такую же крутизну, как и функция H'(w) при N=100 и практически плоскую вершину в интервале спектра сигнала. Качество работы фильтра для сигнала, приведенного на рис. 4.2.1, можно видеть на рис. 4.3.2.

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

p(n) = .

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

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

D = N·D_p/p.

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

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

                                               = 0.9222                     при А<21.

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

                                               = 0                              при А<21.

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

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

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

N = pD/D_p.

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

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

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

  1. А= -20 log d.  А= 34.                                                               2. N= p (A-7.95)/(14.36 D_p).  N= 18.

  3. b= 0.5842(A-21)^0.4 +0.07886(A-21).  b= 2.62.   4. h_о= (w_в-w_н)/p.  h_о= 0.3

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

  6. p_n= J_o{b} / J_o{b}.  p_n = 1.00, 0.997, 0.9882, .......

  7. Oператор фильтра: h_n= h(n)p_n,  n=0,1,2,...,N.  h_-n=h_n.  h_n= 0.3000, 0.04508, -0.2420, ........

  8. Проверка по формуле: H(w) =h_n cos nw,   0 £ w £ p.

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

                (!!!КР7 - Разработка программы расчета сглаживающих НЦФ).

                (!!!КР8 - Разработка программы расчета полосовых НЦФ).

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

Передаточная функция. Из выражения для производной

d(exp(jwt))/dt = jw exp(jwt)

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

h_n [exp(jwn)-exp(-jwn)] = 2j h_n sin nw,

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

H(w) = 2j(h₁ sin w + h₂ sin 2w + ... + h_N sin Nw),

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

y_n = h_n(s_k+n - s_k-n).

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

H_н(w) = w,  w £ w_в,     H_н(w) = 0,  w_в< w £ w_N.

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

                            h(n) = (1/p)H_н(w) sin(npw/w_N) dw,  n = 0,1,2,...                (4.4.1)

Принимая, как обычно, w_N = p (Dt = 1)  и решая (4.4.1) при H_н(w) = w, получаем:

                                    h_n = (1/p)[sin(nw_в)/n² - w_в cos(nw_в)/n],                          (4.4.2)

h_о = 0,   h_-n = -h_n.  

Проверка:                  H_н(w) =h_n sin nw = 2h_n sin nw .                 (4.4.3)

На рис. 4.4.1 приведен пример расчета коэффициентов дифференцирующего фильтра на интервал {0-0.5}p при Dt=1 (w_в = p/2). Операторы дифференцирующих фильтров, как правило, затухают очень медленно и, соответственно, достаточно точная реализация функции (4.4.3) весьма затруднительна.

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

Для оценки возможных погрешностей дифференцирования усеченными операторами произведем расчет фильтра при w_в = p/2. По формулам (4.4.2) определяем: 

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

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

у =h_n  s_o-n = 2 n h_n,

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

Такое существенное расхождение с действительным значением производной объясняется тем, что при w=0 тангенс угла наклона реальных передаточных функций фильтра, как это видно на рисунке 4.4.2, весьма существенно отличается от тангенса угла наклона аппроксимируемой функции H(w)= w. На рис. 4.4.3 приведены частотные графики относительной погрешности дифференцирования s = H_н'(w)/H_н(w) с вычислением значений на нулевой частоте по пределам функций при N ® ¥. На рис. 4.4.4 приведен пример операции дифференцирования s_*h гармоники s с частотой w_o оператором с N=10 в сопоставлении с точным дифференцированием ds/dk.

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

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

h_n = (w_нcos nw_н-w_вcos nw_в)/(np) + (sin nw_в-sin nw_н)/(n2p).

(!!!КР9- Разработка программы расчета НЦФ дифференцирования)               

(!!!КР10- Оценка возможностей усечения операторов НЦФ, умноженных на весовые функции).

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

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

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

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

                                                      H(w) = h_о+ 2h_n cos nw.                                    (4.5.1)

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

                                                     H(w) =g_n (cos w)ⁿ =g_n xⁿ,                           (4.5.2)

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

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

     g(x)= (1+x)^z (1-x)^r,                                        (4.5.3)

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

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

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

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

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

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

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

  Произвести расчет ФНЧ с гладкой частотной характериcтикой с перегибом характеристики в точке p/3.

  За исходную функцию принять функцию (4.5.3).

В лекции "Понятие политического режима в политологии" также много полезной информации.

  1. x= cos(p/3)= 0.5= (z-r)/(z+r). Принято: z=3, r=1. Исходный многочлен: g(x) = (1-x)(1+x)³ = 1+2x-2x³-x⁴.

  2. h(x)=  g(x)dx = C+x+x²-0.5 x⁴-0.2 x⁵.      При х = -1,  h(-1) = 0,  откуда С=0.3.  При х=1,  h(1)=1.6.

  Отсюда:  H(x)= (3+10x+10x²-5x⁴-2x⁵)/16.      g_n = {3/16, 10/16, 10/16, 0, -5/16, -2/16}.

  3. Применяя рекурсивное преобразование, получаем коэффициенты ФНЧ: h_n= {(98, 70, 20, -5, -5, -1)/256}.

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

(!!!КР11- Разработка программы расчета гладких полосовых фильтров)