Всем привет!
Вот есть задача - генерировать квадратурный сигнал для DDC. DDC сам реализован внутри DSP, поэтому надо алгоритм чтобы его можно было-бы реализвать потом на аасемблере.
Используется плавующая точка.
В идеальном случае ф-я должна выглядеть так:
result = GenSinCos(float freq); где result это структура содержащая Cos и Sin, т.е. сигналы сдвинуты на 90 градусов.
freq - 0...1.
Ф-я вызывается для каждого отсчета. Начальная фаза в общем не очень важна. Главное чтобы 90 градусов было. По идее это даже не ф-я а макрос.
Как-то так.
Может кто уже встречал подобное? Или может уже есть готовые решения?

Особенность в том, что наверняка есть решение которое не вызывает ф-ю sin и cos каждый раз.

Спасибо!

PS. Табличный метод не предлагать!

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

sin = sqrt(1 - cos^2) ?

Это прекрасно!! хороший синус. будет всегда не меньше 1. а иногда и побольше.



Что касается начальной задачи, то решение есть!!! явочным порядком можно вызывть синус и косинус например каждый 10-й раз.. или даже 100-й. можно вообще ни разу не вызывать. всё зависит ведь от необходимой точности описания функций

c(n)=cos(2*pi*f0/sf*n)

s(n)=k*c(n)+c(n-1)-k*s(n-1)

k=(sin(2*pi*f0/sf)-1)/(sin(2*pi*f0/sf+pi/2))

Это точно! Вот хороший синус: sin(pi/2 - 2.06345i) = 4.0001.

большое спасибо за правку, торопился и ошибся

Ну, конечно, CORDIC. Какие могут быть сомнения?!

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

Поясню. Есть методы генерирования синусоидально сигнала на основе биквада. Т.е. там просто рекурсивно вычисляется каждое следуюшее значение. Генератор в чистом виде.
Там меняя коэффициенты биквада, можно менять частоту. Там фаза не важна, важна только частота.
Я думаю, что должны быть методы, которые позволяют генерировать синус и косинус вместе, подобным образом. Дополнительным условием для данной ф-и есть разница фаз в 90 градусов.
Т.е либо начальные условия для генератора задавать, либо еще что-то. Но в общем это должно быть рекурсивное вычисление, где каждый следующий отсчет зависит от предыдущих. Это не синус и косинус в чистом виде. Нужен именно генератор (так как он более экономичен к вычислениям).

Школьная тригонометрия рулит..

sin[2*pi*f*(t+Δt)] = sin[2*pi*f*t] * cos[2*pi*f*Δt] + cos[2*pi*f*t] * sin[2*pi*f*Δt];
cos[2*pi*f*(t+Δt)] = cos[2*pi*f*t] * cos[2*pi*f*Δt] - sin[2*pi*f*t] * sin[2*pi*f*Δt];

Полагая, что t+Δt == (n+1)*Δt, находим формулы рекурсии:

sin_n+1 == sin[2*pi*f*(n+1)*Δt] = sin[2*pi*f*n*Δt] * cos[2*pi*f*Δt] + cos[2*pi*f*n*Δt] * sin[2*pi*f*Δt] == sin_n * cos[2*pi*f*Δt] + cos_n * sin[2*pi*f*Δt];
cos_n+1 == cos[2*pi*f*(n+1)*Δt] = cos[2*pi*f*n*Δt] * cos[2*pi*f*Δt] - sin[2*pi*f*n*Δt] * sin[2*pi*f*Δt] == cos_n * cos[2*pi*f*Δt] - sin_n * sin[2*pi*f*Δt];

Понятно, что sin[2*pi*f*Δt] и cos[2*pi*f*Δt] это константы, зависящие от f и Δt и выбираются так, чтобы 1/(f*Δt) = 2^N == n_max+1.

Начальные значения:

sin₀ = 0;
cos₀ = 1;

Значения после переполнения счетчика n = (n_max+1)%2^N == 0, полагаем равными:

sin_2^N = 0;
cos_2^N = 1;

Это предотвратит накопление ошибок..

Соглашусь FatRobot, выбор алгоритма должен определяться тем, сколько надо SFDR и сколько mips можно потратить. Отбрасывать что-то с порога глупо.

Спасибо. Надо попробовать.
Как попробую - напишу.
Пока все выглядет очень даже ничего!

Насчет любых DSP не согласен. Недавно решали задачу, в которой нужно вычислять много арктангенсов, DSP TigerSHARC 101, CORDIC + SIMD дал результат лучше (в смысле вычислительных затрат) чем оптимизированный полиномиальный, хотя отдельно взятый CORDIC проигрывает.

ТС: хочу обратить Ваше внимание, что в DDC с целью подавления спуров при генерации квадратурной гармоники используются некоторые хитрости, например, младшие разряды заменяются на соответствующий отрезок ПСП. Если не ошибаюсь, в Матлабе есть такая модель (или у Xilinx-а в генераторе корок, запамятовал уже).

Тут всё на веру, конечно. Не ставя под сомнение квалификацию. Но вот здесь http://www.opensource.apple.com/source/Lib...ce/Intel/atan.c
за где-то 15 умножений и 15 сложений получается точность double на входе и выходе для atan.

кто/что мешает генерить синус-косинус двумя параллельными биквадратными звеньями?
начальная фаза целиком и прлностью определяется начальными условиями звеньев.
умножений в 2 раза меньше чем в школьнотригонометрическом методе.

Там есть страшная вещь - зависимость по данным, которая рушит конвейер, от нее не уйти и сходу просаживается производительность (в случае TigerSHARC) минимум в 2 раза, возможности SIMD при использовании плавающей точки ограничены, использование фиксированной точки для реализации приведенного алгоритма сомнительно, т.к., скорее всего, очень быстро накопится достаточно большая ошибка. Ну и в приведенной реализации нет операции деления Im / Re, что в случае DSP, как правило, делается каким-нибудь Ньютоном-Рафсоном, CORDIC же позволяет не выполнять деление.

А во сколько раз больше сложений?

в 1.

Как вы оцениваете, при таком подходе через сколько сгенерированных отсчетов существенно изменятся амплитуды квадратур (+/- 10% разность или каждая) и/или разность их фаз (pi/20)? За счет округлений.

Все операции, например, соответствуют 754 float.

Я полагаю, что за пару часов при 10 MSPS и при 1 MHz на выходе от квадратур в таком генераторе уже ничего не останется.

Никак не оцениваю, ибо оценка качества таких генераторов и необходимость периодического сброса сводит все их прелести почти к 0.
Это тс фанатеет рекурсиями, которые на конечной разрядности дают рост ошибок.
А для меня полиномы/кордики рулят.

LUT с линейной интерполяцией - наше все Ну и всякие навороты типа noise shaping или phase dithering, если охота шибко умным побыть.

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

PS. Кстати, а вот и ответ: генератор делается на основе формулы предложенной Blackfinом и таблицей. Т.е. если у меня в потоке обрабатывается по 128 отсчетов, то все что мне нужно это каждый раз инициализировать значение синуса и косинуса, значениями из таблицы, а потом считать рекурсию на 128 отсчетов. И так далее.

грустный смайлик..

Ну, вам виднее.

И нужно помнить, что ошибки амплитуды и фазы будут расти примерно как:

sqrt(N)*0.5 [LSB] = sqrt(128)*0.5 [LSB] ~= 5.7 [LSB].

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

У меня плавучка.
Там немного подругому. Но в любом случае 128 отсчетов это не так много даже для 24 бит фиксированной точки. Но при этом есть возможность увеличить разрешающую способность генератора в 128 раз, или 128*K, где K - может быть и больше. В итоге мы получим разрешающую способность - N*M, где N - длина таблицы, а M - количество отсчетов генератора. Круто короче.
В общем надо пробовать реализвать.
В любом случае - спасибо! Я как-то сразу сам не додумался до этого.
Это тот редкий случай когда обсуждение на форуме приводит к рождению нового эффективного алгоритма.

Умножать на матрицу поворота и нормировать на каждом такте. Для низких частот можно упростить считая, что sin(x) ~ x, cos(x) ~ 1. Можно не упрощать, если частота фиксированная то sin/cos считаются один раз.

так ошибка будет накапливаться.
по сравнению с вариантом от blackfin, можно еще одно умножение сэкономить если Sin[n+1] выразить через Sin[n] и Sin[n-1].
Sin[n+1] = Sin[n]*K - Sin[n-1];
да и со стабильностью вроде будет получше, не?
в целых числах если считать, то вообще ничего никуда не расползается, а с ошибками округления floatов - не знаю.

http://web.archive.org/web/20120308235958/...ors.html#Second Order Oscillator

гы, да ладно...
http://electronix.ru/forum/index.php?s=&am...t&p=1298367

ну и конечно, все расползается.

чукча - не читатель, чукча - писатель.

школьно-тригонометрический - да, а вот осциллятор второго порядка вроде как абсолютно стабильный должен быть с целочисленной арифметикой.

угу. всенепременно.

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

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

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


int a0, a1 = 0, a2 = 100, K=450;
for (int i = 0; i < 1000; i++){
a0 = (K*a1)/256 - a2;
a2 = a1;
a1 = a0;
printf("%d\n",a0);
}
через 765*N шагов, значения a0,a1,a2 опять станут -100, 0, 100, почему же он не развалился раз он не устойчивый?

Я посчитал. Устойчивость действительного генератора определяется только шумами квантования отсчетов,
в отличие от комплексного, где и квантованные отсчеты, и квантованные коэффициенты играют роль.



Хороший пример. Много подтверждает (limit cycles при единственном значении коэффициента).

Попробуйте так:
К = 45 * pi * pi;

а как же: "А для любой частоты при арифметике с конечной разрядностью это, увы, не выполняется." ?

45*pi*pi = 444
арифметика-то целая, в частоту немного не попадёт из-за округления, но работать будет стабильно:
0
100
173
200
173
100
0
-100
-173
-200
-173
-100
0

.

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

Если элементы матрицы поворота посчитать как настоящие sin/cos без приближения, то погрешности будут не намного больше чем в лобовом способе с вызовом sin/cos на каждом такте. Ничего не расползется, нормирование для этого есть.

Если же считать sin(x) ~ x, то будет погрешность задания частоты и как следствие уход фазы. Чем крупнее шаг, тем больше погрешность. По амплитуде погрешность постоянная и небольшая, и тоже зависит от шага.

Да и sin/cos всегда согласованы по фазе в отличии от методов раздельной их генерации.

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

Что-то мы скатились на обсуждение ерунды.

Ситуация с ф-й генератора вот какая.
Когда мы описываем те или иные методы, надо всегода смотреть на то, как данная ф-я будет использоваться.
Вот у меня задача достаточно общая:
есть приложение в котором используется DDC. Для этого DDC мне нужен генератор квадратуры с постоянной частотойк. Значение частоты может меняться, но не в реальном времени.
Т.е. получаем вот что. То что считается постоянно в реал-тайме критично ко времени. То что относится к расчету параметров - не критично. В итоге получается, что расчитать один раз sin/cos/step и т.д. это не проблема, но зато потом просто умножать нару раз и выбирать значения из таблицы.
В общем метод Blackfinа самое то. Сам генератор там состоит из вдух сложений и четырех умножений. А инициализация приращения и шаг таблицы расчитываются неспеша как угодно. Причем это "именно то" при условии что значения генератора реинициализируются через фиксированные промежутки (в моем случае 128 отсчетов).

Методы генерирования на основе расчета синуса и косинуса "налету" с помощью полиномов и т.д. имеют свое применение. Например когда у меня частота меняется в реальном времени и я ее должен постоянно подстраивать на основе ФАПЧ или что-то такое. В этом случае конечно надо делать генератор на основе прямого расчета синуса и косинуса.
То-же самое относится когда генератор должен работать на "бесконечном интервале времени". Тут шумы никуда не уйдут, и генератор просто "зазвенит".

В общем, побеждает как всегда компромис.

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

Компоненты вектора на единичной окружности дают всегда согласованные sin/cos, как они могут расползтись по фазе между собой? Удержание на единичной окружности обеспечивает нормирование. Итого дрейф может быть только по общей фазе/частоте, ну и амплитуда может стабилизироваться с небольшим смещением нормы от единицы.


У меня не звенит на бесконечном времени. Да и частоту тоже перестраивать пробовал в реальном времени.

а он и развалился. сразу. результат никак не гармоника.

а что же это тогда?
возьмите калькулятор и разницу меджу {0, 100, 173, 200, 173, 100, 0} и 200*sin(2*PI*i/12) посчитайте.


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

Конечно. Даже в вашем простом примере в диапазоне частот (К=1..511) амплитуда на выходе меняется на 4 разряда или в 16 раз.
Если отсчеты и коэффициенты принять 16 разрядными, то изменение амплитуды в зависимости от частоты будет уже в пределах 8 разрядов. Попытки как-то ограничить разрядность умножителя будут приводить к описываемым эффектам. То, что при полной разрядности арифметических блоков генератор будет стабилен, мы обсудили.

Зависимость амплитуды от частоты я уже даже за недостаток не считаю.

Амплитуда определяется начальными значениями a1 и a2, а не только K который задаёт частоту. и соответственно амплитуду можно сделать какой угодно, так что замечание вообще мимо.
что такое полная разрядность? вот выше пример при 9х9->(18 с округлением до 10) разряном умножении, он как, нестабилен?
давайте до 4х разрядов ограничим
int a0, a1 = 0, a2 = 10, K=12;
...
a0 = (K*a1)/16 - a2;
работать он от этого не перестанет.

зы с printfом то что не так?

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

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

Т.е. я считаю приращение по таблице для каждых 128 отсчетов, и матрицу поворота для каждого отсчета.
В итоге все работает вот как:
Когда я запускаю обработку 128 отсчетов, я загружаю начальные значения из таблицы (это константа), загружаю приращения синуса и косинуса (то-же константы расчитанные один раз при задании частоты) и потом считаю 128 раз поворот. После этого я делаю щаг по таблице, и сохраняю состояние указателя таблицы. В следующий раз повторяется все то-же самое.
В итоге, для таблицы в 1024 отсчета, и пакета 128 отсчетов, я могу генерировать синус с частотой кратной 128*1024 ( = 131072). Если нужно, то я могу переинициализировать не каждые 128 отсчетов, а каждые 1024, и получим мы кратность частоты 1024*1024, что в общем афигенно!
Соотношением сигнал/шум мы конечно жертвуем в некоторой степени, но "за все надо платить".

Да. у вас крутой генератор.
Снимаю чепчик...