Автоматическое определение интерполянта с наименьшим количеством параметров Текст научной статьи по специальности «Математика»

УДК 517.518.8

Д. А. Никитин, К. В. Сафонов

АВТОМАТИЧЕСКОЕ ОПРЕДЕЛЕНИЕ ИНТЕРПОЛЯНТА С НАИМЕНЬШИМ КОЛИЧЕСТВОМ ПАРАМЕТРОВ*

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

Ключевые слова: интерполяция, идентификация, рекурсивные цифровые фильтры.

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

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

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

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

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

Суть данного метода заключается в следующем. Исходный набор значений принимается за начало импульсной характеристики рекурсивного цифрового фильтра неизвестного порядка с одинаковым количеством коэффициентов ак и Ьк. Этот фильтр рассчитывается таким образом, чтобы минимизировать его порядок. Затем по полученным коэффициентам можно определить вид функциональной зависимости, которой подчиняются члены исходной последовательности. А после этого, опять же только по коэффициентам, можно полностью выписать эту функцию. Выявленная таким способом функция будет иметь не известный заранее вид, а будет наиболее простой из некоторого множества функций и при этом с минимально возможным количеством параметров.

Для синтеза цифрового фильтра (ЦФ) в данной постановке задачи воспользуемся уравнением цифровой фильтрации:

N М

У( п) = Х Ькх(п - к) -^ акУ(п - к). (1)

к=0 к=1

Согласно постановке задачи, М = N + 1. Так как последовательность У = {ук} является импульсной характеристикой и в нашей постановке задачи известна заранее, подставим в (1) вместо х, - единичный импульс, а вместо у, - элементы исходной последовательности. Получим систему линейных уравнений

Г1

- >0 - >1

-Ун-1 - Ун

- Уо

- ун-2

- Ун-1

ун+м-1 ун+м-2

-Уо - У1

V Ь Л

°о

-Уо

Ум-2 Ум-1 у Vам У Ун +м

н

а

(

Уо >1 >2 >> Ун+1

Л

(2)

*Работа выполнена в рамках реализации Федеральной целевой программы «Научные кадры и научно-педагогические кадры инновационной России на 2009-2013 годы» (контракт П1032 от 27.05.2010).

Неизвестными в системе (2) являются коэффициенты ак и Ьк, а также минимально возможные М и Ы, при которых система имеет единственное решение.

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

В рамках данной работы был разработан алгоритм расчета рекурсивных цифровых фильтров. Более подробное его описание приводится в [2-3]. Сначала определяется необходимый порядок фильтра, а затем решается система (2).

Доказано [2; 4], что данный алгоритм синтеза будет иметь решение как раз в случае, когда значения элементов последовательности У являются значениями функции, взятыми на равномерной сетке. Такое доказательство проведено для нескольких видов функций (любых полиномиальных функций, любых экспоненциальных и любых синусоид), отдельно для каждого вида. Эти три вида функций составляют множество функций, среди которых будет отыскиваться интерполянт с наименьшим количеством параметров. Приведенные ниже теоремы также дают и порядок, соответствующий функции конкретного вида. А значит, устанавливается однозначное соответствие между функцией и количеством коэффициентов соответствующего фильтра.

Таким образом, алгоритм расчета фильтра является первым шагом предлагаемого метода автоматической интерполяции. Фильтр, найденный таким образом, предлагается называть фильтром, соответствующим последовательности У.

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

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

Для показательной и синусоидальной функций правила идентификации следующие.

Для показательной функции:

М = 2, N = 1, а1 + а2 = -1, а2 > 0, а2 Ф 1. (3)

Для синусоидальной функции:

М = 3, N = 2, а! = - а2, -1 < а2 < 3, а3 = -1. (4)

После идентификации вида интерполянта необходимо определить все его параметры (коэффициенты в записи функции). Для полинома каждой степени зависимость коэффициентов прямой связи (Ьк) соответствующего фильтра от значений параметров полинома однозначна и определяется в ходе решения системы (2). Данная зависимость линейна, таким образом, для определения обратной зависимости необходимо всего лишь решить систему линейных уравнений. Решить такую систему для полинома степени п необходимо всего один раз - полученные формулы могут быть использованы для всех полиномов такой же степени. Несколько примеров для полиномов младших степеней приведено в табл. 2. Для остальных видов функций - показательной и синусоидальной - формулы также вычислены, хотя зависимость уже и не является линейной.

Для показательной функции у(х) = к-аЬх + с:

к = Ь1 + Ь0 а1,

а = -а1 -1, с = ф0 + ¿і)/(і - аі).

(5)

Для синусоидальной функции у(х) = а-$>т(Ьх + с) + ё:

= Ь2(2 - а2) - Ь0 - Ь1

(3 - а2)8т

аг^

Ь2 (2 - а2) - Ь0 - Ь1 а2 +1

Ь0 - Ь1 - а2Ь2 V 3 - а2

а2 -1

Ь = агссо8—-------,

с = аг^

(6)

Ь2(2 - а2) - Ь0 - Ь1 а2 +1

_ Ьо - Ь1 - а2Ь2 ^ і о4 I 3

а =

Ь0 + ь + Ь2

3 - а2

Значение аргумента арккосинуса в (6) принимает

только допустимые значения

а2 - 1

2

<1 в соответст-

вии с условиями (4) для синусоидальной функции (так как -1 < а2 < 3).

2

Таблица 1

Коэффициенты ОС фильтров, соответствующих полиномам

Интерполирующая функция Формула Порядок соответствующего фильтра Коэффициенты ак

Линейная у = кх + Ь 2 (-2 1)

Квадратичная у = ах2 + Ьх + с 3 (-3 3 -1)

Кубическая у = ах3 + Ьх2 + сх + а 4 (-4 6 -4 1)

Полином 4-й степени у = ах4 + Ьх3 + сх + ах + е 5 (-5 10 -10 5 -1)

Таблица 2

Зависимость параметров полиномов от коэффициентов Ък соответствующего рекурсивного цифрового фильтра

Интерполирующая функция Формула Выражения для вычисления параметров

Линейная у і = ^і + 5, 50 = Ъ0 + Ъ,, 5, = -Ъ,

Квадратичная 2 Уі = Sо■l + .517 + 52 5 Ъ0 + Ъ1 + Ъ2 5 Ъ0 - Ъ1 - 3Ъ2 „ = Ъ 50 = 2 ’ 5 = 2 ’

Кубическая Уі = Е ■ік к=0 Ъ0 + Ъ1 + Ъ2 + Ъ3 Ъ0 - Ъ2 - 2Ъ3 50 = 6 ’ 51 = 2 2Ъ0 - Ъ, + 2Ъ, + 11Ъ3 — 0 1 2 3 , = -Ъ3 2 6 3 3

Полином 4-й степени 4 Vі к Уі = А 5к ■1 к—0 „ Ъ0 + Ъ1 + Ъ2 + Ъ3 + Ъ4 0 24 3Ъ0 + Ъ1 - Ъ2 - 3Ъ3 - 5Ъ4 51 — 1 12 _ 11Ъ0 - Ъ, - Ъ2 + 11Ъ3 + 35Ъ4 2 24 _ 3Ъ0 - - 3Ъз - 25Ъ4 ^ = , оо — , ¿4 і?4 12 4 4

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

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

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

Далее кратко изложим доказанные теоремы.

Теорема 1. Для последовательности У длины Ь > 2(п + 1), члены которой являются равноотстоящими отсчетами полиномиальной функции одного переменного у і = рп(і) = Оо' іп + а1 - г'”-1 + ... + ап, ао Ф 0, и заданы любые (п + 1) подряд идущих члена Ут, Ут+ъ- -, Ут+п, существует единственный рекурсивный цифровой фильтр порядка п + 1 (М = п + 1, N = п), импульсная характеристика которого совпадает с У.

Доказательство. Часть 1. Пусть т = 0, а отсчеты функции берутся с шагом, равным единице. Используя выражение (1), построим начальную систему уравнений для М + N + 1 первых элементов последовательности У, приняв ее за импульсную характеристику искомого ЦФ. Матрица системы:

Р —

і 0 0 • 0 0 0 - 0 0

0 і 0 • 0 -Рп (0) 0 •• 0 0

0 0 і • 0 -Рп (і) -Рп (0) 0 0

0 0 0 • • і -Рп (п - 1) - Рп (п - 2) • - Рп (0) 0

0 0 0 • 0 -Рп (п) - Рп (п - 1) • -Рп (і) - Рп (0)

0 0 0 • 0 -Рп (п + 1) -Рп (п) • - Рп (2) -Рп (і)

0 0 0 • 0 -Рп (2п - 1) -Рп (2п - 2) • -Рп (п) -Рп (п - 1)

0 0 0 • 0 -Рп (2п) -Рп (2п - 1) •• • -Рп (п + 1) -Рп (п)

Элементарными преобразованиями приведем матрицу к верхнетреугольному виду. Если построить эту систему для большего числа уравнений, то она аналогичными преобразованиями будет приводиться к такому же треугольному виду - все строки, начиная с (п + 1)-й, будут нулевыми.

Теперь рассмотрим расширенную матрицу этой системы. При выполнении аналогичных преобразований матрица приводится к аналогичному треугольному виду - дополнительный столбец оказывается нулевым.

Кроме того, все элементы на главной диагонали этих треугольных матриц оказываются равными гп(ак) = -п! ао. Таким образом, система совместна (так как ранги матриц системы равны) и имеет единственное решение (так как система обладает полным рангом при а0 Ф 0, что соответствует условию теоремы).

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

У, = X а1 ' ' 51еР)П-1 = X (а1 ' 51еРП-1)' (г')П-1 =

Доказательство. Рассмотрим сначала случай, когда У является отсчетами показательной функции, взятой с шагом, равным единице - , = 0, 1, ..., Ь - 1. Используя выражение (1), построим начальную систему уравнений для первых четырех элементов последовательности У, приняв ее за импульсную характеристику искомого ЦФ (т. е. х(п) = {1, 0, 0, 0, ...}). Элементарными преобразованиями матрица системы приводится к треугольному виду:

Р =

(1 0 0 0 Л

0 1 к + с 0

0 0 к • аЬ + с к + с

у 0 0 к • а2Ь + с к • аЬ + с

1 0 0 0 Л

0 1 к + с 0

0 0 к • (аЬ -1) к + с

0 0 0 с • (1 - аЬ )

Аналогично для расширенной матрицы системы:

= Ха} ’ ,п 1, , = 0, Ь -1, ] = 0, п.

1

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

Часть 2. Пусть т > 0. Для любых п + 1 фиксированных точек существует единственный полином, проходящий через них. Следовательно, зная у,, , = т, т + п , можно построить систему уравнений:

Р =

(1 0 0 1

0

к + с

0 0 к • а + с

0

0

к + с

к + с Л

; Ь , к • а + с

, 2Ь .

к • а + с

у0 0 к • а2Ь + с к • аЬ + с к • азЬ + с у

(1 0 0 1

0

(к + с)

0 0 к • (аЬ -1)

0 0 к + с

0 Л

0

0

0 0

0

с • (1 - аЬ ) 0

Xа1 • -1 = Ут ,

1

XаІ • (т + 1)И-1 = Ут+1>

Xа, • (т + «)И 1 = Ут

Данная система имеет единственное решение относительно а для фиксированного т > 0. Значит, если фиксированы у,, , = т, т + п, то фиксированы коэффициенты полинома а. А если фиксированы коэффициенты полинома, то фиксированы у, , = 0, п , так как они однозначно выражаются через коэффициенты полинома. Следовательно, при т > 0 справедливо доказательство из части 1 данной теоремы.

Теорема 2. Для последовательности У = у(,) длины Ь > 4, члены которой являются равноотстоящими отсчетами показательной функции у, = _Д/) = каъ' + с, а Ф 0, а Ф 1, Ъ Ф 0, к Ф 0, и заданы первые четыре члена у0, уь у2, у3, существует единственный рекурсивный цифровой фильтр второго порядка (М = 2, Ж = 1), импульсная характеристика которого совпадает с У.

Таким образом, ранги матриц Р и Р равны при любых параметрах к, а, Ь и с, а ф 0, а ф 1, Ь ф 0, к ф 0. Значит, система определена и совместна, т. е. имеет единственное решение.

Построим теперь те же матрицы, но с Ь строками, Ь > М + N + 1. После аналогичных элементарных преобразований матрицы приводятся к следующему треугольному виду:

Р-

1 0 0

0 1 к + с

0 0 к • (аЬ -1)

0 0

0 0 0 0

0 0 (1 0 0 1

Р

0

0

к + с

0 0 к • (аЬ -1)

0 0 0

0 0 0

0 0 0

0 0

к + с с • (1 - аЬ) 0

0

0 0

к + с

0 Л 0 0

с • (1 - аЬ) 0

0

0

Следовательно, га^(Р) = га^( Р ) = 4, и добавленные Ь - М - N - 1 уравнений линейно зависимы от уравнений начальной системы. Поэтому решение, найденное для начальной системы, состоящей из четырех уравнений, будет удовлетворять последовательности отсчетов показательной функции, любой длины большей либо равной четырем. Это решение представляет собой четыре коэффициента рекурсивного цифрового фильтра: Ь0, Ь1, а1, а2.

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

у. = к • аь(г"‘ер) + с = к • а(ь"‘ер>г + с = к • аь'-г + с, г = 0, ..., Ь - 1.

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

коэффициентом в показателе степени. Следовательно, доказательство справедливо при любом постоянном шаге.

Теорема 3. Для последовательности У = у(г) длины Ь >6, члены которой являются равноотстоящими отсчетами синусоидальной функции Уг = /(') = а^ш(Ь-/ + с) + ё, а ф 0, Ь ф пк, к е2, ёф 0, и заданы первые шесть членов у, у1, ., у5, существует единственный рекурсивный цифровой фильтр третьего порядка (М = 3, N = 2), импульсная характеристика которого совпадает с У.

Доказательство. Рассмотрим сначала случай, когда У является отсчетами синусоидальной функции, взятой с шагом, равным единице - г = 0, 1, . , Ь - 1. Используя выражение (1), построим начальную систему уравнений для первых шести элементов последовательности У, приняв ее за импульсную характеристику искомого ЦФ (т. е. х(п) = {1, 0, 0, 0, ...}). Элементарными преобразованиями матрица системы приводится к следующему треугольному виду:

P =

(1 0 0 1 0 0 00 00

(1 0 0 0 0 0

0 1 0 a sin c + d 0 0

0 0 1 a sin(b + c) + d a sin c + d 0

0 0 0 a sin(2b + c) + d a sin(b + c) + d a sin c + d

0 0 0 a sin(3b + c) + d a sin(2b + c) + d a sin(b + c) + d

V 0 0 0 a sin(4b + c) + d a sin(3b + c) + d a sin(2b + c) + d y

0 0 0 0

0 a sin c + d 0 0

a(sin(b + c) + d a sin c + d 0

2d(1 - cosb) a sinb cos c + d(1 - cosb) a sinc + d

0 a sinb(cos(b + c) - cos c) a(sin(b + c) - sinc)

4 a sin b sin2 —

0 0 0

cos(b + c) - cos c

т. е. га^(Р) = 6 при любых а ф 0, Ь ф пк, к е2, ёф 0.

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

Построив те же матрицы, но с Ь строками, Ь > М + N + 1, убеждаемся, что га^(Р) = га^( Р ) = 6, и добавленные Ь - М - N - 1 уравнений линейно зависимы от уравнений начальной системы. Поэтому решение, найденное для начальной системы, состоящей из шести уравнений, будет удовлетворять последовательности отсчетов синусоидальной функции, любой длины большей либо равной шести. Это решение

представляет собой шесть коэффициентов рекурсивного цифрового фильтра: b0, b1, b2, a1, a2, b3.

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

yi = a • sin(b • (i • step) + c) + d =

= a • sin((b • step) • i + c) + d = a • sin(b' • i + c) + d,

I = 0,..., L - 1.

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

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

Библиографические ссылки

1. Марпл-мл. С. Л. Цифровой спектральный анализ и его приложения. М. : Мир, 1990.

2. Никитин Д. А., Ханов В. Х. Синтез рекурсивных цифровых фильтров по импульсной характеристике, определяемой элементарной математической функцией // Цифровая обработка сигналов. 2008. № 3. С. 10-14.

3. Ханов В. Х., Никитин Д. А. Алгоритм анализа числовых последовательностей // Вестник СибГАУ 2006. Вып. 6(13). С. 11-15.

4. Никитин Д. А. Теоремы о существовании и порядках цифровых рекурсивных фильтров с импульсными характеристиками определенной формы // Ин-форм. технологии и мат. моделирование (ИТММ-2009) : материалы VIII Всерос. науч.-практ. конф. с междунар. участием (13-14 нояб. 2009 г.). Ч. 2. Томск, 2009. С. 144-146.

D. A. Nikitin, K. V. Safonov

AUTOMATIC DETECTION OF INTERPOLANT WITH THE FEWEST PARAMETERS

A method that allows, for the initial set of points on a uniform grid, to determine automatically an interpolant with the fewest parameters in the following set offunctions: polynomials, exponential functions, sine, any linear combination of the above functions.

Keywords: interpolation, identification, digital recursive filters.

© Никитин Д. А., Сафонов К. В., 2011

УДК 004.932.2

Н. Ю. Петухов

РАСПОЗНАВАНИЕ ТЕКСТУРНЫХ ИЗОБРАЖЕНИЙ НА ОСНОВЕ СТАТИСТИЧЕСКИХ И ФРАКТАЛЬНЫХ ПРИЗНАКОВ

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

Ключевые слова: текстура, фрактальная размерность, статистические характеристики, нейронная сеть.

В последнее время все большую актуальность приобретают междисциплинарные исследования,

в частности, в современных науках о Земле исследования все чаще проводятся на стыке различных научных направлений [1]. Для комплексной оценки современного состояния различных природных структур необходимо проводить исследования с использованием не только классических научных методов, разработанных и апробированных в системе наук о Земле, но и новейших физических, математических, компьютерных знаний и технологий, которые позволяют моделировать и прогнозировать возможные тенденции в изменении структуры и свойств природных объектов. Одним из таких инструментов, позволяющих анализировать современное состояние природных объектов, является фрактальный анализ. Данный метод позволяет оценить характер самоподобия природного объекта, раскрыть его фрактальные свойства.

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

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

Алгоритмическая реализация фрактальных методов оценки текстур. На практике для измерения фрактальной размерности обычно используют три алгоритма: метод покрытия поверхности эталонами; дисперсионное масштабирование, основанное на