[Allregia 9]

Привет.

Генерирую синус DDS-ом. Не по таблице а счетом (хотя потом может перенесу и в таблицу, но пока быстордействия хватает и посчитать).

На каждый семпл вызывется:

	        data = (int32_t)( 0x7FFFFFFF * arm_sin_f32((PI*pDDS)/0x7FFFFFFF));
        pDDS+=dDDS;
	

 

Счетчик в аргументе (pDDS)  и приращение (dDDS)- uint32_t, data - int32_t.

dDDS - расчитывается заранее, в зависимости от частоты и семплрета.

Все прекрасно работает, за исключением одного - в сигнале присутствует зеркалка, частотой Fsample-Fsignal, и чем Fsignal ближе к Fsample/2, тем выше ее амплитуда (с приближениек к половине семлрейта - они уже чуть ли не сравниваются).

пример: Fs=48kHz, частота сигнала 20кгц - алиас на 28кгц амплитудой примерно -100дб от основной, частоа 23кгц - алиас на 25кшц с уровнем около -27дБ.

Вопрос - что я делаю "не так", или оно так и должно быть?

 

P.S. Хм, похоже так и должно быть - проверил на чужих программах, генерирующи сигналы, там тоже самое....

 

 

 

[MrBearManul 0]

2 часа назад, Allregia сказал:

Вопрос - что я делаю "не так", или оно так и должно быть?

Банальный вопрос: у вас фильтр низких частот после ЦАПа имеется с частотой среза 48/2=24 кГц? Наверное, будет лучше, если этот фильтр активный.

[jcxz 184]

4 часа назад, Allregia сказал:

На каждый семпл вызывется:

	        data = (int32_t)( 0x7FFFFFFF * arm_sin_f32((PI*pDDS)/0x7FFFFFFF));
        pDDS+=dDDS;
	

 

Зачем так громоздко? Да ещё с float-ом. Для вычисления синуса float не нужен. 

Кстати и выражение у вас неправильное. Должно быть то, что в скобках: (PI*pDDS)/0x80000000u. 

Возможно из-за этого также могут быть косяки в спектре.

А ещё чуть оптимальнее (возможно) будет: pDDS * (float)(PI / 0x80000000u).

 

 

4 часа назад, Allregia сказал:

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

Вообще-то генерация синусоиды на ARM/DSP по таблице как правило медленнее чем вычислением. Конечно - если с умом делать.

Причём - может даже кратно медленнее.

[jcxz 184]

4 часа назад, Allregia сказал:

P.S. Хм, похоже так и должно быть - проверил на чужих программах, генерирующи сигналы, там тоже самое....

https://nag.ru/articles/article/103332/teorema-kotelnikova-dlya-chaynikov-prostyimi-slovami.html   

[V_G 8]

8 часов назад, Allregia сказал:

Вопрос - что я делаю "не так", или оно так и должно быть?

Так и должно быть: DDS генерирует лес палок, нужную для себя надо отфильтровывать. Для оценки этого леса существует софт, например, AD предлагает ADIsimDDS. Он, правда, нацелен на синтезаторы от AD, но если понять принцип, можно попробовать увидеть свой лес...

[Allregia 9]

6 hours ago, MrBearManul said:

Банальный вопрос: у вас фильтр низких частот после ЦАПа имеется с частотой среза 48/2=24 кГц? Наверное, будет лучше, если этот фильтр активный.

В звуковых картах они имеются?   Сейчас они для проверки исользуются - одна для генерации на 48кгц, другая для анализа, на 192кгц.

Могу попробовать без железа, без преобращзования в аналог и обратно, на компе.

 

5 hours ago, jcxz said:

Зачем так громоздко? Да ещё с float-ом. Для вычисления синуса float не нужен. 

Кстати и выражение у вас неправильное. Должно быть то, что в скобках: (PI*pDDS)/0x80000000u. 

Возможно из-за этого также могут быть косяки в спектре.

А ещё чуть оптимальнее (возможно) будет: pDDS * (float)(PI / 0x80000000u).

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

Флоатом оказалось самое оптимальное, осенно при последующем умножении на амплитуду.

Можно и не флоатом, если исползовать другую функцию - arm_sin_q31(pDDS/2), но тогда придется амплитуду менять не одним умножением, а переводом в int64 и делением на коэффициент.

Про pDDS * (float)(PI / 0x80000000u) не очень понял. Зачем результат деления PI на константу переводить во флоат, если PI и так флоат и результат деления естественно тоже?

 

 

Quote

 

 

Вообще-то генерация синусоиды на ARM/DSP по таблице как правило медленнее чем вычислением. Конечно - если с умом делать.

Причём - может даже кратно медленнее.

Как одна выборка из таблицы может быть быстрее вычисления синуса, даже скоростными методами (в которых внутри тоже выборка из таблицы есть)?

 

4 hours ago, jcxz said:

https://nag.ru/articles/article/103332/teorema-kotelnikova-dlya-chaynikov-prostyimi-slovami.html   

А причем тут ТК? Я-же  выше Fs/2 не подымаюсь.

[_sda 0]

2 часа назад, Allregia сказал:

А причем тут ТК? Я-же  выше Fs/2 не подымаюсь.

Ваша частота находится в первой зоне Найквиста, а её образ во второй. Его вы и видите на Fsample-Fsignal. Для вещественного сигнала это нормально.

Ставьте аналоговый фильтр на Fs/2 для удаления образов.

[andyp 9]

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

[V_G 8]

1 час назад, _sda сказал:

Ваша частота находится в первой зоне Найквиста, а её образ во второй. Для вещественного сигнала это нормально.

Кстати, да: я использую AD9859 c максимальной тактовой 400 МГц для генерации частот порядка 190 МГц. И вместо фильтрации зеркальной компоненты 210 МГц при тактовой 400 МГц предпочел вторую зону с тактовой 240 МГц: основная частота получается 50 МГц, нужная мне 190 МГц и далее палка на тактовой 240. 190 и 240 при этом разделить значительно проще, чем 190 и 210.

Топикстартеру, если ему нужно генерить звуковой диапазон частот, можно посоветовать лишь повысить тактовую, например, до 96 кГц, и использовать хороший ФНЧ на выходе. Если же нужно генерить узкий спектр в районе 20 кГц, то наоборот, стоит понизить тактовую и выфильтровывать полезный сигнал из второй зоны полосовым фильтром.

[jcxz 184]

2 часа назад, Allregia сказал:

Флоатом оказалось самое оптимальное, осенно при последующем умножении на амплитуду.

Можно и не флоатом, если исползовать другую функцию - arm_sin_q31(pDDS/2), но тогда придется амплитуду менять не одним умножением, а переводом в int64 и делением на коэффициент.

Не совсем понял, что Вы тут имели в виду....  

Если Вы про умножение на 0x7FFFFFFF результата функции синуса, то нет никаких проблем сделать то же самое на fixed point:

result = (int32)((int64)0x7FFFFFFF * (int64)arm_sin_q31(...) >> 32) << 1;

На ARM это должно скомпилиться всего в 2 команды: SMMUL, LSL.

 

Цитата

Про pDDS * (float)(PI / 0x80000000u) не очень понял. Зачем результат деления PI на константу переводить во флоат, если PI и так флоат и результат деления естественно тоже?

Я не знаю что такое у Вас PI, так как Вы не привели его определение. Поэтому привёл явно.

У меня например аналогичная константа определена как: #define M_PI 3.14159265358979323846

Это лучше чем float. Сами подумайте почему.

 

Цитата

Как одна выборка из таблицы может быть быстрее вычисления синуса, даже скоростными методами (в которых внутри тоже выборка из таблицы есть)?

Во-первых: внутри функции вычисления синуса не факт что есть какие-то "выборки из таблицы". Обычно внутри (имхо) вычисление идёт полиномом. Для этого таблица не нужна.

Во-вторых: А кто говорил про "вычисление синуса"? Для генерации синусоиды с постоянным фазовым шагом достаточно знания школьного курса тригонометрии и нужна всего одна MAC-команда на один генерируемый сэмпл. Которая выполняется за 1 такт CPU. А даже одна выборка из таблицы (во флешь или SDRAM) по произвольному смещению может занимать множество тактов CPU. И не нужно считать синус на каждом шаге.

[Allregia 9]

1 hour ago, jcxz said:

Не совсем понял, что Вы тут имели в виду....  

Если Вы про умножение на 0x7FFFFFFF результата функции синуса, то нет никаких проблем сделать то же самое на fixed point:

result = (int32)((int64)0x7FFFFFFF * (int64)arm_sin_q31(...) >> 32) << 1;

На ARM это должно скомпилиться всего в 2 команды: SMMUL, LSL.

Я с q31 так и делал.

 

1 hour ago, jcxz said:

 

Я не знаю что такое у Вас PI, так как Вы не привели его определение. Поэтому привёл явно.

У меня например аналогичная константа определена как: #define M_PI 3.14159265358979323846

Это лучше чем float. Сами подумайте почему.

Это дабл, у меня проц толко с сингл  FPU.

 

1 hour ago, jcxz said:

 

Во-первых: внутри функции вычисления синуса не факт что есть какие-то "выборки из таблицы". Обычно внутри (имхо) вычисление идёт полиномом. Для этого таблица не нужна.

Ну достаточно посмотреть в исходники arm_sin_XXX:

 * Computes the trigonometric sine function using a combination of table lookup
 * and linear interpolation.  There are separate functions for
 * Q15, Q31, and floating-point data types.

	float32_t arm_sin_f32(
  float32_t x)
{
  float32_t sinVal, fract, in;                           /* Temporary variables for input, output */
  uint16_t index;                                        /* Index variable */
  float32_t a, b;                                        /* Two nearest output values */
  int32_t n;
  float32_t findex;
	  /* Special case for small negative inputs */
  if ((x < 0.0f) && (x >= -1.9e-7f)) {
     return x;
  }
	  /* input x is in radians */
  /* Scale the input to [0 1] range from [0 2*PI] , divide input by 2*pi */
  in = x * 0.159154943092f;
	  /* Calculation of floor value of input */
  n = (int32_t) in;
	  /* Make negative values towards -infinity */
  if (x < 0.0f)
  {
    n--;
  }
	  /* Map input value to [0 1] */
  in = in - (float32_t) n;
	  /* Calculation of index of the table */
  findex = (float32_t) FAST_MATH_TABLE_SIZE * in;
	  index = ((uint16_t)findex) & 0x1ff;
	  /* fractional value calculation */
  fract = findex - (float32_t) index;
	  /* Read two nearest values of input value from the sin table */
  a = sinTable_f32[index];
  b = sinTable_f32[index+1];
	  /* Linear interpolation process */
  sinVal = (1.0f-fract)*a + fract*b;
	  /* Return the output value */
  return (sinVal);
}
	q31_t arm_sin_q31(
  q31_t x)
{
  q31_t sinVal;                                  /* Temporary variables for input, output */
  int32_t index;                                 /* Index variables */
  q31_t a, b;                                    /* Four nearest output values */
  q31_t fract;                                   /* Temporary values for fractional values */
	  /* Calculate the nearest index */
  index = (uint32_t)x >> FAST_MATH_Q31_SHIFT;
	  /* Calculation of fractional value */
  fract = (x - (index << FAST_MATH_Q31_SHIFT)) << 9;
	  /* Read two nearest values of input value from the sin table */
  a = sinTable_q31[index];
  b = sinTable_q31[index+1];
	  /* Linear interpolation process */
  sinVal = (q63_t)(0x80000000-fract)*a >> 32;
  sinVal = (q31_t)((((q63_t)sinVal << 32) + ((q63_t)fract*b)) >> 32);
	  return sinVal << 1;
}
	

 

1 hour ago, jcxz said:

Во-вторых: А кто говорил про "вычисление синуса"? Для генерации синусоиды с постоянным фазовым шагом достаточно знания школьного курса тригонометрии и нужна всего одна MAC-команда на один генерируемый сэмпл. Которая выполняется за 1 такт CPU. А даже одна выборка из таблицы (во флешь или SDRAM) по произвольному смещению может занимать множество тактов CPU. И не нужно считать синус на каждом шаге.

Надо подумать, Я пока не очень понял.

 

Но пока интереснее другое - беру программу, генерирующую 20кгц при 48кгц семплрейте, и сам делаю такое-же. Подаю каждый синал на свой ЦАП. Оцифровываю одновременно 2-х канальным АЦП с семплрейтом 96 или 192 (это не влияет) и смотрю спектр

У них зеркалка на 28кгц с уровнем -90дБ а у меня  -20дБ.  Железо практичски то-же самое. Уровень 2-й гармоники примерно одинаковый, -110дБ, т.е. это не искажения ЦАП/АЦП. Чего они могут делать по другуму?

[jcxz 184]

1 час назад, Allregia сказал:

Это дабл, у меня проц толко с сингл  FPU.

Ну и что? Компилятор ведь у вас с double наверное?  

1 час назад, Allregia сказал:

Ну достаточно посмотреть в исходники arm_sin_XXX:

Как-то очень громоздко.....  Полиномом гораздо быстрее будет.

1 час назад, Allregia сказал:

Надо подумать, Я пока не очень понял.

https://ru.wikipedia.org/wiki/Тригонометрические_тождества

формула 6.2. 

1 час назад, Allregia сказал:

Чего они могут делать по другуму?

Может у вас ещё где-то баг в программе...

[Allregia 9]

1 hour ago, jcxz said:

Ну и что? Компилятор ведь у вас с double наверное?  

 

Разумеется :)

Но это толко если этот расчет попадет не в ран-0тайм а в компайл-тайм.

 

1 hour ago, jcxz said:

Как-то очень громоздко.....  Полиномом гораздо быстрее будет.

https://ru.wikipedia.org/wiki/Тригонометрические_тождества

формула 6.2. 

 

Произведение синуса на косинус? И как это поможет?

 

[jcxz 184]

23 минуты назад, Allregia сказал:

Но это толко если этот расчет попадет не в ран-0тайм а в компайл-тайм.

(PI / 0x80000000u) - это выражение уже считается компилятором в double с соответствующей точностью, и только результат конвертируется в float.

А ведь могут быть и гораздо более сложные выражения. Их лучше считать в double. для бОльшей точности. Тем более что это ничего не стоит. А уже результат (который идёт в МК) только переводить во float. Поэтому я и написал такое выражение и поэтому у меня такое M_PI.

 

Цитата

Произведение синуса на косинус? И как это поможет?

Ну подумайте же немного!

 

Выразите sin(a+b) из той формулы:

sin(a+b) = 2*sin(a)*cos(b)-sin(a-b)

где: a - некий начальный угол; b - угловое приращение на каждом шаге (на сэмпл).

Так как b - константа, то 2*cos(b) - тоже. Вычисляем его заранее и присваиваем K: K=2*cos(b)

получаем:

sin(a+b) = K*sin(a)-sin(a-b)

sin(a+b) - это значение выражения, которое надо вычислить на некоем шаге i. Назовём его Y(i).

А значит sin(a) - это значение того же самого выражения на предыдущем шаге (i-1). Или: Y(i-1).

А sin(a-b) - это то же значение на шаге (i-2): y(i-2).

Итого: Y(i)=K*Y(i-1)-Y(i-2)

Теперь считаем стартовые значения: Y(i-1), Y(i-2) обычным способом (функцией синуса), запоминаем.

И теперь для генерации каждого последующего значения: Y(i), Y(i+1), Y(i+2), ... нужна только одна MAC-команда на отсчёт.

 

Элементарно же!  

[andyp 9]

47 minutes ago, jcxz said:

Элементарно же! 

Угу. Что там со стабильностью у такого генератора? Слыхал, что у него полюса прямиком на единичной окружности ;)