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

УДК 517.518.85

В. Н. Исаков

Оптимальная регулярная локальная сплайновая интерполяция сигналов

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

Ключевые слова: интерполяция сигнала, интерполирующая функция, интервал дискретизации, сплайн-интерполяция, сплайн-базис.

о см

■ч-

О!

<

I

(0 та

0 ^

СО та

1

о.

<и

3

и <и со

см ■ч-ю

с?

см ■ч-ю см

Введение

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

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

Классические математические подходы к задачам интерполяции подробно изложены в работе [1].

Одно из перспективных современных направлений в теории интерполирования -сплайновая интерполяция. Математические свойства сплайнов подробно исследованы в работах [2, 3]. Практические аспекты сплай-новой интерполяции рассмотрены в работе [4]. Однако основная часть указанных работ по сплайнам посвящена глобальной интерполяции. Локальная сплайновая интерполяция применительно к характеристикам летательных аппаратов рассмотрена в работе [5], применительно к сигналам - в источниках [6-12].

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

При регулярной интерполяции сетка дискретизации равномерна и бесконечна, а на каждом частном интервале дискретизации формирование интерполирующей функции осуществляется по одним и тем же правилам. Обозначим координаты узлов интерполяции ^ = пТ, sn = s(nT), где п = 0, ±1, ±2, ...; Т -период дискретизации сигнала. Запишем выражение для интерполирующей функции в виде

) = X 5пФоО - К ) =

= X $(пТ)фс(* - пТ),

(1)

где порождающая функция ф0(7) удовлетво ряет условию

[1, к = 0,

10, к = ±1, ±2, ±3, ....

Фо(^) =

(2)

- © Исаков В. Н., 2016

Сходимость интерполяционного процесса при неограниченном уменьшении периода

п=-со

дискретизации обеспечивается при выполнении условия [13]

Фг

2п k

\T, k = 0,

10, k = ±1, ± 2, ± 3, ...,

(3)

где Ф 0(ю) - спектральная плотность ф0(/).

Из условий (2), (3) следуют нормировки

|фо(/)dt = Т,

—со

+<х>

| Фю = 2п.

-оэ

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

Локальные фундаментальные интерполяционные обобщенные сплайн-базисы

При сплайновой интерполяции фрагменты интерполирующей функции уп (), описывающие ее на каждом частном интервале дискретизации / е \}п; tn+l], стыкуются в узлах интерполяции так, чтобы вместе с условием согласованности (2) обеспечить непрерывность ее первых г производных (значение г характеризует степень гладкости):

v„ (tn) = sn, vи On+O = i.

я+1'

Vn (tn ) ^ dm )sn+m ,

m=-M

M

Vn ОП+l) = ^ dm1 Sn+l+m ,

m=-M

• • • i

M

Wn^i^n ) _ ^ dm ) Sn+m ,

m=-M

M

v = £

(4)

n+1+m'

m=-M

где {^ту}М=-м - коэффициенты линейного алгоритма оценки 1-й производной сигнала;

М - порядок алгоритма оценки производных.

На рис. 1 показан пример интерполяционного фрагмента.

Фрагмент у п () интерполирующей функции определяется значениями М предшествующих и М последующих относительно интервала г е [^; отсчетов sn_м, ..., эп_ь s„,

Sn+2, ..., \+И+1 (М > 0 - порядок интерполяции) и представляет собой обобщенный многочлен

N-1

у„ (t) = £ Cn,kwk (t - tn).

(5)

к=0

Здесь [wk(t^ -стема функций (первичный базис);

линейно независимая си-

" sn-3

*п-3

с sn-2

<> S,

п-1

«Л+1

sn+2

" sn+3

2 *п-\ ^п 1 ги+2

Рис. 1. Пример интерполяционного фрагмента

те

X (U

ч

те

Q.

те

О О.

Ё ф

ц

(Ч

о см

■ч-

О!

<

I

о та

г

о ^

со та г о.

<и

3

и <и со

см ■ч-ю

с?

см ■ч-ю см

(П (П

СпЛ - коэффициенты, которые ввиду выражения (1) линейно связаны с отсчетами сигнала:

Б =

м+1

Сп, к = X ат,

к^п+т'

(6)

т=-М

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

N = 2г + 2.

(7)

N-1

X а(М+1),кЩ *(-(М+1))Х * е [*(-(М+1)); *(-М)],

к=0

N-1

Xа0,к™к(* - *о), * е [Ч;

(8)

к=0

N-1

А = W

(9)

в котором

А =

а

- М ,0

а

- М+1,0

а

М+1,0

а_

М ,1

а_

М+1,1

а

М+1,1

а- М, N -1 а- М+1, N-1 ■■■ аМ+1, N -1

W =

<0,0 40 . . <N-1,0

Ч0,1 <1,1 . . <N-1,1

4,0 4,0 . . <N-1,0

<1 Ч) . . Ч-ц

w0:o) 40? Ч0 . Ч,1 . . 4-1,0 . ч£-1Д

0 . . 0 1 0 0 . .. 0

0 . . 0 0 1 0 . .. 0

а(1) М-м . а« а0(1) а(1) а2(1) . .0

0. . а-2 а® а0О) а« . а (1) м

а (г) М-м . а-? 4Г) а(г) а(г) . . 0

0. . а-? ) а0г) а1г) . а (г) м

(10)

Рассматривая восстановление единичного отсчета, из формул (5), (6) для порождающей функции получим

Фо(*) =

где = П\™Т);

ат,к = [ А]к,т+М .

Таким образом, порождающая функция (8) может быть найдена из уравнения (9). Тогда в виду (1) для каждого интерполирующего фрагмента сможем записать

X а( - М ),кЩ *М ), * е [*М; *М+1], I к=0

где коэффициенты определяются из решения матричного уравнения

М+1

Vп(?) = Х Эп+мЪ^ - (п + т)Т). (11)

Как видно из выражения (11), вычислительные затраты на расчет значения интерполирующей функции в некоторый момент времени составляют 2М +1 умножений при условии, что значение порождающей функции, соответствующее этому моменту времени, хранится в памяти компьютера. Последнее обычно имеет место, поскольку моменты времени, в которых будет рассчитываться интерполирующая функция, известны заранее. При этом вычислительные затраты определяются только порядком интерполяции и не зависят от степени гладкости и выбора первичного базиса. Однако если указанное условие не выполнено, то на каждое из 2М +1 умножений добавляются N операций умножения для расчета (8) и умножения , необходимые для расчета функций первичного базиса. При этом в соответствии с уравнением (7) чем больше степень гладкости, тем больше вычислительные затраты, они составляют (2М + 1)((2г + 2) + пм,). Анализ сходимости при локальной регулярной сплайновой интерполяции сигналов

В статье [13] показано, что сходимость регулярного метода интерполяции имеет место, когда он обеспечивает восстановление единицы. Рассматривая восстановление единицы,

из системы уравнений (4) получим условия сходимости при локальной сплайновой интерполяции:

м

£ dm=о, I=1,..., г. (12)

Согласно теореме о дифференцировании сигнала для преобразования Фурье, дополнительно укажем, что последовательность коэффициентов алгоритма оценки производной нечетного порядка должна быть нечетно-симметричной, а в случае производной четного порядка - четно-симметричной:

ит — а(-т) ) Ы0 _

I _ 0, 1, 2,..., т _ -М, ..., М;

= ^й, I = 1 , 2 , . . . , т = -М,..., М. (13)

Анализ искажений при локальной регулярной интерполяции сигналов

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

Т(ю) = Ф0(ю) X s(nT)e-j^anT.

Поскольку спектр дискретного сигнала

(ш) = X <пТ)е-™т = - X ^(ш - шд),

где £(ш) - спектральная плотность сигнала 2п

юд = —— частота дискретизации, последнее выражение перепишем в виде

1 +1»

Т(ю) = 5д(ю)Ф0(ю) = Фо(ю)- X 5(ю-юд).

Спектральная плотность интерполирующей функции представляет собой результат выделения спектра исходного сигнала £ (ш) из спектра дискретного сигнала £д(ш) путем умножения порождающей функции Ф 0(ю) (рис. 2) на спектральную плотность. Точным такое выделение является в случае, когда Ф0(ю) представляет собой прямоугольную функцию, что при локальной интерполяции исключено, поскольку ввиду ограниченности во времени порождающей функции ф0(7) ее спектральная плотность Ф0(ю) является аналитической функцией и, не являясь тождественным нулем, не может обращаться в нуль ни на каком конечном интервале. В связи с этим при условии шд > 2шт интерполирующую функцию у(7) можно формально рассматривать как результат искажений исходного сигнала э^), при которых для каждой гармоники сигнала частоты ш независимо от других при умножении на Ф0(ш) происходит изменение ее амплитуды и наложение на сигнал колебаний частот кюд + ю, к = 0, 1, 2, ..., с амплитудами, пропорциональными Ф0(кшд + ш). При этом комбинационного взаимодействия между

о см

01

<

I

о та

0 ^

СО та

1

о.

ф

£

и

V

со

см ■ci-io

с?

см ■ci-io см

(П (П

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

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

Ки (W) = T х

(14)

х>/(Т - Ф0(ш))2 + Фо(шд + ш) + Ф2(2®д + + Ф2(3<Ид + ш) +...

При практических расчетах учитывается конечное число нежелательных гармоник

KN (ш) = - х

(15)

- Ф0(ш))2 + Ф2(шд + ш) + Ф2(2шд + ш) +... + Ф2( + ш)

Анализ зависимости Ки(ш) может быть положен в основу определения требуемого значения частоты дискретизации и/или выбора метода интерполяции при заданном уровне искажений. Выбор допустимого уровня искажений Ки,доп зависит от специфики решаемой задачи и определяет граничную частоту шгр для метода интерполяции, такую, что Ки(ш < Шгр ) < Ки,доп. Для восстановления сигнала с максимальной частотой спектра шт с заданным уровнем искажений требуется обеспечить шт < Шгр , при этом каждая гармоника сигнала частоты ш восстанавливается с Ки(ш) <

и,доп-

< к

Постановка задачи синтеза оптимального интерполяционного сплайн-базиса

При заданных допустимом уровне искажений

Ки,доп, первичном базисе {к (О}}-1, порядке интерполяции М и степени гладкости г требуется отыскать набор элементов матрицы Б (10), удовлетворяющий условиям сходимости (12), который обеспечивает максимальную граничную частоту Ш гр коэффициента искажений (14):

юГр ((}М-м ) max. (16)

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

Задача синтеза оптимального интерполяционного сплайн-базиса (16) - это многопараметрическая однокритериальная задача оптимизации. Замкнутый вид целевой функции задачи получить затруднительно, что обусловило отыскание приближенных (квазиоптимальных) решений численными методами, включая методы случайного поиска и генетический алгоритм [14, 15], и потребовало разработки программного обеспечения большого объема.

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

Отметим также, что дополнительное привлечение взаимосвязей (13) позволяет сократить количество параметров задачи. Некоторые квазиоптимальные локальные полиномиальные сплайн-базисы

При полиномиальной сплайновой интерполяции первичный базис образован степенными функциями wk(t) = tk, k = 0, ..., N -1, где N определено (7). Порождающие функции оптимальных интерполяционных сплайн-базисов при Ки,доп = 0,001 для интерполяции порядка M = 1, 2, 3 и степеней гладкости r = 1, 2 приведены в таблице. Там же даны значения граничных частот в нормированных единицах. Графики коэффициентов искажения показаны на рис. 3. Номера линий на рис. 3 соответствуют номерам, указанным в крайнем правом столбце таблицы.

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

r M Фо(0 (Ю^)гр №

1 1 fl - 2,4586т2 + 1,4586т3, т = 11 / T |, 111 < T, Pn(t) = 1 [-0,5414т + 1,0829т2 - 0,5414т3, т = 11 / T | -1, T < 111 < 2T 0,75 1

2 Фо(г) = 1 - 2,2343т2 +1,2343т3, т = 11 / T |, 111 < T, -0,7657т + 1,3826т2 - 0,6169т3, т = 11 / T | -1, T < 111 < 2T, 0,1489т - 0,2978т2 + 0,1489т3, т = 11 / T | -2, 2T < 111 < 3T 1,02 2

3 Фо(г) = 1 - 2,1173т2 +1,1173т3, т = 11 / T |, 111 < T, -0,8827т + 1,4979т2 - 0,6152т3, т = 11 / T | -1, T < 111 < 2T, 0,2675т- 0,4840т2 + 0,2165т3, т = 11 / T | -2, 2T < 111 < 3T, -0,0511т + 0,1021т2 - 0,0511т3, т = 11 / T | -3, 3T < 111 < 4T 1,02 3

2 1 Фо(г) = 1 - 1,0182т2 - 4,3122т3 + 7,2100т4 - 2,8796т5, т = 11 / T |, | t| < T, -0,5310т + 0,5091т2 + 1,6588т3 - 2,7208т4 + 1,0839т5, т = 11 / T| -1, T < 111 < 2T 0,70 4

2 Фо(г) = 1 - 1,1779т2 - 2,8917т3 + 5,0617т4 - 1,9922т5, т = 11 / T |, 111 < T, -0,7448т + 0,5957т2 + 2,1489т3 - 3,2373т4 + 1,2375т5, т = 11 / T | -1, T < 111 < 2T, 0,1315т - 0,0068т2 - 0,7685т3 +1,0314т4 - 0,3876т5, т = 11 / T | -2, 2T <111 < 3T 1,36 5

3 Фо(г) = 1 - 1,36т2 - 1,7925т3 + 3,5101т4 - 1,3576т5, т = 11 / T |, 111 < T, -0,8451т + 0,7473т2 +1,7872т3 - 2,6821т4 + 0,9927т5, т = 11 / T | -1, T < 111 < 2T, 0,2461т - 0,0570т2 -1,1009т3 + 1,4422т4 - 0,5304т5, т = 11 / T | -2, 2T <| 11< 3T, -0,0537т - 0,0103т2 + 0,3533т3 - 0,4607т4 + 0,1715т5, т = 11 / T | -3, 3T < 111 < 4T 1,75 6

Ки(а>)

0,5

А //у/

3 f

1 \ 2 / ^4

\

5

0,2

0,4 0,6 0,8 1,0 1,2 1,4 1,6 1,8 2,0 2,2 2,4 2,6

Рис. 3. Коэффициент искажений для методов интерполяции из таблицы:

1 - М = 1, г = 1; 2 - М = 2, г = 1; 3 - М = 3, г = 1; 4 - М = 1, г = 2; 5 - М = 2, г = 2; 6 - М = 3, г = 2

2,8

(оТ

те

х

(U

ч

те Q.

те

о

о.

£

V

ц

(Ч

о см

01

<

I

со та

0 ^

СО та

1

.

ф

£

и ф со

см ■ci-io

с?

см ■ci-io см

(П (П

Заключение

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

Численно-аналитическими методами получены квазиоптимальные решения поставленной задачи синтеза оптимального интерполяционного сплайн-базиса в виде порождающих функций фундаментальных интерполяционных полиномиальных сплайн-базисов при интерполяции порядка M = 1, 2, 3 и степенях гладкости r = 1, 2. Даны практические рекомендации по их использованию: допустимые искажения обеспечиваются, если максимальная частота в спектре сигнала не превышает граничной частоты коэффициента искажений.

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

За рамками статьи остался поиск оптимальных локальных полиномиальных сплайн-базисов при более высоких порядках локальной интерполяции и степенях гладкости. Список литературы

1. Гончаров В. Л. Теория интерполирования и приближения функций. М.: Гос. изд-во техн.-теор. лит., 1954. 330 с.

2. Ahlberg J. H, Nilson E. N, Walsh J. L. The theory of splines and their applications. New York: Academic Press, 1967. 292 p.

3. Завьялов Ю. С., Квасов Б. И., Мирошниченко В. Л. Методы сплайн-функций. М.: Наука, 1980. 355 с.

4. De Boor C. A practical guide to splines // Applied Mathematical Sciences. 1978. Vol. 27. 292 s.

5. Сухорученков Б. И., Меньшиков В. А. Методы анализа характеристик летательных аппаратов. М.: Машиностроение, 1995. 365 с.

6. Денисенко А. Н., Исаков В. Н. Применение различных методов восстановления непрерывных сигналов по их дискретным значениям // Радиотехника. 2001. № 10. С. 16-20.

7. Денисенко А. Н., Исаков В. Н. Сравнительный анализ методов восстановления непрерывных сигналов по их дискретным значениям, основанных на локальных сплайнах различной степени гладкости // Вопросы повышения эффективности радиоэлектронных систем: межвуз. сб. науч. трудов. М.: МИРЭА, 2001. С. 101-106.

8. Денисенко А. Н., Исаков В. Н. Методы сплайновой интерполяции сигналов // 52 на-уч.-техн. конф. МИРЭА. Сб. трудов. Ч. 2. Физико-математические науки. М.: МИРЭА, 2003. С. 79-84

9. Денисенко А. Н. Сигналы. Теоретическая радиотехника. М.: Горячая линия-Телеком, 2005. 704 с.

10. Crochiere R E., Rabiner L. R Interpolation and decimation of digital signals: a tutorial review // Proceedings of the IEEE. 1981. Vol. 69. No. 3. Pp. 300-331.

11. Schafer R. W., Rabiner L. R. A digital signal processing approach to interpolation // Proceed-ingsofthe IEEE. 1973. Vol. 61. No. 6. Pp. 692-702.

12. Unser M. Splines: a perfect fit for signal and imageprocessing // IEEE Signal Processing Magazine. 1999. Vol. 16. № 6. Pp. 22-38.

13. Исаков В. Н. Сходимость при регулярной интерполяции и локальные интерполяционные базисы // Наукоемкие технологии. 2013. № 4. С. 40-46.

14. Исаков В. Н. Алгоритмы случайного поиска в задачах численной оптимизации // Сб. науч. трудов II междунар. науч.-практ. конф. «Актуальные проблемы и перспективы развития радиотехнических и инфокоммуника-ционных систем» (РадиоИнфоКом-2015). М.: МГТУ, МИРЭА, 2015. С. 146-151.

15. Рутковская Д., Пилиньский М., Рутков-ский Л. Нейронные сети, генетические ал-

горитмы и нечеткие системы: пер. с польск. И. Д. Рудинского. М.: Горячая линия-Телеком, 2007. 452 с.

16. Исаков В. Н. Фундаментальные интерполяционные базисы и спектральный анализ

при локальной интерполяции обобщенными сплайнами // Вестник МГТУ МИРЭА. 2015. № 1 (6). С. 144-154.

Поступила 30.01.17

Исаков Владимир Николаевич - старший преподаватель кафедры теоретической радиотехники и радиофизики Федерального государственного бюджетного образовательного учреждения высшего профессионального образования «Московский технологический университет», г. Москва.

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

Optimum regular local spline interpolation of signals

The study deals with an optimum approach to regular local signal interpolation by means of generalised splines. For the special case of local regular polynomial spline interpolation we derive quasi-optimal interpolation bases and provide corresponding recommendations dealing with selecting interpolation order and order of smoothness. Keywords: signal interpolation, interpolating function, quantization interval, spline interpolation, spline basis.

Isakov Vladimir Nikolaevich - Assist. Professor, Department of Theoretical Radio Engineering and Radiophysics of Federal State Budget Educational Institution of Higher Education "Moscow Technological University", Moscow. Science research interests: theory of signals, digital processing of signals, signal interpolation, extrapolation and approximation, optimum signal processing.