Генерация синуса

Мне показалось не все знают такие простые способы. Можно генерировать последовательные значения sin(t) без таблиц и каких либо тяжелых вычислений. Суть именно в том, что генерируется последовательность на регулярной сетке по времени, для примера sin(0.1), sin(0.2), и тд. Для этого случая есть простой способ.


sin(t) можно получить как решение дифференциального уравнения.



Которое можно пребразовать к системе уравнений первого порядка, обозначив производную от x по времени за новую переменную.



Или тоже самое с использованием матриц.



Здесь, x — значение sin(t), v — скорость изменения sin(t) или cos(t).

Решение линейных стационарных систем, такого простого вида,



может быть получено с помощью матричной экспоненты,



Где, x0 — начальное состояние системы, t — время.

Матричная экспонента определяется через ряд.



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





Здесь, T — период дискретизации.

Сходим в octave, проверим.


A = [0 1; -1 0];
T = 0.1;

t = 0:T:10;
x = [1; 0];
xl = [];

Ad = expm(A * T);

for i=t
	xl(:,end+1) = x;
	x = Ad * x;
end

plot(t, xl(1,:));
grid on

pause




Как можно видеть в цикле осталось только умножение A на x. Это четыре умножения в общем случае. Но можно использовать приближенное вычисление Ad.




%Ad = expm(A * T)
Ad = eye(2) + (A * T)




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

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


for i=t
	xl(:,end+1) = x;
	x = Ad * x;

	len = sqrt(x' * x);
	x *= 1/len;
end




Чудесно, амплитуда вернулась на место. И надо заметить никуда теперь не денется. Как долго бы генерация не длилась. Однако, была использована функция sqrt. Не самая легкая в вычислении.

Попробуем заменить её на, что нибудь более простое. Ведь точное нормирование не требуется, достаточно того, чтобы процесс сходился.

Возмем вот такое грубое приближение в точке 1.




	%len = sqrt(x' * x);
	%x *= 1/len;
	x *= (3 - (x' * x)) / 2;




Это финальный результат, 5 умножений.


x = Ad * x;
x *= (3 - (x' * x)) / 2;


Для задания частоты меняем Ad, ни одного умножения.


Ad = eye(2) + (A * T);


Замечу, что при крупном шаге будет заметно отклонение частоты.

У этого способа есть геометрическая интерпретация, x — вектор в двухмерном пространстве, Ad — матрица поворота на угол T в радианах. Можно использовать оба элемента вектора, если нужен синус одновременно в разных фазах.

Извиняюсь если кому не понятна матлабо-октавская запись, по компонентам расписывать лениво. Поиграйтесь с кодом в octave, все должно стать ясно.

Да и наверно можно ещё ускорить расчет, если очень хочется. Например делать нормирование не на каждом шаге.

Добавка

alekseyb и _pv показали на другой способ, требующий меньше вычислений.





Формальный вывод понятен (ссылки в комментах) но до интуитивного понимания пока не расковыривал.

p.s. пока сохранял в черновики, получил в ответ «Hacking attempt».
  • +11
  • 15 октября 2012, 20:42
  • amaora

Комментарии (100)

RSS свернуть / развернуть
«p.s. пока сохранял в черновики, получил в ответ «Hacking attempt».»
нуцк выж синус хакнули :)
+3
Перенес в алгоритмы.
0
Не, всё круто…
Но вы что, серьёзно предлагаете, чтобы засунуть это в микроконтроллер, ставить Матлаб? Да ещё и изучать его? Оригинальное решение. :-)
А, забыл! Ещё и выяснять, что такое octave. :-) Я в первый раз слышу.
0
  • avatar
  • Alfa
  • 15 октября 2012, 23:23
Если оно лично Вам неизвестно — значит это лично Вам ненужно (или не судьба). Бывает.
Но здесь достаточно людей, владеющих матлабом и октавой.

Чтобы «засунуть» — ненужно, конкретный алгоритм уже разжёван и в рот положен. Можно просто сказать «спасибо».

С уважением
-1
За статью действительно спасибо. Но…

x *
— это шо?
x'
— это шо?
eye(2)
— мда…

Вывод: алгоритм — да. Реализация — нет.
0
x *

Умножение.

x'

Транспонирование.

eye(2)

Матрица с единицами по главной диагонали и нулями в остальных элементах размерности 2x2.

Можете ещё посмотреть g*gle-> «octave basics», там это все есть.
0
Было бы неплохо, если бы ты еще на С код привел.
+1
  • avatar
  • Vga
  • 16 октября 2012, 00:03
Если не ошибаюсь Matlab Compiler позволяет конвертировать M-файлы автоматически в C. Интересно было бы взглянуть на этот сконвертированный код — возможно автор захочет написать продолжение статьи.
0
Это ж вроде Octave, а не Matlab?
0
Окатава работает на языке матлаба. Совместимость очень высокая. Если не используются специфические фичи матлаба — часто бывает достаточно небольшой правки. А то и вовсе без неё.
0
Да, матлаб умеет генерить сишный код, но вот разобраться в нём, мягко скажем, сложно.
0
Пардон муа, но эти несколько строк — даже без бумажки, в уме можно перевести. Ну ок — с бумажкой.
Достаточно помнить как делается перемножение и транспонирование матриц.

При матрице 2x2 это вообще выходит код простой и небольшой.
0
Пардон муа, но эти несколько строк написаны на языке, который я не знаю.
+1
Вы точно уверены, что не знаете? ;) Обычно в институте такое проходят.
0
Матлаб? Ой ли? Хотя у меня был курс матлаба, но я его уже забыл основательно, и он читался вместое курса Delphi или C++.
0
Какой матлаб?! =))) Я про запись матричных операций НА ДОСКЕ =) Почти оди в один как тут.
Перемножение и транспонирование. Ну ок, транспонирование на доске могут указывать надчёркиванием вместо апострофа. Но апострофом тоже часто обозначают.
0
Перемножение и транспонирование я знаю. Мне непонятен текст на Октаве. И мне казалось, матлаб куда больше на С походил.
0
Синтаксис у матлаба и октавы один и тот же. Сугубо IMHO это может быть похожим на C не больше чем perl или tcl.
0
Тогда я тем более оный язык не знаю. Мне в листингах разве что +-*/ понятны.
0
Пардон муа за мой французский, но вот лично от тебя таких слов ну никак не ожидал увидеть. Вроде опытный разработчик, хороший программист, а начинает непонятный холивар не тему того, что он не понимает простой скрипт на M (особенно с учётом того, что он его когда-то изучал).
Синтаксис M в целом не сложен и разобраться в нём с нуля среднестатистическому императивщику гораздо проще, чем, например, в лиспе, хаскеле или даже луа. Нужно только осознавать, что бывают нетипизированные языки и понимать изначальную специфику матлаба (работа с матрицами и векторами). ИМХО, единственная конструкция, сложная для понимания далёкому от M программисту — «xl(:,end+1) = x;», но и её смысл можно интуитивно понять с учётом стоящей выше строчки «xl = [];». Ну ещё нетипичная для ЯВУ общего назначения операция транспонирования — но это не повод устраивать перепалку =).
ИМХО, единственное, в чём неправ автор — не указал явно, что октав — это, грубо говоря, опенсоурсная имплементация матлаба, потому что слово Matlab знают все современные инженеры, даже если они его ни разу в глаза не видели, а вот слово Octave — только те, кто увлекается опенсорсом, линуксом, или просто предпочитает честно юзать софт, а покупать матлаб не хочет =).
0
Да, не понимаю. И что с того? Изучить руки пока не доходят, тем более что я вообще математическими пакетами не особо интересуюсь.
0
P.S. Я бы даже сказал, что суть метода я понял не из формул и объяснений автора, а из короткой приписки «геометрический смысл» в конце поста.
0
Статья так красиво начиналась: "… простые способы.… генерировать… без таблиц и… тяжелых вычислений.… на примере sin(0.1), sin(0.2), и тд. Для этого случая есть простой способ."
А в итоге вылилось в сложную «сложную» математическую формулу, да еще с "… Теперь будем подставлять костыли, чтобы избавиться от проблем ...."

и тепреь не знаешь как воспринимать: толи наука, толи шаманство.
помоему для ряда sin(0.1), sin(0.2) проще составить таблицу значений и осуществлять выборку, чем заниматься подобного рода вычислений. А где есть матлаб, сопроцессор, и т.д. там зачастую проще вычислить
+1
… вычислить синус.
0
Простите пожалуйста, но в каком конкретном месте оно сложно ДЛЯ ПРОЦЕССОРА?

Вот здесь, где нет ни одного умножения (A*T = константа) и это вычисляется ОДНОКРАТНО?
Ad = eye(2) + (A * T);

Или вот здесь, где несколько умножений за шаг и нет ни одного деления (x/2 = сдвиг вправо)?
x = Ad * x;
x *= (3 — (x' * x)) / 2;

При этом, в этом ЧАСТНОМ случае, транспонирование делается «в лоб», с импользованием смещений. Прямо в операции перемножения матриц. Оно даже не идёт за отдельную операцию.
0
Если я правильно понимаю, x' * x — получение квадрата длины вектора. В обычных языках программирования это обычно несколько иначе делается. Так что транспонирование и не нужно.
0
слишком много вычислений, что является не самым оптимальным вариантом.
Вот здесь всё намного проще сделано Генерация DTMF или любого синуса
0
нормальная статья. идея тож неплохая.
а то что некоторые не знают что матлаб удобен для моделирования и на МК его никуда ставить не надо — не вина автора. а код для МК можете сами написать на любом яп
+1
  • avatar
  • kest
  • 16 октября 2012, 08:52
Интересно конечно, но я сам преподаю матлаб в универе, и пару вещей не понял, в частности про octave и про expm (вернее отличие от «просто» exp). Такие вещи в статье конечно надо было описать.
0
вот вроде читаем эту статью не в газете… до соседнего таба в броузере дотянуться — ровно 5 мм… гугл там открыть — пару секунд… делаем это ежедневно и помногу… википедию все знаем…

Вот никак не могу понять — как такое может являться ЗНАЧИМЫМ препятствием для человека с ТЕХНИЧЕСКИМ складом ума?! Любопытство и пытливость ума уже отменены? Про октаву узнать из первых строк статьи вики — несколько секунд. Сравнить в матлабе help по exp/expm — не больше минуты.

Мне немного странно, что автор должен был описать каждую используюмую функцию. Инструмент указан, имена даны — всё есть.
0
Сравнить в матлабе help по exp/expm — не больше минуты.
Который предварительно надо найти, скачать и поставить.
0
Это был ответ человеку, который ПРЕПОДАЁТ матлаб… Не находите, что сам матлаб у него уже скорее всего есть?
А скачать октаву (когда от гугла уже известно что это матлаб) — ещё проще.

А если погуглить три простых слова «octave exp expm» — в ПЕРВОЙ ЖЕ строке находится искомое. Где сразу же объясняется что это «Compute the exponential OF A MATRIX.»

Я уже не говорю о том, что такую разницу exp/expm ПРЕПОДАВАТЕЛЮ неплохо бы знать. Ну ок — всё знать нельзя, но догадаться уж можно.
0
По вашей логике, чтобы понять вашу статью, мне надо еще полчаса рыть в интернете. Тогда в чем ценность вашей статьи? если она не дает полное представление об ее основной идее?
0
1) статья НЕ МОЯ
2) даже с учётом того, что она не моя — у меня не возникло таких АЗБУЧНЫХ вопросов, как у ПРЕПОДАВАТЕЛЯ
3) мне более чем странно, что найденное мною по строке «octave expm» в гугле ЗА ПЯТЬ СЕКУНД, у Вас занимает пол-часа упорных поисков
4) ещё более странным мне кажется Ваше обобщение о том, что если статья непонятна ЛИЧНО ВАМ, то она не имеет общей ценности

Она даёт очень даже хорошее представление об основной идее. Но всё это она даёт лишь тому, кто ХОТЯ БЫ имеет представление об соответствующих инструментах.

По Вашей логике выходит, что точно так же, человек не знакомый с Verilog, может заявить о том, что статьи на тему FPGA не имеют никакой ценности тут. А человек, не знакомый с микроконтроллерами, может заявить о том, что всё, что сложнее дискретной логики — не имеет здесь ценности.

Это всего лишь вопрос УРОВНЯ знаний — Ваших, моих, прочих коллег. Но что Важнее — это гораздо больше вопрос того, насколько человек ГОТОВ эти знания получить.

Нестрашно не знать, страшно не стремиться знать.
+1
Вопрос в тематике сайта на котором размещается статья. Если работа в MatLab и тому подобных пакетах не является основной тематикой сайта, то при написании статьи необходимо, на мой взгляд, хотя бы написать комментарии к коду, из уважения к читающим. Я понимаю, если бы на сайте Exponenta разместить такую статью, то конечно можно кидать код и ничего не объяснять. Но если я там размещу пример программирования FPGA на Verilog, без объяснения кода и начну говорить, что каждый человек с инженерным образованием при помощи гугла сам разберется, что написано в коде, главное идея. Как вы думаете, какая будет реакция?
+2
Я понял Вашу мысль, отчасти это справедливо. И я бы с Вами согласился, если бы не одно «но».
Не электроника является прикладным инструментом в математических работах, а математика в элеткронных.

С другой стороны — тематика сообщества, это не нечто закостенелое, типа «паяльник, транзистор и тестер». Время идёт, электроника развивается, электроника доступная любителям — тоже. Чтобы не остаться с динозаврами в прошлом, нужно постоянно осваивать новые знания.

На мой сугубо личный взгляд, в современной любительской электронике, решающей интересные задачи, без матлаба уже трудно. Ради любопытства — посмотрите результат этого поиска we.easyelectronics.ru/search/topics/?q=dsp

Пусть DI меня поправит, но несмотря на название EASYelectronics, здесь есть место не только начинающим. И по-моему это прекрасно.

Сегодня про матлаб не знает только ленивый. И слышать апелляции к тому, что это сложно и непонятно — очень странно. Потому что когда лично мне сложно и непонятно — я воспринимаю это, как новую ступеньку для освоения. Начинаю разбираться и осваивать. А не жаловаться.
0
Потому что когда лично мне сложно и непонятно — я воспринимаю это, как новую ступеньку для освоения. Начинаю разбираться и осваивать. А не жаловаться.
Так давайте просто писать заголовки статей и делать ссылки на гугле. Не знаешь, что обозначает заголовок и начинаешь разбираться и осваивать. Тоже вариант :)
0
Вас не коробит пользоваться демагогическими приёмами в обсуждении? Доведение идеи до абсурда — один из них…
0
Нет. Вообще рассмотрение критических случаев хорошо позволяет понять некоторые вещи. Так же на предельных случаях объясняется много разных способов решения, вроде даже некоторых теорем. Так почему тогда я не могу использовать граничные случаи.
По поводу обязательных знаний по матлабу, я с вами не согласен.
Во первых, вы скинули ссылку, замечательно. Обратите внимание, что большая часть статей, по тому как именно сделать в матлаб, они обучающие. А в данном случае использование среды без объяснений — очень большая разница!
Во вторых, на сегодняшний день все среды разработки делаются таким способом, что бы весь цикл проектирования можно было провести в одной среде, не прибегая к помощи других программных продуктов. По тому же пути пытается идти и Матлаб. По-этому если я использую другую среду, которая мне позволяет реализовывать все функции, которые мне необходимы, зачем мне еще изучать Матлаб?
+1
Это всё замечательно, но только Ваше предложение оставить одни заголовки и ссылки не является критическим случаем по отношению к этой статье. Это упрощение до степени потери свойств. Нерелевантная абстракция. Во-первых в силу конкретности статьи, последовательности изложения материала и собственной авторской подачи идеи.

А во-вторых — в силу контекста. Который хоть и не предполагает обязательных знаний инструмента, но ВПОЛНЕ ОБОСНОВАННО предполагает его хорошую известность. Эта статья имеет право не быть обучающей ровно в силу того факта, что обучающих (и заметьте — непрофильных этому ресурсу) статей по матлабу — много. Эта же статья ценна именно в силу профильности.
И любые апелляции к неизвестности инструмента в этом контексте работают плохо. 2012-й год на дворе.

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

Поэтому, если Вы используете другую среду, и Вам непонятно написанное тут, и незачем изучать матлаб, то… То я всего лишь могу выразить самое искрене сожаление по поводу Вашего, глубоко частного, случая. Но глобальных выводов делать не стану.
0
Думаю, больше не стоит разводить базар :). Мысль я вашу понял. Чем больше знаешь, тем круче. Согласен.
0
_pv дал ниже ссылку на статью, которая рассказывает то же самое, но короче и понятнее.
0
Не сказал бы, что там изложено поятнее.
0
Там, кстати, сразу же ошибка в уравнении колебаний: в нем и X и Y одновременно. Опечатка, конечно…
0
и про начальные условия еще скромно умолчали, при x[-1] и x[-2] == 0 не заведётся.
вот более правильный аппноут:
www.ti.com/lit/an/spra096a/spra096a.pdf
ссылку на его перевод выше уже приводили:
0
Статья хорошая, но на сайте есть и школьники и студенты младших курсов. Если у меня, преподавателя есть «некоторые затруднения», то у них могут возникнуть «трудности». На данный ресурс я захожу, чтобы отдохнуть, пообщаться, узнать новое. Мой комментарий имеет рекомендательный характер автору (коим вы не являетесь) на будущее.

Хорошее должно стремиться к лучшему, а эта статья правда хорошая. Надеюсь понятно к чему я?
+1
:) Мне это чем то напомнило сцену из фильма О чем говорят мужчины :)

«А гренка в нашем ресторане называется croûton. Это точно такой же кусочек поджаренного хлеба, но гренка не может стоит 8 долларов, а croûton — может». А дальше ты начинаешь искать хоть какой-то вкус, отличающий этот крутон от гренки. И находишь!
0
0
Эээ… Я сегодня не выспавшийся неверное, а какое это имеет отношение к моему комментарию? Я просто не знаю, как реагировать :)
0
Если посмотреть начало видео, станет понятно. Пошли такие комментарии, что можно провести аналогию данной статьи с крутым рестораном.
0
Это вы к тому, что некоторые смотрят с высока на полтора десятка человек, которым непонятна пара слов из этой статьи? Или я опять не в теме? :)
0
Наверное так. Просто пришла такая ассоциация в голову, когда прочитал все коменты. :)
0
дожили =))) «послать в гугл на пять секунд» = «прослыть снобом»
0
Глубоко не вникал, но мне кажется, что вряд ли предложенный вариант будет проще, а главное — быстрее, чем табличный, со значениями на четверть периода… Особенно на контроллерах без аппаратного умножения.

А про интуитивную понятность табличного метода практически для всех — и упоминать как — то вроде бы незачем, и так ясно…

И табличный способ более универсален, — можно задать не только синус, но и любую другую форму, лишь сменив значения в таблице, не меняя самого алгоритма формирования выходного напряжения в программе, а не подстраиваться к каждому случаю заново, использую математические пакеты и «костыли»…
+2
  • avatar
  • SWG
  • 16 октября 2012, 10:05
IMHO, автор не предлагал это как «лучшую замену табличного». Он предложил ЕЩЁ ОДИН вполне интересный способ. Лично мне — приятно его знать теперь. Где-то этот способ может оказаться выгоднее, где-то нет. Это нормально.

Простой пример — если мне понадобится аппаратный синус в FPGA, то метод автора запросто может стать более выигрышным.
0
На МК — пожалуй, а вот в статье по ссылке ниже, ориентирующейся на DSP, табличный метод получается одним из худших, уступая по производительности (и остальным параметрам тоже) даже полиномиальному (через ряд Тейлора).
0
Не понимаю, почему все так прицепились к синтаксису а ля Матлаб? Все операции вполне очевидны из текста поста. Нет, чтобы сказать, какой автор молодец, оформил формулки в ТеХе… Замечания автору: 1) зря не сказали, что исходная система первого порядка эквивалентна уравнению колебаний (второго порядка); 2) можно было бы упомянуть, что компоненты вектора икс — ортогональные sin и cos; 3) немного покоробила фраза про определение через ряд — из нее можно вынести, что и скалярная экспонента вводится также; 4) формулу изменения частоты, имхо, надо подправить.
Есть и вопрос! :) Что известно про соотношение частоты дискретизации (шаг по времени) и желаемой частоты гармоники, когда начинает сильно проявляться ошибка разложения матричной экспоненты в два первых члена рада? В плане отклонения от желаемой частоты. И через сколько периодов проявится — ошибка-то будет накапливаться. Хотя, я так думаю, что эту проблему можно спрятать за правильным выбором сетки: чтобы N шагов по времени точно укладывались на периоде. Его же мы знаем изначально?
Лично мне, правда не специалисту в программировании и электронике, было очень интересно прочитать о таком способе. И т.н. костыль, вообще говоря, костылем и не является.
+1
Не откажите в любезности пожалуйста — не могли бы Вы развернуть значительно подробнее пп.1,3 и 4? И про накопление ошибки разложения тоже… Буду искренне признателен. Спасибо.
0
Отчего же отказать? Пожалуйста!

1) Автор топика первой формулой приврдит систему ОДУ первого порядка:

Простой подстановкой одного в другое получаем одно ОДУ, но уже второго порядка:

Это общее правило, что ОДУ N-ого порядка эквивалентен с некоторой системой из N ОДУ 1-ого порядка. Но полученное уравнение ни что иное, как обычное уравнение колебаний. К примеру, механического маятника с единичной частой. В общем же виде для циклческой частоты оно выглядит так:

Ну а честное решение этого уравнения:

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

2)

3) Касательно разложения матричной экспоненты в ряд:

по типу ряда Тейлора для обычной скалярной функции — на сколько мне известно, это именно определение функции от оператора (матрицы). То есть, если надо определить , то раскладываем скалярную около нуля и затем «иксы» заменяем на оператор (матрицу). А, вот, ряд Тейлора для скалярной функции — это не определение, а теорема. То есть, например, вводится не через ряд, а, ЕМНИП, обобщением понятия возведения в степень на множество действительных чисел.

Поэтому если «икс» — косинус, то автоматом получаем и синус при таком матричном подходе.

4) Про изменение частоты. Для единичной частоты матрица системы имеет вид, как показал автор:

Однако для случая произвольной частоты должно быть:

Чтобы получалось приведенное мною выше уравнение колебаний.
Автор же для изменения частоты предлагает такой код:
Ad = eye(2) + (A * T);

А надо изменять антидиагональные элементы (которые пока по модулю равны единичке). Сделать это можно умножением:
Ad = eye(2) + (A * T)*omega_0;

Ну, или прибавлять матрицу, у которой на диагонали нули, а на антидиагонали —
Вот, как-то так…

*) Что касается накопления ошибки: так дело в том, что когда матрчиную экспоненту аппроксимируем только первыми двумя членами, то получаем некоторую ошибку в частоте. Следовательно, если строим длинный во времени синус, то из-за наличия этой будет все больше проявлять сдвиг по времени двух гармоник: генерируемой численным методом и желаемой . Если же строить только первый период, а потом пользоваться периодичностью (что будет особенно удобно, если отсчетов по времени укладывается точно в период), эта ошибка будет нивелирована.

И не надо никакого Матлаба в контроллер запихивать! Что может быть проще линейной алгебры в действительных числах? Да еще в пространстве размерности 2.
0
Тогда от табличного метода эта метода отличается только тем, что четверть периода хранится во флеше в случае табличного метода, а в случае рекуррентной генерации мы эту четверть формируем алгоритмически. Да же знаю, это быстрее?
0
Ну, примерно… Только если Вам надо получить значения синуса с малым шагом, то и таблица будет расти. Признаться, не знаю, какова каноническая реализации такой таблицы… Например в свей программе на CUDA я делал её из «половинных углов». То есть хранил значения синуса и косинуса от a, 2a, 4a, 8a,… до угла в 2*pi. А угол «а» — самый мелкий, определяемый шагом. А дальше по формулам типа «синус от суммы углов» получал нужное мне значение функции. Правильно подобрав сетку по отношению к периоду, достаточно анализировать каждый бит значения угла, от которого надо найти функцию, отдельно. Битик в «1» — прибавляем n*a, в ноль — не прибавляем.
0
Только все эти аппроксимации иногда дают такую вычислительную нагрузку, что полиномиальный вариант оказывается быстрее.
0
Я думаю, что, как это часто бывает, для каждой конкретной задачи подходит своя вариация решения :) Я алгоритмов специально не ботал — придумал велосипед. Но ездит у меня достаточно быстро. Потому, что квадратность его колес очень хорошо повторяет рельеф дороги :) Мне нужно постоянно генерировать синус и косинус на одной и той же сетке. Подогнав параметры сетки, под требуемые значения периодов я вообще избавился от каких бы то ни было аппроксимаций с момента вычисления cos и sin в узлах сетки, которые лежат в таблице. А этим вычислением занимается библиотечная функция из C++, которая оказалась весьма точнее полиномиального разложения в штук 7-8 членов! Дальше работают только тригонометрические формулы, а их точность ограничивается лишь точностью машинной арифметики.
0
Потому, что квадратность его колес очень хорошо повторяет рельеф дороги :)
Хорошая фраза)
0
Исправил. За исключением введения частоты, мне кажется так проще, хотя и не правильно с какой-то стороны.

Если думать про точность частоты, я прихожу либо обратно к вычислению sin() либо к медленной подстройке под опорный сигнал.
0
Прошу прощения! Формула с v(t) и последнее предложение в п.3 убежало из п.2 :(
0
а если решать дифф уравнение второго порядка Y'' = -w^2 * X всё получается гораздо красиввее:
Y[i+1] = K*Y[i] — Y[i-1], где K = 2 * cos(2*PI*F),
то есть всего одно умножение и вычитание для вычисления очередного значения синуса и никаких проблем со стабильностью, K считается один раз заранее для выбранной частоты.
этот осциллятор используется в алгоритме Герцеля.

www.ied.com/~petr/Oscillators.html
+2
  • avatar
  • _pv
  • 16 октября 2012, 12:56
+1
мега-спасибо =)
0
Не зря я думал, что можно лучше.
0
Только ничего не придумал, и поискать не догадался :)
0
А по делу. Как я пока понял, стабильность здесь такая же как и с поворотом, только в целых числах процесс не расходится из-за немного неточно посчитанной матрица Ad т.к. здесь нет такой избыточности коэффициентов.
0
Со стабильностью понятно, но мне интересно, не могли бы владельцы аналога Matlab'а взять табличный метод генерации, потом рекуррентный второго порядка и, скажем, предъявить разность между выходами для какого-то t после непрерывной генерации n=10, n=100, n=1000, n=10000 периодов синусоиды? Что касается табличного метода, то тут у меня сомнений нет, что будет то одно и тоже, а вот для второго случая я чуток сомневаюсь.
0
  • avatar
  • uni
  • 16 октября 2012, 14:17
Для разрядности: 8, 16, 32 естественно. Думаю, что это можно в octave организовать.
0
По фазе может уползти, особенно если Ad (или K в методе _pv) считать приближенно. Но можно синхронизироваться с опорой, если это важно.
0
Между прочим, в цифровой технике всё приближенно, иногда даже до 8 бит в целом исчислении. Поэтому для этих случаев в Matlab'е добавили специальные целочисленные типы, чтобы адекватно рассматривать задачу не в double точности, которая всегда будет давать симпатичный результат.

Надо бы автору изображать не православый отрезок [0, 1], а более типичный, местный [0, 255].
+1
Да, кстати, если вернуться к школьному курсу, а не использовать диффуры, то, общем, ясно, что синус (и косинус сразу) получается элементарным повотором вектора (1,0) вокруг центра координат. Используя комплексную запись для поворота вектора на элементарный угол, известную формулу Эйлера, мы и получаем эту самую экспоненту.

А вообще, в курсе ЦОС проходится такая штука как метод БПФ. Там же вводится такое понятие как «поворачивающий множитель» (может кто помнит). Так вот это и есть эта экспонента, она же матрица поворота для тех, кто изучал работу с графикой в 2D.
0
  • avatar
  • uni
  • 16 октября 2012, 15:21
А вот так не проще?

Выводим рекурсивный фильтр на границу устойчивости, подаем на вход импульс для запуска — и получаем цифровой генератор.
0
Есть вариант этого метода, в котором потеря точности за счет итераций намного меньше.
Синус и косинус суммы углов:
cos(θ+δ) = cosθ - [αcosθ + βsinθ]
sin(θ+δ) = sinθ - [αsinθ - βcosθ]

Здесь α≡2sin²(δ/2), β≡sinδ.

Начиная с θ₀=0, s₀=0, с₀=1, получаем:
s' = s - [αs - βc]
c' = c - [αc + βs]

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

При вычислениях с плавающей точкой эта схема дает ошибку в нескольких младших битах за тысячи шагов и часто используется в пакетах для вычисления БПФ.
0
  • avatar
  • ztkl
  • 16 октября 2012, 18:20
В конце статьи хорошо бы добавить пример кода на каком-нибудь популярном языке программирования или псевдокоде. Я в институте матлаб не изучал, конструкция типа
xl(:,end+1) = x;
мне не особо понятна, не говоря уж о том, что такие листинги частично определяются как комментарии.
0
Проверил. Работает!!!
0
Небольшой ликбез из комплексной арифметики:

Вывод формулы

Пояснения по поводу поворачивающего множителя:

Поворачивающий множитель
0
У этого способа есть геометрическая интерпретация, x — вектор в двухмерном пространстве, Ad — матрица поворота на угол T в радианах.
Посоветуйте пожалуйста книжку по вышке, такую с геометрическими интерпретациями матричных вычислений и т.п. Конечно, матрицы в институте проходил и все нормально решал, но как-то проникнутся матрицами не получилось. Если например с интегралами или производными, ну или с рядами Тейлора и Фурье в голове отложилась интерпретация, а вот с матрицами никак.
0
  • avatar
  • mux
  • 22 октября 2012, 14:29
Нинай насчет книжек по вышке, но мне это понятно благодаря тому, что я работал с 3D графикой. Там как раз линейная алгебра в основе всего и книжки по 3D обычно содержат тот или иной ликбез по ней.
0
Хорошие ВУЗовские книги:
Ильин В.А., Позняк Э.Г., «Линейная алгебра»
Свешников А.Г. Тихонов А.Н., «Теория функций комплексной переменной»

Касательно матриц приложенных к графике, их «геометрический смысл» (типа матрицы поворота), наверно, можно почерпнуть из
Порев В.Н., «Компьютерная графика»

Все это это в djvu на просторах Сети.
0
Насколько я понял, статья интересна тем, что:
1) Ещё один способ генерирования синуса, достаточно медленный;
2) Который позволит табличному методу не хранить во флэше кучу (1/4 периода) значений синуса.
Тогда вопрос, для какого МК (ЦАП) будет приемлем данный алгоритм.
0
Для мощных числомолотилок или же там, где очень не хватает флеша. В плане первого — подобные методы способны на DSP и других процессорах с быстрым MAC обогнать даже табличный по скорости.
0
Я пришел к этому способу не от задачи «надо генерировать синус» а от «надо управлять вращающейся системой».

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

Но все измениться если взять сервомашинку без ограничений на количество полных оборотов. Надо будет добавлять «костыли» чтобы вращение происходило по кратчайшей траектории, т.к. топологии пространств теперь не совпадают. Если в действительности двигаться в одну сторону (вращать механизм сервы) то можно оказаться в исходной точке. Но если делать движение «виртуально» прибавляя к угловой координате приращение то в исходную точку попасть нельзя.

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

Примерно такие же дела в контроллере BLDC, если связываться с синусоидальным управлением.
0
Синус за одно умножение и вычитание
0
Не очень точно выходит. Пара синус-косинус гораздо точнее
0
Ну зато амплитуда стабильная, да и частота довольно точная. Только небольшая ошибка по фазе. А при генерации вращением вектора амплитуду нужно стабилизировать.
0
Кстати, нашел причину фазовой ошибки. У тебя сгенерированный синус имеет 0 в точке x[1], а образцовый — x[0]. Если исправить инициализацию на x[0]:=0; x[1]:=sin(δ) (либо x[1]=δ), то графики полностью совпадают.
0
Было же уже. Правда тут наглядно. Заодно вбил в маткад, поигрался — в частности, мне было интересно, насколько оно критично к начальным условиям. Как выяснилось, они влияют исключительно на амплитуду получающегося сигнала. Ну и при X0:=-δ оно вполне неплохо работает.
0
Голова! А я уж было забраковал такой синус.
Можно еще добавить модуляцию (здесь разность коэффициэнтов 0.01 приводит к возникновению биений из-за разности частот)
0
вот из за таких «ценителей прекрасного» уже и не хочется писать статей для сообщества. То octave не торт, то таблицы все же лучше, а автор почти имбецил и проч. Автору респект! Даже за донесение знаний об Matlab и Octave в массы (сам пользуюсь сабжем очень часто)
0
Спасибо и за математику и за Октаву. А то такая усталость от стандартных, общепринятых зазубренных методов. Без математики вообще труба, когда долго нет задачи где нужно соображать, когда всё долго слишком просто и не надо учиться вроде — потом наступает момент истины.
0
Я вот хоть убейте, не пойму, зачем, было с матрицами лезть в Октаву???
Ведь приведённое автором матричное уравнение легко решается в уме и приводит к БИХ фильтру вида Н(z)=1/(1+z^{-2}), т.е. всё сводится к стандартным рекуррентным формулам БИХ -фильтра, которые уже были здесь, к счастью, приведены.
Автор попугал всех дифференциальными уравнениями, матрицами, но в конце концов все его рассуждения свелись к элементарной формуле, приведённой AlexX1810: D[n]=X[n]+aD[n-1]-D[n-2]
0
Я этот алгоритм пытался в ПЛИС реализовать. Да, работает, но на длительном периоде (десятки секунд при 19000 Гц частоте синуса и 1 МГц дискретизации) рекурсивный алгоритм «разваливается», сигнал просто исчезает. Нужно после целого большого числа периодов перезагружать исходными значениями.
0
Который конкретно алгоритм? Тут их несколько.
0
Есть вариант этого метода, в котором потеря точности за счет итераций намного меньше.
Синус и косинус суммы углов:
cos(θ+δ) = cosθ — [αcosθ + βsinθ]
sin(θ+δ) = sinθ — [αsinθ — βcosθ]

Здесь α≡2sin²(δ/2), β≡sinδ.

Начиная с θ₀=0, s₀=0, с₀=1, получаем:
s' = s — [αs — βc]
c' = c — [αc + βs]

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

В таком виде:
module gen
(input clk, input init,
output wire signed [15:0] data);

parameter signed [31:0] cosw=32'd2147478986;
parameter signed [31:0] sinw=32'd8948915;
reg signed [63:0] cost;
reg signed [63:0] sint;
reg signed [63:0] cc;
reg signed [63:0] ss;
reg signed [63:0] cs;
reg signed [63:0] sc;
reg signed [63:0] tmp;
reg signed [63:0] tmp2;

initial
begin
cost<=32'h7FFFFFFF;
sint<=32'h0;
end

always @(posedge clk)
if (init)
begin
cost<=32'h7FFFFFFF;
sint<=32'h0;
end
else
begin
cc<= cost * cosw;
ss<= sint * sinw;
cs<= cost * sinw;
sc<= sint * cosw;
tmp<=(cc — ss);
tmp2<=(cs + sc);
cost<=$signed(tmp[62:31]);
sint<=$signed(tmp2[62:31]);
// cost <= $signed({cost[63:32]});
// sint <= $signed({sint[63:32]});
end

assign data=$signed(sint[31:16]);

endmodule
0
Только зарегистрированные и авторизованные пользователи могут оставлять комментарии.