#помощь — Синусоида в C++


Содержание

Описание функций языка Си

All | _ | A | B | C | D | E | F | G | H | I | J | K | L | M | N | O | P | Q | R | S | T | U | V | W | X | Y | Z

sin, sinf, sinl – расчет синуса

double sin (double a);
float sinf (float a);
long double sinl (long double a);

a – значение угла в радианах.

Синус угла, заданного аргументом.

Функции рассчитывают значение синуса угла, заданного в радианах.

Аргумент и возвращаемое значение функции sin являются числами с плавающей точкой двойной точности (тип double, точность не менее десяти значащих десятичных цифр, разрядность — 64).

Аргумент и возвращаемое значение функции sinf являются числами с плавающей точкой (тип float, точность не менее шести значащих десятичных цифр, разрядность — 32).

Аргумент и возвращаемое значение функции sinl являются числами с плавающей точкой повышенной точности (тип long double, точность не менее десяти значащих десятичных цифр, разрядность — 80).

В примере рассчитывается синус от 1.23 радиан с помощью функций sin, sinf и sinl, и результат выводится на консоль. Обратите внимание на точность полученных результатов. У синуса, рассчитанного с помощью функции sinf, будет самая маленькая точность, а у рассчитанного с помощью функции sinl – самая большая.

Аргумент: 1.23 рад.
sinf : 0.94248878955841064453
sin : 0.94248880193169748409
sinl : 0.94248880193169750410

График функции sin(x)

Задача Для x∈[0;7] построить график функции y=sin(x).

График функции y=sin(x) симметричен относительно горизонтальной оси, поэтому должен быть построен в окне со смещением относительно левого верхнего угла окна, принятого за начало отсчета.

Координата X меняется в пределах [0;7] , координата Y — [-1;1] . Величины MAX_X и MAX_Y представляют собой область допустимых значений по осям координат:

  • MAX_X = maxX — minX = 7 ;
  • MAX_Y = maxY — minY = 2 .

Смещение графика функции представляет собой положение первой точки относительно начала отсчета левого верхнего угла:

  • смещение по оси X : OffsetX = minX*w >;
  • смещение по оси Y : OffsetY = maxY *height/MAX_Y

Значение minX=0 представляет собой минимальное значение координаты x из области допустимых значений.

Значение maxY=1 представляет собой максимальное значение координаты y из области допустимых значений. Рассматривается именно максимальное значение, поскольку за начало отсчета принят левый верхний угол окна, и координата y увеличивается по направлению вниз.

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

  • масштабный коэффициент X : ScaleX = width / MAX_X ;
  • масштабный коэффициент Y : ScaleY = height / MAX_Y .

Вычисление координат следующей точки (x;y) графика в окне будет осуществляться по формулам:

  • координата X=OffsetX + x*ScaleX ;
  • координата Y=OffsetY + y*ScaleY ,

где (x;y) — координаты точки, полученные из функции y=sin(x) .

Реализация на C++

Примечание : Для корректной сборки приложения используется многобайтовая кодировка.

Документация

Сгенерируйте синусоиду, с помощью внешнего сигнала в качестве источника времени

Simulink / Математические операции

Описание

Функциональный блок Синусоиды выводит синусоидальную форму волны. Блок может действовать в основанном на времени или основанном на выборке режиме.

Примечание

Этот блок совпадает с блоком Sine Wave , который появляется в библиотеке Math Operations. Если вы выбираете Use simulation time для параметра Time в диалоговом окне блока, вы получаете Функциональный блок Синусоиды .

Основанный на времени режим


Блок вычисляет выходную форму волны.

В основанном на времени режиме значение параметра Sample time определяет, действует ли блок в непрерывном режиме или дискретном режиме.

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

>0 заставляет блок действовать в дискретном режиме.

Для получения дополнительной информации см. Настройку времени выборки.

Блокируйте поведение в непрерывном режиме

При работе в непрерывном режиме блок Sine Wave может стать неточным из-за потери точности, как время становится очень большим.

Блокируйте поведение в дискретном режиме

Значение параметров Sample time, больше, чем нуль, заставляет блок вести себя, как будто это управляло блоком Zero-Order Hold , шаг расчета которого установлен в то значение.

Таким образом, можно создать модели с источниками синусоиды, которые чисто дискретны, а не модели, которые являются гибридными непрерывными/дискретными системами. Гибридные системы являются по сути более комплексными и в результате занимают больше времени, чтобы моделировать.

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

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

грех ( t +Δt ) =sin (t) cos (Δt) +sin (Δt) cos (t) cos ( t +Δt ) =cos (t) cos (Δt) −sin (t) грех (Δt)

В матричной форме эти тождества:

[ грех ( t +Δt ) cos ( t +Δt ) ] = [ cos (Δt) грех (Δt) −sin (Δt) cos (Δt) ] [ грех (t) cos (t) ]

Поскольку Δ t является постоянным, следующее выражение является константой:

[ cos (Δt) грех (Δt) −sin (Δt) cos (Δt) ]

Поэтому проблема становится одним из умножения матриц значения греха (t) постоянной матрицей, чтобы получить грех ( t +Δt ) .

Дискретный режим уменьшает, но не устраняет накопление ошибок округления, например, (4*eps) . Это накопление может произойти, потому что вычисление блока вывод на каждом временном шаге зависит от значения вывода на предыдущем временном шаге.

Методы, чтобы обработать ошибки округления в дискретном режиме

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

Вставьте блок Saturation , непосредственно нисходящий из блока Sine Wave.

Путем установления пределов насыщенности для блока Sine Wave вывод можно удалить проскакивание из-за накопления ошибок округления.

Настройте блок Sine Wave, чтобы использовать математическую библиотечную функцию sin() , чтобы вычислить, блокируют вывод.

На диалоговом окне блока Sine Wave, набор Time к Use external signal так, чтобы входной порт появился на значке блока.

Соедините сигнал часов с этим входным портом с помощью блока Digital Clock .

Установите шаг расчета сигнала часов к шагу расчета блока Sine Wave.

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

Основанный на выборке режим

Основанный на выборке режим использует эту формулу, чтобы вычислить вывод блока Sine Wave .

y=Asin ( 2π (k+o) /p ) +b

A является амплитудой синусоиды.

p является количеством выборок времени на период синусоиды.

k является повторяющимся целочисленным значением, которое колеблется от 0 до p –1.

o является смещением (сдвиг фазы) сигнала.

b является смещением сигнала.

В этом режиме Simulink ® устанавливает k, равный 0 на первом временном шаге, и вычисляет блок вывод, с помощью формулы. На следующем временном шаге Simulink постепенно увеличивает k и повторно вычисляет вывод блока. Когда k достигает p, Simulink сбрасывает k к 0 прежде, чем вычислить блок вывод. Этот процесс продолжается до конца симуляции.

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


математика — Самая быстрая реализация синуса, косинуса и квадратного корня в C ++ (не очень точная)

Я гуглю вопрос в течение прошедшего часа, но есть только указания на ряд Тейлора или некоторый пример кода, который либо слишком медленный, либо не компилируется вообще. Ну, большинство ответов, которые я нашел в Google: «Google, он уже задан», но, к сожалению, это не

Я профилирую свою игру на недорогом Pentium 4 и обнаружил, что

85% времени выполнения тратится на вычисление синуса, косинуса и квадратного корня (из стандартной библиотеки C ++ в Visual Studio), и это, по-видимому, сильно зависит от процессора ( на моем I7 те же функции получили только 5% времени выполнения, а игра waaaaaaaaaay Быстрее). Я не могу ни оптимизировать эти три функции, ни вычислить синус и косинус за один проход (они взаимозависимы), но мне не нужны слишком точные результаты для симуляции, поэтому я могу жить с более быстрым приближением.

Итак, вопрос: Какой самый быстрый способ вычислить синус, косинус и квадратный корень для числа с плавающей точкой в ​​C ++?

РЕДАКТИРОВАТЬ
Таблица поиска более болезненна, так как в результате Cache Miss на современных процессорах обходится гораздо дороже, чем Taylor Series. В наши дни процессоры работают очень быстро, а кеш — нет.

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

Итак, обновленный вопрос: есть ли быстрая оптимизация для квадратного корня?

EDIT2

Я использую квадратный корень для вычисления расстояния, а не нормализации — не могу использовать быстрый алгоритм обратного квадратного корня (как указано в комментарии: http://en.wikipedia.org/wiki/Fast_inverse_square_root

EDIT3

Я также не могу работать на квадрате расстояний, мне нужно точное расстояние для расчетов

Решение

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

НО, если вы настаиваете на вычислениях во время выполнения, вы можете использовать ряд синус или косинус Тейлора …

Для получения дополнительной информации о серии Тейлор … http://en.wikipedia.org/wiki/Taylor_series

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

Кроме того … Не умножайте свой x ^ n с начала каждый раз … например. умножьте x ^ 3 на x еще два раза, а затем еще на два, чтобы вычислить показатели.

Другие решения

Во-первых, серия Taylor НЕ является лучшим / самым быстрым способом реализации синуса / cos. Это также не то, как профессиональные библиотеки реализуют эти тригонометрические функции, а знание наилучшей числовой реализации позволяет настроить точность, чтобы повысить скорость более эффективно. Кроме того, эта проблема уже широко обсуждалась в StackOverflow. Вот только один пример .

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

Наконец, давайте поговорим о коде на вашем старом ПК. Проверьте научная библиотека gsl gnu (или числовые рецепты) реализации, и вы увидите, что они в основном используют формулу аппроксимации Чебышева.

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

Кстати: правило 0 для такого рода проблем: если вам нужны подробности реализации какой-либо специальной функции или числового метода, вам следует взглянуть на код GSL перед любыми дальнейшими действиями — GSL является СТАНДАРТНОЙ числовой библиотекой.

РЕДАКТИРОВАТЬ: вы можете улучшить время выполнения, включив агрессивные флаги оптимизации с плавающей запятой в gcc / icc. Это снизит точность, но кажется, что это именно то, что вы хотите.

РЕДАКТИРОВАТЬ 2: Вы можете попытаться создать сетку грубых грехов и использовать процедуру gsl (gsl_interp_cspline_periodic для сплайна с периодическими условиями), чтобы сплайновать эту таблицу (сплайн уменьшит ошибки по сравнению с линейной интерполяцией => вам нужно меньше точек на вашей таблице = > меньше кеша не хватает)!

Вот гарантированная максимально быстрая функция синуса в C ++:

О, вы хотели лучшую точность, чем | 1.0 |? Хорошо читайте дальше.

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

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

Большинство людей, которые задают этот вопрос, не понимают важности компромисс между производительностью и точностью. В частности, вам придется сделать некоторые выборы в отношении точности вычислений, прежде чем делать что-либо еще. Сколько ошибок вы можете терпеть в результате? 10 ^ -4? 10 ^ -16?

Если вы не можете измерить ошибку любым методом, не используйте ее.

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

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

За эти годы было много постепенных улучшений оригинальных алгоритмов CORDIC. Например, в прошлом году некоторые исследователи в Японии опубликовали статья на улучшенном CORDIC с лучшими углами поворота, что уменьшает требуемые операции.

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

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

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

Квадратные корни, как и сейчас, обычно рассчитываются по некоторому варианту Алгоритм Ньютона-Рафсона , обычно с фиксированным числом итераций. Обычно, когда кто-то провернул «удивительно новое» Алгоритм создания квадратных корней, это просто замаскированный Ньютон-Рафсон.


Полиномы Ньютона-Рафсона, CORDIC и Чебышева позволяют вам поменять скорость на точность, поэтому ответ может быть столь же неточным, как вы хотите.

Наконец, когда вы закончили все свои модные тесты и микрооптимизацию, убедитесь, что ваша «быстрая» версия на самом деле быстрее, чем версия библиотеки. Вот типичная реализация библиотеки fsin () ограничен доменом от -pi / 4 до pi / 4. И это чертовски медленно.

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

ТЛ: д-р; зайти в Google «аппроксимация по синусу» или «аппроксимация по косинусу» или «аппроксимация квадратного корня» или «теория приближений .»

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

Число с плавающей запятой, определенное IEEE-754, использует некоторый определенный бит, представляющий времена описания множественных основанных 2. Некоторые биты предназначены для представления базового значения.

Это постоянное время расчета квадратного корня

Основано на идее http://forum.devmaster.net/t/fast-and-accurate-sine-cosine/9648 и некоторое ручное переписывание для улучшения производительности в микро-бенчмарке. В итоге я получил следующую косинусную реализацию, которая используется в физическом моделировании высокопроизводительных вычислений, что является узким местом при повторных вызовах из-за большого числа номеров. Это достаточно точно и намного быстрее, чем справочная таблица, в особенности деление не требуется.

Компилятор Intel по крайней мере также достаточно умен для векторизации этой функции при использовании в цикле.

Если определено EXTRA_PRECISION, максимальная ошибка составляет около 0,00109 для диапазона от -π до π, при условии, что T является double как это обычно определяется в большинстве реализаций C ++. В противном случае максимальная ошибка составляет около 0,056 для того же диапазона.

Для x86 аппаратные инструкции FP с квадратным корнем выполняются быстро ( sqrtss является скалярной одинарной точностью). Одинарная точность быстрее, чем двойная, поэтому обязательно используйте float вместо double для кода, где вы можете позволить себе использовать меньшую точность.

Для 32-битного кода вам обычно нужны опции компилятора, чтобы заставить его выполнять математику FP с инструкциями SSE, а не x87. (например. -mfpmath=sse )

Чтобы получить C sqrt() или же sqrtf() функции, чтобы встроить как раз sqrtsd или же sqrtss , вам нужно скомпилировать с -fno-math-errno . Наличие математических функций errno на NaN обычно считается ошибкой проектирования, но стандарт требует этого. Без этой опции gcc вставляет ее, но затем выполняет сравнение + ветвь, чтобы увидеть, был ли результат равен NaN, и, если это так, вызывает библиотечную функцию, чтобы он мог установить errno , Если ваша программа не проверяет errno после математических функций нет опасности в использовании -fno-math-errno ,

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

код gcc для случая NaN кажется слишком сложным; он даже не использует sqrtf возвращаемое значение! Во всяком случае, это тот беспорядок, который вы на самом деле получите без -fno-math-errno , для каждого sqrtf() Звоните на сайт в вашей программе. В основном это просто код-блат, и никто из .L8 Блок когда-либо запускается, если взять sqrt числа> = 0.0, но в быстром пути еще есть несколько дополнительных инструкций.

Если вы знаете, что ваш вклад в sqrt не равен нулю, вы можете использовать быструю, но очень приблизительную обратную инструкцию sqrt, rsqrtps (или же rsqrtss для скалярной версии). Один Итерация Ньютона-Рафсона обеспечивает почти ту же точность, что и аппаратная одинарная точность sqrt инструкция, но не совсем.

sqrt(x) = x * 1/sqrt(x) , за x!=0 , так что вы можете вычислить sqrt с помощью rsqrt и умножения. Они оба быстрые, даже на P4 (было ли это актуально в 2013 году)?

На P4, возможно, стоит использовать rsqrt + Ньютоновская итерация для замены одного sqrt даже если вам не нужно ничего делить на это.

Смотрите также ответ, который я недавно написал об обработке нулей при расчете sqrt(x) как x*rsqrt(x) , с итерацией Ньютона. Я включил некоторые обсуждения ошибки округления, если вы хотите преобразовать значение FP в целое число, и ссылки на другие соответствующие вопросы.

  • sqrtss : Задержка 23c, без конвейера
  • sqrtsd : Задержка 38c, без конвейера
  • fsqrt (x87): задержка 43 с, без конвейера

rsqrtss / mulss : Задержка 4с + 6с. Возможно один на 3c пропускной способности, так как они, очевидно, не нуждаются в одном и том же модуле выполнения (mmx против fp).

SIMD упакованные версии несколько медленнее

  • sqrtss / sqrtps : Задержка 12c, один на пропускную способность 3c
  • sqrtsd / sqrtpd : Задержка 15-16c, по одному на пропускную способность 4-6c
  • fsqrt (x87): задержка 14-21cc, по одному на пропускную способность 4-7c
  • rsqrtss / mulss : Задержка 4с + 4с. Один на 1с пропускной способности.
  • Векторные версии SIMD 128b имеют одинаковую скорость. Векторные версии 256b имеют немного большую задержку, почти половину пропускной способности. rsqrtss версия имеет полную производительность для 256b векторов.

С итерацией Ньютона rsqrt Версия не намного, если вообще быстрее.

Номера от Экспериментальные испытания Агнера Фога . Посмотрите его руководства по микроархам, чтобы понять, что делает код быстрым или медленным Также смотрите ссылки на x86 тег вики.

ИДК, как лучше всего рассчитать грех / cos. Я читал, что оборудование fsin / fcos (и только чуть медленнее fsincos это делает оба сразу) не самый быстрый способ, но ИДК, что есть.

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

Алгоритм предварительно вычисляет значения sin и cos для 0.5 ^ i, i = 1, …, N. Затем мы можем объединить эти предварительно вычисленные значения, чтобы вычислить sin и cos для любого угла вплоть до разрешения 0,5 ^ N

Это реализация метода ряда Тейлора, ранее приведенного в ответ Акеллехе .

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

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

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

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

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

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

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

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


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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Добавка

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

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

p.s. пока сохранял в черновики, получил в ответ «Hacking attempt».

  • easymath,
  • генератор
  • +11
  • 15 октября 2012, 20:42
  • amaora

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

  • GYUR22
  • 15 октября 2012, 22:55
  • DIHALT
  • 15 октября 2012, 23:23

  • Alfa
  • 15 октября 2012, 23:23

Если оно лично Вам неизвестно — значит это лично Вам ненужно (или не судьба). Бывает.
Но здесь достаточно людей, владеющих матлабом и октавой.

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

  • dipspb
  • 15 октября 2012, 23:30

За статью действительно спасибо. Но…

— это шо?
— это шо?
— мда…

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

  • Alfa
  • 15 октября 2012, 23:44

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

Можете ещё посмотреть g*gle-> «octave basics», там это все есть.

  • amaora
  • 15 октября 2012, 23:52
  • Vga
  • 16 октября 2012, 00:03
  • xeccentricx
  • 16 октября 2012, 00:07

  • Vga
  • 16 октября 2012, 00:20
  • dipspb
  • 16 октября 2012, 00:30
  • zombie
  • 16 октября 2012, 14:16

Пардон муа, но эти несколько строк — даже без бумажки, в уме можно перевести. Ну ок — с бумажкой.
Достаточно помнить как делается перемножение и транспонирование матриц.

При матрице 2×2 это вообще выходит код простой и небольшой.

  • dipspb
  • 16 октября 2012, 00:13
  • Vga
  • 16 октября 2012, 00:19
  • dipspb
  • 16 октября 2012, 00:28

  • Vga
  • 16 октября 2012, 00:32
  • dipspb
  • 16 октября 2012, 00:36
  • Vga
  • 16 октября 2012, 00:52
  • dipspb
  • 16 октября 2012, 01:01
  • Vga
  • 16 октября 2012, 01:06
  • Alatar
  • 16 октября 2012, 08:58

  • Vga
  • 16 октября 2012, 11:09
  • Vga
  • 16 октября 2012, 11:11

Статья так красиво начиналась: «… простые способы.… генерировать… без таблиц и… тяжелых вычислений.… на примере sin(0.1), sin(0.2), и тд. Для этого случая есть простой способ.»
А в итоге вылилось в сложную «сложную» математическую формулу, да еще с «… Теперь будем подставлять костыли, чтобы избавиться от проблем . «

и тепреь не знаешь как воспринимать: толи наука, толи шаманство.
помоему для ряда sin(0.1), sin(0.2) проще составить таблицу значений и осуществлять выборку, чем заниматься подобного рода вычислений. А где есть матлаб, сопроцессор, и т.д. там зачастую проще вычислить

  • evgeniy
  • 16 октября 2012, 05:14
  • evgeniy
  • 16 октября 2012, 05:16

Простите пожалуйста, но в каком конкретном месте оно сложно ДЛЯ ПРОЦЕССОРА?

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

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

При этом, в этом ЧАСТНОМ случае, транспонирование делается «в лоб», с импользованием смещений. Прямо в операции перемножения матриц. Оно даже не идёт за отдельную операцию.

  • dipspb
  • 16 октября 2012, 12:21
  • Vga
  • 16 октября 2012, 12:42

  • alekseyb
  • 16 октября 2012, 06:37
  • kest
  • 16 октября 2012, 08:52
  • geovas
  • 16 октября 2012, 09:55

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

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

Мне немного странно, что автор должен был описать каждую используюмую функцию. Инструмент указан, имена даны — всё есть.

  • dipspb
  • 16 октября 2012, 11:51

Сравнить в матлабе help по exp/expm — не больше минуты.

  • Vga
  • 16 октября 2012, 11:55

Это был ответ человеку, который ПРЕПОДАЁТ матлаб… Не находите, что сам матлаб у него уже скорее всего есть?
А скачать октаву (когда от гугла уже известно что это матлаб) — ещё проще.

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

Я уже не говорю о том, что такую разницу exp/expm ПРЕПОДАВАТЕЛЮ неплохо бы знать. Ну ок — всё знать нельзя, но догадаться уж можно.

  • dipspb
  • 16 октября 2012, 12:01
  • opolo84
  • 16 октября 2012, 13:01

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

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

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

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

Нестрашно не знать, страшно не стремиться знать.

  • dipspb
  • 16 октября 2012, 13:24
  • opolo84
  • 16 октября 2012, 13:43

Я понял Вашу мысль, отчасти это справедливо. И я бы с Вами согласился, если бы не одно «но».
Не электроника является прикладным инструментом в математических работах, а математика в элеткронных.

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

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

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

Сегодня про матлаб не знает только ленивый. И слышать апелляции к тому, что это сложно и непонятно — очень странно. Потому что когда лично мне сложно и непонятно — я воспринимаю это, как новую ступеньку для освоения. Начинаю разбираться и осваивать. А не жаловаться.

  • dipspb
  • 16 октября 2012, 13:59

Потому что когда лично мне сложно и непонятно — я воспринимаю это, как новую ступеньку для освоения. Начинаю разбираться и осваивать. А не жаловаться.

  • opolo84
  • 16 октября 2012, 14:15
  • dipspb
  • 16 октября 2012, 14:34
  • opolo84
  • 16 октября 2012, 14:43

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

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

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

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

  • dipspb
  • 16 октября 2012, 15:02
  • opolo84
  • 16 октября 2012, 15:14
  • Vga
  • 16 октября 2012, 13:50
  • dipspb
  • 16 октября 2012, 14:05
  • V_oron
  • 16 октября 2012, 14:37
  • _pv
  • 16 октября 2012, 14:47

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

Хорошее должно стремиться к лучшему, а эта статья правда хорошая. Надеюсь понятно к чему я?

  • geovas
  • 16 октября 2012, 13:52

:) Мне это чем то напомнило сцену из фильма О чем говорят мужчины :)

«А гренка в нашем ресторане называется croûton. Это точно такой же кусочек поджаренного хлеба, но гренка не может стоит 8 долларов, а croûton — может». А дальше ты начинаешь искать хоть какой-то вкус, отличающий этот крутон от гренки. И находишь!

  • opolo84
  • 16 октября 2012, 14:00
  • opolo84
  • 16 октября 2012, 14:07
  • geovas
  • 16 октября 2012, 14:12
  • opolo84
  • 16 октября 2012, 14:24
  • geovas
  • 16 октября 2012, 14:32
  • opolo84
  • 16 октября 2012, 14:35
  • dipspb
  • 16 октября 2012, 14:37

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

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

И табличный способ более универсален, — можно задать не только синус, но и любую другую форму, лишь сменив значения в таблице, не меняя самого алгоритма формирования выходного напряжения в программе, а не подстраиваться к каждому случаю заново, использую математические пакеты и «костыли»…

  • SWG
  • 16 октября 2012, 10:05

IMHO, автор не предлагал это как «лучшую замену табличного». Он предложил ЕЩЁ ОДИН вполне интересный способ. Лично мне — приятно его знать теперь. Где-то этот способ может оказаться выгоднее, где-то нет. Это нормально.

Простой пример — если мне понадобится аппаратный синус в FPGA, то метод автора запросто может стать более выигрышным.

  • dipspb
  • 16 октября 2012, 11:56
  • Vga
  • 16 октября 2012, 13:53
  • V_oron
  • 16 октября 2012, 12:34
  • dipspb
  • 16 октября 2012, 13:28

Отчего же отказать? Пожалуйста!

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

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

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

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

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

2)

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

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

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

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

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

Чтобы получалось приведенное мною выше уравнение колебаний.
Автор же для изменения частоты предлагает такой код:

А надо изменять антидиагональные элементы (которые пока по модулю равны единичке). Сделать это можно умножением:

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

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

И не надо никакого Матлаба в контроллер запихивать! Что может быть проще линейной алгебры в действительных числах? Да еще в пространстве размерности 2.

Задача «Вычисление sin(x)». Пример решения

Постановка задачи. Вычислить значение тригонометрической функции sin(x) от произвольного значения аргумента x.

Вариант решения 1. Поищем готовый вариант решения. Тригонометрия – раздел математики. Предположим, что в библиотеке System реализован класс, связанный с математическими функциями (в библиотеках прежних языков всегда была функция извлечения квадратного корня из числа – sqrt()). Как же узнать название класса?

Воспользуемся интеллектуальной подсказкой. Внутри метода Main() консольного приложения наберем одну букву M в английской раскладке. В выпадающем меню выберем класс Math (по-английски математика – Mathematics или maths). Наведя мышкой на это слово, прочитаем:
«class System.Math. Предоставляет константы и статические методы для тригонометрических, логарифмических и иных общих математических функций».
Мы почти у цели.

После слова Math ставим точку и снова выбираем что-либо, например метод Sin. Прочтем подсказку:
« double Math.Sin(double a). Возвращает синус указанного угла».
Выбрав Sin, откроем скобку « ( ». Снова появится подсказка:
«double Math.Sin(double a) Возвращает синус указанного угла. a: — Угол, измеряемый в радианах)».
Итак, если нам необходимо вычислить y=sin(x), где переменные x и y должны быть объявлены типом double, то это можно сделать одним оператором:
y = Math.Sin(x);

Если вы обратили внимание на третью подсказку, то вы знаете, или хотите узнать, что за единица измерения угла 1 радиан. В школе чаще всего вас приучили к измерению углов в градусах (прямой угол — 90°, развернутый — 180°, и т.д.), а при необходимости еще и в минутах и секундах.
Разработчики библиотеки не стали вводить еще один новый тип данных («градусный»), но обеспечили возможность для перевода градусов в радианы, добавив константу – число π (double Math.PI).
По определению, 1 радиан — центральный угол, длина дуги которого равна радиусу окружности, 2π радиан = 360°. Вспомните, что длина окружности радиусом R равна 2πR.
Составив пропорцию, вы сможете переводить угол a, заданный в градусах, в угол x, заданный в радианах: x = a· π /180.
Проверка: a=360°, x=2π; a=30°, x=0,523598775598299.

Программа, реализующая вычисление sin(a), угол а задан в градусах, представлена ниже:

Результат выполнения программы:

Вариант решения 2.

Понятно, что метод Sin( ) есть некоторая подпрограмма класса Math, которая вызывается по имени. Однако может быть вам будет интересно узнать, как она устроена. Для этого мы дополним наш проект методом Sin2(x), в котором самостоятельно реализуем некоторый алгоритм вычисления sin(x). Тогда библиотечную функцию Sin(x) мы сможем использовать для задания ожидаемого результата при тестировании, а обсуждая алгоритм вычисления Sin2(x), мы освоим элементы языка C# и приобретем некоторые навыки программирования.

Заглянем в справочник по высшей математике, раздел «Ряды». Из него мы узнаем, что функция sin(x) может быть представлена бесконечным рядом:
В теории, вы скажете все красиво, но бесконечный ряд – это бесконечное время вычислений, кроме того, возможно переполнение при возведении в степень и вычислении знаменателя.

Справка. Через n! обозначается функция ФАКТОРИАЛ — произведение целых чисел от 1 до n. По определению 0!=1. Например, 3!=1·2·3; 7!=1·2·…·6·7.

Отметим, что ряд знакопеременный: + , — , + , — , …
Предположим, что |x| 2 / ((n+2)(n+1)) > 1 (много больше 1), ведь в этом случае наш критерий работать не будет? Тут пригодится знание такого свойства функции sin(x) как периодичность. Для угла x, заданного в радианах:
sin(x+2πk) = sin(x), где k – любое целое число, 0, ±1, ±2, ±3, … .
Для угла а, заданного в градусах:
sin(a+360°·k) = sin(a), где k – любое целое число.
Таким образом, если угол х по модулю не будет больше 2π, |x| 2 / ((n+2)(n+1))
т.е., начиная с n=7 убывание гарантировано.

Перейдем к программированию функции Sin2(x). Используем принцип 7 – проектирование сверху-вниз. Объявим в классе Program метод:
static double Sin2(double x)
а перед оператором Console.ReadKey(); в методе Main() вставим:
y = Sin2(x);
Console.WriteLine(«sin(<0>)=<1>», x, y);

Наша функция, также как и библиотечная, получает на вход значение угла в радианах и возвращает значение sin(x). Тип аргумента и самой функции задан как double. Приступим к разработке начинки метода Sin2() – блока в фигурных скобках.

Реализация 2.
Очевидный, но не лучший алгоритм состоит в непосредственной реализации формулы ряда: начальное значение суммы ряда s = x, затем в цикле for добавлять каждый вычисленный член ряда с чередованием знака rn=x n /n!, для чего будет необходимо написать еще две функции: возведение в степень и вычисление факториала.
Не совсем понятно, сколько раз выполнять цикл. Хотелось бы увязать необходимое число членов ряда с требуемой точностью вычислений eps.
Добавим в метод Sin2() локальную константу double eps = 1E-15.
Вместо цикла for применим цикл while (условие) < >.
Введем две переменные: double r – активный член ряда, int n – степень.
Для учета смены знака используем функцию модуля числа – Math.Abs( ).

Тогда метод Sin2() может быть записан так:

Выполним тестирование программы:

Метод Объяснение
Угол a, градусы Угол x, радианы Sin(x) Sin2(x)
30 0,523598775598299 0,5 0,5
60 1,0471975511966 0,8660254037844439 0,8660254037844438
90 1,5707963267949 1 1
270 4,71238898038469 -1 -1
360 6,28318530717959 -2,45E-16 -4,16E-15
10000 174,532925199433 -0,984807753012209 -3,07952784012437E+58

При углах, не превышающих 360°, оба метода дают практически одинаковый результат, погрешность почти не превышает eps=1E-15. При больших значениях углов Sin2() иногда выдает неверный результат (|sin(x)|

#помощь — Синусоида в C++

PS
пробовал сложение ряда sin(x)=x-x 3 /3!+x 5 /5!-x 7 /7.
но там после 6-ого шага значение только увеличивается.

Всего записей: 29 | Зарегистр. 04-01-2008 | Отправлено: 19:47 04-04-2008
Cheery

.:МордератоР:.

Редактировать | Профиль | Сообщение | Цитировать | Сообщить модератору kein2 3

Цитата:

но там после 6-ого шага значение только увеличивается.

да потому что x должен быть в радианах и в пределах от -pi до pi
http://en.wikipedia.org/wiki/Taylor_series
7 членов хватает для аппроксимации

———-
Away/DND
Всего записей: 52737 | Зарегистр. 04-04-2002 | Отправлено: 19:49 04-04-2008 | Исправлено: Cheery, 19:55 04-04-2008
kein2 3

Newbie

Редактировать | Профиль | Сообщение | Цитировать | Сообщить модератору я проверял на значении pi с 20 цифрами после запятой, самое меньшее значение функции: на 6-ом шаге -0.00044, на 7-ом +0.001

это нормально?

Всего записей: 29 | Зарегистр. 04-01-2008 | Отправлено: 21:31 04-04-2008
Cheery

.:МордератоР:.

Редактировать | Профиль | Сообщение | Цитировать | Сообщить модератору kein2 3

Цитата:

я проверял на значении pi с 20 цифрами после запятой, самое меньшее значение функции: на 6-ом шаге -0.00044, на 7-ом +0.001

что вы проверяли??

по y — значение, по x — число членов ряда

———-
Away/DND
Всего записей: 52737 | Зарегистр. 04-04-2002 | Отправлено: 02:32 05-04-2008 | Исправлено: Cheery, 02:32 05-04-2008
kein2 3

Newbie

Редактировать | Профиль | Сообщение | Цитировать | Сообщить модератору ааа, все понял. я тут ошибку небольшую допустил, сейчас все вроде бы нормально спасибо Cheery за помощь
Всего записей: 29 | Зарегистр. 04-01-2008 | Отправлено: 16:54 05-04-2008
susuman

Newbie

Редактировать | Профиль | Сообщение | Цитировать | Сообщить модератору Держите синус из FreeBSD (/usr/src/lib/msun/src/k_sin.c) :

/* __kernel_sin( x, y, iy)
* kernel sin function on

[-pi/4, pi/4] (except on -0), pi/4

0.7854
* Input x is assumed to be bounded by

pi/4 in magnitude.
* Input y is the tail of x.
* Input iy indicates whether y is 0. (if iy=0, y assume to be 0).
*
* Algorithm
* 1. Since sin(-x) = -sin(x), we need only to consider positive x.
* 2. Callers must return sin(-0) = -0 without calling here since our
* odd polynomial is not evaluated in a way that preserves -0.
* Callers may do the optimization sin(x)

x for tiny x.
* 3. sin(x) is approximated by a polynomial of degree 13 on
* [0,pi/4]
* 3 13
* sin(x)

#помощь — Синусоида в C++

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

Варианты реализации DDS на MCU

Когда используется прямой цифровой синтез синусоидального сигнала (DDS), требуется вычислять величину сигнала в моменты времени, определяемые тактовым сигналом. Это означает, что необходимо с высокой скоростью вычислять значения синуса. Для чего можно использовать несколько методов, предъявляющих различные требования к производительности процессора и объёму используемой памяти.

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

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

С помощью формул приведения, тригонометрическая функция может быть выражена через тригонометрическую функцию другого аргумента, в определённых ситуациях меньшего (по модулю), чем исходный. Благодаря этому, полная таблица синусов, расcчитанная для значений аргументов из полного периода [0, 2π), может быть заменена таблицей для четвёртой части периода, т.е. может быть уменьшена по размеру в 4 раза. Тем самым мы экономим значительный объём памяти, но теряем в быстродействии из-за необходимости выполнять дополнительные вычисления.

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

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

Быстрое вычисление синуса суммированием ряда

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

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

Данные ряды сходятся и имеют в качестве суммы функцию sin x, cos x соответственно при любом x. Однако, чем меньше x, тем сходимость более быстрая. Допустим, если |x| * Значит, пользуясь формулами (3), (4), можем вычислять sin с точностью не хуже 4e-5 и cos с точностью до 4e-6. А это даже более высокая точность, чем требуется при управлении 12-разрядным DAC с его 1LSB=1/4096 (порядка 2e-4).

* На самом деле тригонометрические функции любого аргумента можно выразить через sin, cos ещё меньших углов, но в данном случае это не приведёт к упрощению вычислений.

Теперь о том, как выразить тригонометрическую функцию произвольного аргумента \(\phi\), через sin, cos величины, по модулю не превосходящей π/4. Представим аргумент в виде $$ \phi=\frac2+\alpha, \\ -\pi/4 \le \alpha \lt \pi/4, $$ где n — целое число. Тогда, например, функция sin в соответствии с формулой синуса от суммы может быть выражена следующим образом $$ \sin \phi=\sin\frac2 \cos\alpha+\cos\frac2 \sin\alpha, $$ $$ \sin\phi=\sin\alpha, \text < если >n=4z \text < или >(n\&3==0), \\ \sin\phi=\cos\alpha, \text < если >n=4z+1 \text < или >(n\&3==1), \\ \sin\phi=-\sin\alpha, \text < если >n=4z+2 \text < или >(n\&3==2), \\ \sin\phi=-\cos\alpha, \text < если >n=4z+3 \text < или >(n\&3==3), $$ здесь z — некоторое целое число.

Осталось для данного \(\phi\) подобрать подходящие значения \(\alpha\) и n. Заметим, что $$ \phi=\frac2+\alpha=\frac<\pi>2\left(n+\frac<\alpha><\pi/2>\right), \tag <5>\\ -1/2 \le \frac<\alpha> <\pi/2>\lt 1/2, $$ тогда $$ 0 \le \frac<\alpha><\pi/2>+\frac 1 2 \lt 1. $$ Следовательно, в соответствии со свойствами функции — целой части числа и с учётом того, что n — целое, можно записать: $$ floor\left(\frac<\alpha><\pi/2>+\frac 1 2 \right)=0 \Rightarrow \\ floor\left(n+\frac<\alpha><\pi/2>+\frac 1 2 \right)=n. $$ Но из (5) следует, что $$ n+\frac<\alpha><\pi/2>=\frac<\phi><\pi/2>, $$ значит, $$ n=floor\left(\frac<\phi><\pi/2>+\frac 1 2 \right), \\ \alpha=\phi-n\frac<\pi>2=\phi-\frac<\pi>2 floor\left(\frac<\phi><\pi/2>+\frac 1 2 \right), $$ где \(-\pi/4 \le \alpha \lt \pi/4\).

Практическая реализация

Ниже приведён пример реализации быстрого вычисления тригонометрических функций с использованием арифметики с фиксированной точкой на C++ . Предполагается, что числа с фиксированной точкой представлены целыми числами со знаком размером 32 бит. На целую и дробную части выделяем по 16 бит.

В связи с весьма ограниченной точностью чисел, ошибка вычислений будет превышать значения, которые дают формулы (3), (4). Это связано с погрешностями округления и накоплением ошибки при выполнении последовательности операций. Численный эксперимент с перебором всех возможных значений аргумента с фиксированной точкой от -π/4 до π/4 показывает, что наибольшая ошибка по модулю не превышает 6e-5 при вычислении значений sin и не превышает 4e-5 при вычислении значений cos. Числу 6e-5 в формате fix16 соответствует целое 2 16 *6e-5=65536*6e-5, что приблизительно равно 4, а значит 2..3 младших бита полученного значения не являются достоверными, в то время как 13..14 старших битов дробной части являются верными. Точности полученного результата вполне достаточно для работы с DAC с разрядностью до 12. Ошибка вычислений не превышает 1/4 LSB 12-разрядного DAC.

Вычисление синуса как функции от индекса фазы

Теперь мы имеем функции, которые дают нам возможность быстрого вычисления значений sin, cos, при условии что аргумент x отвечает условию |x| k=0..2^<32>-1, \> 0 \le \phi \lt 2\pi. $$ Но реализованные ранее функции _fast_sin, _fast_cos требуют, чтобы аргумент по модулю не превышал значения π/4. Воспользуемся тем, что как мы уже установили, \(\phi\) можно представить в виде $$ \phi=n\frac<\pi>2+\alpha, -\pi/4 \le \alpha \lt \pi/4, $$ где n — целое число, вычисляется следующим образом: $$ n=floor\left(\frac<\phi><\pi/2>+\frac 1 2\right), $$ тогда $$ \alpha=\phi-n\frac<\pi>2. $$

В нашем случае \(\phi=2\pi k/2^<32>\), тогда $$ n=floor\left(\frac<2\pi k><2^<32>>\frac 2<\pi>+\frac 1 2 \right)= floor\left(\frac<4k><2^<32>>+\frac 1 2\right)= \\ floor\frac<4k+2^<31>><2^<32>>=floor\frac<4k+4\cdot 2^<29>><2^<32>>, \\ n=floor\frac><2^<30>>, $$ что с использованием средств языка C++ можно записать в виде

В зависимости от значения n, \(\sin\phi=\pm\sin\alpha\) или \(\sin\phi=\pm\cos\alpha\), где $$ \alpha=\phi-n\frac<\pi>2. $$ Бит [0] числа n определяет функцию, через которую выражается \(\sin\phi\) (0 — sin, 1 — cos), бит [1] определяет знак, с которым берётся эта функция (0 — «+», 1 — «-«).

Удобнее сначала вычислить соответствующий \(\alpha\) индекс фазы: $$ k_<\alpha>=\frac<2^<32>><2\pi>\alpha=\frac<2^<32>><2\pi>\frac<2\pi><2^<32>>k-\frac<2^<32>><2\pi>\frac<\pi>2 n, \\ k_<\alpha>=k-n\cdot 2^ <30>$$ или

Тогда $$ \alpha=2\pi\frac><2^<32>>. $$ Здесь \(k_<\alpha>/2^<32>\) можно рассматривать как число \(k_<\alpha>\) в формате fix32 (целая часть равна 0, все 32 бита представляют дробную часть числа).

Функция, возвращающая значение синуса для данного индекса фазы может быть реализована следующим образом (использует определённые ранее _fast_sin, _fast_cos ):

Погрешность вычисления не превышает 7e-5.

Пример программы для микроконтроллера

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

Микроконтроллер STM32F100RB при работе с тактовой частотой 24 МГц, способен осуществлять синтез сигнала с опорной частотой до 250 кГц, что позволяет получать сигнал с частотой от 0 до 60..100 кГц (при компиляции проекта с оптимизацией по быстродействию). Но это при практически полном задействовании вычислительных ресурсов микроконтроллера. Результаты получены при запуске из RAM, при запуске из ROM архитектура ARM обеспечивает большее быстродействие.

Если используется компиляция без оптимизации (например, при отладке), можно достичь опорных частот до 90 кГц.

Максимально достижимая частота оказалась ниже, чем при использовании табличного метода, где поддерживаются значения опорной частоты до 400 кГц даже в случае применения «тяжёлых» вариантов метода.

Следует очень аккуратно выбирать опорную частоту — реализуемы только значения, которые можно получить из тактовой частоты используемого таймера путём деления на целочисленный делитель (из диапазона 2..65536). Чтобы избежать ошибок, если требуется нецелое значение опорной частоты, лучше задать требуемый коэффициент деления непосредственно в коде, соответствующим образом сконфигурировав используемый таймер.

Как на C++ сгенерировать синусойду и вывести её на колонки?

Надо сгенерировать синусойдный звук (с определённой частотой и амплитудой) и вывести его на колонки, Подскажите плиз, как это сделать?

Link90
может маркером?

DirectSound
заполняешь буфер своим синусом и поешь.
просто до безобразия.

Link90
OpenAL — «заполняешь буфер своим синусом и поешь» (с) MOD

Link90
Audiere — создаеш поток и поешь

FruityLoops — ставишь 3xOSC и поешь :)

asig oscil amp, freq, ftnum
out asig

(ну и еще парочку строчек) и поешь :-)))

записываешь на плеер то, что вывели все предыдущие ораторы, подключаешь к колонкам и поешь. =)

Раньше была хорошая статья «ПРОГРАММНАЯ ГЕНЕРАЦИЯ ЗВУКОВЫХ СЭМПЛОВ»
на http://www.programme.ru/ , но теперь этот домен стал свободен(

Цитирую:
Например, для 16-килобайтного сэмпла синусоидального звука, с частотой 1000Hz и качеством сэмпла 44100Hz — наши параметры должны иметь следующие значения: samplerate=44100, wavefrequency=1000, samplelength=16384.

Особенного пояснения требуют параметр wavevolume и формат представления звука.
Для дискретного представления звука записывают значение громкости (амплитуды волны) в определенный момент времени; мы будет пользоваться для этого переменной vol.
Естественно, качество воспроизведения звука пропорционально зависит от битности его представления (8-bit, 16-bit, 24-bit и т.д.). Для 8-bit – значения vol от 0..255, для 16-ти – 0. 65536, для 24-х – 0. 16777216. Какой вариант выбрать?

Смотря какая у вас задача; но меньше 16-ти я бы не советовал (хотя бывают и исключения – когда надо сократить объем сэмпла взамен качеству). Все приведенные здесь примеры расчитаны на 16-битный звук (поэтому и vol, и buffer[] – типа short, 16-bit)

Вот отрывок кода из статьи как сгенерировать синусойду.

Написать программу, которая выводит синусы и косинусы. На языке Си. Заранее спасибо

Ответ

#include /* Подключаем библиотеку, чтоб работали функции cout, cin */
#include /*Подкючаем библиотеку, чтоб работали функции cos, pi */
#define PI 3.14159265 /*Задаем сколько у нас будет значит число ПИ */
using namespace std;

main () <
float a, result; /*Вводим 2 переменные — действительные числа: одно — число градусов угла, второе — полученный результат*/
cout «Vvedite znachenie ugla» > a; /* Вводим значение угла с клавиатуры*/
result = cos (a*PI/180); /*высчитываем значение косинуса, переводя значениезаданного угла в радианы */
cout «Kosinus » » raven » /*Выводим полученное на экран */
return 0; /*Возвращает нулевое значение */
>

Цукерберг рекомендует:  Слой с аннотацией на CSS3
Понравилась статья? Поделиться с друзьями:
Все языки программирования для начинающих