Требуется сформировать 2 синуса (естественно в цифре). Можно пойти табличным методом - но для 8 бит и 1к хватит, а для 18 бит и метра может мало оказатся. Есть другой метод http://www.gaw.ru/html.cgi/txt/doc/...v/max2000_5.htm .
Но интересные эфекты: даже при инт64 вылазеют гармоники на уровне -70 дб + к этому и частота не точно заданная.
Попробывал так:
X(n)=k*X(n-1)-X(n-2), где k=2*cos(2*Pi*F/Fd);
для 0 градусов X(-1)=0, X(-2)=-A*sin(2*Pi*F/Fd);
для 90 градусов X(-1)=A, X(-2)=A*sin(pi/2+2*Pi*F/Fd)
90 гр. частота на сотые доли процентов выше и уползает вперёд по отношению к 0 гр.
Сейчас хочу через приблежённое вычеслени по Тейлору попробывать. Но явно по вычислениям более накладно будет.
Может ещё есть способ с минимум вычислений сгенерить синус?
(Всё должно в итоге оказатся в FPGA Спартан3).

[quote name='aal' date='Jan 24 2006, 12:29' post='80720']
Требуется сформировать 2 синуса (естественно в цифре). Можно пойти табличным методом - но для 8 бит и 1к хватит, а для 18 бит и метра может мало оказатся. Есть другой метод http://www.gaw.ru/html.cgi/txt/doc/...v/max2000_5.htm .
[\quote]

Дайте точную ссылку на метод

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

Для примера: Нам требовалось иметь синусы 425 Гц, 2100 Гц, 2600 Гц на выходе
ИКМ кодека (8000 выборок/сек по 8 бит).

Сначала думали придеться ставить ПЗУ, но поигравшись с подбором
кратных частот нашли замены: 421,053 Гц (19 выборок), 2105,26 (19 выборок),
2600 (40 выборок).

Причем Synplify очень хорошо оптимизирует таблицы, и вышеприведенные синусы
заняли в ПЛИС Actel Ex256 не больше 10 %.

http://www.gaw.ru/html.cgi/txt/doc/micros/...v/max2000_5.htm

Требуется генерить синусоиды с частотами от 300 Гц до 2,5 МГц с шагом 1 Гц.
Желательно с нелиненостями -90 дБ. Проект - селективный вольтметр, с полосами 30 и 100 Гц, желательно и со спектроанализатором.
Основная идея - сигнал с ацп перемножаем с синусом и косинусом по отдельности, фильтруем и децимируем. перемножаем каждый из них на самого себя, складываем и из постоянной состовляющей извлекаем корень квадратный - значение уровня сигнала в выбранной полосе на выбранной частоте.
Т.к. в спартане умножители 18 битные, то и синус желательно сформировать 18 битный. Частота опорного генератора 10 МГц.

Что-то в постановке не то: -90дБ до 2.5МГц с шагом 1 Гц с тактовой 10МГц, это не описка?
Попробуйте посмотреть такой метод:
http://archive.chipcenter.com/dsp/DSP080201F1.html
и смежные статьи на том же сервере.

Когда-то от аналоговых времен помню, что пробовал генерить синусоиду при кусочно-линейной аппрокс. по методу треугольников.
Возможно, если удастся подобрать что-то в этом духе, ( делал я это для реализации в альтере лет 8 назад и деталей уже не помню ) то из таблицы надо будет брать не сам сигнал, а только поправку. Следовательно объем таблицы значительно уменьшится.
Удачи!

Есть такая штука как резонантный фильтр, подробности по ссылке
http://www.research.scea.com/research/pdfs...rmath_GDC02.pdf

Стандартный способ вычисления значения функции заключается в полиномиальном приближении ее на отрезке и отображении любого значения функции на этот отрезок. Например, для косинуса и синуса:
1 Разбиваем один период косинуса на четыре части, учитывая, что cos(-x)=cos(x); cos(2*pi*n+x)=cos(x); cos(pi+x)= -cos(x); sin(x)= -cos(pi/2+x). Пользуясь этими соотношениями, значение аргумента приводится к отрезку [0; pi/2].
2.Участок функции с аргументом [0; pi/2] аппроксимируем полиномом: sin(x) ~ A0+A1*x+A2*x^2+...An*x^n, где А0...Аn - коэффициенты полинома. Для наилучшего приближения (в смысле наименьших средних квадратов), если склероз не изменяет, это к-ты полинома Лежандра, и ни в коем случае не Тейлора! Для Ваших требований хватит n=5-6. О вычислении к-тов нужно посмотреть в справочнике по математике. Значение функции вычисляется рекуррентно, на Спартане должно реализоваваться в пол-пинка...
Можно попробовать аппроксимировать косинус и на отрезке [-pi/2; pi/2], если это будет удобнее для реализации на ПЛИС. Тогда останутся только четные степени раздожения, но вычислений будет несколько больше.

Вообще-то такую штуку можно реализовать и на недорогом DSP...

Я в своих изысканиях пришёл к выводу что лучше всего получается если сводить всё к отрезку [-pi/4;+pi/4] и считать либо синус либо косинус как разложение в ряд в окрестности нуля. На таком интервале они быстрее сходятся.
В данной конкретной задаче придётся считать в любом случае обе функции, и синус и косинус.

Ещё одно интересное разбиение - на четыре (на самом деле получается три) разных полинома:
sin(x) и cos(x) при x=[-pi/8;+pi/8] - как обычно - в окрестности нуля, а на интервале
x=[+pi/8;+3pi/8] считать разложением в окрестности pi/4, коэффициенты для sin и cos получаются одинаковыми с точностью до знака.

Чем ряды считать, лучше уж по фомуле sin(a+b) - быстрее будет это точно ( это в принципе и есть резонантный фильтр, я просто настаиваю на том, чтобы вы ознакомились с предложенным мной документом ). Особый шик - при некоторых частотах умножений для вычисления синусов/косинусов вообще не потребуется.

Спасибо Всем откликнувшимся!
Изучаю Вашу инфу. Проведу эксперименты - обязательно отчетаюсь.

Все-таки, попробуйте использовать для разложения многочлены Лежандра - будете немало удивлены. Дело в том, что Тэйлор наилучшим образом приближает ф-цию только в бесконечно малой (дельта-) окрестности конкретной _точки_, а пол-м Лежандра - на _отрезке_.

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

Вопрос таким образом сводится к выбору наилучшего сочетания: объём таблицы - степень интерполяционного полинома.

Сейчас посмотрел в исходники IAR-овского sin() - там над вопросами эффективности не особо задумывались, считают Тейлором суммируя до 9-й степени для float и аж до 17-й для double.
То ли дело фортрановские библиотеки для PDP - когда-то они были для меня источником алгоритмов (интернет тогда уже был, только до нас ещё не добрался, а гугл ещё не выдумали - за каждой фигнёй приходилось тащится в библиотеку, за N километров!, по морозу! не то что нынешняя молодёжь! )

Нет-нет, табличным способом задаются только коэффициенты при степенях аппроксимирующего полинома, которые находятся аналитически, а узловые значения функции задавать вовсе не нужно. Точно так же, как в тейлоре, но, в отличие от него, эти к-ты будут зависеть от степени полинома. Наверное, и в ИАРе используется не тейлоровское разложение, здесь может сбить с толку абсолютно идентичный способ вычисления значения полинома. Кстати, представление лежандром оптимально только в смысле наименьших средних квадратов на отрезке, поэтому я и предложил его для генератора - будет давать минимум гармоник. Можно задаться другим критерием, напр., минимизацией масимальной ошибки на отрезке, тогда к-ты поинома будут слегка другими. Если мне опять склероз не изменяет, приближение по Лежандру дает для синуса на отрезке [0;pi/2] по 3 верных двоичных разряда за итерацию, что хорошо сочетается с цифрами, приведенными Вами. Тейлор же такой точности (24р за 9 итераций) не даст. Если выложите исходный текст ИАРовского синуса, постараюсь выяснить это - самомУ стало интересно.

Точно, от к-тов Тейлора немного отличаются в младших разрядах, наверняка какой-нибудь полином.
Уступив настоятельной просьбе ознакомился с предложенным KPAH документом, вставляет.
Исходник синуса прилагается.
#define __TWOOPI 0.63661977236758134308
Нажмите для просмотра прикрепленного файла

SSerge, так вам надо постоянно генерировать синусоиду во времени или всё-таки синусы произвольных аргументов считать?

Для первого случая ничего лучше sin(a + b) не придумаешь.

Для второго случая: я писал библиотеку для TI DSP 64XX. Использовал свой формат с плавающей точкой мантисса - 23 бита, только нормированные числа так что значащих бит вообще 22. Синус считал тейлором на отрезке 0 - Пи/4 (для синуса на этом отрезке полином был такой x + A*x^3 + B*x^5 + C*x^7) в итоге для всего допустимого диапазона значений аргумента (что-то типа до 2^24) получается погрешность 4 ULP от floata, то есть 2 ULP от моего представления. A, B и C похожи на тейлоровские но слегка модифицированные чтобы минимизировать относительную погрешность. По любому больше чем 18 бит точности (не забываем ещё и про argument reduction, который у меня есть и соответственно он тоже привносит ненулевую погрешность). Без argument reduction на 0 - Пи/4 ошибка 1 ULP флоата (в принципе это неустранимая ошибка так как мантисса на 1 бит меньше).

вот на всякий случай кодд:

CAFloat sin1 = -0.16666666666666666666666666666667f;
CAFloat sin2 = 0.0083333333333333333333333333333333f;
CAFloat sin3 = -0.0001984126984126984126984126984127f*0.9925f;

inline CAFloat CF_sin_reduced(const CAFloat &x)
{
if(_ext(x.x, 0, 24) < -11)
return x;

CAFloat arg = x * x;
return x + x*(sin1 + arg*(sin2 + arg*sin3))*arg;
}

Посмотрите эту статью.
Можно зайти на сайт GlobalDSP
Стастья за декабрь 2003 года -> Algorithms (Fast, Continuous, Sine Wave Generator)
Судя по статье, результаты не плохие. Там же есть пример кода на DSP.

Успехов.Нажмите для просмотра прикрепленного файла

Попробуйте применить DDS (прямой цифровой синтез) или, по крайней мере, применить те же алгоритмы

Это же чистая табличная реализация, с нее-то и начинали. А в Spartan'е, между прочим, память свободная не завалялась? А то можно взять и просто сделать аккумулятор фазы бита на 32 и иметь ~0,01Гц разешение (см. теорию DDS лучше по datasheet на AD9854) . Хотя -90dB получишь где-то в районе 32к-точечной таблицы...

А микросхему DDS поставить не проще?

Я реализовал DDS на сигнальном процессоре. Для сокращения объема выборок (оставил 1024 отсчета) использовал линейную интерполяцию между соседними отсчетами (старшие 10 бит фазового аккумулятора адресуют выборку из таблицы, а младшие (32-10) бит используются для интерполяции). На интерполяцию ушло 2 умножения и 3 сложения. Проверял в MathCad, отличие от чистой синусоиды не более 5*10^-4 % (я, правда, не учитывал, что выборки 16-битные).

А что, линейная интерполяция... не так уж и плохо... А хранить можно только pi/2, усложнив адресацию и добавив (+/-).
Кстати, где-то видел науку о кусочно-линейной интерполяции синуса, причем с нелинейным разбиением по времени (в зависимости от скорости нарастания). Приводило это все к радикальному уменьшению размера памяти выборок.
5e-4% - это не так и много, около -66dB, а как получено это значение (среднее квадратичное отклонение или...)? Посчитать бы еще Кг.

Думаю, что ещё раз в 10 улучшить интерполяцию синуса можно если линейно интерполировать производную, а потом брать от неё интеграл. Например если нарисовать производную в виде треугольного сигнала, то первообразная будет очень похожа на синус.

Дык производная ж - тот самый синус: sin'(x) = cos(x) = sin(x+pi/2).
И про интегрирование не понял, как его делать? Каждый новый отсчет в копилку класть? Тогда быстро за пару десятков периодов набежит ошибка.

Есть две табличные точки А и Б. Б-А = интеграл от производной. При линейной интерполяции подразумевается производная = const = (Б-А)/шаг таблицы. Однако если производную наклонить так, чтобы "определённый интеграл" не изменился, то получится гораздо большая точность. При этом производная будет прямой наклонённой линией. Никаких накоплений для этого делать не надо. Немного больше арифметических действий и фсё.

Готовый DDC, например AD6620, не хотите попробовать?

А так, для этой задачи не столько нелинейность синуса важна, сколько SFDR (хотя они взаимосвязаны). При этом спектральные "палки" можно с помощью дизеринга "размазать". У нас, например, -110dBc SFDR получалось на таблице 256x25 бит с линейной интерполяцией и дизерингом фазы (около -90dBc без дизеринга).

Подскажите кто-нибудь сколько я насчитал децибелл.
Составил прогу на паскале с линейной интерполяцией синуса по таблице 256*16 бит. Максимальное отклонение по центру между колонками составило 3.1e-6 относительно самого синуса (+-1.0).

Аномалия какая-то!!!
Таблица 128*13 бит с линейной интерполяцией даёт ещё большее качество - 0.9e-6 !!!


Кто-нибудь может мне это объяснить?

Странно. При ширине выхода таблицы 16бит, максимальное отклонение от "истинного" синуса должно было быть порядка 1.5e-5. Или в программе что-то не то, или я не правильно понял, что такое "Максимальное отклонение по центру между колонками" .

При отклонении от "идеального" синуса на 3.1e-6, SFDR даже в самом поганом случае (вся "энергия" отклонения "ушла" в один спур), думаю, будет ниже -110dBc.

Вот прога на паскале. Если кто хочет меня проверить.


Правда? А для 0.9e-6 сколько децибелл будет?

Я вообще не врубаюсь. Кто-нибудь, ткните меня носом, где я ошибся.
_____________

Да. Поспешил, людей...
Там в проге в конце надо было поставить
WriteLn(er).

Теперь правильный результат:
для 256*16 бит ошибка составит 8.2e-5

Если я правильно понял, то ты смотришь отклонение от "идеального" синуса толко по середине между табличными значениями и тебе просто "повезло" (что само по себе забавно). Попробуй вычислить отклонение для, скажем, 8 точек на шаг таблицы (всего 1К точек), включая значения в "узлах" таблицы.

Для 0.9e-6, думаю, SFDR будет ниже -120. Для оценки: SFDR <= 20*lg(d/A), где A - амплитуда тона,
а d - максимальное отклонение от "идеального" синуса.

Еще можно sin cos от фазового аккумулятора CORDICом посчитать, если тактовая позволяет то вообще последовательным, на FPGA вс ё это очень хорошо ложится.

Уважаемые коллеги! Для сведения. В статье Никулин Э. С., Пахоменков Ю.М. Метод построения низкочастотных генераторов сигналов сложной формы // Судостроительная промышленность. Автоматика и телемеханика. – 1992. – Вып. 14. – С. 37-48 описан метод воспроизведения временных зависимостей, путём решения автономного дифференциального уравнения, которое, в частности, для функции x(t)=Asin(wt) на четверти периода имеет вид dx/dt=w [A·(1-x^2)^0.5]. В этом уравнении множитель в скобках [] соответствует функции y(t)=Acos(wt). Метод неплохо работает для низкочастотных сигналов и удобен там, где для других задач, например для вычисления расхода рабочей среды, уже имеются средства извлечения корня.

К сожалению до первоисточника не добрался. Попробовал в лоб посчитать в MatLab. Получил комплексные числа. Потом подумал... Четверть периода - это pi/2? т.е. ~1.57 радиан. Таким образом, под корнем имеем отрицательные числа. Да и сигнал на синус не похож абсолютно, втыкается под прямым углом в ось абсцисс, по науке д.б. 45 градусов.

KPAH
Большое спасибо за ссылку. Очень документ порадовал (хотя такой метод генерации синуса я знал и раньше).
0FF/2 А есть ли ещё что-нибудь в этом же духе, т.е. в виде справочника по разным методам, алгоритмам и прочее?

Нажмите для просмотра прикрепленного файлаПояснение
В автономном уравнении dx/dt=w [A·(1-x^2)^0.5]. Четверть периода - это abs(wt)< pi/2?.
Для большей наглядности прикрепляю рисунок