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

Таким образом, разработан алгоритм индивидуальной градуировки МПП различных «спектральных» типов систематической погрешности, позволяющий получить градуировочную характеристику МПП, удовлетворяющую заданной точности и требующую минимальный объём памяти; разработано ПО оптимизации процесса градуировки МПП, позволяющее во время градуировки определить необходимые «спектральные» свойства систематической погрешности градуируемого МПП и выбрать оптимальный для данного «спектрального» типа статической характеристики способ градуировки; исследована зависимость эффективности индивидуальной градуировки МПП от «спектральных» свойств систематической погрешности.

Литература

1. Надеев А.И. Интеллектуальные магнитострикционные преобразователи информации // Датчики и системы. 2002. № 5. С. 16 - 20.

2. Надеев А.И., Радов М.Ю., Вдовин А.Ю. Способ нормирования статической погрешности магнитострикционных преобразователей перемещений // Датчики и системы 2002. № 5. С. 25 - 26.

3. Свидетельство об официальной регистрации программы для ЭВМ № 2004612143 от 17. 09. 2004 г.

4. Надеев А.И., Радов М.Ю., Вдовин А.Ю. Автоматизация многофакторных экспериментов при исследовании маг-нитострикционных преобразователей // Материалы науч.-техн. конф. с участием зарубежных специалистов. «Дат-чик-2003». 2003.

5. Радов М. Ю. Автоматизация экспериментальных исследований магнитострикционных преобразователей перемещений // Тез. докл. VIII Междунар. конф. «Образование. Экология. Экономика. Информатика». (Астрахань 15 - 20 сент. 2003 г.) Астрахань, 2003.

Астраханский государственный технический университет 18 февраля 2005 г.

УДК 681.3(088.8): 518.5

ЦЕЛОЧИСЛЕННЫЕ АЛГОРИТМЫ ГЕНЕРАЦИИ ГАРМОНИЧЕСКИХ СИГНАЛОВ

© 2005 г. И.Н. Булатникова

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

равного ~ 1,11). 242

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

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

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

Цифровая интерполяция состоит в замене исходной линии (кривой) ломаной линией, проходящей через узлы решётки с единичным отстоянием узлов друг от друга по вертикали и по горизонтали.

Основными параметрами при цифровой интерполяции прямой являются конечные приращения по

осям X и У интерполируемого отрезка прямой АХ и АУ . Их отношение задаёт крутизну прямой.

При цифровой интерполяции произвольных кривых это отношение задаёт угловой коэффициент касательной к кривой. Он различен в каждом 1-м узле и поэтому появляется индекс у АХг и АУг. Для произвольных кривых АХг выбирается произвольным целым числом, определяющим быстродействие и точность цифровой интерполяции.

Суть цифровой линейной интерполяции можно видеть на рис. 1.

АХ = 5, АУ = 3, Е 0 =-1

\

/ и 1

0 0 1 2 3 4 5

Рис. 1. Цифровая интерполяция прямой

Как и в случае линейной интерполяции, вводится оценочная функция, учитывающая отклонение в каждом г-м узле (при |АХг | > |АУ-1)

Е г =K к 2 |Д7.| .

(1)

3

2

1

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

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

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

Приведём правило выбора очередного шага в обобщённом алгоритме (при |ДХ,-

I > N I).

Если Ег > 0 , то х,+1 = х + 8^пДХ,, у,+1 = Уi;

иначе

xM = хг + signAX, Уг+1 = Уг + Sign AYi •

(2)

Корректировку оценочной функции запишем следующим образом

Ег+1 = Ег + |АХг+1 | - 2 |AY+ | - |AXj | Sign Е,

(3)

где ДХг, ДХг +1, ДУг+1 - значения величин ДХ и ДУ в текущем (г-м) и очередном узлах интерполяции.

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

Генерацию синусоиды проведём с использованием вспомогательной функции - ее производной -косинусоиды, т. е. одновременно будут генерироваться

хх синусоида М 81и(—) и косинусоида М со8(—), где

М М

М - масштабный коэффициент.

Введем обозначения = М зт(—) и

М

г

сг = М соб(—). Тогда приращения величин Si и Сг

М

на интервале (г, г+1) будут следующими:

= ± с + (сг + Дсг ) ; М 2

Ac = _ Sj + (Sj +Asj)

г M 2

(4)

(5)

Формулы (4) и (5) получены на основе усреднения производных в начале и конце единичного отрезка между соседними узлами интерполяции. При этом

учтено, что производная от синуса равна косинусу, а производная от косинуса - минус синусу. Упрощая, имеем систему уравнений:

2С, + ДС, = 2МД5,;

-(25, + Д5,) = 2МДСг.

Решая её относительно Д5, и ДСг, получаем

4Мсг - 24Мз{ + 2сг Дsi =-—--; Дсг =--—--.

4M 2 +1

4M 2 +1

Отсюда с учётом того, что средняя крутизна на интервале равна отношению приращения функции на этом интервале к длине (она равна единице) интервла, получим:

. с 4Мсг - 2si . С 4М5г + 2сг ДУ,5 = г--'-; ДУ,С = -г--- ;

4M2 +1

4M2 +1

AXf = AXC =AX = r •

Чтобы исключить деление на (4М + 1), увеличим ДУ5, ДУС, ДХ в (4М2 + 1) раз, т. е. выберем г равным (4М 2 + 1).

Окончательно запишем:

AYiS = 4Mc г _ 2s г; A Y;C = _(4Ms г + 2с г); AXS = 4M 2 +1.

(6)

Формулы (6) даны в терминах обобщённого алгоритма [1].

Следует отметить, что усреднение крутизны кривых и с■ проведено по формулам, эквивалентным интегрированию по методу трапеций. Фактические приращения больше в q раз = 2М tg (2М) 1).

Компенсируем уменьшение ДХ4 и ДУс, обусловленное применением формулы интегрирования по методу трапеций. Для этого увеличим их в q раз. Однако более удобно уменьшить в q раз дх . Итак,

4M2 +1 „ , ,

AX =-т. Заменив tg(2M)

его разложе-

2М tg (2М )-

нием в степенной ряд, при М ^^ получим ДХ = 4М 2 + 2/3 .

При больших значениях М величиной 2/3 можно пренебречь, иначе умножаем на 3/2, чтобы получить целые величины

ДУ,5 = 6Мс, - 34,; ДУ,С = -(6Муг + 3с,); (7)

ДХ,5 = 6М 2 +1.

2

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

путём прибавления к старому значению 6Мя, -1 (или 6Мс,-1) величины ±6М, если было приращения на ±1 синуса (или косинуса).

Другим источником погрешности определения

АУг и АУС является округление текущих значений гг и сг до ближайших целых значений при использовании их в (6). Попытаемся устранить его.

Относительно оптимального алгоритма цифровой линейной интерполяции, лежащего в основе предложенного алгоритма, следует заметить, что вводимая в [2] величина аг, связанная с Ег соотношением

Ег + 2 |АУг — I АХ 1—11

а г =-!—:-—г--, равна отстоянию по гори-

2 |А7г—1|

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

АУг —1 1 Ег +2АУг —1 у = V —а -= V +----

^точн ^ г АХ, —1 ^ 2 2АХг —1

где АУг—1 значение АУг для синуса (или АУС для косинуса) в (г - 1)-м узле. Однако в этом случае необходимо провести операцию деления, которую можно заменить сдвигом, выбрав АХ равной целой степени двух. Кроме того, в выражениях (7) появится умножение 6Мгг, 6Мсг, которое раньше заменялось прибавлением ±6М к старому значению произведения, если было приращение ±1 у синуса (косинуса). Теперь эту операцию осуществить невозможно, так как приращение синуса (косинуса) не равно ±1 .

Обеих трудностей можно избежать одновременно, выбрав М равным целой (положительной) степени двух, а АХ - равным 4М2, что всегда допустимо, так как 4М2 >> 2/3 . Однако в этом случае период синусоиды, равный 2пМ, будет строго фиксированным, что неудобно при необходимости изменения частоты гармонических сигналов.

Конечно, при использовании формулы (6) потребуются дополнительные разряды вычислителя для представления дробной части отсчетов у точн , промежуточных и окончательных результатов при расчётах АУгг, АУгС и Ег+1, что также крайне нежелательно.

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

Применим алгоритмы (2), (3), (6) для случая, когда амплитуда синусоиды (косинусоиды) равна N (|N < М), т.е. для случая генерации

у = N бш—, г = N соб-— . Определив их производ-

ММ

/ г , у ные и представив их как у =--, г =—, можно

ММ

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

Таким образом, алгоритм (2), (3), (6) сохраняет достоверность и при N < М. Отличие алгоритма для

случая N1 < М от алгоритма для N = М заключается в том, что начальные значения синуса и косинуса будут равны 0 и N (• 0 = 0, с0 = К). В итоге амплитуда сигнала равна N а период равен 2пМ.

Итак, в окончательном виде алгоритм генерации гармонических сигналов (синуса и косинуса) таков (в нижеприведенных формулах •г и верхний индекс • относятся к синусу, а сi и верхний индекс с - к косинусу):

• 0 = 0, с 0 = М; х 0 = 0;

АХ 0 = 6М 2 +1, АУо = 6М 2, А У0С = —3М;

Е• =|АХ0| — 2|АУо |, ЕС =|АХ^ — 2|АУ0С|.

Для г = 0,1,2,...

хг+1 = хг + 1, ах,+1 = АХг .

Если Ег > 0, то гг+1 = гг, иначе гг+1 = гг + signАУis Если Ei > 0, то Сг+1 = с{, иначе Сг+1 = Сг + signАУic

АУг+1 = 6Мс г+1 — 3гг+1, АУг+1 = —(6Мг г+1 + 3с г+1),

Е*+1 = Е• + |АХг +11 — 2 |АУг+1 — |АХг| ^пЕ*,

ЕгС+1 = ЕС + |АХг+11 — 2|АУг+^ — |АХг|signEcl .

Возврат цикла г = г +1.

Результаты цифрового моделирования предложенного алгоритма приведены в таблице. Для первой четверти синусоиды и косинусоиды в числителе даны максимальные ошибки (по модулю) генерации при разных М и N а в знаменателе - число отсчётов, для которых абсолютная ошибка (по модулю) превышает 0,5.

Таблица

Результаты моделирования алгоритма генерации гармонических сигналов

Функция Параметры

N=10 M=10 N=75 M=100 N=100 M=100 N=375 M=500 N=500 M=500 N=1000 M=1000

x N sin — M 0,5/0 0,5/0 0,510/2 0,504/2 0,526/8 0,521/9

x N cos — M 0,5/0 0,514/1 0,534/2 0,507/1 0,511/8 0,511/16

Генерация синуса и косинуса в последующих четвертях при необходимости осуществляется цифровой интерполяцией фрагмента кривой синуса и косинуса для первой четверти в обратном направлении (сменой знаков АХ, АУг, АУС ) с соответствующей модификацией алгоритма. Возможен и другой способ генерации цифровых отсчётов гармонических сигналов путем предварительного насчёта приращений цифровых отсчетов уг по (2), причём хг является номером отсчёта, а сам факт приращения (или отсутствия тако-

вого) кодируется одним битом (1 или 0). Кортеж этих битов для последовательных отсчётов запоминается в ОЗУ (удобно ОЗУ последовательного типа) для одной четверти периода синусоиды или (и) косинусоиды. Затем в режиме генерации оно выдаёт этот кортеж битов на реверсивный счётчик, который, подсчитывая их, при достижении величины ±M переключается в режим вычитания (сложения), а ОЗУ - в режим обратного (прямого) пересчета адресов. Преимущество второго способа (с использованием быстродействующего ОЗУ) состоит в возможности генерации синусоиды (косинусоиды) с большей частотой, чем быстродействие вычислителя, работающего в режиме интерполяции.

Очевидно, что с помощью алгоритма одновременного получения синусоиды и косинусоиды можно осуществить цифровую генерацию окружности. Однако известны более простые алгоритмы круговой интерполяции [2], поэтому более интересной представляется генерация эллипса с произвольным эксцентриситетом e. Для этого генерируются две гармонические кривые разной амплитуды с постоянным сдвигом фазы П 2.

Проведя масштабирование одной кривой, например, косинусоиды, на величину е, получаем выражения для АХ, AYs, AYc

AY? = 6Mct /e - 3sj; AYic = -(6Mest + 3сг);

АХ = 6M 2 + 1, (8)

которые должны быть входными параметрами для алгоритма (1), (2), (3), генерующего синусоиду с амплитудой M и косинусоиду амплитудой eM , где M -большая полуось эллипса. Для сохранения целочис-ленности алгоритма константы 6M/e и 6Me должно быть целыми, либо округлёнными до целого. Конечно, во втором случае получается определённая погрешность, величина которой зависит обратно пропорционально от M.

Из сравнения предложенного метода с традиционными, например с вычислением синуса через степенной ряд или по рекуррентной формуле si+1 = si cosh + ci sinh (где h = 1/M - шаг интерполяции), вытекает, что он превосходит традиционные по быстродействию и аппаратурным затратам.

Действительно, разложение синуса в степенной ряд при методической погрешности вычислений менее 10-5 требует наличия, как показано в работе [3], пяти членов. Машинная реализация потребует 10 умножений и 5 сложений. Хотя во втором случае (рекуррентная формула) умножений и сложений меньше (соответственно 4 и 2), но разрядность представления констант cos h и sin h, а также текущих значений si и ci должна быть в два-три раза (обозначим символом n) больше, чем это требуется для представления номеров отсчётов функции, так как погрешность округления произведений накапливается и передаётся в последующие отсчёты. Таким образом, в вычислительном плане второй метод эквивалентен 4n2 умножениям и 2n сложениям обычной разрядности.

Умножение проводится минимум в 10 раз медленнее сложения (например, быстродействие однокристального микропроцессора при умножении и сложении [4]), поэтому на каждый отсчёт первый способ (степенной ряд) требует временных затрат, эквивалентных 105, а второй (при n = 2) - 164 сложениям. Предложенный способ использует всего 6 сложений на один отсчёт. Причём этот выигрыш увеличивается с возрастанием требований по точности генерации.

Для подтверждения этого выигрыша приведём микроассемблерную программу, которая весьма короткая:

LOR: SUB R4, R0; вычитание BPL KOR; усл. переход SUB # T, R4; вычитание DOL: DEC R4; приращение на 1 DEC R5; приращение на 1 SUB # P, R1; вычитание ADD # Q, R0; сложение SAL: ADD # T, R4; сложение

HALT; остановка KOR: SUB R1, R3(вход); вычитание BPL LOR; усл. переход DEC R1; приращение на 1 INC R2; уменьшение на 1 ADD # Q, R3; сложение SUB R4, R0; вычитание BPL SAR; усл. переход

BR DOL; безусл. переход.

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

Соответственно, координаты точек эллипса будут изменяться по закону синуса и косинуса.

Исходные значения содержимых регистров (для

-1) следующие:

R0 = M 2 +1/6 = M 2; R1 = 2(M /e)y0 - х0;

R2 = х0; R3 = M2 +1/6 = M 2;

R4 = 2Mex0 + y0 ; R5 = y0 ,

где (х0, y0) - начальный узел интерполяции. Текущие координаты Х y¿) - в регистре R2 и R5 (после стопа). В качестве констант используются следующие величины: P = 2M/e , Q = 2М 2 , T = 2eM .

Вывод

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

Одним из основных подходов при исследовании устойчивости нелинейных динамических систем является исследование устойчивости по первому линейному приближению, т.е. изучение устойчивости линейной части этих систем.

Часто на практике коэффициенты характеристических полиномов исходной системы зависят от к технических параметров а, моделируемого объекта, которые могут быть модифицированы в процессе проектирования, изготовления и эксплуатации и, следовательно, образуют целое множество возможных значений коэффициентов характеристических полиномов а, (а 1,..., ак), , = 0,..., п .

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

Теорема Харитонова, в сильной степени, помогает решить эту проблему робастной устойчивости [5], однако она касается только семейств интервальных полиномов, т. е. того случая, когда выпуклое множество коэффициентов полиномов Гурвица задано в виде п-мерного прямоугольного параллелепипеда.

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

Литература

1. Булатникова И.Н. и др. Обобщенный алгоритм цифровой интерполяции // Новые информационные технологии. М., 2004. С. 6-8.

2. Анишин Н.С., Тивков А.М. Оптимальный алгоритм цифровой линейной информации// Изв. вузов СССР. Приборостроение. 1983. № 8. С. 56-59.

3. Байков В.Д., Смолов В.Б. Аппаратурная реализация элементарных функций в ЦВМ. Л., 1974.

4. Однокристальные микропроцессоры комплекта БИС серии К1801 / В.Л. Дшхунян и др. // Микропроцессорные средства и системы. 1984. № 4. С. 12-18.

2004 г.

Определение 1. Интервальным полиномом называется множество полиномов вида

Р = {/(А) = а0 + а{К + а2А 2 +... + ап-1А п-1 + А п },(1)

где а, - действительные числа, удовлетворяющие неравенствам а, < а, < а, , 0 <, < п .

Определение 2. Интервальный полином является полиномом Гурвица, если любой полином, принадлежащий этому множеству, является полиномом Гурвица.

Введем в рассмотрение четыре «угловых» полинома:

/1 (X) = а0 + а 1Х + а2А2 + а3А3 + а4А4 +... + А п,

/2(А) = а0 + + а2А2 + а3А3 + а4А 4 +... + А п,

Г _ _ Г (2)

/3(А) = а 0 + + а 2А2 + а 3А3 + а 4А 4 +... + А п },

/4(А) = а 0 + + а 2А2 + а 3А3 + а 4А 4 +... + А п }.

Теорема Харитонова. Для того чтобы интервальный полином был полиномом Гурвица необходимо и достаточно, чтобы четыре «угловых» полинома (2) являлись полиномами Гурвица [5].

Доказательство. Необходимость. Так как четыре «угловых» полинома (2) принадлежат множеству полиномов (1), то они являются полиномами Гурви-ца, если интервальный полином (1) является полиномом Гурвица.

Достаточность. Покажем, что доказательство непосредственно вытекает из критерия Михайлова.

Кубанский государственный технический университет, г. Краснодар 6 декабря

УДК 517.929

НЕОБХОДИМЫЕ И ДОСТАТОЧНЫЕ УСЛОВИЯ СУЩЕСТВОВАНИЯ ВЫПУКЛОЙ ОБЛАСТИ УСТОЙЧИВОСТИ В ПРОСТРАНСТВЕ КОЭФФИЦИЕНТОВ ХАРАКТЕРИСТИЧЕСКОГО МНОГОЧЛЕНА

© 2005 г. Л.Д. Блистанова, Г.А. Зеленков, Е.В. Стрюк