CC BY
48
о
m
УДК 519.216 ББК 22.3
О ЕДИНСТВЕННОСТИ КВАЗИГАРМОНИЧЕСКОГО ПРЕДСТАВЛЕНИЯ 1
В.К. Игнатьев, А.В. Никитин, С.В. Юшанов
Показано, что если для сигнала существует квазигармоническое представление с ограниченной медленно меняющейся частотой, то оно единственно. Показано, что при условии медленности изменения параметров системы, представимой в виде дифференциального уравнения второго порядка с переменными коэффициентами, решение уравнения существует в виде квазигармонического представления с медленно меняющейся амплитудой и частотой. Предложен регуляризированный алгоритм оценивания параметров сигнала с квазигармоническим представлением и метод оценивания медленно меняющихся коэффициентов дифференциального уравнения второго порядка. Приведены результаты статистического моделирования.
Ключевые слова: квазигармоническое представление, медленно меняющиеся параметры, мгновенная частота, обратная задача, условная корректность, регуляризация Тихонова.
Введение
Во многих радиофизических и радиотехнических задачах широко используется представление сигнала в виде квазигармонического колебания
х(о = л(0сов[е(0], (!)
причем производная полной фазы е(0 называется мгновенной частотой ю(^) ю(/) = е(/) [19]. Для о однозначного описания сигнала х(0 с помощью двух функций - амплитуды (огибающей) А (0 и г полной фазы е(0 необходимо то или иное дополнительное условие. Например, если полная фаза ^ определяется как аргумент аналитического сигнала z(t) = А(0ехр[/е(0] = х(0 + /у(0, построенного о на основе исследуемого колебания х(0 с помощью преобразования Гильберта [23]:
5 у(0 ==н[x{t)]==i-)^ct==-^\x^^x^^-cн==л{t)sm[е{t)], (2)
03 n-)t -Т 2Я-) Т
<С и
В то мгновенная частота ю(0 сигнала (1) равна
х^) у^) - Х^) у^)
o(t ) = 0(t) =-
х 2(t) + f(t)
Применяются и другие определения полной фазы [22], которые, как и определение Га-бора [23], являются формально-математическими и могут приводить к физическим проти-© воречиям [3].
Отметим, что вопрос о единственности представления колебания в форме (1) существенен, если амплитуде и (или) мгновенной частоте ставятся в соответствие измеряемые физические величины, например, амплитуда определяет энергию процесса, а мгновенная частота в задаче доплеровской локации связана со скоростью объекта. В этом случае нужны физические условия однозначного определения амплитуды и фазы. При естественных предположениях, что амплитуда А(0 и мгновенная частота ю(0 - положительные и непрерывные функции, а ю(0 ограничена как сверху, так и снизу
0 < ш т < <а(?) < шм, (3)
одним из таких условий может быть ограниченность носителя спектра мгновенной частоты ш(0.
1. Единственность квазигармонического представления
Пусть функция х(0 имеет нули в моменты времени tk, образующие монотонную последовательность без конечных предельных точек. Из непрерывности и положительности функции А(0 следует, что
9$) = 9к = тк - т/2. (4)
Здесь принято, что в момент времени t = t0 = 0 функция х(0 проходит через ноль с положительной производной. Таким образом, нахождение полной фазы 9(0 сводится к решению интерполяционной задачи для неэквидистантной последовательности точек (узлов интерполяции) tk.
Задача интерполяции (4) всегда имеет решение, не обязательно в виде ряда Лагранжа, в классе целых аналитических функций, какими бы ни были действительные числа tk и 9k [5]. Это решение не является единственным. Пусть целые функции 9(г) и 9х(г), где z = t + т, являются решением интерполяционной задачи (4). Из условия (3) следует, что в точках z = tk функция cos[9(z)] имеет простые нули. Тогда функция г) = (91( г) -9 (г) )/ ео8[9( г)] тоже является целой аналитической, так как нули знаменателя совпадают с нулями числителя. Поэтому общее решение интерполяционной задачи (4) в классе целых функций имеет вид
ВД = 9(2) + х(г)соБ[9(г)], (5)
где х(г) - произвольная целая функция.
Покажем, что решение интерполяционной задачи (4) единственно в классе целых функций с финитным спектром. Пусть существует конечный предел
юс = lim
C T ^œ
1 t
— fœ(t)dt 2T
= lim
k^œ
2%k
h " t-k
(6)
Естественно, что шт < шС < шм. Из ограниченности мгновенной частоты следует, что полная фаза 9(0 растет при ^ да не быстрее, чем t. Если носителем спектра функции 9(0 является отрезок [-шВ, шВ], причем шВ < шс, то в силу теоремы Винера - Палея - Шварца [4; 21] ее аналитическое продолжение 9(г) является целой функцией конечной степени, меньшей, чем шс.
Целая функция конечной степени имеет первый порядок. Из теоремы о порядке произведения целых функций [12] следует, что для того, чтобы функция 9х(г), не равная тождественно 9(г), также имела конечную степень, меньшую шс, необходимо, чтобы и функция cos[9(z)] имела порядок не выше первого. Целая функция cos(z) имеет первый порядок, поэтому сложная функция cos[9(z)] может иметь первый порядок только в случае, если 9(г) - многочлен [24]. Но многочлен степени выше первой растет на действительной оси быстрее, чем что противоречит условию (3) ограниченности мгновенной частоты.
Таким образом, различные решения интерполяционной задачи (4) возможны лишь при 9(0 = ю^ - л/2. Но спектр функции вида (5)
0^0 = - л/2 + х(0соэ(Ю(/ - л/2) (7)
ограничен сверху частотой ю0 + ю^ > юс > юВ, где ю^ - верхняя граничная частота спектра функции Х(0, так как для сигнала (1) с полной фазой (7) из условия (6) следует юс = ю0, что и доказывает единственность решения интерполяционной задачи (4) в классе целых функций с финитным спектром. При этом финитности спектра амплитуды А(0 и самого сигнала х(0 не требуется.
Если предел (6) не существует, класс единственности решения интерполяционной задачи (4) определяется условием юВ < ют. Условие положительности амплитуды А(0 также можно ослабить. Достаточно того, чтобы функция А(0 имела нули в счетном множестве точек tm, которые полагаются известными. Тогда из множества узлов интерполяции tk следует исключить простые нули сигнала х(0, совпадающие с числами tm.
Поскольку производная целой функции имеет тот же порядок и степень, что и сама функция [12], мгновенная частота в силу теоремы Винера - Палея - Шварца [4; 21] также должна быть целой функцией с финитным спектром и может быть представлена интегралом Фурье
)= 1 ю(р )ехр(р Мр. Тогда для производных мгновенной частоты справедлива оценка
|ю (и)(0| < Рюпв, Р = Цю( р )|йр .
(8)
Сигнал, для которого существуют константы Р и юВ такие, что условия (3) и (8) выполняются в любой момент времени, естественно назвать сигналом медленно меняющейся частоты (ММЧ).
С другой стороны, при выполнении условия (8) коэффициенты разложения с аналитической
I I < юпвР
функции ю^)/юс в ряд Тейлора в произвольной точке удовлетворяют условию |с„| < ^ю . Тогда
для степени а этой функции получаем [12] а = 1;т ( ^п! сп ) < ю в и по теореме Винера -
уп! сп i < ю ,
п ^ да ^
Палея - Шварца [4; 21] мгновенная частота ю(0 имеет финитный спектр, ограниченный частотой юв < юс. Таким образом, если для сигнала существует представление (1) с медленно меняющейся частотой, удовлетворяющей условиям (3) и (8), то оно единственно. В общем случае это представление может не являться аналитическим сигналом.
Вопрос об условиях, которым должен удовлетворять сигнал, чтобы для него существовало квазигармоническое представление (1) с медленно меняющейся частотой, относится скорее не к радиофизике, а к функциональному анализу. Он сводится к задаче об условиях, которым должна удовлетворять последовательность tk нулей функции х(0, чтобы числа 9к = лк - л/ 2 могли быть значениями функции с финитным спектром, взятыми в моменты tk [7]. Поскольку последовательность 9к бесконечно большая, эти условия могут быть излишне ограничительными или не иметь ясного физического смысла [11]. Например, может быть усилено условие (6) [7]:
к =л, 5к < * , к = 0,1,... к юс к 2к+1
На практике существование решения в классе единственности можно обосновать, например, связью параметров сигнала с медленно меняющимися параметрами породившей его физической системы. К примеру, в ЯМР известны уравнения Блоха для поперечных проекций магнитного момента [9]:
-ю
в
Мх = уМВ - Мх/Г2, Му = -уИхВ -Му/Г2, Вх = 0, Ву = 0, Вг = B(t),
из которых легко получается уравнение для свободной прецессии поперечной намагниченности в медленно меняющемся продольном магнитном поле В(^:
х^) + р^ )х(0 + q(t)x (t) = 0, (9)
где х(0 = Мх или Му, p(t) = 2/ Т2 - В/В, q(t) = у2 В2 - В/ (ВТ2) + 1/ Г/. Если возможно по измеренному, например, спиновым детектором [9] сигналу х{() восстановить коэффициенты р(^) и q(t) уравнения (9), то можно проследить за динамикой изменения продольного магнитного поля В(0 и времени поперечной релаксации Т2(0, что важно в задачах спектроскопии, томографии, исследовании геомагнитного поля и т. д. Для этого необходима аналитическая связь параметров квазигармонического представления (1) с коэффициентами уравнения (9).
2. Асимптотическое решение
В общем случае будем рассматривать колебательную систему с одной степенью свободы, у которой некоторые параметры медленно (по отношению к периоду собственных колебаний) изменяются. Тогда мы приходим к рассмотрению дифференциального уравнения второго порядка с медленно меняющимися коэффициентами (9), которое подстановкой
x(t ) = y(t)exp
1 t
- - J P (t')dt '
2 о
(10)
приводится к виду
y(t ) + ^ 2(t) y (t ) = 0, (11)
где Q2(t) = q(t) __ Ip2(t). (12)
2 dt 4 2
В случае медленного изменения параметров осциллятора функцию Q2(t) удобно представить в виде
Q 2(t) = Q If (t / т),
где Q0 - параметр размерностью частоты;
f - безразмерная функция величиной порядка единицы; т - характерное время ее изменения.
Условие медленности изменения параметров осциллятора можно выразить соотношением 1/Q0t = ц, ц << 1, или записать следующим образом:
Q(t)~ цQ2(t), ц <<1. (13)
Если уравнение (9) описывает осциллятор с малыми потерями, то с учетом физических допущений p2(t) ~ цд(0. При таких условиях с математической точки зрения уравнение (11) с медленно меняющимися коэффициентами эквивалентно уравнению с малым параметром ц при старшей производной. Такие уравнения решают с помощью асимптотических методов разложения по малому параметру [16; 19; 20]. Для построения приближенного решения удобно перейти к безразмерной переменной = t/т [19].
Асимптотическое решение уравнения (11) на некотором интервале I = ^2] при выполнении условия (13), а также f е C" (I) можно представить рядом [20]:
у © = ехр| — " Е а* к.
>=о
Удобно представить этот ряд в виде
у © = ехр{}|,—к ак (ОС}.
Первые члены этой асимптотики имеют вид:
а
_ /'(£)
,£) = ±/), а о(^) = -4/(|) =
аД) =
1 / "ф + . 5 (/ '(О)2
/ф 32 /2ф
а для последующих членов верны рекуррентные соотношения:
а.Д) =
2/)
ак © + Еа, ©а * -. ф 1.
]=0
(14)
Будем считать, что /£) > 0, и это соответствует колебаниям осциллятора, а не апериодическому затуханию [19]. Таким образом, имеем два формально асимптотических решения
у^, —) =^17^X3 К - + * {а *
/ © I — к=
(15)
Использованная процедура приводит к асимптотическому ряду, поэтому ограничиваются конечным отрезком ряда, например, при использовании приближения ВКБ [20] двумя его членами:
у.2 & —) = — { Щ"Ш"} + О—).
(16)
же].
Оценка остаточного члена равномерна по на I, то есть |О(—)| < С—, где С не зависит от [там При заданных начальных условиях можно представить решение в тригонометрической форме:
у (0 =
cos
jQ(t')dt ' + ф0
л,
4Щ
)], е(0 = | Q(t " )Л"
+ ф0
Здесь Л0 и ф0 - константы, определяемые из начальных условий. Возвращаясь к х(^, получим:
х^) = Л— ехр
',/ОД
1 t
- ^ { р (t'
cos
{О (Г )dt " + ф0
= Л(t)cos[е(t)],
где Л(Г) = 7Ще>Ф
- ^ { р с №'
е(t) = ' "+ф 0
(17)
(18)
0
Видно, что введенная ранее как производная полной фазы 9(0 частота ю(1) [19] совпадает с 0(1). Кроме того, как следует из условия (13), рассмотренное приближение неприменимо вблизи точек поворота 0(1) = 0 [10] - это согласуется с условием, что ю(1) положительная и ограниченная снизу функция.
Для применимости приближенного решения (16) необходимо выполнение условия (13) [19], однако само по себе оно может оказаться недостаточным, поскольку получено путем оценки по порядку величины различных членов в уравнении (11). В действительности надо требовать малости последовательных членов разложения в решении этого уравнения, а она может не обеспечиваться малостью отбрасываемого члена в уравнении (15). Так, если в экспоненте решения содержится член, возрастающий со временем по близкому к линейному закону, то малость первой производной в уравнении (15) не мешает тому, что на достаточно больших временах этот член может «набрать» большую величину; ВКБ-приближение оказывается тогда неприменимым для анализа сигнала на больших временных интервалах [10].
Рассмотрим для примера еще один член асимптотического разложения цаД^). Подставим выражение для этого члена из (14) в(15). Интегрируя, получим решение в виде
j± i - imШ
f 1MG) I I
32J f5'2(e)
de
(19)
Поправка, вносимая третьим членом разложения, является мнимой величиной и входит в экспоненту решения, что эквивалентно появлению вещественной поправки в фазе выражения (17). Члены в квадратной скобке в (19) должны быть малы по сравнению с 1 [там же].
При рассмотрении члена асимптотического ряда с номером к возникает слагаемое вида I(к VI32, которое также должно быть порядка единицы. Для него условие (13) будет иметь вид
0(к) (0| ~ Ца 0 а(к-1) (о|, ц << 1. (20)
В выражение (18) для амплитуды входит множитель 1Д^0(1), а из условия медленности изменения 0(0 и вдали от точек поворота можно говорить, что амплитуда тоже является медленно меняющейся функцией времени. Следует также заметить, что условия (8) и (20) идентичны.
3. Метод решения обратной задачи
Мы показали, что при условии медленности изменения параметров системы, представи-мой в виде дифференциального уравнения второго порядка с переменными коэффициентами, можно говорить, что решение системы существует в виде квазигармонического представления (1) с медленно меняющейся амплитудой и частотой. Мы будем решать обратную задачу -вычисление коэффициентов в уравнении вида (9) при известном х(0. Такая задача относится к типу обратных коэффициентных задач и в общем случае является некорректной [17], но возможно построение корректного решения на некоторых множествах [6; 18]. Медленно меняющиеся функции входят в класс функций с ограниченной вариацией, на множестве которых можно построить корректное решение. По определению, вариация непрерывно дифференцируемой функции на интервале равна интегралу нормы ее производной на этом интервале [14]. Для медленно меняющегося коэффициента 0(1) вариация за время его характерного изменения т ~ 1/ ю(0 при выполнении условия (20) равна
V(Q) =J||Ö(0||dt ~ i4|Q2(t)|d
I,
то есть вариация медленно меняющегося коэффициента на интервале т имеет порядок О(ц). Нетрудно показать, что таким же соотношениям удовлетворяют и исходные коэффициенты р(1)
и q(t), то есть они являются функциями с ограниченной вариацией. Но для применения регуляри-зирующих алгоритмов необходимо получить соотношения для самих коэффициентов. Воспользуемся представлением сигнала в виде (17)
x(t) = a(t) sin {e(t )}= b(t )exp{s(t)}sin{e(t)},
и обозначим e(t) = ®(t), S (t) = у (t). Условия медленности изменения параметров:
e(t) = ¿>(t) ~ цю2 (t), s(t) = y(t) ~ |ao(t), s(t) = y(t) ~ ц2ю2 (t),
b(t) ~ цш(0, 0 <ц<< 1. (21)
Из соотношений (12) и (18) следует связь между p(t), q(t) и параметрами сигнала ю(0, y(t):
q(t) = 1 p 2(t) + 2 p(t) + ю2(t), p(t) = -2y(t). (22)
Таким образом, для нахождения p(t) и q(t) достаточно получить выражения для ю(0 и y(t). Построим по аналогии с выражениями для коэффициентов авторегрессии метода Прони [8; 13; 15] второго порядка функции следующего вида:
) = x 2(t - 2А) - x (t - 4A)x(t) d (t) = X (t -A) - x(t - 2A )x(t) x 2(t - 2A) - x(t - 3A)x (t -А У x2(t - 2A) - x(t - 3A) x(t -A),
x 2(t - 2A) - x(t - 3A) x(t -A) * 0, (23)
где A - временной интервал, такой, что A << т или ю(^ << 1. Сначала рассмотрим знаменатель этих соотношений. Представим параметры сигнала разложением в ряд Тейлора около центральной точки t - 2A
х(1 - А) = [а(1 - 2А) + а(г - 2А)>т{9(! - 2А) + ю(! - 2А)А + ю(! - 2А)А2 /2}+ 5х,
х(1 - 2А) = а (I - 2 А)8т{0(! - 2 А)}, х(г - 3А) = [а(1 - 2А) - а(I - 2А)^1п{9(1 - 2А) - ю(! - 2А)А + ю(! - 2А)А2/2}+ 5х.
Здесь 5х отвечает за величину отброшенных членов разложения (5х = О(ц)). Числители также разложим вокруг центральной точки - для е(1) это I - 2А:
х(1) = [а(1 - 2А) + 2а(1 - 2А)]зт{9(1 - 2А) + 2ю(! - 2 А) А + 2ю(1 - 2А )А2}+ 5х,
x(t - 2A) = a(t - 2A)sin{e(t - 2A)}, x(t - 3A) = [a(t - 2A) - 2a(t - 2A)]sin{e(t - 2A) - 2ra(t - 2A)A + 26(t - 2A)A2}+ 5 Для d(t) центральной точкой будет t - A:
x(t) = [a(t - A) + a(t - A)]sin{e(t - A) + o(t - A)A + Ю(t - A)A2/2}+ 5x, x(t - A) = a(t - A) sin{e(t - A)},
хЦ - 2А) = [а^ -А) - аЦ - А)^т{е (} - А) - ю^ - А )А + со ^ - А) А2 / 2}+ 5 х.
Подставим найденные для отсчетов сигнала разложения в (23) и, применяя известные тригонометрические формулы, ограничиваясь порядком малости О(ц), получим
e(t) « 4cos2 {co(t - 2А)А} + 5e, 5d = О»,
d(t) « exp^ 2y(t - - А)А \ + Sd, 5e = О(ц).
(24)
Из (24) легко получить соотношения для параметров сигнала, считая остаточные члены 5а и 5 малыми величинами:
, ч 1 иe(t + 2А) I 1 , ,
c(t ) = а arccosj^^-J-j, y(t) = — ln{d (t + 2А)}.
Из (22) получаем выражения для коэффициентов уравнения (9):
(25)
q(t) = -4 p 2(t) + 2 p (t ) +
1 U e(t + 2А) ~
—arccos s --
А I 2
, p(t) = - lln{d(t + 2А)},
(26)
которые, как показывают численные эксперименты, неустойчивы к шумам. Для удовлетворения условиям корректности [6; 17; 18] мы должны рассматривать соотношения на некотором интервале LА ^ > 4), на котором выполняются условия (21). Из соотношений (26) и определения вариации функции как интеграла нормы ее производной, можно утверждать, что в случае ограниченности вариации коэффициентов р(0 и q(t) на интервале LА, промежуточные функции и с1^) также являются функциями с ограниченной вариацией на интервале LА. Запишем выражения (23) на интервале LА в виде
e(t - ХА+ А) =
х2(t - LА - А) - x(t - LА - 3А)x(t - LА + А) x2(t - LА - А) - x(t - LА- 2А)х^ - LА) '
e(t) =
x 2(t - 2 А) - x(t - 4 А^(0 x2 (t - 2А) - x(t - 3А) x(t - А) '
d(t - LА +А ) =
x 2(t - LА) - x (t - LА- А) x(t - LА +А) x2 (t - LА - А) - x(t - LА - 2А)x(t - LА) '
d(t) =
x2(t - А) - x(t - 2А) x(t ) x2 (t - 2А) - x(t - 3А) x(t - А):
(27)
и будем регуляризировать эти выражения на множестве функций с ограниченной вариацией, а затем по соотношениям (26) переходить уже к искомым коэффициентам уравнения (9). Каждую из систем (27) можно рассматривать как операторное уравнение 1-го рода вида Av = u, причем вместо точных значений u = Av и A нам известны такие их приближения ug и непрерывный оператор Ar, что ||us - Щ <5, ||Av - Av|| < y(h,| |v||). Для решения уравнения при известных Ar и ug необязательно знать величину изменения v. Можно ограничиться более широким классом изменения для первого приближения, а затем решать задачу методом градиента, минимизируя функци-
2
онал невязки Ф(^) =||Д^ - н5|| на множестве функций, где V не превосходит заданного приращения на рассматриваемом интервале [18].
Недостатком решения поставленной задачи на множестве функций ограниченной вариации, связанным с тем, что ограниченность вариации не предполагает медленного изменения функции, является ступенчатое изменение значений, что дает погрешности определения производной в выражении (26). Априорная информация о медленном изменении функций е(1) и d(t) позволяет провести дополнительную регуляризацию путем сглаживания окном длиной М, по аналогии с методом усреднения для медленно меняющихся функций [2], и после этого вычислять коэффициенты р(1) и q(t). Задавшись условием минимума невязки, получаемого после сглаживания решения, можно построить адаптивный алгоритм подбора оптимального значения М.
Функции е(1) и d(t) имеют разный характерный интервал изменения, поэтому в общем случае стоит рассматривать их на интервалах разной длины. Обозначим интервалы для е(1) и d(t) как LeА и LЙА, тогда выражения (27) следует записать в виде:
е(! - LeА + А) =
х2(! - Le А - А) - х(1 - Le А - 3А) х(1 - Ц А + А) х2 (г - Ье А - А) - х(1 - Ь А - 2 А)х(1 - Ье А) 5
е(0 =
х2(! - 2 А) - х(1 - 4А)х(t) х2 (t - 2А) - х(1 - 3А )х(1 - А)'
d(t - ЬЛА + А) =
х2 (t - А) - х(1 - А - А) х(1 - А + А) х2(! - А - А) - х(1 - Ьа А - 2А)х (! - А)'
d(t) =
х 2(! -А) - х(1 - 2 А) х(1 ) х2 (t - 2А) - х(1 - 3А) х(1 - А)'
(28)
Соответственно, и сглаживать полученные решения также следует окнами различной длины. Таким образом, коэффициенты уравнения (9) р(!) и q(t) можно оценить по соотношениям (26), решая системы (28) методом регуляризации и сглаживанием полученных зависимостей некоторым окном. Покажем это с помощью численного моделирования на примере уравнения Матье.
4. Численное моделирование
Уравнение Матье является частным случаем уравнения (11), в котором коэффициент О2(0 изменяется по периодическому закону. В размерных переменных оно имеет вид
у + а02(г - 2h008(20„О)у = 0, (29)
то есть в уравнении (11) обозначено
а 2(0 = а 2 (? - 2Ь соэ(2а 0 о).
Для уравнения Матье (29) известны периодические решения, которые носят название функций Матье. Они представляются асимптотическими рядами, которые сходятся тем быстрее, чем меньше коэффициент h в (29). Коэффициент 2 зависит от h и тоже рассчитывается с помощью асимптотических формул [1]. Асимптотическое решение уравнения (29) может быть найдено также с помощью ВКБ-приближения в квазигармонической форме (17).
Численное моделирование проводилось с помощью программы, написанной на языке FORTRAN, подпрограмм решения некорректных задач, описанных в работе [18] (PTIGR1, PTIGR2), а также библиотеки IMSL, которая содержит функции для вычисления асимптотических решений уравнения Матье и коэффициентов h и z. По заданному h находилось соответствующее собственное значение z, после чего для заданного Q0 вычислялось асимптотическое решение y(t) необходимого порядка. Для перехода от Q2(t) к соотношениям для p(t) и q(t) предполагалось, что p(t) меняется следующим образом:
p(t) = p0 + p, cos(Q,t ), (30)
а q(t), как следует из (12):
q(t) = Q0г - 2Q2hcos(2Qot) + 4(Po + Pi cos^t))2 - 4дЦ sm(Q,t).
В этом выражении второе слагаемое не является медленно меняющимся в смысле условия (20), а его введение обусловлено структурой уравнения Матье (29) для того, чтобы получить аналитическое решение уравнения (9) с переменным коэффициентомр(1). На динамику реальной инерционной системы это быстро осциллирующее слагаемое не влияет. Выделим как имеющую физический смысл медленно меняющуюся часть q(t)
q(t) = а22 +1 (р0 + р соз^О)2 - 2АО зт^) = q(t) + 200Лсоз^Щ). (31)
При этом должны выполняться условия на медленность изменения параметров:
~(0~ ц~2(/), р(!)~ цр 2(0, 0 <ц<< 1, а исследуемый сигнал из соотношения (10) должен иметь вид:
x(t ) = y (t)expJ- 2
p ot +-^sin(Qit)
К исследуемому сигналу добавлялся аддитивный белый шум с нулевым средним и среднеквадратичным отклонением (СКО) с^ и решались системы уравнений (27) методом условного градиента. Численное моделирование производилось с дискретными аналогами всех полученных выражений с шагом дискретизации А, а производная в выражении (26) вычислялась по разностной схеме. В качестве сглаживающего применялось окно Хэмминга [13]. Определялась систематическая погрешность как среднее значение относительного отклонения требуемой величины по 100 реализациям, а затем полученные значения усреднялись по интервалу наблюдения. Также рассчитывалась случайная погрешность как значения СКО на реализациях, которые потом усреднялись по интервалу наблюдения. На рисунках 1 и 2 приведены зависимости систематических погрешностей (Ар и А^) и случайных погрешностей (ср и с^) оценок Р(1) и q(t) соответственно, от значения СКО с^ при А = 0,21; Ь = 653; h = 0,1; 00 = 1,3; = 0,027; р0 = 0,035; рх = 0,01. Значения функций получены сглаживанием регуляризованного решения окном длиной М = 137.
Видно, что как систематическая, так и случайная погрешности оценки коэффициента р(0 увеличиваются с ростом с^ значительно быстрее, чем для ~ (/). Это связано с тем, что порядок величин этих коэффициентов разный:р2(0 ~ ^(0. Уменьшить погрешности возможно подбором интервала ЬА, а также при наличии какой-либо дополнительной априорной информации (например, об уровне шума, величине или виде самих коэффициентов), выбором
алгоритма регуляризации [18]. В погрешность оценивания р(0 значительный вклад вносят погрешности вычисления из-за малости самого коэффициента. В целом можно сказать, что рассмотренный метод дает состоятельную оценку коэффициентов, но оценка ) смещена, хотя и очень незначительно. Видимо, это смещение также обусловлено погрешностью вычислений и отброшенным остаточным членом в формуле (16). Смещение и случайная погрешность оценки медленно меняющейся частоты ю(0 из выражения (25) практически совпадает с графиком на рисунке 2.
Рис. 1. Систематическая Ар (1) и случайная стр (2) погрешности оценивания коэффициентар(0
в зависимости от уровня шума ст,
Рис. 2. Систематическая Аq (1) и случайная ст^ (2) погрешности оценивания коэффициента q(t)
в зависимости от уровня шума ст
Следует отметить, что условие медленности изменения частоты не предполагает ограниченности спектра сигнала (17), поэтому при моделировании использовали аддитивный белый шум, который имеет неограниченный спектр. Если же известна априорная информация о полосе частот исходного сигнала, на практике будет использован фильтр, ограничивающий спектр шума. В этом случае при моделировании следует использовать шум с ограниченным спектром, что уменьшит погрешности.
Заключение
Помимо задачи определения коэффициентов р(0 и q(t) рассмотренный алгоритм возможно применять для оценки и других параметров. К примеру, в общем случае при оцифровке сигнала шаг дискретизации А не является фиксированным. Если А задается, например, кварцевым резонатором, то при изменении температуры будет изменяться и частота резонатора, а следовательно, и шаг дискретизации, то есть, в общем случае, будет зависимость от времени А(0. При дискретизации гармонического сигнала с известной фиксированной частотой ю, но не обязательно постоянной амплитуды а(0, например сигнал квантового эталона частоты, можно получить соотношения для оценки А(0 при условии, что А(0 является медленно меняющейся функцией. Диск-ретизированный сигнал представим в виде
х = а + ф)
п п V п т /
и обозначим
п п
tn =ЕА,, А, =А0 +а,, tn = пА0 + = пА0 +5п.
1=1 1=1
Здесь А7 - значение шага дискретизации в 7-й момент времени; А0 - заданный шаг дискретизации; а7 - отклонение шага дискретизации А7 от А0 (изменение этой величины принимаем за малый параметр ц). Составим соотношение вида (23) для отсчетов сигнала в момент времени tn
х , + х, Л
с =-п1-п+1, х Ф 0,
тогда, повторяя процедуру по разложению отсчетов сигнала в ряд Тейлора вокруг центральной точки (в данном случае хп) и пренебрегая слагаемыми больше первого порядка малости по ц, получим соотношение вида (25)
~ 1 Г с 1 А = — агссо8| — I
ю I 2 У
п
Рассматривая это соотношение в другие моменты времени и решая систему, аналогичную (27), методом регуляризации, можно оценить медленное изменение шага дискретизации. Этот алгоритм позволяет определять нестабильность частоты кварцевых генераторов или сличать вторичные или рабочие эталоны частоты с первичным.
Также возможно предложенным методом решить задачу о разделении двух близких частот в спектре. В этом случае модель сигнала запишем в виде
х($) = Л^ )с08(ю/ + ф1) + АДОсОЗ^ t + ф2).
Введем следующие обозначения:
ю1 +ю2 ю1 -ю2 ,,
ю0 = 1 2 2 , М«^
X
П
0(0 =
А1 81П(0/ + ф1) - А2 81П(0/ - ф2) А1 008(0/ + ф1) + А2 008(0/ - ф2)
А(!) = ^А12 + А22 + 2А1А2 008(20/ + ф1 - ф2).
Тогда исходный сигнал преобразуется к виду (1)
х(/) = а(/)008( ю0/ + э(/)).
Так как А(/) и Э(/) зависят от разностной частоты, то они являются медленно меняющимися функциями, а следовательно, применяя рассмотренный алгоритм, можно оценить исходные частоты и ю2, даже если они спектрально не разрешены при ширине спектра огибающих АДО и А2(/), большей чем О.
Подобных примеров можно привести много, но несомненно, что априорная информация о медленности изменения параметров системы позволяет получать достаточно простые алгоритмы и решения задачи их оценивания.
ПРИМЕЧАНИЕ
1 Работа частично поддержана грантом РФФИ №2 10-07-9713 р_а.
СПИСОК ЛИТЕРАТУРЫ
1. Абрамовиц, М. Справочник по специальным функциям с формулами, графиками и математическими таблицами / М. Абрамовиц, И. Стиган. - М. : Наука, 1979. - 832 с.
2. Боголюбов, Н. Н. Асимптотические методы в теории нелинейных колебаний / Н. Н. Боголюбов, Ю. А. Митропольский. - М. : Наука, 1974. - 408 с.
3. Вайнштейн, Л. А. Разделение частот в теории колебаний и волн / Л. А. Вайнштейн, Д. Е. Вакман. -М. : Наука, 1983. - 288 с.
4. Гельфанд, И. М. Обобщенные функции и действия над ними / И. М. Гельфанд, Г. Е. Шилов. - М. : Физматгиз, 1958. - 276 с.
5. Гельфонд, А. О. Исчисление конечных разностей / А. О. Гельфонд. - М. : Физматгиз, 1959. - 212 с.
6. Денисов, А. М. Введение втеорию обратных задач / А. М. Денисов. - М. : Изд-во МГУ, 1994. - 208 с.
7. Евграфов, М. А. Асимптотические оценки ицелые функции / М. А. Евграфов. - М. : Наука, 1979. - 320 с.
8. Игнатьев, В. К. / В. К. Игнатьев, А. В. Никитин, С. В. Юшанов // Изв. вузов. Электромеханика. -2009.- №> 2. - С. 28-32.
9. Квантовая радиофизика / под ред. В. И. Чижика. - СПб. : Изд-во СПбГУ, 2004. - 689 с.
10. Ландау, Л. Д. Квантовая механика. Нерелятивистская теория / Л. Д. Ландау, Е. М. Лифшиц. - М. : Физматлит, 2004. - 800 с.
11. Левин, Б. Я. Распределение корней целых функций / Б. Я. Левин. - М. : Гостехтеориздат, 1956. - 632 с.
12. Леонтьев, А. Ф. Целые функции. Ряды экспонент / А. Ф. Леонтьев. - М. : Наука, 1983. - 176 с.
13. Марпл-мл., С. Л. Цифровой спектральный анализ и его приложения : пер. с англ. / С. Л. Марпл-мл. -М. : Мир, 1990. - 584 с.
14. Натансон, И. П. Теория функций вещественной переменной / И. П. Натансон. - М. : Наука, 1974. - 480 с.
15. Никитин, А. В. / А. В. Никитин, С. В. Юшанов // Физика волновых процессов и радиотехнические системы. - 2006. - Т. 9, №° 2. - С. 57- 63.
16. Островский, Л. А. Введение в теорию модулированных волн / Л. А. Островский, А. И. Потапов. -М. : Физматлит, 2003. - 400 с.
17. Тихонов, А. Н. Методы решения некорректных задач / А. Н. Тихонов, В. Я. Арсенин. - М. : Наука, 1979. - 286 с.
18. Тихонов, А. Н. Численные методы решения некорректных задач / А. Н. Тихонов, А. В. Гончарский, В. В. Степанов, А. Г. Ягола. - М. : Наука, 1990. - 232 с.
19. Трубецков, Д. И. Линейные колебания и волны : учеб. пособие / Д. И. Трубецков, А. Г. Рожнев. -М. : Изд-во физ.-мат. лит., 2001. - 416 с.
20. Федорюк, М. В. Асимптотические методы для линейных обыкновенных дифференциальных уравнений / М. В. Федорюк. - М. : Наука, 1983. - 352 с.
21. Хургин, Я. И. Финитные функции в физике и технике / Я. И. Хургин, В. П. Яковлев. - М. : Наука, 1971.- 408 с.
22. Cohen, L. // Proc. IEEE. - 1989. - V 46, №№ 7. - P. 941-981.
23. Gabor, D. // JIEE. - 1946. - V 93. - P. 429-457.
24. Polya, G. // Journal London Math. Soc. - 1926. - V. 1. - P. 12-15.
ABOUT THE UNIQUENESS OF QUASIHARMONIOUS REPRESENTATION
VK. Ignatjev, A. V Nikitin, S. V Yushanov
It is shown that if for a signal exists quasiharmonious representation with the limited slowly varying frequency, it is unique. It is shown that under condition of a slowness of change of parameters of system, representable in the form of the differential equation of the second order with variable factors, the equation decision exists in the form of quasiharmonious representation with slowly varying amplitude and frequency. It is proposed regularizing algorithm of estimation of parameters of a signal with quasiharmonious representation and a method of estimation of slowly varying factors of the differential equation of the second order. Results of statistical modeling are shown.
Key words: quasiharmonic representation, slowly varying parameters, instantaneous frequency, inverse problem, condition of well-posed problem, Tikhonov regularization.
