Сглаживающие фильтры и фильтры аппроксимации /л24/.

Предположим, что требуется осуществить сглаживание (регуляризацию, аппроксимацию) по методу наименьших квадратов (МНК) равномерного по аргументу массива данных.

Фильтры МНК 1-го порядка (МНК-1). Простейший способ аппроксимации по МНК произвольной функции s(t) - с помощью полинома первой степени, т.е. функции вида y(t) = A+Bt (метод скользящих средних). В качестве примера произведем расчет симметричного фильтра на (2N+1) точек с окном от -N до N.

Для определения коэффициентов полинома найдем минимум функции приближения (функцию остаточных ошибок). С учетом дискретности данных по точкам t_n = nDt и принимая Dt = 1 для симметричного НЦФ с нумерацией отсчетов по n от центра окна фильтра (в системе координат фильтра), для функции остаточных ошибок имеем:

s(A,B) = [s_n - (A+B·n)]².

Дифференцируем функцию остаточных ошибок по аргументам 'А, В' и, приравнивая полученные уравнения нулю, формируем 2 нормальных уравнения:

(s_n-(A+B·n)) º s_n - A 1 - B n = 0,

(s_n-(A+B·n))·n º n×s_n - A n - B n² = 0,

С учетом очевидного равенства n = 0, результат решения данных уравнений относительно значений А и В:

А = s_n , B = n×s_n / n².

Подставляем значения коэффициентов в уравнение аппроксимирующего полинома, переходим в систему координат по точкам k массива y(k+t) = A+B·t, где отсчет t производится от точки k массива, против которой находится точка n = 0 фильтра, и получаем в общей форме уравнение фильтра аппроксимации:

y(k+t) = s_k-n + t n×s_k-n / n².

Для сглаживающего НЦФ вычисления производятся непосредственно для точки k в центре окна фильтра (t = 0), при этом:

y_k = s_k-n. (2.1.1)

Импульсная реакция фильтра соответственно определяется (2N+1) значениями коэффициентов b_n = 1/(2N+1). Так, для 5-ти точечного НЦФ:

h(n) = {0.2, 0.2, 0.2, 0.2, 0.2}.

Передаточная функция фильтра в z-области:

H(z) = 0.2(z^-2+z^-1+1+z¹+z²).

Рис. 2.1.1.

Коэффициент усиления дисперсии шумов:

K_q = S_n h²(n) = 1/(2N+1),

т.е. обратно пропорционален ширине окна фильтра. Зависимость значения K_q от ширины окна приведена на рис. 2.1.1.

Частотная характеристика фильтра (передаточная функция фильтра в частотной области) находится преобразованием Фурье импульсной реакции h(n) (фильтр симметричный, начало координат в центре фильтра) или подстановкой z = exp(-jw) в выражение передаточной функции H(z). И в том, и в другом случае получаем:

H(w) = 0.2[exp(2jw)+exp(jw)+1+exp(-jw)+exp(-2jw)]. (2.1.2)

Можно использовать и непосредственно уравнение фильтра, в данном случае уравнение (2.1.1). Подадим на вход фильтра гармонический сигнал вида s_k = exp(jwk). Так как сигнальная функция относится к числу собственных, на выходе фильтра будем иметь сигнал y_k = H(w)exp(jwk). Подставляя выражения входного и выходного сигналов в уравнение (2.1.1), получаем:

H(w) exp(jwk) = 0.2 exp(jw(k-n))= 0.2 exp(jwk) exp(-jwn).

Отсюда, выражение для передаточной функции:

H(w) = 0.2 exp(-jwn) = 0.2[exp(2jw)+exp(jw)+1+exp(-jw)+exp(-2jw)],

что полностью идентично выражению (2.1.2).

Следует запомнить: если оператор фильтра известен, то для получения его частотной характеристики достаточно подставить сигнал exp(jwn) непосредственно в линейное уравнение фильтра. Тем самым выполняются сразу 2 операции: производится z- преобразование h(n) и подставляется z = exp(-jwn), т.е. осуществляется трансформация h(n)→ h(z) → H(w).

Так как импульсная реакция фильтра МНК симметрична (функция h(n) четная), частотное представление передаточной функции должно быть вещественным, в чем нетрудно убедиться, объединив комплексно сопряженные члены выражения (2.1.2):

H(w) = 0.2(1+2 cos w+2 cos 2w).

Альтернативное представление передаточной функции H(w) для фильтра с произвольным количеством коэффициентов 2N+1 нам достаточно хорошо известно, как нормированный фурье-образ прямоугольной функции, каковой по существу и является селектирующее окно фильтра (2.1.1):

H(w) = sin((N+1/2)w)/[(N+1/2)w] = sinc((N+1/2)w). (2.1.3)

Рис. 2.1.2. Сглаживающие фильтры МНК.

Графики передаточных функций (2.1.3) приведены на рисунке 2.1.2. По графикам можно видеть коэффициент передачи сигнала с входа на выход фильтра на любой частоте. Без ослабления (с коэффициентом передачи 1) сглаживающим фильтром пропускается (и должен пропускаться по физическому смыслу сглаживания данных) только сигнал постоянного уровня (нулевой частоты). Этим же определяется и тот фактор (который стоит запомнить), что сумма коэффициентов сглаживающего НЦФ всегда должна быть равна 1 (отсчет ненормированного дискретного фурье-преобразования на частоте w = 0 равен сумме значений входной функции).

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

Рис. 2.1.3. Фильтрация шумов фильтрами МНК 1-го порядка.

Частотное представление передаточных функций позволяет наглядно видеть особенности фильтров и целенаправленно улучшать их характеристики. Так, если в рассмотренном нами фильтре с однородной импульсной реакцией h_n = 1/(2N+1) уменьшить два крайних члена в 2 раза и заново нормировать к сумме S h_n = 1, то частотные характеристики фильтра заметно улучшаются. Для нахождения передаточной функции модифицированного фильтра снимем в выражении (2.1.3) нормировку (умножим на 2N+1), вычтем значение 1/2 крайних членов (exp(-jwN)+exp(jwN))/2 = cos(wN) и заново пронормируем полученное выражение (разделим на 2N). Пример новой передаточной функции при N=3 также приведен на рисунке 2.1.2. Передаточные функции модифицированных таким образом фильтров приводятся к нулю на частоте Найквиста, при этом несколько расширяется полоса пропускания низких частот и уменьшается амплитуда осцилляций в области подавления высоких частот. Если смотреть на сглаживание, как на операцию подавления высокочастотных помех, то модифицированные фильтры без сомнения больше соответствует своему целевому назначению.

Рис. 2.1.4.

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

Допустим, что нужно обеспечить максимальное подавление дисперсии шумов при минимальном искажении верхней граничной частоты сигнала f_в, на которой мощность шумов равна мощности сигнальной гармоники f_в. Значение f_в равно 0.08 частоты Найквиста дискретизации данных, т.е. f_в = 0.04 при Dt=1. Относительные значения мощности (дисперсии) гармоники и шума принимаем равными 1. Спектр модели сигнала + шума в сопоставлении с передаточными функциями фильтров приведен на рис. 2.1.4.

По формуле (2.1.3) вычисляем коэффициенты К_у(f_в) усиления фильтров с N от 0 до 6 на частоте f_в (см. таблицу 2.1.1). При мощности гармоники W_u = 1 амплитудное значение гармоники на входе фильтра равно U = = 1.41. Мощности гармоник на выходе фильтров в зависимости от N:

W_u(N)= 0.5·[U· К_у(f_в)]².

Соответственно, при мощности входного шума W_q=1 мощности шумов на выходе фильтров будут численно равны коэффициентам усиления дисперсии шумов W_q(N) = W_q·K_q(N).

Рис. 2.1.5.

Максимум отношения

К_с/ш(N) = W_q(N)/W_u(N)

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

Рис. 2.1.6.

При К_у(f_в) > 0.5 и W_u(N) = W_q(N) = 1 численные значения величины d²(N) = 1/ К_с/ш(N) в первом приближении могут служить оценкой s²(N) квадрата среднего квадратического отклонения выходных сигналов от "чистой" гармоники f_в, заданной на входе. Свидетельством этому служат последние строки таблицы 2.1.1, где приведены результаты математического моделирования фильтрации по данным условиям на выборке 10000 точек. На рис. 2.1.6 приведены результаты сопоставления расчетных d²(N) и модельных s²(N) значений данных коэффициентов. Эффект фильтрации можно видеть на рис. 2.1.7, где приведен пример сигналов моделирования на ограниченном отрезке данных.

Рис. 2.1.7. Сигналы на входе и выходе фильтра МНК 1-го порядка.

Фильтры МНК 2-го порядка (МНК-2) рассчитываются и анализируются аналогично. Рассмотрим квадратный многочлен вида y(t)=A+B·t+C·t². Для упрощения анализа ограничимся симметричным сглаживающим НЦФ с интервалом дискретизации данных Dt=1.

Минимум суммы квадратов остаточных ошибок:

s(A,B,C) = [s_n-(A+B·n+C·n²)]². (2.1.4)

Система уравнений после дифференцирования выражения (2.1.4) по А, В,С и приравнивания полученных выражений нулю:

A 1 + B n + С n² = s_n.

A n + B n² + С n³ = n·s_n.

A n² + B n³ + С n⁴ = n²·s_n.

При вычислении значения квадратного многочлена только для центральной точки (t=0) необходимости в значениях коэффициентов В и С не имеется. Решая систему уравнений относительно А, получаем:

n⁴ s_n - n² n²s_n

A = -----------------------­ (2.1.5)

1 n⁴ - [ n²]²

При развертывании выражения (2.1.5) для 5-ти точечного НЦФ:

y_o = (17 s_n - 5 n²s_n) /35 = (-3·s_-2+12·s_-1+17·s_o+12·s₁-3·s₂) /35. (2.1.6)

Импульсная реакция: h_n = {(-3, 12, 17, 12, -3)/35}.

Передаточная функция фильтра:

H(z)= (-3z^-2+12z^-1+17+12z¹-3z²)/35. (2.1.7)

Рис. 2.1.8. Сглаживающие фильтры МНК.

Аналогичным образом выражение (2.1.5) позволяет получить импульсную реакцию для 7, 9, 11 и т.д. точек фильтра:

³h_n = {(-2,3,6,7,6,3,-2)/21}.

⁴h_n = {(-21,14,39,54,59,54,39,14,-21)/231}.

⁵h_n={(-36,9,44,69,84,89,84,69,44,9,-21)/459}.

Подставляя значение z = exp(-jw) в (2.1.7) или непосредственно в (2.1.6) сигнал s_n = exp(jwn) и объединяя комплексно сопряженные члены, получаем частотную характеристику 5-ти точечного сглаживающего фильтра МНК второго порядка:

H(w) = (17+24 cos(w)-6 cos(2w))/35.

Вывод формул передаточных функций для 7, 9, 11-ти точечных фильтров МНК предлагается для самостоятельной работы.

Рис. 2.1.9.

Вид частотных характеристик фильтров при N=3 и N=5 приводится на рис. 2.1.8. При сравнении характеристик с характеристиками фильтров МНК-1 можно видеть, что повышение степени полинома расширяет низкочастотную полосу пропускания фильтра и увеличивает крутизну ее среза. За счет расширения полосы пропускания главного частотного диапазона при тех же значениях N коэффициенты усиления дисперсии шумов фильтров МНК-2 выше, чем фильтров 1-го порядка, что можно видеть на рис. 2.1.9.

Рис. 2.1.10.

Методика выбора окна фильтра под частотные характеристики входных сигналов не отличается от фильтров МНК 1-го порядка. На рис. 2.1.10 приведены значения d²(N) и s²(N) фильтров МНК-2 в сопоставлении со значениями фильтров МНК-1 для частоты f_в = 0.08 Гц при Dt=1. Из сопоставления видно, что по своим характеристикам подавления шумов фильтры МНК-2 примерно соответствуют фильтрам МНК-1 при в 2 раза большей ширине окна. Об этом же свидетельствует и пример моделирования фильтрации, приведенный на рис. 2.1.11.

Рис. 2.1.11.

Модификация фильтров. Фильтры МНК второго порядка (равно как и другие фильтры подобного назначения) также можно модифицировать по условию H(w) → 0 при w → p. Один из простейших методов модификации заключается в следующем. В выражение передаточной функции (со всеми коэффициентами фильтра, вида (2.1.7)) подставляем z = exp(-jw), заменяем значения концевых коэффициентов фильтра на параметры, принимаем w = p, и, приравняв полученное выражение нулю, находим новые значения концевых коэффициентов, после чего сумму всех коэффициентов нормируем к 1 при w = 0.