Сглаживающие кубические сплайны (СКС) находят широкое применение для фильтрации шумов измерений экспериментальных данных или, другими словами, для сглаживания зашумленных данных [1, 2]. К основному преимуществу СКС по сравнению с другими алгоритмами фильтрации можно отнести тот факт, что сглаживающий сплайн и его первые две производные являются непрерывными функциями на всем интервале построения сплайна. Это позволяет использовать СКС для устойчивого вычисления производных по зашумленным измерениям обрабатываемой функции, в частности в задачах непараметрической идентификации как линейных [3–5], так и нелинейных [6, 7] динамических систем. Для последнего класса систем СКС используется для вычисления смешанных производных [8]. СКС имеет некоторый управляющий параметр – параметр сглаживания, от величины которого существенно зависит ошибка сглаживания зашумленных данных [2, 9]. Выбор этого параметра является основной проблемой при использовании сглаживающих сплайнов на практике. Предложены несколько алгоритмов выбора, позволяющие с приемлемой для практики точностью оценить оптимальное значение параметра сглаживания (ошибка сглаживания минимальна) при различной априорной информации о дисперсии шума измерений [2, 8]. Хотя СКС по принципу своей работы является низкочастотным фильтром (удаляющим высокочастотные – шумовые составляющие обрабатываемого сигнала) традиционный подход к построению СКС не позволяет интерпретировать СКС как такой фильтр и, как следствие, решать задачи анализа и синтеза СКС в терминах низкочастотного фильтра и его характеристик, которые естественны и понятны экспериментатору.

Поэтому в данной работе предлагается частотная модель СКС и вводятся динамические характеристики СКС, позволяющие решать задачи анализа и синтеза СКС как низкочастотного фильтра, а именно:

− строится амплитудно-частотная характеристика СКС, включающая все основные параметры СКС;

− вводятся понятия аппаратной функции сплайна и эквивалентной ширины аппаратной функции сплайна;

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

Материалы и методы исследования

В начале приведём некоторые определения, нужные для построения частотной модели СКС [1, 9]. Предположим, что на интервале [T1,T2] заданы n узлов и в этих узлах измерены значения функции:

, (1)

где ηi – случайный шум (погрешность) измерений с нулевым средним и дисперсией (равноточные измерения). Функция Sn,α(t) называется сглаживающим кубическим сплайном дефекта единица, если: а) на каждом отрезке функция Sn,α(t) является кубическим многочленом вида

(2)

б) функция Sn,α(t) дважды непрерывно дифференцируема на всём интервале [T1,T2].

Сглаживающий кубический сплайн Sn,α(t) в общем случае не проходит через точки , а проходит в некоторой окрестности (ее размеры определяются величиной параметра сглаживания) от этих точек, обеспечивая тем самым сглаживание (фильтрацию) шума измерений. Для однозначного определения коэффициентов СКС (2) в узлах t1, tn задаются краевые условия сплайна, из которых наиболее часто используют:

− условия на значения первой производной вида

; (3)

− условия на нулевые значения второй производной

; (4)

− или сочетание условий (3), (4) (например, в узле t1 задается (3), а в tn – условие (4)).

Доказано, что СКС, построенный при краевых условиях (3) или (4), доставляет минимум функционалу [1]:

, (5)

где pj – весовые множители, отражающие точность j-го измерения (в случае равноточных измерений задаются одинаковыми, например ), α – параметр сглаживания.

Перейдем к построению частотной модели СКС. Для этого представим СКС в виде «черного» ящика, на вход которого поступает последовательность измеренных значений , а выходом является функция Sn,α(t) с непрерывными производными до второй включительно и определенная на интервале [T1,T2]. При этом будем полагать, что шаг сетки равный и равен Δt и весовые множители pj одинаковы и равны p.

Первоначально установим связь между значениями , и вторыми производными Mj сплайна Sn,α(t) в узлах tj. За исключением двух первых и двух последних узлов эта связь выражается следующим разностным уравнением [1, 2]:

. (6)

В формуле (6)

. (7)

Уравнение (6) отражает структуру дискретного фильтра, преобразующего входную последовательность {f̃j} в последовательность {Mj}. Для построения частотной характеристики такого фильтра обратимся к так называемому D-преобразованию бесконечной последовательности , определяемому соотношением

.

где . Заметим, что HD(ω) является периодической функцией с периодом 2π / Δt. Обозначим через – D-преобразования дискретных последовательностей {Mj}, {f̃j}. Тогда частотная характеристика рассматриваемого дискретного фильтра определяется как отношение . С учетом свойств D-преобразова-ния для сдвинутых последовательностей [ ], получаем частотную характеристику дискретного преобразования {f̃j} ⇒ {Mj}:

(8)

Далее рассмотрим преобразование дискретных значений второй производной в непрерывную вторую производную сплайна:

, (9)

что соответствует линейной интерполяции с частотной характеристикой:

. (10)

Для нахождения Sn,α(t) по второй производной выполним двойное интегрирование, что соответствует преобразованию с частотной характеристикой:

. (11)

Очевидно, что частотная характеристика СКС равна произведению частотных характеристик (8), (10), (11) и определяется выражением

. (12)

Характеристика (12), по сути, является амплитудо-частотной характеристикой динамического звена {f̃j} ⇒ {Mj}. Введем понятие аппаратной функции СКС kα(t) как обратное преобразование Фурье от частотной характеристики (12):

. (13)

Возникает вопрос: как использовать понятия частотной характеристики и аппаратной функции СКС для его анализа и синтеза? Предположим, что функция f(t) имеет Фурье-образ, равный ФF(ω). Тогда Фурье-образ ФSα(ω) сплайна Sn,α(t), построенного по точным значениям {fj = f(tj)}, будет определяться выражением

ФSα(ω) = Kα(ω)ФF(ω), – ∞< ω < ∞. (14)

В соответствии с теоремой о свертке [10, с. 527–529], сплайн Sn,α(t) можно представить как результат интегрирования:

(15)

и его можно интерпретировать как «усреднение» функции f(t) аппаратной функцией сплайна. Очевидно, что чем меньше «ширина» аппаратной функции, тем меньше Sn,α(t) отличается от f(t) и тем меньше систематическая ошибка сглаживания. Для количественной оценки ширины аппаратной функции введем характеристику, которую назовем эквивалентной шириной аппаратной функции (ЭШАФ) W(α):

. (16)

Величина W(α) уже позволяет говорить о разрешающей способности СКС, а именно: если функция имеет две δ-образных структуры, то они проявятся в СКС (возможно, в сглаженном виде), если «расстояние» (по шкале аргумента t) между ними будет не меньше W(α) (см. результаты вычислительного эксперимента).

Замечание. В силу конечной длительности функций f(t), k(τ) пределы интегрирования ta, tb, tc, td в (15), (16) имеют конечные значения, что позволяет вычислить (15), (16), применяя подходящие квадратурные формулы.

Таким образом, введение характеристики W(α) позволяет говорить об анализе разрешающей способности СКС – когда по заданному параметру сглаживания вычисляется значение W(α) и о задаче синтеза СКС – когда по заданному максимальному значению W(α) вычисляется величина параметра сглаживания. Так как сглаживающая способность сплайна увеличивается по мере возрастания параметра сглаживания и W(α) при этом также увеличивается, то, решая задачу синтеза сплайна по заданной максимальной величине W(α), можно говорить о построении сплайна с наименьшим уровнем случайной ошибки среди всех сплайнов с ЭШАФ не более W(α).

Результаты исследования и их обсуждение

Первоначально рассмотрим результаты вычисления частотной характеристики (12) и аппаратной функции (13) сплайна с параметрами Δt = 0.05, p = 1. На рис. 1 изображены графики частотной характеристики Kα(ω) при разных параметрах сглаживания: сплошная кривая – α = 0; штриховая кривая – α = 0.0001. точечная кривая – α = 0.1. На рис. 2 показаны графики аппаратных функций сплайна при тех же параметрах сглаживания (соответствие типов кривых те же, что и на рис. 1). Хорошо видно, что при увеличении параметра сглаживания уменьшается ширина полосы пропускания СКС (как низкочастотного фильтра), а следовательно, увеличивается подавление шума измерений, но при этом увеличивается ширина аппаратной функции, что вызывает ухудшение разрешающей способности сплайна – известное противоречие между систематической и случайной составляющей ошибки фильтрации любого линейного алгоритма фильтрации.

Перейдем к рассмотрению эквивалентной ширины аппаратной функции СКС, определяемой соотношением (16). На рис. 3 представлен график характеристики W(α) от параметра сглаживания. Используя этот график, попробуем определить параметр сглаживания по заданной величине W(α) = WS = 0.2. Для этого по оси ординат (рис. 3) откладываем величину 0.2 и проводим штриховую прямую до пересечения с графиком W(α) (точка пересечения обозначена буквой А), а затем из этой точки опускаем перпендикуляр на ось абсцисс – и точку пересечения с этой осью обозначим как В. Значение α в точке В примерно равно α0.2 = 3 ∙ 10–4, и, по сути, это значение является «графическим» решением нелинейного уравнения W(α) = WS.

Рис. 1. https://s.top-technologies.ru/pic/2022/5-1ная характеристика СКС при различных значениях параметра α

Рис. 2. Аппаратная функция СКС при различных значениях параметра α

Рис. 3. Зависимость характеристики W(α) от параметра сглаживания α

Рис. 4. СКС при заданной разрешающей способности

На рис. 4 показаны графики четырех δ-образных функций. Расстояние между левой парой кривых равно 0,2, т.е. равно заданной величине WS, а расстояние между правой парой равно 0.5. Точечной кривой на рис. 4 показаны значения сглаживающего сплайна, построенного при α0.2 = 3 ∙ 10–4.

Видно, что в сплайне четко «проявились» две левых δ-образных функций, расстояние между которыми равно заданной величине ЭШАФ построенного СКС. Следовательно, по величине ЭШАФ сплайна можно определить расстояние между «тонкими» деталями функции, сохраненные в построенном СКС, т.е. можно говорить о разрешающей способности СКС, которую можно характеризовать ЭШАФ. Так как расстояние между правыми δ-образных функций в 2,5 раза больше заданной ЭШАФ, то эти функции полностью сохранились в построенном СКС.

Заключение

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

− задачу анализа СКС – когда по заданному параметру сглаживания определяют введенные частотные и временные характеристики сплайна;

− задачу синтеза СКС – когда по заданной ЭШАФ WS находят значение параметра сглаживания СКС путем решения нелинейного уравнения W(α) = WS численными методами.

Исследование выполнено при финансовой поддержке РНФ в рамках научного проекта № 22-21-00409.