КОМПЬЮТЕРНОЕ МОДЕЛИРОВАНИЕ COMPUTER SIMULATION HISTORY
УДК 519.677
DOI: 10.18413/2411-3808-2018-45-2-298-311
ВЫЧИСЛЕНИЕ ЭЛЛИПТИЧЕСКИХ ФУНКЦИЙ ЯКОБИ ДЛЯ РАСЧЕТА ХАРАКТЕРИСТИК ФИЛЬТРА КАУЭРА
CALCULATING THE JACOBI ELLIPTIC FUNCTIONS TO CALCULATE THE CHARACTERISTICS OF CAUER FILTER
В.П. Волчков, А.В. Мирошниченко V.P. Volchkov, A.V. Miroshnichenko
Московский технический университет связи и информатики, Россия, 111024, г. Москва, ул. Авиамоторная, 8а
Moscow Technical University of Communications and Informatics, 8а Aviamotornaya St., 111024, Moscow, Russia
E-mail: volchkovvalery@mail.ru, Mirosh.A.V@yandex.ru
Аннотация
Рассматривается проблема расчета и аналитики эллиптических фильтров (фильтров Кауэра), заключающаяся в вычислении эллиптических функций Якоби, применяемых для аппроксимации передаточной функции фильтров Кауэра. Приводятся интегральные выражения для определения прямых и обратных эллиптических функций Якоби, зависящие от двух параметров u, k . Анализируются различные алгоритмы вычисления прямых, обратных по параметру u и обратных по параметру k эллиптических функций Якоби. Приводятся и сравниваются результаты, полученные различными алгоритмами. Даются рекомендации по применению разных методов вычисления эллиптических функций Якоби в задачах синтеза и анализа фильтров Кауэра.
Abstract
The problem of calculation and analysis of elliptic filter (Cauer filter), which consists in calculating the Jacobi elliptic functions, used for mathematical approximation Cauer filter transfer function. Provides integral for determining the direct and inverse Jacobi elliptic functions that depend on two parameters u, k. For direct Jacobi elliptic functions were developed five algorithms of calculation: the arithmetic-geometric mean method with Landen's transformation, approximation using trigonometric and hyperbolic functions, Taylor series and Fourier series expansion. For inverse by u Jacobi elliptic functions analyzed four algorithms of calculation: the arithmetic-geometric mean method, Cotes formulas, Simpson's method and the Monte Carlo method. Discusses the problem of calculation the inverse by parameter k Jacobi elliptic function. Offered possible solutions to this problem, among which a complete enumeration and attempt approximation through trigonometric and hyperbolic functions. Describes and compares the results obtained in different ways. Recommendations on the use of different methods of calculating the Jacobi elliptic functions in problems of synthesis and analysis of Cauer filters.
Ключевые слова: Фильтр Кауэра, эллиптические функции Якоби, метод арифметико-геометрического среднего, преобразование Ландена, ряд Тейлора, ряд Фурье, формулы Котекса, метод Симпсона, метод Монте-Карло.
Keywords: Cauer filter, Jacobi elliptic functions, the arithmetic-geometric mean method, Landen's transformation, Taylor series, Fourier series, Cotes formulas, Simpson's method, Monte Carlo method.
Введение
По результатам работ, связанных с аналитикой аналоговых фильтров [Мирошниченко, Волчков, 2017; Мирошниченко, 2017], был получен интересный вывод: наибольшую сложность анализа среди классических аналоговых фильтров представляет фильтр Кауэра. На первый взгляд кажется, что фильтр, открытый еще в 1930-х годах и активно применяющийся в радиотехнике, должен быть изучен вдоль и поперек, однако даже вычисление его амплитудно-частотной характеристики (АЧХ) вызывает огромные проблемы.
В целом передаточная функция фильтра Кауэра идентична передаточной функции фильтра Чебышева 1-го рода, за исключением функций со^(х) и агссо^(х), вместо них используются функции сё(х,к) и агссё(х,к), являющиеся эллиптическими функциями Якоби. В частности квадрат АЧХ фильтра Кауэра имеет вид:
Н 0)Г =-5—5-1-•
1 1 1 + есО (N ■ агссО(а, к), к)
Главная проблема работы с эллиптическими функциями заключается не только в их сложной аналитике, а в их эффективном вычислении. Даже в учебниках, рассчитанных на механико-математические факультеты таких вузов, как МГУ [Смирнов, 1958], данные функции даются крайне бегло и форменным образом не понятно, как их считать. В профильной литературе по аналоговым фильтрам, например, в [Лем, 1982], вопрос синтеза эллиптических фильтров заканчивается его упоминанием и приведением нескольких формул. На просторах интернета можно найти различные формулы, взятые «с потолка», без ссылок на какие-то источники, верность которых вызывает большие вопросы.
В данной работе мы рассмотрим и сравним различные способов вычисления эллиптических функций Якоби, которые можно применить для анализа фильтра Кауэра, например, для расчета его импульсной характеристики [Мирошниченко, 2017]. Для этого потребуется воспользоваться справочными монографиями [Абрамовиц и др., 1979; Березин, Жидков, 1962], где анализируются различные формы представления функций Якоби и их вычислительные аспекты. Различные приложения, связанные с синтезом цифровых алгоритмов обработки сигналов отражены в авторских публикациях [Волчков, Санников, 2017; Волчков, Безруков, 2017; Санников, Волчков, 2017; Тактакишвили, Волчков, 2017; Волчков, Асирян, 2017; Волчков, Тактакишвили, 2017; Sannikov, Volchkov, 2017; Безруков, Волчков, 2017; Волчков, Санников, 2016; Волчков, Шурахов, 2016; Волчков, Санников, 2016; Волчков, Шлома, Поборчая, 2015; Волчков, Шлома, 2014; Волчков, Уваров, 2015; Волчков, Поборчая, Шлома, 2013; Волчков, Шурахов, 2014; Волчков, 2009; Волчков, Петров, 2009; Волчков, Петров, 2009а; Черноморец, Волчков, 2012; Бурданова, Мама-тов, Волчков, 2011; Волчков, Петров, 2010; Волчков, Шурахов, 2012; Волчков, Шурахов, 2012а].
Определение эллиптических функций Якоби
Введем эллиптический интеграл вида
Р О©
и =
0
который интегрально связывает параметры и , р и к . Тогда можно определить следующие эллиптические функции [Абрамовиц и др., 1979]:
^п(и, к) = з1п(р), сп(и, к) = ео8(р), ёп(и, к) = ^ 1 — к2 з1п2(р) (2)
Через них остальные эллиптические функции можно выразить по формулам:
,, 1Ч сп(и, к) ёп(и, к) 1
са(и,к) =-, ёс(и, к) =-, т(и, к) =-, (3)
ёп(и, к) сп(и, к) ^п(и, к)
sd(и, к) = к) , пс(и, к) =-1-, О^м, к) = ап(и к),
ёп(и, к) сп(и, к) $п(и, к)
ь а ©, • (1)
— к2 81П2(©)
пё(и, к) =-1-, 8с{и, к) = к), С8(и, к) = сп(и к).
ёп(и, к) сп(и, к) 8п(и, к)
Период эллиптических функций 8п(и, к) и сп(и, к) отличен от привычного 2я , и
равен 4К (к) (у функции ёп(и, к) период равен 2К (к) ), где функция К (к) определяется
выражением [Абрамовиц и др., 1979]:
пП ё©
K (к) = f -
J0 у! 1 - к2 sin2 (0)
с параметром к е [0,1].
Следует заметить, что каждая эллиптическая функция имеет две обратные, но, к сожалению, даже в справочниках мирового уровня [Абрамовиц и др., 1979] они не рассмотрены.
Представим эллиптическую функцию в виде
У = f(u, к), (4)
где f - некая эллиптическая функция, например: sn, cn, dn и т.д.
Тогда определим обратную функцию по аргументу u:
u = arcf (y, к) . (5)
А также обратную функцию по параметру k:
к = агск[(u, у). (6)
Для большего понимания обозначений приведем примеры написания прямых и обратных эллиптических функций: у = cd(u, к), u = arccd(у, к), к = arc^d(u, у) .
Различные способы вычисления прямых эллиптических функций Якоби
В этой главе будет рассмотрена проблема вычисления прямых эллиптических функций Якоби, которые можно описать выражением (4).
Основной способ нахождения эллиптических функций - это метод арифметико-геометрического среднего (АГС) [Абрамовиц и др., 1979]:
1 шаг: задаем начальные условия для вспомогательных параметров an, bn, cn:
а = 1, b0 = V1 — к2, c0 = к — для sn(u, k), cn(u, k), dn(u, k) a = 1, b0 = cos^), c = sin^) - для K(k)
2 шаг: вычисляем an, bn,cn по рекуррентной формуле до заданной точности:
an = (an—1 + bn—1) 1 2 bn =4aZK~l, cn = (an—1 — bn—1) 1 2 .
3 шаг: находим % или вычисляем K(к) :
Т
K (к) =— aN , (N = 2NaNu
4 шаг: находим % по следующей рекуррентной формуле:
Г c 1 (n—1 = 0.5arcsin — sin((n) + 0.5%.
I an J
В результате данного шага мы получаем угол %, который будет с заданной точностью равен углу % из выражения (1).
5 шаг: находим искомые значения:
sn(u, к) = sin(%0), cn(u, к) = cos(%), dn(u, к) = cos(%) I cos(% — %).
Как можно видеть, метод арифметико-геометрического среднего хоть и позволяет получить результат, но из-за того, что он требует разложения в бесконечный ряд, для получения точных аналитических выражений он не пригоден.
Следующие 2 метода основаны на интересном свойстве эллиптических функций Якоби. При достаточно малом к « 0 эллиптические функции вырождаются в тригонометрические, а при к «1 эллиптические функции вырождаются в гиперболические.
На основе данных свойств можно допустить аппроксимацию эллиптических функций через тригонометрические или гиперболические [Абрамовиц и др., 1979]. При k — 0 :
sn(u, k) - sin(u) + 0.25k2 (u — sin(u) cos(u)) cos(u), cn(u, k) — cos(u) + 0.25k2 (u - sin(u) cos(u)) sin(u) , dn(u, k) — 1- 0.5k2 sin2(u) .
При k — 1:
sn(u, k) — th(u) — 0.25(1 — k 2)(sh(u)ch(u) — u) sech2(u), cn(u, k) — sech(u) — 0.25(1 — k2 )(sh(u)ch(u) — u)th(u) sech(u) , dn(u, k) — sech(u) + 0.25(1 — k2 )(sh(u)ch(u) + u)th(u) sech(u) . Приведенные формулы также являются громоздкими, однако позволяют привести нерекуррентную аналитику.
Главной проблемой данных способов являются указанные требования к параметру k . Однако есть возможность перерасчета эллиптических функций с заданным значением k, через эллиптические функции с меньшим или большим значениями k' = je [0,1] с помощью преобразования Ландена [Абрамовиц и др., 1979]. А именно, идея преобразования Ландена заключается в замене k ^ j и u ^ v, для уменьшения или увеличения параметра k .
Понижающее преобразование Ландена позволяет понизить параметр k в 5-20 раз. Проводим замену:
, 1—V1—k2 х
к ^ л =-. , u ^v =-,
1+ V1 — k2 1 +u
после этого расчет осуществляется по следующим выражениям:
sn(u,k) = (1 — ¡)sn(v,u), cn(u,k) = cn(v,¡)dn(v,u), dn(u,k) = dn2(v,u) —2(1 — u) . 1 + jusn (v, ju) 1 + jusn (v, ju) (1 + л) + dn (v, ju)
Повышающее преобразование Ландена позволяет повысить параметр k в 1.5-5 раз. Проводим замену:
k = 41 — ¡Л U = 7-7, u ^v = -—,
1 + k 1 + j
после чего расчет осуществляется по следующим выражениям:
snuk) = (1 + u>n(v,u)cn(v,u), cn(u,k) = (1 + j1)(dn2(v,u) — j) ,
dn(v, u) judn(v, ju)
dn(u, k) = (1 - u)(dn2(v,u) + j) judn(v, ju)
Четвертый способ расчета эллиптических функций основан на разложении их в ряд Тейлора [Абрамовиц и др., 1979]:
, 1Ч 1 + k2 3 1 + 14k2 + k4 5 1 + 135k2 + 135k4 + k6 7
sn(u,k) = u--u н--u--u + ...
3! 5! 7!
, , 1 2 1 + 4к2 4 1 + 44к2 + 16к4 6 сп(и,к) = 1--и л--и--и +...
2! 4! 6!
, , п , к2 2 4к2 + к4 4 16к2 + 44к4 + к6 6 ап(и,к) = 1--и л--и--и +...
2! 4! 6!
Однако рекуррентные формулы для вычисления коэффициентов полиномов не известны, а табличные коэффициенты в литературе [Абрамовиц и др., 1979] даны только для первых четырех членов.
Пятый способ вычисления заключается в разложении эллиптических функций в ряд Фурье [Абрамовиц и др., 1979].
sn(u, к) = 2ж Т -V ' кК (к) ¿01
q
2 n+1
cn(u, к ) =
dn(u,к)=
q
2ж ^ qn+05
sin(2nv + v),
т-
и=П 1
кК (к) П=01 + q
cos(2nv + v),
ж
2ж
и+0.5
2K (к) К (к) „
-жК (ф-к2)
т 1
n=1 1
+ q
-cos(2nv),
где q = e
К ( к )
V = -
жи
(к)
На рис. 1 представлены графики функции её(х, к) при трех значениях (к =0.1 -сплошная линия, к =0.5 - пунктирная линия, к =0.9 - штрихпунктирная линия), рассчитанные по всем 5 методами. Метод АГС был применен для 20 итераций, ряд Тейлора, в виду отсутствия рекуррентной формулы, представлен всего 4 членами, ряд Фурье представлен 10 членами (дальнейшее увеличение членов практически не влияет на результат, погрешность меньше Ю-4).
/7 ч\
if у\
if i J 1 1 v \ \ \ I I
i / \
/ / \ j
/ J \ /
1 J s
l 4, I / V / \ \ v \J: у
\ f ,\ /7 \\ '/ \ \ if 1Л ч 1 - \ 1J i \ ¡1 \ \Jr \/j
\ / \ / \ 1 \ ;
V V
А
B
C
у \ / \ f\j \ / X
A \ / Л ' Л
V v /'
\ \ fi Л / /
v. \ tj Л I /
V \ /' ■Л ¡7
V. \ /7 V\ / /
V \j! \v /
■ V/v' ■ЧХУ
D
E
Рис. 1. Графики функции cd(u,k), рассчитанные для параметра: k=0.1 - сплошная линия, k=0.5 - пунктирная линия, k=0.9 - штриховая линия. А - метод арифметико-геометрического
среднего, B - аппроксимация тригонометрическими функциями, C - аппроксимация гиперболическими функциями, D - разложение в ряд Тейлора, E - разложение в ряд Фурье Fig. 1. cd(u k) charts calculated for the parameter k=0.1 - solid line, k=0.5 - dotted line, k=0.9 - dashed line.
A - method of arithmetic-geometric mean, B - approximation of trigonometric functions, C - approximation of hyperbolic functions, D - Taylor series expansion, E - decomposition in Fourier series
Как можно видеть из рис. 1, метод АГС дает адекватные результаты при любом значении х. Методы аппроксимации тригонометрическими и гиперболическими функциями дает адекватный результат только на первом периоде х = ±2К (к). Разложение в ряд Тейлора дает адекватный результат только на полупериоде х = ±К(к) . Разложение в ряд Фурье хоть и дает адекватный результат изменения её(х, к), аналогичный методу АГС, но требует нормировки (нормировочный коэффициент равен 1/ её (0, к) и учтен при построении графиков на рис. 1).
Различные способы вычисления обратных по аргументу и эллиптических функций Якоби
В этой главе будет рассмотрена проблема вычисления обратных по аргументу и эллиптических функций Якоби, которые можно описать выражением (5).
Обратные эллиптические функции, как и тригонометрические, можно представить через интеграл. Согласно выражению (1) функция (5) представляет собой определенный интеграл с приделами 0 и р, и параметром к .
Заметим, что интеграл (1) имеет два входных аргумента р и к, а искомая функция (5) имеет аргументы у и к . Поэтому нам надо преобразовать у в р .
Данное преобразование получается из выражений (2):
у = ^п(и, к) = ъ1п(ф) ^ ф = агс8т(у), у = еп(и, к) = соъ(ф) ^ ф = агссоъ(у) (7)
У =
ёп(и, к) = ^ 1 - к2 ъ1п2(ф) ^ ф = агсъ1п
у2 -1
к
2
у = её (и, к) = . С0Ъ(Ф) ^ у 2(1 - к2 ъ1п2(ф)) = соъ2(ф) ^ у2 = соъ2(ф) + к2 ъ1п2(ф) у2 ^ - к2 ъ1п2(ф)
2 1(21 ^ у2 = 1 - ъ1п2(р) + к2ъ1п2(р)у2 ^—-—-т- = ъ1п2(р) ^ р = агсъ1п ' у
(-1 + к2у2) ^ " ]1(-1 + к2у2)
По аналогии можно вывести формулы для любой эллиптической функции.
Теперь, зная р, наша задача упирается во взятие не берущегося интеграла (1):
р
и =
о
- к2 ъ1п2(0) '
Этот интеграл носит название нормального эллиптического интеграла Лежандра 1-го рода. Данный интеграл имеет два возможных подхода к его решению.
Первый подход - применение метода АГС. Слегка изменив алгоритм АГС из предыдущего раздела, мы сможем его применить к вычислению данного интеграла:
1 шаг: задаем начальные условия для вспомогательных параметров аи, К, еи:
а = 1, К =>/1 - к2, е = к.
2 шаг: вычисляем ап, К, еи по рекуррентной формуле до заданной точности:
а = (а 1 + К 0/2, К =\1аА7, е = (а 1 - К 0/2.
п V п-1 п-1 ? п у п-1 п-1 ? п V п-1 п-1/
3 шаг: находим рп по следующей рекуррентной формуле (р0 находим через (7)):
ри_! = 0.5агсъ1п
и. ^
—п(рп) + О.р, =>2рп--рп = штат -^тОп)
V ап
V ап у
^81п(2р„-1 -рп) = —ът(рп) =>81п(2р„-1)с08(р„)-С08(2р„-1)81п(р„) ^ът^п)
ап ап
^ 81п(2р,) С08(рп) = (^ + соъ(2р,)) 81п(рп) ^ 51П(2рп-1) = ^^ ^
ап ^ + С0ъ(2р С08(рп)
ап
^р. = агсг (/"^ ). (8)
— + С0ъ(2рп-1) ап
Выражение (8) не совсем строгое. Так как в формуле стоит аг^, являющийся периодической функцией, то выражение (8) верно будет представить как:
рп = а^(с 51п(2рп-1) ) + пп,
п- + С0ъ(2рп-1)
где п е [0,1,...,да] некоторая константа. Определить п аналитически не представляется возможным, однако вычислить п можно из условия:
(е Л
рп_, = 0.5агсъ1п ъ1п(рп) + 0.5рп. (9)
V ап у
Требуется перебрать все п е [_рй_! /2],...,2\_р^ /2] + 2] и найти то единственное, для
которого выполняется условие (9).
Отметим, что при реализации данного алгоритма на ЭВМ следует учесть, что тригонометрические функции вычисляются до определенного знака после запятой, и по этому данное условие может не выполниться из-за расхождения на (10-2...10-4), вызванного приближенным вычислением тригонометрических функций. Поэтому следует указывать не точное равенство, а условие, при котором разность по модулю была бы меньше некого небольшого значения.
4 шаг: находим и :
ры
u =
2N а
2 aN
Второй подход - применение численных алгоритмов для вычисления определенных интегралов.
Мы рассмотрим три наиболее интересных алгоритма. Для упрощения формул введем обозначение:
f (x) = . 1 .
ф - к2 sin2 (x)
Тогда (1) можно представить как:
u = f , d 0 = = f f (0)d0 . 1^1 - к2 sin2 (0) о
Первый алгоритм - применение формул Котекса (средних прямоугольников) [Березин, Жидков, 1962]. Он описывается следующим выражением:
р-0 р-0
Г П и-1 (к +1)--Ьк--„-1 „0^,1
u = f f (0)d0 = lim р-0 f f (-n-) = lim р f f P —) .
о nn к^о 2 nt~0 n 2
При программной реализации n берут достаточно большим, чтобы погрешность вычисления была незначительной.
Второй алгоритм - метод Симпсона [Березин, Жидков, 1962]: р n-i n u = J f (0)d0 = lim P (f (0) + f (р) + 2f f (P 2к) + 4f f (P (2k -1))) . 0 6n к 2n кк 2n
где также n берут достаточно большим, чтобы погрешность вычисления была незначительной.
Третий алгоритм - метод Монте-Карло, включающий 5 шагов: 1 шаг - задаем начальные условия:
p = 0, q = 0.
an
2 шаг - случайным образом выбираем 2 числа:
х = Я(0,ф), у = Я(0,/(ф)), где Я(а, Ь) - равномерное распределение от а до Ь.
3 шаг - делаем перерасчет р и ц:
р = р + 1, если f (х) > у q = д + 1, если f (х) < у
4 шаг - повторяем 2 и 3 шаг достаточное большое число раз.
5 шаг - находим значение ^
и = ф/ ф)— . р+д
К сожалению, все 4 метода являются численными, что значительно затрудняет их аналитику. Метод Монте-Карло, в отличие от трех других, требует наличия генератора случайных чисел. Поэтому выбор метода заключается в выборе оптимального соотношения скорости схождения к трудоемкости вычисления.
Различные способы вычисления обратных по параметру к эллиптических функций
Якоби
В этой главе будет рассмотрена проблема вычисления обратных по параметру к эллиптических функций Якоби, которые можно описать выражением (6).
Параметр к связан с и и ф через (1)
ф
0
и =
- к2 б1П2 (0) '
Аналитически вывести из выражения (1) параметр к крайне затруднительно, поэтому для решения данной задачи применим методы из предыдущих разделов.
Метод АГС в теории позволит вычислить данный параметр наиболее точно, однако для этого надо подобрать коэффициенты а, Ь и с так, чтобы выполнилось условие:
фп_х = 0.5агс8т
^ е ^
и =
81и(ф„) + 0.5ф„,
,ап ) (10)
фм
2ыа 2 аы
Усложняет задачу то, что коэффициенты а, Ь и с связаны выражениями:
а = (а 1 + Ь 0/2, Ь =\1а А7, е = (а 1 -Ь 0/2 .
п V п—1 п—1 ? п у п—1 п—1? п V п—1 п—1/
При этом эти коэффициенты имеют начальные условия:
а0 = 1, Ь = >/ 1 — к2, е0 = к .
Видно, что если удастся выразить аналитически функции
ап = /1 (u, ф), Ьп = ЛОфХ еп = /3(и,ф) ,
тогда легко можно будет найти параметр к, однако нахождение данных функций крайне затруднительно.
Наиболее эффективным будет способ перебора параметра к во всем диапазоне от 0 до 1 и проверка его условием (10).
Главную сложность в улучшении сходимости полного перебора составляет зависимость периода эллиптических функций от параметра к. Все методы вычисления эллиптических функций, кроме метода АГС, позволяют адекватно вычислить функцию только на первом полупериоде. Не зная параметра к, мы не можем утверждать, что попадем в первый полупериод, поэтому невозможно гарантировать точное вычисление любым другим способом, кроме полного перебора с применением метода АГС.
Существует способ, позволяющий применить аналитические формулы - это аппроксимация тригонометрическими и гиперболическими функциями:
sn(u, k) = y w sin(u) + 0.25k2(u - sin(u) cos(u)) cos(u) ^ к = cn(u, к) = y w cos(u) + 0.25k2 (u - sin(u) cos(u)) sin(u) ^ к =
4( y - sin(u))
V
(u - sin(u) cos(u)) cos(u)
4( y - cos(u))
(u - sin(u) cos(u)) sin(u)
dn(u, k) = y w1 - 0.5k2 sin2 (u) ^ k =
2(1 - y)
sin (u)
cd (u, k) = cn(u,k) = y,
dn(u, k) ^ k = sn(u, k) = y ^ k =
cos(u) + 0.25k2 (u - sin(u) cos(u)) sin(u) 1 - 0.5k2 sin2 (u)
V
y - cos(u)
0.25(u - sin(u)cos(u))sin(u) + 0.5 y sin2(u) w th(u) - 0.25(1 -k22(sh(u)ch(u) -u) sech2(w) :
(11)
\
y - th (u) + 0.25(sh(u)ch (u) - u) sech2 (u)
0.25(sh(u)ch (u) - u)sech2(u) cn(u, k) = y w sech(u) - 0.25(1 - k2 )(sh(u)ch(u) - u)th(u) sech(u) ^
k =
V
y - sech(u) + 0.25(sh(u)ch(u) - u)th(u) sech(u)
0.25(sh(u)ch(u) - u)th(u) sech(u) dn(u, k) = y w sech(u) + 0.25(1 -k2)(sh(u)ch(u) + u)th(u) sech(u) ^
k =
y - sech(u) - 0.25(sh(u)ch(u) + u)th(u) sech(u)
V
cd(u, k) = cn(u k) = y,
dn(u, k)
-0.25(sh(u)ch(u) + u)th(u) sech(u)
sech(u) - 0.25(1 - k 2)(sh(u)ch(u) - u)th(u) sech(u)
sech(u) + 0.25(1 - k2 )(sh(u)ch(u) + u)th(u) sech(u)
k =
"V
(y -1) sech(u) + 0.25 sech(u)(y(sh(u)ch(u) + u)th(u) + (sh(u)ch(u) - u)th(u))
(12)
0.25(^к(и)ек(и) - и)Хк(и) ъесЬ(и) + у(^(и)е^(и) + и)й(и) ъесЬ(и)) На рис. 2 представлены графики функций (11) и (12) при параметре у, рассчитанном
для к = 0.5 .
-0.5
A B
Рис. 2. Графики нахождения параметра k по предварительно рассчитанному методом арифметико-геометрического среднего параметру y: A - аппроксимация тригонометрическими функциями, B - аппроксимация гиперболическими функциями. Fig. 2. Graphic find parameter k pre-determined method of arithmetic-geometric mean of the y parameter: A - approximation of trigonometric functions, B - approximation hyperbolic functions
Построив график данных функций, мы обнаруживаем очень плачевную ситуацию, данные функции имеют вертикальные разрывы и позволяют получить результат только на очень ограниченном участке периода:
Г (0,1 - 0,9) K (k) „
r \ i^K(к) - для тригонометрической аппроксимации (13)
r е (0,1 - 0,8)K(k) - для гиперболической аппроксимации
Так как функция К(к) не поддаётся предварительному расчету, без знания искомого параметра к мы не сможем предугадать, будет соблюдено условие (13) или нет. Поэтому применение алгоритма (11-12) невозможно.
Из этого можно сделать вывод, что на данный момент единственным алгоритмом, позволяющим получить гарантированный результат, является полный перебор параметра к в интервале [0,1] и проверка его методом АГС.
Сравнение алгоритмов вычислений
Сравним алгоритмы вычисления прямых эллиптических функций, взяв в качестве критерия точность вычисления. Так как метод АГС позволяет максимально приблизиться к эталонному значению, то все другие методы будем сравнивать с ним.
В табл. 1 приведены значения среднеквадратичных отклонений а каждого из методов. В качестве эталона был выбран метод АГС использующий 20 рекуррентных итераций. Расчет производился по формуле:
* = «(!/л)Е (X - У )2, ^
Гхг - значение методом А. Г С [у - значение сравниваемым методом Причем значения ряда Фурье были нормированы с коэффициентом 1/ cn(0, к), ряд Тейлора рассчитывался только по полупериоду 0 ± K(к), а для тригонометрических и гиперболических расчетов не применялось преобразование Ладена.
Таблица 1 Table 1
Точность различных методов вычисления Accuracy of different methods of calculation
k 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
Триг. 8.110-6 1.3-10-4 6.8-10-4 0.0023 0.006 0.0138 0.0306 0.0693 0.1782
Гипер. 0.0571 0.0551 0.0517 0.0467 0.0402 0.0319 0.0222 0.0117 0.0026
Тейлор 0.0011 0.004 0.0093 0.0182 0.0323 0.0548 0.0925 0.1703 0.1837
Фурье 0.0447 0.0452 0.0602 0.0965 0.1532 0.2314 0.3366 0.4825 0.7112
Из таблицы 1 видно, что при к < 0,7 лучшую точность показал метод аппроксимации тригонометрическими функциями, а при к > 0,7 - метод аппроксимации гиперболическими функциями. Однократное применение понижающего преобразование Ладена позволит уменьшить значение к в 5.. .20 раз, тем самым получить даже для к « 0,9 точность эквивалентную к « 0,1 ^ 0,2, а именно 10—4. ..10—6. Применение повышающего преобразования Ладена позволит увеличить к в 1,5.5,0 раз и для гиперболической аппроксимации при к « 0,3 получить точность эквивалентную к « 0,6...0,7, а именно 0,02...0,03. Из этого можно сделать вывод, что для тригонометрической аппроксимации преобразование Ладена дает очень хороший результат - улучшает точность на несколько порядков, в то время как для гиперболической аппроксимации удается улучшить точность только в 2-3 раза.
Разложения в ряды показали заметно худшую точность вычисления (для ряда Тейлора были взяты первые 4 члена, для ряда Фурье первые 10 членов).
Теперь сравним алгоритмы вычисления обратных по параметру и эллиптических функций Якоби.
Так как все методы являются численными, то в качестве критерия сравнения выберем скорость их сходимости.
В табл. 2 приведены значения количества итераций для достижения заданной точности. В качестве исходных параметров для расчета служили: к = 0.75, <р = 1.9728. Так же заметим, что количество итераций весьма условное значение, и может сильно варьироваться в зависимости от исходных данных.
Из таблицы (см. табл. 2) видно, что метод АГС хоть и является самым громоздким, но имеет лучшую сходимость, и, применяя 5-10 его итераций, можно добиться феноменальной точности порядка 10—5...10—8 практически при любых исходных параметрах.
Таблица 2 Table 2
Количество итераций для достижения заданной точности
T ie number of iterations to achieve specified accuracy
Погрешность 0.1 0.01 0.001 10-4 10-5
АГС 1 2 2 3 3
Средних прям. 2 4 16 32 64
Симпсон 2 4 6 8 10
Монте-Карло 200 2000 20000 200000 2000000
Методы численного интегрирования показали себя примерно одинаково по критерию сходимость/сложность. Метод средних прямоугольников обладает худшей сходимостью, но при этом значительно проще метода Симпсона. Метод Симпсона же, наоборот, имеет лучшую сходимость, но в плане вычислительной сложности он заметно проигрывает. Метод Монте-Карло приведен в данной работе как некая альтернатива «классическим» численным методам, и данный эксперимент показал, что он может работать. Однако он включает трудно реализуемый генератор случаемых чисел и обладает худшей сходимостью, поэтому его применение не целесообразно.
Сравним алгоритмы вычисления обратных по параметру к эллиптических функций Якоби. В качестве критерия возьмем среднеквадратическое отклонение полученного к с исходным значением. Для этого при разных параметрах к методом АГС рассчитаем значения у = cd(и,k), после этого по полученным значениям у найдем k = arckcd(и, у) . Для улучшения точности, принудительно ограничим функции в диапазоне [0,1].
Сравним величину среднеквадратического отклонения параметра к на всем периоде функции её. Полученные результаты занесем в табл. 3.
Гх,. - исходное значение k
*=J(i/(x-y)2, где
i=1
y - значение k полученное функцией arckcd
Таблица 3 Table 3
Среднеквадратическая погрешности вычисления параметра k
к 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
Триг. 8.4-10-5 3.8-10-5 0.0013 0.004 0.0078 0.0172 0.0252 0.0325 0.0421
Гипер. 0.4832 0.3726 0.2822 0.199 0.134 0.0845 0.0511 0.0269 0.0163
Как видно из таблицы 3, среднеквадратическая погрешность довольно небольшая и позволяет вычислять параметр к с точностью порядка 10-3...10-4, однако данный показатель это «средняя температура по больнице». Как уже было сказано в разделе 5, адекватный результат с очень хорошей точностью получается только при и е [0.05К(к) - 0.95К(к)]. Как только переменная и выходит из этого диапазона появляются разрывы и отклонения. Параметр К (к) может принимать любые значения в диапазоне [0.5, да), следовательно, диапазон адекватного расчета будет зависеть от неизвестного параметра к и принимать любые значения в диапазоне: [а, Ь], где a е[0.01, да), Ь е[0.51, да).
Выводы
1. В рамках данной статьи представлен практически весь объем алгоритмов, доступных в источниках, позволяющих вычислять эллиптические функции Якоби. Были рассмотрены обратные функции Якоби и алгоритмы их вычислений. Данный анализ охватывает практически весь математический аппарат, который требуется для синтеза эллиптических фильтров.
2. Различные вариации метода АГС показали себя лучше всего во всех рассмотренных выше разделах. Поэтому можно рекомендовать его в качестве основного для вычисления эллиптических функций Якоби.
3. Метод АГС требует довольно большого количества математических операций, поэтому в качестве более «компактного» аналога можно предложить аппроксимацию тригонометрическими и гиперболическими функциями. Следует заметить, что данная аппроксимация требовательна к параметру к. Уменьшить данное требование можно, применив понижающие или повышающие преобразования Ландена.
Список литературы References
1. Мирошниченко А.В., Волчков В.П., 2017. Три канонические формы представления импульсной характеристики аналогового фильтра Чебышева 1-го рода. Телекоммуникации и информационные технологии. 1: 89-93.
Miroshnichenko A.V., Volchkov V.P., 2017. Three canonical forms for representing the impulse response of a Chebyshev analog filter of the first kind. Telecommunications and information technologies. 1: 89-93.
2. Мирошниченко А.В., 2017. Обобщенный метод вывода канонических форм импульсных характеристик аналоговых фильтров. Фундаментальные проблемы радиоэлектронного приборостроения. М., МИРЭА, 4: 1194-1197.
Miroshnichenko A.V., 2017. Generalized method of output of canonical forms of impulse characteristics of analog filters. Fundamental problems of electronic instrument engineering. Moscow, MIREA, 4: 1194-1197.
3. Смирнов В.И., 1958. Курс высшей математики. Том 3, часть 2. М., Государственное издательство физико-математической литературы, 674.
Smirnov V.I., 1958. Course of higher mathematics. Vol. 3, Part 2. Moscow, State Publishing House of Physical and Mathematical Literature, 674.
4. Лем Г., 1982. Аналоговые и цифровые фильтры. Расчет и реализация. Под ред. к. т. н. Теплюка И.Н. М., Мир, 592.
Lem. G., 1982. Analog and digital filters. Calculation and implementation. Ed. k. t. n. Teplyuka I.N. Moscow, Mir, 592.
5. Абрамовиц М., Стиган И., Липман Д., Мак Ниш А. и д.р., 1979. Справочник по специальным функциям с формулами, графиками и математическими таблицами. Под ред. Абрамовиц М., Стиган И. М., Наука, 832.
Abramovic M., Stigan I., Lipman D., Mak Nish A. and other, 1979. A handbook on special functions with formulas, graphs and mathematical tables. Ed. Abramovic M., Stigan I. Moscow, Nauka, 832.
6. Березин И.С., Жидков Н.П., 1962. Методы вычислений. Том 1. М., Государственное издательство физико-математической литературы, 464.
Berezin I.S., Zhidkov N., 1962. Computational methods. Vol. 1. Moscow, State Publishing House of Physical and Mathematical Literature, 464.
7. Волчков В.П., Санников В.Г., 2017. Синтез предыскаженных финитных сигнальных базисов для борьбы с межсимвольной интерференцией. Системы синхронизации, формирования и обработки сигналов». М., Издательский дом Медиа паблишер. 8(3): 28-30.
Volchkov V.P., Sannikov V.G., 2017. Synthesis of preemphasis finite signal bases for reduction of intersymbol interference. Synchronization system, shaping and signal processing. M., Publishing House Media Publisher. 8(3): 28-30.
8. Волчков В.П., Безруков И.М., 2017. Сравнительный анализ помехоустойчивости OTDM систем с цифровым канальным прекодером и эквалайзером. Системы синхронизации, формирования и обработки сигналов». М., Издательский дом Медиа паблишер. 8(3): 25-27.
Volchkov V.P., Bezrukov I.M., 2017. Comparative analysis of noise immunity of OTDM systems with digital ^annel precoder and equalizer. Synchronization system, shaping and signal processing. M., Publishing House Media Publisher. 8(3): 25-27.
9. Санников В.Г., Волчков В.П., 2017. Синтез цифрового когерентного модема для одной модели канала связи с двойной случайностью. Материалы 27-й Международной Крымской конференции «СВЧ-техника и телекоммуникационные технологии». (КрыМиКо-2017), Севастополь, Сев.ГУ. 401-407.
Sannikov V.G., Volchkov V.P., 2017. Synthesis of digital coherent modem for one model of channel of communication with the double randomness. Materials of 27-th international Crimean Conference «Microwave equipment and telecommunication technologies». (KriMiKo-2017), Sevastopol, Sev.GU. 401-407.
10. Тактакишвили В.Г., Волчков В.П., 2017. Обработка траекторных данных в условиях априорной неопределенности. Материалы 27-й Международной Крымской конференции «СВЧ-техника и телекоммуникационные технологии». (КрыМиКо-2017), Севастополь, Сев.ГУ, 359-363.
Taktakishvili V.G., Volchkov V.P., 2017. Trajectory Processing of data in terms of a priori uncertainty. Materials 27-th international Crimean Conference "Microwave equipment and telecommunication technologies. KriMiKo-2017, Sevastopol, Sev.GU, 359-363.
11. Волчков В.П., Асирян В.М., 2017. Вычислительно эффективный алгоритм формирования оптимального базиса Вейля-Гейзенберга. Фундаментальные проблемы радиоэлектронного приборостроения. М., МИРЭА. Часть 4, 1151-1154.
Volchkov V.P., Asiryan V.M., 2017. Computationally efficient algorithm of forming an optimal basis for the Weyl-Heisenberg. Fundamental problems of electronic instrument engineering. Moscow, MIREA. Part 4,1151-1154.
12. Волчков В.П., Тактакишвили В.Г., 2017. Многогипотезный алгоритм радиолокационного слежения с адаптацией по порогу обнаружения. Фундаментальные проблемы радиоэлектронного приборостроения. М., МИРЭА. Часть 4, 977-980.
Volchkov V.P., Taktakishvili V.G., 2017. Multihypothesis radar tracking algorithm with the adaptation on the threshold of detection. Fundamental problems of electronic instrument engineering. Moscow, MIREA. Part 4, 977-980.
13. Sannikov V.G., Volchkov V.P., 2017. Synthesis of preemphasis finite signal bases for reduction of intersymbol interference. 2017 Systems of Signal Synchronization, Generating and Processing in Telecommunications, SINKHR0INF0-2017, 799751, Institute of Electrical and Electronics Engineers In.
14. Безруков И.М., Волчков В.П., 2017. Применение сингулярного разложения для улучшения стойкости алгоритма Прони к аддитивным шумам. Телекоммуникации и информационные технологии». 4(1): 57-62.
Bezrukov I.M., Volchkov V.P., 2017. Use singular decomposition to improve resistance to Proni algorithm additive noise. Telecommunications and information technologies. 4(1): 57-62.
15. Волчков В.П., Санников В.Г., 2016. Адаптивное канальное прекодирование в телекоммуникационных системах. Материалы 26-й Международной Крымской конференции «СВЧ-техника и телекоммуникационные технологии». (КрыМиКо-2016), Севастополь, Сев.ГУ. 665-671.
Volchkov V.P., Sannikov V.G., 2016 Adaptive channel precoding in telecommunication systems. Materials of 26-th international Crimean Conference «Microwave equipment and telecommunication technologies». (KriMiKo-2016), Sevastopol, Sev.GU. 665-671
16. Волчков В.П., Шурахов А.А., 2016. Применение технологии MIMO с линейным прекоди-рованием в радиорелейных системах прямой видимости. Фундаментальные проблемы радиоэлектронного приборостроения. М., МИРЭА. 16(5): 210-216.
Volchkov V.P., Shurahov A.A., 2016. MIMO technology with linear precoding in microwave line-of-sight systems. Fundamental problems of electronic instrument engineering. M., MIREA. 16(5): 210-216.
17. Волчков В.П., Санников В.Г., 2016. Синтез канальных прекодеров для цифровых систем связи с финитным сигнальным базисом. Электросвязь, 4: 41-45.
Volchkov V.P., Sannikov V.G., 2016. Synthesis of channel precoder for digital communication systems with finite signal basis. Electrosvyaz. 4: 41-45.
18. Волчков В.П., Шлома А.М., Поборчая Н.Е., 2015. Параметрический спектральный анализ случайных сигналов с использованием рекуррентных циркулянтных моделей скользящего окна. Электросвязь. 5: 20-23.
Volchkov V.P., Shloma A.M., Poborchaja N.E., 2015. Parametric spectral analysis of random signals with use of recurrent circulant models of the moving window. Electrosvyaz. 5: 20-23.
19. Волчков В.П., Шлома А.М. 2014. Синтез рекуррентных фильтров скользящего окна в базисах функций Виленкина-Крестенсона. Электросвязь. 5: 17-19.
Volchkov V.P., Shloma A.M. 2014. Synthesis of recurrent filters of the moving window in bases of Vilenkin-Krestenson functions. Electrosvyaz. 5: 17-19.
20. Волчков В.П., Уваров C.C. 2015. Аппроксимация узкополосных случайных процессов с помощью комплексной рекурентной m-модели скользящего окна второго порядка. T-Comm: Телекоммуникации и транспорт. 9(3): 54-61.
Volchkov V.P., Uvarov S.S. 2015. Approximation of narrowband random signal using vectotial recurent circulant m-model of sliding window of the second order. T-Comm: telecommunications and transport. 9(3): 54-61.
21. Волчков В.П., Поборчая Н.Е., Шлома А.М., 2013. Аппроксимация случайных сигналов с помощью рекуррентных циркулянтных моделей скользящего окна. Научные ведомости БелГУ. Сер. История. Политология. Экономика. Информатика. 22(165): 135-143.
Volchkov V.P., Poborchaya N.E., Shloma A.M., 2013.The random process representation with the help of vectorial recurrent circulant model of the second order. Nauchnye vedomosti BelGU. Istoriya. Politologiya. Ekonomika. Informatika. [Belgorod State University Scientific Bulletin. History Political science Economics Information technologies]. 22(165): 135-143.
22. Волчков В.П., Шурахов А.А. 2014. Возможности линейных прекодеров по управлению ресурсами и характеристиками систем MIMO. Научные ведомости БелГУ. Сер. История. Политология. Экономика. Информатика. 8(179): 172-181.
Volchkov V.P., Shurahov A.A. 2014. Posibilities of liner precoders to manage resoureses and performance of mimo systems. Nauchnye vedomosti BelGU. Istoriya. Politologiya. Ekonomika. Informatika. [Belgorod State University Scientific Bulletin. History Political science Economics Information technologies]. 8(179): 172-181.
23. Волчков В.П. 2009. Новые технологии передачи и обработки информации на основе хорошо локализованных сигнальных базисов. Научные ведомости БелГУ. Сер. История. Политология. Экономика. Информатика. 15(70): 181-189.
Volchkov V P. 2009. A new technology of transmitting and processing of jnformation based on well-localized signal basis. Nauchnye vedomosti BelGU. Istoriya. Politologiya. Ekonomika. Informatika. [Belgorod State University Scientific Bulletin. History Political science Economics Information technologies]. 15(70): 181-189.
24. Волчков В.П., Петров Д.А., 2009. Условия ортогональности обобщенных базисов Вейля-Гейзенберга для OFTDM сигналов. Научные ведомости БелГУ. Сер. История. Политология. Экономика. Информатика. 15(70): 190-199.
Volchkov V.P., Petrov D.A. 2009. Generalized Weyl-Heysenbeg bases orthogonality conditions for OFTDM sygnals. Nauchnye vedomosti BelGU. Istoriya. Politologiya. Ekonomika. Informatika. [Belgorod State University Scientific Bulletin. History Political science Economics Information technologies]. 15(70): 190-199.
25. Волчков В.П., Петров Д.А., 2009a. Оптимизация ортогонального базиса Вейля-Гейзенберга для цифровых систем связи, использующих принцип OFDM/OQAM передачи. Научные ведомости БелГУ. Сер. История. Политология. Экономика. Информатика.1(56): 102-112.
Volchkov V.P., Petrov D.A., 2009a. Orthogonal Weyl-Heisenberg basis optimisation for digital communication systems based on OFDM/OQAM. Nauchnye vedomosti BelGU. Istoriya. Politologiya. Ekonomika. Informatika. [Belgorod State University Scientific Bulletin. History Political science Economics Information technologies]. 1(56): 102-112.
26. Черноморец А.А., Волчков В.П., 2012. О свойствах квазисубполосных и G-субполосных матриц. Научные ведомости БелГУ. Сер. История. Политология. Экономика. Информатика. 21(120): 127-134.
Chernomorets A.A., Volchkov V.P., 2012. About properties of quasisubband and G-subband matrixes. Nauchnye vedomosti BelGU. Istoriya. Politologiya. Ekonomika. Informatika. [Belgorod State University Scientific Bulletin. History Political science Economics Information technologies]. 21(120): 127-134.
27. Бурданова Е.В., Маматов Р.А., Волчков В.П., 2011. Техническая реализация алгоритмов обработки радиолокационных изображений при обнаружении объектов наземной поверхности. Научные ведомости БелГУ. Сер. История. Политология. Экономика. Информатика. 19(108): 196-203.
Burdanova E.V., Mamatov R.A., Volchkov V.P., 2011. Technical realization of algorithms of processing radiolo-katsionnyh of images at detection of objects on the terrestrial surface. Nauchnye vedomosti BelGU. Istoriya. Politologiya. Ekonomika. Informatika. [Belgorod State University Scientific Bulletin. History Political science Economics Information technologies]. 19(108): 196-203.
28. Волчков В.П., Петров Д.А., 2010. Обобщенная теорема Найквиста для OFTDM сигналов. Системы синхронизации, формирования и обработки сигналов. М., Издательский дом Медиа паблишер. 1(1): 28-32.
Volchkov V.P., Petrov D.A., 2010. The generalized Nyquist's theorem for OFTDM of signals. Synchronization system, shaping and signal processing. M., Publishing House Media Publisher. 1(1): 28-32.
29. Волчков В.П., Шурахов А.А., 2012. Синтез двухкритериальных линейных прекодеров для системы MIMO. Системы синхронизации, формирования и обработки сигналов. М., Издательский дом Медиа паблишер. 3(1): 35-39.
Volchkov V.P., Shurahov A.A., 2012. Synthesis of two-criteria linear precoder for the MIMO system. Synchronization system, shaping and signal processing. M., Publishing House Media Publisher. 3(1): 35-39.
30. Волчков В.П., Шурахов А.А., 2012a. Исследование эффективности алгоритмов линейного прекодирования в системах MIMO. Электросвязь. 5: 15-16.
Volchkov V.P., Shurahov A.A., 2012a. Research of efficiency of algorithms of a linear precoding in the MIMO systems. Electrosvyaz. 5: 15-16.