CC BY
44
РЕАЛИЗАЦИЯ РЕКУРСИВНОГО ЦИФРОВОГО ФИЛЬТРА НА ОСНОВЕ ШТРАФНОГО Р-СПЛАЙНА
Е.А. Кочегурова1, Д. У1 1Национальный исследовательский Томский политехнический университет, Томск, Россия
Аннотация
В работе показана возможность проектирования рекурсивного цифрового фильтра с использованием штрафного Р-сплайна. Аналитически получены и исследованы частотные и временные характеристики сплайн-фильтра для данных, поступающих в режиме реального времени. Исследовано влияние параметров штрафного Р-сплайна на показатели эффективности интерпретации входной измерительной информации. Закономерности, полученные в процессе частотного анализа сплайн-фильтра, подтверждены примером восстановления до-плеровской функции.
Ключевые слова: Р-сплайн, сглаживающий сплайн, рекурсивный цифровой фильтр, аппаратная функция, системная функция.
Цитирование: Кочегурова, Е.А. Реализация рекурсивного цифрового фильтра на основе штрафного Р-сплайна / Е.А. Кочегурова, Д. У // Компьютерная оптика. - 2018. - Т. 42, № 6. -С. 1083-1092. - БСЯ: 10.18287/2412-6179-2018-42-6-1083-1092.
Введение
Эффективность систем обработки информации во многом определяется качеством входной информации, поступающей от реальных объектов наблюдения и управления. Поэтому в составе современных вычислительных устройств нередко присутствуют цифровые фильтры (ЦФ), реализованные аппаратно или программно. При аппаратной реализации системная (аппаратная) функция измерительного преобразователя корректируется динамическим звеном с заданной передаточной функцией [1]. Алгоритмическое или программное сглаживание основано на использовании математических методов цифровой фильтрации. Основными методами цифровой фильтрации для данной задачи являются частотная селекция сигналов и оптимальная (адаптивная) фильтрация [2, 3].
Требования к характеристикам современных ЦФ весьма высоки и определяются селективными свойствами фильтров [4 - 7]. Так, во временной области эти требования определяются точностью и сохранением формы полезного сигнала. В частотной - требованиями к полосе пропускания сигнала и подавлению помехи. Стремление к сокращению переходной зоны между этими полосами приводит к увеличению окна и, соответственно, порядка фильтра [8].
Альтернативное решение основано на применении рекурсивных фильтров, имеющих бесконечную импульсную характеристику (БИХ-фильтр) при невысоком порядке ЦФ.
Классические методы создания БИХ-фильтров базируются на их аналоговом прототипе с последующей аппроксимацией требуемых свойств фильтра на основе физически реализуемой дискретной системы [2, 3, 9]. При таком подходе требования к ЦФ полностью определены требованиями к АЧХ (реже к ФЧХ) фильтра. Это снижает эффективность фильтрации с позиции последующей интерпретации измерительной информации. Подход к проектированию ЦФ, позволяющий объединить разнообразные требования к качеству фильтрации, возможен на основе численных методов математического программирования.
В задачах интерпретации данных целесообразно использовать модельные функции на основе простого математического базиса. Именно этим объясняется популярность и эффективность сплайн-функций в прикладных задачах интерпретации (прогноз, фильтрация, дифференцирование, выделение тренда и др.) [10 - 12].
Для создания рекурсивного сплайн-фильтра необходима сплайн-функция, допускающая реализацию в реальном времени, т.е. с рекуррентной формой вычисления. Базисные Р-сплайны [13 - 17] и штрафные Р-сплайны [18 - 25] имеют высокое быстродействие только при оптимальном выборе узлов. Выбор узлов является одним из способов адаптации сглаживающего сплайна [18, 19, 24 - 26]. Адаптация значительно повышает эффективность сплайнов.
Другое направление адаптации связано с оценкой штрафного параметра на основе регрессионных и вероятностных подходов [18, 19, 21, 27 - 29].
Предлагаемая здесь модификация Р-сплайна реализует вариационный подход для объединенных в группы данных [30 - 32]. Создание группы для входных данных автоматически решает проблему выбора узлов. А компактная форма оценки параметров рекуррентной сплайн-функции позволяет реализовать сплайн в реальном масштабе времени (РМВ). Кроме того, анализ частотных и временных характеристик сплайна с позиций ЦФ позволяет более наглядно оценить прикладные свойства сплайна для задач интерпретации данных.
Штрафной P-сплайн для группируемых данных
Получение классического штрафного сплайна S(t) по отсчетам y(ti), i = 1, n основано на оптимизации специального вида функционала, определенного на всем интервале наблюдения [a, b] (B.W. Silverman, 1994)
J (S) =X-j [S'(t)] 2dt + J [S (ti) - y(ti)]2
(1)
где первое слагаемое вместе со сглаживающим параметром определяет штраф кривизны; второе слагае-
мое реализует обычный метод наименьших квадратов.
Для реализации режима реального времени функционал (1) модифицирован [23] и определен для каждого 1-го звена сплайна на основании входных данных, собранных в группы по к измерений между отсчетами /0, 4 с шагом дискретизации Д/
з (£) = (1 -Р)(кд )21 [ я" (/)]Ч/-
+Р£[ я $) - у $ )]2.
1 = 0
(2)
В (2) введен сглаживающий сомножитель р, нормированный по аналогии с [13], но только для каждого звена сплайна. Это сужает диапазон выбора ре[0, 1] и придает сглаживающему параметру физический смысл.
Для реализации штрафного сплайна для 1-го звена в реальном времени
Я,- (т) = а'0 + а{ -х + а2 -х2 + а33 -х3, - q <х< к - q, (3)
нами были получены [32] коэффициенты кубического сплайна в рекуррентном виде. При этом q - номер отсчета Ц внутри '-го звена (] = 0,к; q = 0,к -1), в котором сопряжены непрерывные производные сплайна Я(кЧ^;!) + = Я(к44)-, к = 0,1 для дефекта 2 и к = 0 для дефекта 1. Разрывные коэффициенты а'2, а'3 найдены из условия минимизации функционала (2):
дЗЯ = 0дЗЯ = 0
да'2 да\
Соотношение моментов вычисления (х) и сопряжения (д) '-го звена сплайна порождает несколько вычислительных схем рекуррентного сплайна, подробно исследованных в [33].
Описание штрафного сплайна в виде разностного уравнения
Для исследования рекуррентного сплайна методами теории цифровой фильтрации требуется представить его в виде разностного уравнения. Входной информацией такого фильтра является последовательность измерений '-го звена {у(/})}, 1 = 0,к (решетчатая функция). Выходная информация - значение сплайна Я(/х), т = q +1, q = 0, к-1. Я(/х), q = 1, т = 1, к . Обозначив Я, = Я (/Х), получим я-т= Я (/Х-1) = а0, Я- т= Я '(/Х-1) = а/.
и
Суммы измерений ^ у(/'к), входящие в функцио-
к =0
нал (2), обозначены для рекуррентного случая как
'+к
^ук-х . Решетчатые функции {Я,}, {у,}, ' = 1, 2, 3...
к ='
имеют одинаковый интервал дискретизации, совпадающий с интервалом дискретизации входных данных Д/.
Для введенных обозначений перепишем уравнение сплайна (3) как разностное уравнение для решетчатой функции {Я'}, ' = 1, 2, 3.
У0Я1 - У1 (х)Я'-х - у2 (х)Я/-х =
''+к г+к (4)
где
= У3 (х)Х ук-х (к - х)2 +у 4 (х)Х ук-х (к - х)3,
к =' к =''
у 0 = ВС - А2;
У1 (х) = У0 + рх2 (АН3 -СН2) + рх3 (АН2 - ВН3); У2(х) = У0х + рх2 (АН4 - СН3) + рх3 (АН3 - ВН4); У 3(х) = р( х2С -х3 А); У 4 (х) = р(х3 В -х2 А). А = 6(1 -р)к4 +рН5; В = 4(1 -р)к3 +рН 4; С = 12(1 -р)к5 +рН6;
Н„ = ±к".
Разностное уравнение (4) соответствует уравнению с переменными коэффициентами у 1 (х), ] = 0,4. Порядок уравнения соответствует значению х, хе[1, к], где х - момент вычисления значения сплайна. Однако если момент вычисления сплайна х не изменяется на разных звеньях, то коэффициенты
У1, ] = 0,4, будут постоянными относительно звена
сплайна и параметр х опускается. В (4) и далее q = 0, т.е. сопряжение звеньев сплайна осуществляется в начале звена.
Последовательно раскладывая в выражении (4)
Я/-х,Я-х и далее Я/-2х,Я/-2х,...,Я0,Я0 (что соответствует разложению непрерывных коэффициентов
-1, а0-2,а;-2,..., а00, а10), получаем разност-
сплайна а
ное уравнение в виде
Я=1
01 (р, к, q)£ ук-х (к -х)2 +
к=' 1+к
+02 (р, к, q)£ ук-х (к -х)3
+ С.
(5)
Вид разностного уравнения (5) указывает на линейность уравнения относительно измерений у(4). Причем коэффициенты 01 (р,к, q), 02 (р,к, q) полностью определены настроечными параметрами сплайна р, к, q, а константа С зависит только от начальных условий сплайна.
Для подтверждения линейности сплайн-фильтра проведены эксперименты по оценке свойств линейных преобразований: сохранение частотного состава входного гармонического сигнала и сохранение закона распределения плотности вероятности входного случайного сигнала [2, 3]. Было установлено, что в спектре выходного сигнала после сплайн-
к=0
I =0
преобразования не появляются дополнительные значимые гармоники. Их амплитуда составляет (0,00001 - 3) % от амплитуд существующих гармоник во входном сигнале. Также установлено, что после сплайн-преобразования случайного сигнала типа белый шум в ограниченной полосе частот с нормальным законом распределения выходной сигнал также подчинен нормальному закону распределения. Согласованность эмпирического закона распределения с теоретическим, рассчитанная по %2- критерию Пирсона, составляет 96 % (при объеме выборки 1000 значений).
Линейность сплайн-фильтра делает возможным использовать для его исследования аппарат линейных динамических систем.
Исследование частотных и временных
характеристик рекурсивного сплайн-фильтра
Традиционно проектирование программного ЦФ [5, 8 - 12] включает задание требуемой передаточной функции фильтра в спектральной области, ее аппроксимацию с заданной точностью, переход в 7-область и, наконец, получение импульсной весовой функции на основании обратного 7-преобразования. И таким образом, основным математическим инструментом цифровой фильтрации является уравнение свертки аппаратной функции фильтра и входного сигнала [34].
При анализе сплайн-фильтра уравнением преобразования во временной области является уравнение сплайна вида (3) или разностное уравнение (4). Особый интерес для понимания сути и наглядности сплайн-фильтра имеет его описание в частотной области [16, 33, 35 - 38]. Кроме этого, появляется альтернативное описание сплайн-фильтра во временной и частотной области.
Цифровые фильтры по аналогии с линейными динамическими системами могут быть представлены следующими основными характеристиками [4 - 8]:
- системная функция (частотная передаточная или комплексный коэффициент передачи для непрерывных систем);
- аппаратная функция (импульсная весовая для непрерывных систем).
Для решения разностного уравнения (4) использовано дискретное преобразование Лапласа, в данном случае в форме 7-преобразования. Получив изображения левой и правой части уравнения (4) S(z) = Si и
Y(z) = уЩ ), можно найти передаточную функцию W(z) = S(z) / Y(z) сплайн-фильтра. Для практического использования целесообразно перейти от системной функции фильтра к комплексному частотному коэффициенту передачи Н(-) = W(ejaЛt), полагая z = ejши М = 1, используя свойства сдвинутых решетчатых функций h
£е(1-т)- • {(1 - т)2 у3 + (к - т)3 у4)
H (-) =-----. (6)
У0 -ye
-у 2юе
Числитель частотной характеристики определен нерекурсивной частью сплайн-фильтра, знаменатель - рекурсивной.
Периодичность функции e)aAt приводит к периодичности частотной характеристики Н(ю) [2, 3, 8]. Поэтому численное исследование амплитудных и фазовых частотных характеристик сплайн-фильтра проведено в области главных частот [-л / At, л / At], где частотные характеристики дискретного и непрерывного фильтра совпадают, а также исключается эффект алиасинга.
Представим частотную характеристику Н(ю) в виде
H(ю) = А(ю) • e-1ф(ю),
где А(ю) = |Г(ю)| - амплитудно-частотная характеристика (АЧХ), ф(ю) = arg [W(ra)] - фазо-частотная характеристика (ФЧХ) фильтра.
Основные параметры сплайна значительно изменяют фильтрующие свойства сплайна. Общую картину этих свойств представляют АЧХ и ФЧХ сплайн-фильтра, приведенные на рис. 1. При исследовании параметры сплайна варьировались в диапазоне их определения: сглаживающий множитель ре[0, 1], момент вычисления те[1, h], длина звена he[3, 20] (h = 3 соответствует использованию 4 измерений для построения полинома третьего порядка, а при h > 20 сплайн сложно использовать в системах реального времени). Вид всех АЧХ на рис. 1а соответствует низкочастотному фильтру с разной степенью влияния параметров. Если рассматривать сплайн-фильтр как динамическое звено второго порядка, то его коэффициент передачи
(1/ю0) J|H (ю)|2 d ю
растет с увеличением параметров т и р. Однако резонансные пики АЧХ возникают при небольших значениях р или И и при приближении т^И.
в) (о г) со
Рис. 1. АЧХ и ФЧХ рекуррентного фильтра (P-сплайн): а) h = 10, т = 1, ре[0, 1]; б) h = 10, р = 1, he[3, 20]; в) h = 10, р = 0,5, те [1, h] г) h = 10, р = 0,5, те[1, h]
На характер ФЧХ (рис. 1г) наибольшее влияние оказывает момент вычисления сплайна х, соответствует ФЧХ системы с запаздыванием. В данной реализации сплайн-фильтра имеем дело с обработкой информации группами [32]. В силу этого ФЧХ имеет опережающий характер, и лишь при х^-к появляется отставание фазовой характеристики.
На полосу пропускания сигнала не столь существенно влияют параметры сплайна х и р, как длина звена сплайна к. И это очевидно, так как при увеличении к сплайн становится более гладким и полоса пропускания уменьшается. Это свойство сплайн-фильтра подтверждено ниже на рис. 2а при анализе ширины аппаратной функции. При малых значениях к инерционность сплайна снижается и в большей степени соответствует АЧХ интерполяционного сплайна. Для идеального фильтра А(ю) = 1, ф(ю) = 0 в заданной полосе частот.
Аналогом импульсной весовой функции линейных динамических систем для дискретных фильтров является аппаратная функция (АФ) g(t)
g(/) = — [ Н(ю) 1аю .
2л
Аппаратная функция часто применяется в уравнении дискретной свертки для аналитической связи сигналов на входе и выходе дискретного фильтра.
Для предлагаемого сплайн-фильтра АФ получена аналитически на основе g(t) и системной функции (6). На рис. 2 продемонстрированы семейства АФ при изменении параметров р, к или х. Во всем диапазоне изменения параметров сплайна АФ имеет затухающий характер. Хорошо прослеживается связь между шириной АФ (Д(р, к)) и частотными характеристиками сплайн-фильтра: полоса пропускания сигнала обратно пропорциональна ширине АФ. Как и при анализе частотных характеристик, ширина АФ в большей степени зависит от длины звена сплайна к и мало зависит от сглаживающего параметра р (рис. 2а). Более гладкий сигнал на выходе фильтра соответствует использованию большого значения к, т.е. узкой полосе пропускания рис. 16 и большой ширине АФ рис. 2а, б.
Момент вычисления х определяет свойство причинности сплайн-фильтра рис. 2в. В данном случае к(/) = 0, / < 0 только для значения х = к, т.е. когда нет запаздывания в вычислениях. В остальных случаях к(/) = 0 при х < к - х, что соответствует обработке данных группами или в системах с накопленными данными.
Кроме требований к установившимся процессам, эффективность и работоспособность ЦФ характеризуется показателями в переходных процессах и условиями устойчивости.
Устойчивость рекурсивного сплайн-фильтра как линейной дискретной системы была оценена аналитически на основе характеристического уравнения [2,
8], соответствующего знаменателю уравнения (6) для х = 1
У0 -У12-1 - У2 1п( 2 ) 2-1 = 0. А(А,р)
5(0
0,3 1 ^^ 1=к
1=6
/ -с=3
2 0
Рис. 2. Аппаратная функция рекуррентного фильтра (Р-сплайн): а) к е[3, 20], ре[0,1]; б) х= 1, р = 0,6; в) к = 10, р = 0,6
Преобразуем данное уравнение, умножая его на 2 и используя разложение 1п(2) в степенной ряд, ограничиваясь одним членом разложения:
2 2-1 0
У02-У1 -У2 2-2-- = 0.
2 +1
Или
С0 22 + С12 + С2 = 0 .
Далее, используя билинейное w-преобразование (Мизеса) 2 = (1+ ю) / (1 - ю), получим
и0+ и^ + и2 = 0 ,
где
С0 = У 0 ; и0 = С0 - С1 + С2 ; С = У0-У1 -2У2; и = 2С0 -2С2; С2 = 2У 2 -У 0; и2 = С0 + С1 + С2.
Замена переменной 2 позволяет свести условия устойчивости 2| < 1 к критерию Гурвица и для уравнения второго порядка имеют вид
с0 - с1 + с2 > 0; < 2с0 - 2с2 > 0; с0 + с1 + с2 > 0.
А с учетом введенных в (4) обозначений для / = 1 условия устойчивости сплайн-фильтра
У0 +р-Р2 > 0; -р-(Р1 -2Р2) > 0; (7)
-р-Р1 > 0,
где
Гр = А(Н2 + Нз) -СН2 -ВНз; [р = А(Нз + Н4) - СНз - ВН4.
Система неравенств (7) определяет условия устойчивости фильтра. Разрешая (7) относительно сглаживающего параметра р, были найдены области устойчивости сплайн-фильтра (табл. 1) при изменении параметров q и к.
Табл. 1. Области устойчивости сплайн-фильтра
h/q 0 1 2 3 4 5 6 7
3 0 - 1 0 - 0,99 0,99- 1
6 0 - 1 0 - 1 0 - 0,98 0,90 - 1 0,99 -1
8 0 - 1 0 - 1 0 - 1 0 - 0,97 0,87 -1 0,97 -1 н н
10 0 - 1 0 - 1 0 - 1 0 - 0,99 0 -0,95 0,84 -1 0,95 -1 н
н - неустойчив
Основная особенность состоит в сужении области устойчивости от момента сопряжения звеньев сплайна q. И при стремлении q к правому концу звена сплайна наблюдается неустойчивость при любом значении параметра р. При сопряжении в начале звена сплайн всегда устойчив и при т < 1. В практических приложениях чаще всего и используется сопряжение в начале звена.
Аналитически оценить устойчивость сплайн-фильтра при т > 1 затруднительно. Однако косвенно устойчивость преобразования можно оценить на основе зависимости погрешностей ст восстановления некоторого модельного сигнала при значительном изменении его длины. Так, на рис. 3а, б представлена зависимость MSPE (Mean Squared Percentage Error) [33] для функции Доплера (8) при изменении числа отсчетов сигнала ne[102, 105]. Аддитивный шум №;(0, 0, 22), т.е. 20 % от максимального значения полезного сигнала.
а)
20 16 12 8 4 0
\ Ч =0,1
I
Г"- рч
...... и," х'лИ
2
MSPE, %
10
L \
\ \ Р=0,2 р=0,8
1
б) 3,0 3,5 4,0 4,5 Рис. 3. Зависимость погрешности восстановления функции от числа отсчетов сигнала, к = 10, р= 0, т = к
Изменение погрешности наиболее различимо при изменении сглаживающего множителя р. И на рис. 3 приведено для максимального значения т = к.
Качество сплайн-преобразования как программного ЦФ невысокого порядка было оценено прямыми методами по реакции на единичное ступенчатое возмущение. Проведенные исследования показали, что графики переходного процесса (1111) Н(пТ) имеют колебательный характер. Все показатели ПП: длительность, перерегулирование, колебательность - ухудшаются при т ^ к, т.е. при уменьшении запаздывания вычисления сплайна к - т.
Для рекуррентного ЦФ важнейшей характеристикой является длительность ПП, которая увеличивается с ростом длины звена сплайна к. И в области изменения параметров сплайна длительность ПП составляет (15 - 20)Д/, что вполне допустимо в системах реального времени.
Результаты вычислительного эксперимента
Для оценки исследованного Р-сплайна реального времени был выбран часто используемый в сравнительном анализе пример доплеровской функции [18, 19]
' (t)=5^. s,„ (« )
(8)
заданной на интервале [0, 1] п дискретными отсчетами с аддитивной помехой №;(0, 0, 22), т.е.
№ = /(/,)+ад.
Данная функция довольно сложна для математической обработки, имеет значительное изменение частоты, амплитуды и является моделью многих процессов в медицине, оптике и других приложениях, основанных на эффекте Доплера (изменение частоты волн в зависимости от движения источника этих волн относительно наблюдателя).
Гистограммы остатков текущего восстановления доплеровской функции в условиях помех №;(0, 0, 22) приведены на рис. 4 при изменении параметров сплайна р, к.
0,Л т
0
-0,1
0,15-
-0,15
а)
0,2 0,6 0,8 б) 3 6 10 15 20 Рис. 4. Гистограммы остатков текущего восстановления функции (8): а) к = 10, б) р = 0,5
Характер и величины остатков подтверждают влияние параметров сплайна на результирующую погрешность, выявленное при частотном анализе сплайн-фильтра. Так, систематическая погрешность, обусловленная шириной аппаратной функции (рис. 2а), растет с увеличением длины звена к или с повышением гладкости сплайна, т.е. с уменьшением р.
Ниже на рис. 5 приведены примеры восстановления доплеровской функции предлагаемым Р-сплайном в режиме реального времени при изменении параметров сплайна р, к. Исходная функция, так же, как на рис. 4 / [19], задана на интервале [0, 1] в 400 дискретных отсчетах.
На рис. 56, г, е представлены зуммированные части функции, соответствующие высокочастотным фрагментам входного сигнала в начале интервала обработки.
/
ч _ JLJ—Ш_ h=io - 1 И I р =0,8
0 0,2 0,4 0,6 0,8 1,0
О 0,2 0,4 0,6 0,8 1,0
s A л
щ h=3 р=0,8 1=1 t
\h=10; p=0,8; i=1\tA— /К
1—IHH—r t
0,1 0,2 0,3
0 0,2 0,4 0,6 0,8 1,0 0,1 0,2 0,3
Рис. 5. Восстановление P-сплайном в режиме реального времени (серая пунктирная линия).
Доплеровская функция (чёрная сплошная линия)
Сопоставляя результаты исследования предлагаемого P-сплайна и приведенные в [18, 19], можно увидеть аналогичное качество восстановления с лучшим из приведенных на рис. 4/ в [19] adaptive P-spline. Следует отметить, что adaptive P-spline в [19] использован с учетом адаптации сглаживающего параметра. Визуально качество восстановления доплеровской функции на рис. 4/ в [19] соответствует варианту, представленному на рис. 5д, е. Гораздо лучшее восстановление высокочастотного фрагмента входного сигнала продемонстрировано на рис. 4г. В этом и состоит преимущество предлагаемого P-сплайна: доступность изменения параметров сплайна на разных звеньях. Кроме того, P-сплайн допускает реализацию в реальном времени, что является явным достоинством в прикладных задачах интерпретации данных.
В качестве другого модельного примера, графически и численно иллюстрирующего свойства рекуррентного сплайн-фильтра, была выбрана функция
) = (1 -12) • exp 1 j. (9)
В вейвлет-анализе данная функция (рис. 6, красная линия) носит название «мексиканская шляпа» и используется в качестве материнского вейвлета в задачах, требующих хорошего пространственного раз-
решения и не требовательных к спектральному разрешению.
f(f)
1 i
- - 7олезныь
сигнал ............n=100
----и =200 =400 Nt
i-wi ¿V» ST. /л s?
..ТЧз^ /в 0 i\ *♦»
A\ d /
4
Рис. 6. Восстановление Р-сплайном в режиме реального времени функции «мексиканская шляпа» при изменении числа отсчетов функции п, к = 10, р = 0,5
На интервале [-10, 10] функция (9) задана п отсчетами. Аддитивная смесь сигнала и помехи §(/) показаны на рис. 6 точками для п = 500.
Следует отметить, что число отсчетов п является самым влиятельным параметром восстановления данной функции (9). И, как видно из рис. 6, удовлетворительное качество фильтрации в РМВ возможно для п = 400 (серая пунктирная линия, графически совпадающая с полезным сигналом). Погрешность МЕРЕ (%) для этого случая составляет 1,29 % при ст^ = 0,2 (20 % от максимального значения полезного сигнала). Для других значений параметров погрешности приведены в табл. 2 для к = 10 и р = 0,5.
Табл. 2. Погрешности MSPE
n | 0 5 10 15 20
100 11,9 11,97 12,0 12,12 12,6
200 4,06 4,07 4,11 4,29 4,31
400 0,85 0,87 1,01 1,19 1,29
500 0,5 0,57 0,78 0,89 1,05
1000 0,1 0,28 0,5 0,86 1,1
При нулевой погрешности удовлетворительные данные достигнуты при п > 400, и в диапазоне п е 4 • 103 -105 погрешность составляет 10-2 - 10-5.
Анализ доплеровской функции позволил графически иллюстрировать исследованные частотные свойства Р-сплайна реального времени. При уменьшении длины сплайна к увеличивается его пропускная способность (рис. 16), что особенно хорошо проявляется при восстановлении высокочастотных сигналов (рис. 4г, е). Однако при увеличении момента вычисления х (рис. 4в) появляются нежелательные осцилляции (всплески) восстановленной функции, обусловленные появлением пиков АЧХ, рис. 1в при х^к. Именно в возможности варьировать параметры сплайна на разных звеньях в зависимости от изменения свойств поступающей информации состоит потенциал Р-сплайна реального времени.
Заключение
В работе показана возможность представления рекуррентного Р-сплайна в виде цифрового фильтра с бесконечной импульсной характеристикой и переменными параметрами. Аналитически найденная системная функция фильтра позволила исследовать частотные характеристики сплайн-фильтра и установить, что во всем диапазоне параметров сплайна фильтр остается низкочастотным. Селективные свойства сплайн-фильтра наиболее ясно иллюстрируются шириной аппаратной функции, которая значительно зависит от параметров сплайна. Исследование влияния основных параметров Р-сплайна показало, что на полосу пропускания в большей степени влияет число отсчетов в звене сплайна h (число нерекурсивных членов разностного уравнения). Закономерности, установленные при частотном анализе сплайн-фильтра, подтверждены восстановлением в реальном времени доплеровской функции.
Благодарности Работа выполнена при поддержке гранта РФФИ (№ 18-07-01007).
Литература
1. Ланге, П.К. Коррекция динамической погрешности инерционного измерительного преобразователя / П.К. Ланге // Вестник Самарского государственного технического университета. Серия: Технические науки. - 2017. - № 2(54). - С. 58-64.
2. Сергиенко, А.Б. Цифровая обработка сигналов: учебн. пособие для вузов / А.Б. Сергиенко. - 4-е изд. - СПб.: БХВ-Петербург, 2011. - 758 с. - ISBN: 978-5-9775-0606-9.
3. Оппенгейм, А. Цифровая обработка сигналов / А. Оппенгейм, Р. Шафер; пер. с англ. - М.: Техносфера, 2007. - 856 с. - ISBN: 978-5-94836-135-2.
4. Бугров, В.Н Синтез целочисленных рекурсивных фильтров с произвольно заданными селективными требованиями / В.Н. Бугров // Цифровая обработка сигналов. - 2016. - № 2. - С. 35-43.
5. Manolakis, D. Applied digital signal processing: Theory and practice / D. Manolakis, V. Ingle. - Cambridge: Cambridge University Press, 2011. - 1009 p. - ISBN: 978-0521-11002-0.
6. Карпенков, А. С. Методика расчёта целочисленного цифрового селекторного нерекурсивного фильтра c заданными добротностью и уровнем подавления / А.С. Карпенков, Ю.В. Гришанович, Д.С. Потехин, Е.П. Тетерин // Вестник Нижегородского университета им. Н.И. Лобачевского. - 2009. - № 6-1. - С. 79-85.
7. Никитин, Д.А. Приложения алгоритма синтеза рекурсивных цифровых фильтров по импульсной характеристике / Д.А. Никитин // Цифровая обработка сигналов. -2009. - № 4. - С. 8-15.
8. Давыдов, А.В. Цифровая обработка сигналов: Тематические лекции / А.В. Давыдов. - Екатеринбург: ИГиГ, ГИН, Фонд электронных документов, 2005. - 185 с.
9. Турулин, И.И. Методика синтеза управляемых цифровых фильтров на основе аналоговых прототипов / И.И. Турулин, Ю.И. Булгакова // Известия ЮФУ. Технические науки. - 2011. - № 2(115). - С. 88-92.
10. Гетманов, В.Г. Методы вычисления аппроксимацион-ных сплайновых функций для задач цифровой фильтра-
ции / В.Г. Гетманов // Известия Российской академии наук. Теория и системы управления. - 2016. - № 5. -С. 47-57. - DOI: 10.7868/S0002338816040077.
11. Майстренко, А. В. Применение методов цифрового дифференцирования сигналов для определения стационарности процессов / А.В. Майстренко, А.А. Светлаков // Научный вестник Новосибирского государственного технического университета. - 2015. - № 2(59). - С. 7-19. - DOI: 10.17212/1814-1196-2015-2-7-19.
12. Zjavka, L. Multi-site post-processing of numerical forecasts using a polynomial network substitution for the general differential equation based on operational calculus / L. Zjavka // Applied Soft Computing Journal. - 2018. -Vol. 73. - P. 192-202. - DOI: 10.1016/j.asoc.2018.08.040.
13. de Boor, C. A practical guide to splines / C. de Boor. -New York: Springer-Verlag, 2001. - 348 p. - ISBN: 978-0387-95366-3.
14. Redd, A.A. Comment on the orthogonalization of B-spline basis function and their derivatives / A.A. Redd // Statistics and Computing. - 2012. - Vol. 22, Issue 1. - P. 251-257. -DOI: 10.1007/s11222-010-9221-0.
15. Шумилов, Б.М. Алгоритмы с расщеплением вейвлет-преобразования сплайнов первой степени на неравномерных сетках / Б.М. Шумилов // Журнал вычислительной математики и математической физики. - 2016. -Т. 56, № 7. - С. 1236-1247. -10.7868/S0044466916070176.
16. Berenguer-Vidal, R. Design of B-spline multidimensional deformable models in the frequency domain / R. Berenguer-Vidal, R. Verdu-Monedero, J. Morales-Sanchez // Mathematical and Computer Modelling. - 2013. - Vol. 57, Issue 7-8. -P. 1942-1949. - DOI: 10.1016/j.mcm.2012.01.011.
17. Galvez, A. Efficient particle swarm optimization approach for data fitting with free knot B-splines / A. Galvez, A. Iglesias // Computer-Aided Design. - 2011. - Vol. 43, Issue 12. - P. 1683-1692. - DOI: 10.1016/j.cad.2011.07.010.
18. Krivobokova, T. Fast adaptive penalized splines / T. Krivobokova, C.M. Crainiceanu, G. Kauermann // Journal of Computational and Graphical Statistics. - 2008. -Vol. 17, Issue 1. - P. 1-20. - DOI: 10.1198/106186008X287328.
19. Yang, L. Adaptive penalized splines for data smoothing / L. Yang, Y. Hong // Computational Statistics and Data Analysis. - 2017. - Vol. 108. - P. 70-83. - DOI: 10.1016/j.csda.2016.10.022.
20. Kouibia, A. Approximation by discrete variational splines /
A. Kouibia, M. Pasadas // Journal of Computational and Applied Mathematics. - 2000. - Vol. 116, Issue 1. - P. 145156. - DOI: 10.1016/S0377-0427(99)00316-7.
21. Tharmaratnam, K. S-estimation for penalized regression splines / K. Tharmaratnam, G. Claeskens, C. Croux, M. Sa-libian-Barrera // Journal of Computational and Graphical Statistics. - 2010. - Vol. 19, Issue 3. - P. 609-625. - DOI: 10.1198/jcgs.2010.08149.
22. Cao, J. Estimating curves and derivatives with parametric penalized spline smoothing / J. Cao, J. Cai, L. Wang // Statistics and Computing. - 2012. - Vol. 22, Issue 5. - P. 10591067. - DOI: 10.1007/s11222-011-9278-4.
23. Kochegurova, E.A. Real-time recovery of functions and their derivatives by variation splines / E.A. Kochegurova, E.S. Gorokhova // Key Engineering Materials. - 2016. -Vol. 685. - P. 920-924. - DOI: 10.4028/www.scientific.net/KEM.685.920.
24. Eilers, P.H.C. Splines, knots, and penalties / P.H.C. Eilers,
B.D. Marx // WIREs Computational Statistics. - 2010. -Vol. 2, Issue 6. - P. 637-653. - DOI: 10.1002/wics.125.
25. Гетманов, В.Г. Алгоритмы вычисления аппроксимаци-онных сплайновых функций с учётом оптимизации расположения сплайновых узлов / В.Г. Гетманов // Автометрия. - 2013. - Т. 49, № 1. - С. 26-41.
26. Денисов, В.И. Исследование алгоритмов выбора оптимальных координат узловых точек в полупараметрических моделях штрафных сплайнов / В.И. Денисов,
B.С. Тимофеев, А.В. Фаддеенков // Научный вестник Новосибирского государственного технического университета. - 2013. - Т. 51, № 2. - С. 35-44.
27. Aydin, D. Optimum smoothing parameter selection for penalized least squares in form of linear mixed effect models / D. Aydin, M. Memmedli // Optimization. - 2012. - Vol. 61, Issue 4. - P. 459-476. - DOI: 10.1080/02331934.2011.574698.
28. Crainiceanu, C.M. Spatially adaptive Bayesian penalized splines with heteroscedastic errors / C.M. Crainiceanu, D. Ruppert, R.J. Carrol, A. Joshi, B. Goodner // Journal of Computational and Graphical Statistics. - 2007. - Vol. 16, Issue 2. - P. 265-288. - DOI: 10.1198/106186007X208768.
29. Walker, C.G. SALSA - a spatially adaptive local smoothing algorithm / C.G. Walker, M.L. MacKenzie,
C.R. Donovan, M.J. O'Sullivan // Journal of Statistical Computation and Simulation. - 2011. - Vol. 81, Issue 2. -P. 179-191. - DOI: 10.1080/00949650903229041.
30. Чичерин, И.В. Сплайн-алгоритмы обработки сигналов измерительной информации в системах автоматизации технологических процессов : дис. ... канд. техн. наук : 05.13.06 / Чичерин Иван Владимирович. - Новокузнецк, 2006. - 173 с.
31. Хуторцев, В.В. Использование метода сплайн-функций при синтезе цифровых алгоритмов фильтрации с группированием наблюдений / В.В. Хуторцев, О.С. Федо-ренко // Радиотехника. - 2010. - № 2. - С. 4-8.
Сведения об авторах
Кочегурова Елена Алексеевна, 1958 года рождения, в 1980 году окончила Томский политехнический институт по специальности «Прикладная математика», доцент кафедры автоматики и компьютерных систем. Область научных интересов: обработка информации в реальном времени, идентификация клавиатурного почерка, биоинсперированные методы. E-mail: kocheg@tpu.ru .
Даньни У, в 2019 заканчивает обучение в магистратуре Томского политехнического университета по направлению «Управление в технических системах». Область научных интересов: прикладная теория фильтрации, биоинсперированные алгоритмы оптимизации, разработка web-приложений. E-mail: danni_0815@J63.com .
ГРНТИ: 27.41.17 28.15.15 Поступила в редакцию 21 февраля 2018 г. Окончательный вариант - 19 октября 2018 г.
REALIZATION OF A RECURSIVE DIGITAL FILTER BASED ON PENALIZED SPLINES
E.A. Kochegurova1, D. Wu1
1 Tomsk Polytechnic University, Tomsk, Russia
Abstract
In this paper the possibility of development of the recursive digital filter using a P-spline is considered. The frequency and time response of the spline filter for real-time data are analytically obtained and investigated. The influence of the P-spline parameters on effectiveness of interpretation is explored with input metrical data. The patterns obtained during spline filter frequency analysis are confirmed by an example of Doppler function restoration.
Keywords: penalized spline, smoothing spline, digital filter, impulse infinite response (IIR filter), instrumental function, amplitude and phase-frequency response.
Citation: Kochegurova EA, D Wu. Realization of a recursive digital filter based on penalized splines. Computer Optics 2018; 42(6): 1083-1092. DOI: 10.18287/2412-6179-2018-42-6-1083-1092.
32. Кочегурова, Е.А. Текущее оценивание производной нестационарного процесса на основе рекуррентного сглаживающего сплайна / Е.А. Кочегурова, Е.С. Горохова // Автометрия. - 2016. - Т. 52, № 3. - С. 79-85. - DOI: 10.3103/S8756699016030109
33. Кочегурова, Е.А. Частотный анализ рекуррентных вариационных P-сплайнов / Е.А. Кочегурова, А.И. Кочегу-ров, Н.Е. Рожкова // Автометрия. - 2017. - Т. 53, № 6. -С. 67-76. - DOI: 10.15372/AUT20170608.
34. Мясников, В.В. Сплайны как средство построения эффективных алгоритмов локального линейного преобразования / В.В. Мясников // Компьютерная оптика. -2007. - Т. 31, № 2. - С. 52-68.
35. Mihajlovic, Z. Frequency domain analysis of B-spline interpolation / Z. Mihajlovic, A. Goluban, M. Zagar // ISIE '99. Proceedings of the IEEE International Symposium on Industrial Electronics. - 1999. - Vol. 1. - P. 193-198. -DOI: 10.1109/ISIE.1999.801783.
36. Guo, W. Smoothing spline ANOVA for time-dependent spectral analysis / W. Guo, M. Dai, H.C. Ombao, R. Von Sachs // Journal of the American Statistical Association. -2003. - Vol. 98, Issue 463. - P. 643-652. - DOI: 10.1198/016214503000000549.
37. Raposo-Sánchez, M.Á. Analog and digital filters with asplines / M.Á. Raposo-Sánchez, J. Sáez-Landete, F. Cruz-Roldán // Digital Signal Processing: A Review Journal -2017. - Vol. 66. - P. 1-9. - DOI: 10.1016/j.dsp.2017.03.003.
38. Kukushkin, Y.A. Rhythmocardiogram approximation methods for calculation of spectral parameters of cardiac rhythm variability / Y.A. Kukushkin, A.I. Maistrov, A.V. Bogomolov // Biomedical Engineering. - 2010. -Vol. 44, Issue 3. - P. 92-103. - DOI: 10.1007/s10527-010-9165-x.
Acknowledgements: The work was funded by the Russian Foundation for Basic Research Grants under grant No. 18-07-01007.
References
[1] Lange PK. Correction of the dynamic error of inertial sensors [In Russian]. Vestnik of Samara State Technical University. Technical Sciences Series 2017; 2(54): 58-64.
[2] Sergienko AB. Digital signal processing [In Russian]. Saint-Petersburg: "BHV-Petersburg" Publisher; 2013. ISBN: 978-5-9775-0606-9.
[3] Oppenheim AV, Schafer RW. Discrete-time signal processing. 3rd ed. Pearson; 2009. ISBN: 978-0-13-198842-2.
[4] Bugrov VN. Integer design of iir filters with difficult selective requirements [In Russian]. Digital signal processing 2016; 2: 35-43.
[5] Manolakis D, Ingle V. Applied digital signal processing: Theory and practice. Cambridge: Cambridge University Press; 2011. ISBN: 978-0-521-11002-0.
[6] Karpenkov AS, Grishanovich YuV, Pothekhin DS, Teterin EP. Design procedure of integer digital selection nonrecur-sive filter with the given Q-factor and suppression level [In Russian]. Vestnik of Lobachevsky University of Nizhni Novgorod 2009; 6-1: 79-85.
[7] Nikitin DA. An algorithm of IIR-filters synthesis on the pulse response applications [In Russian]. Digital signal processing 2009; 4: 8-15.
[8] Davydov AV. Digital signal processing: Thematic lectures [In Russian]. Ekaterinburg: "USMU", "IgiG", "GIN", "Electronic Documents Foundation" Publishers; 2005.
[9] Turulin II, Bulgakova YuI. The technique of synthesis of operated digital filters on the basis of analog prototypes [In Russian]. Izvestiya SFedU. Engineering Sciences 2011; 2(115): 88-92.
[10] Getmanov VG. Evaluation of spline functions for digital filtering problems. Journal of Computer and Systems Sciences International 2016; 55(5): 725-734. DOI: 10.1134/S1064230716040079.
[11] Maystrenko AV, Svetlakov AA. Application of methods of digital signal differentiation to determine a stationary process [In Russian]. Science Bulletin of the Novosibirsk State Technical University 2015; 2(59): 7-19. DOI: 10.17212/1814-1196-2015-2-7-19.
[12] Zjavka L. Multi-site post-processing of numerical forecasts using a polynomial network substitution for the general differential equation based on operational calculus 2018; 73: 192-202. DOI: 10.1016/j.asoc.2018.08.040.
[13] de Boor C. A practical guide to splines. New York: Springer-Verlag; 2001. ISBN: 978-0-387-95366-3
[14] Redd AA. Comment on the orthogonalization of B-spline basis function and their derivatives. Stat Comput 2012; 22(1): 251-257. DOI: 10.1007/s11222-010-9221-0.
[15] Shumilov BM. Splitting algorithms for the wavelet transform of first-degree splines on nonuniform grids. Comput Math and Math Phys 2016; 56(7): 1209-1219. DOI: 10.1134/S0965542516070174.
[16] Berenguer-Vidal R, Verdú-Monedero R, Morales-Sánchez J. Design of B-spline multidimensional deformable models in the frequency domain. Mathematical and Computer Modelling 2013; 57(7-8): 1942-1949. DOI: 10.1016/j.mcm.2012.01.011.
[17] Gálvez A, Iglesias A. Efficient particle swarm optimization approach for data fitting with free knot B-splines. Computer-Aided Design 2011; 43(12): 1683-1692. DOI: 10.1016/j.cad.2011.07.010.
[18] Krivobokova T, Crainiceanu CM, Kauermann G. Fast adaptive penalized splines. Journal of Computational and
Graphical Statistics 2008; 17(1): 1-20. DOI: 10.1198/106186008X287328.
[19] Yang L, Hong Y. Adaptive penalized splines for data smoothing. Computational Statistics and Data Analysis 2017; 108: 70-83. DOI: 10.1016/j.csda.2016.10.022.
[20] Kouibia A, Pasadas M. Approximation by discrete variational splines. Journal of Computational and Applied Mathematics 2000; 116: 145-156. DOI: 10.1016/S0377-0427(99)00316-7.
[21] Tharmaratnam K, Claeskens G, Croux C, Salibian-Barrera M. S-estimation for penalized regression splines. Journal of Computational and Graphical Statistics 2010; 19(3): 609-625. DOI: 10.1198/jcgs.2010.08149.
[22] Cao J, Cai J, Wang L. Estimating curves and derivatives with parametric penalized spline smoothing. Statistics and Computing 2012; 22(5): 1059-1067. DOI: 10.1007/s11222-011-9278-4.
[23] Kochegurova EA, Gorokhova ES. Real-time recovery of functions and their derivatives by variation splines. Key Engineering Materials 2016; 685: 920-924. DOI: 10.4028/www.scientific.net/KEM.685.920.
[24] Eilers PHC, Marx BD. Splines, knots, and penalties. WIREs Computational Statistics 2010; 2(6): 637-653. DOI: 10.1002/wics.125.
[25] Getmanov VG. Algorithms of calculating of approximating spline functions with optimization of the location of spline nodes. Optoelectronics, Instrumentation and Data Processing 2013; 49(1): 21-33. DOI: 10.3103/S8756699013010044.
[26] Denisov VI, Timofeev VS, Faddeenkov AV. Investigation of algorithms for choosing optimal coordinates of nodal points in semiparametric models of penalty splines [In Russian]. Science Bulletin of the Novosibirsk State Technical University 2013; 51 (2): 35-44.
[27] Aydin D, Memmedli M. Optimum smoothing parameter selection for penalized least squares in form of linear mixed effect models. Optimization 2012; 61(4): 459-476. DOI: 10.1080/02331934.2011.574698.
[28] Crainiceanu CM, Ruppert D, Carrol RJ, Joshi A, Goodner B. Spatially adaptive Bayesian penalized splines with het-eroscedastic errors. Journal of Computational and Graphical Statistics 2007; 16(2): 265-288.- DOI: 10.1198/106186007X208768.
[29] Walker CG, MacKenzie ML, Donovan CR, O'Sullivan MJ. SALSA - a spatially adaptive local smoothing algorithm. Journal of Statistical Computation and Simulation 2011; 81(2): 179-191. DOI: 10.1080/00949650903229041.
[30] Chicherin IV. Spline-algorithms for processing signals of measuring information in systems of automation of technological processes [In Russian]. The thesis for the Candidate's degree in Technical Sciences. Novokuznetsk, 2006.
[31] Khutortsev VV, Fedorenko OS. Using of the spline-function method for syntheses of digital filtering algorithms with grouping of the observations [In Russian]. Radioengineering 2010; 2: 4-15.
[32] Kochegurova EA, Gorokhova ES. Current estimation of the derivative of a nonstationary process based on a recurrent smoothing spline. Optoelectronics, Instrumentation and Data Processing 2016; 52(3): 280-285. DOI: 10.3103/S8756699016030109.
[33] Kochegurova EA, Kochegurov AI, Rozhkova NE. Frequency analysis of recurrent variational P-splines. Optoe-
lectronics, Instrumentation and Data Processing 2017; 53(6): 591-598. DOI: 10.3103/S8756699017060085.
[34] Myasnikov VV. Splines as a tool for constructing efficient algorithms of a local linear transform [In Russian]. Computer optics 2007; 31(2): 52-68.
[35] Mihajlovic Z, Goluban A, Zagar M. Frequency domain analysis of B-spline interpolation. Proc IEEE International Symposium on Industrial Electronics 1999; 1: 193-198. DOI: 10.1109/ISIE.1999.801783.
[36] Guo W, Dai M, Ombao HC, Von Sachs R. Smoothing spline ANOVA for time-dependent spectral analysis. Jour-
nal of the American Statistical Association 2003; 98(463): 643-652. DOI: 10.1198/016214503000000549.
[37] Raposo-Sánchez MA, Sáez-Landete J, Cruz-Roldán F. Analog and digital filters with a-splines. Digital Signal Processing: A Review Journal 2017; 66: 1-9. DOI: 10.1016/j.dsp.2017.03.003.
[38] Kukushkin YA, Maistrov AI, Bogomolov AV. Rhythmo-cardiogram approximation methods for calculation of spectral parameters of cardiac rhythm variability. Biomedical Engineering 2010; 44(3): 92-103. DOI: 10.1007/s10527-010-9165-x.
Authors' information
Elena Alekseevna Kochegurova (b. 1958). Candidate of Technical Sciences, associate professor of Information Technology department, Tomsk Polytechnic University. Research interests are information processing in real time, keyboard dynamics, bio-inspected techniques. E-mail: kocheg@mail.ru .
Wu Danni, will graduate fiom Tomsk Polytechnic University in 2019 (master program Control Systems). Her research interests include filter theory, bio-inspected techniques, and web application development. Email: danni_0815@J63.com .
Received February 21, 2018. The final version - October 19, 2018.
