Comments 16
Спасибо, отличная статья! Оптимизация алгоритма на ассемблере в cortex, редко такое встречаю, снимаю шляпу.
Прям DDS на максималках (вспомнился проект генератора DDS на AVR, там был хитрый быстрый алгоритм на ассемблере - использовалась таблица синуса (массив) с выравниванием по 0x100 (+ЦАП R2R), таким образом loop table получалась за счёт автоматического переполнения счётчика, но нужно было дополнять программу NOP-ами, чтобы массив был выровнен)
Табличные способы генерации отсчетов подходят, но не всегда. Все дело в том что начальная фаза сигнала не определена, а в процессе вычисления требуется гладко склеить предыдущий символ со следующим.
Вот этот момент не понял - а в чем здесь проблема? Мы же знаем, в каком месте таблицы мы находимся в данный момент. Нужно сменить частоту - переходим на соответствующее место другой таблицы, у которого то же мгновенное значение.
Это вопрос точности и размера таблицы.
Табличный метод предполагает использование фиксированных сдвигов в таблице. Хорошо, если частота дискретизации кратна частоте синтезируемого сигнала. Тогда, да, Вы просто меняете индекс в таблице с другим инкрементом.
Всегда ли это так?
Я не любитель усложнять себе жизнь, но, к сожалению, так бывает не всегда.
Для подобных случаев требуется аппроксимация между отсчетами таблицы, или вычисление sin/cos напрямую. Для рассчитанных целочисленных sinFixed, cosFixed, Google ИИ привел следующие статистические характеристики погрешности:
Для того чтобы получить подобные характеристики табличным методом, нужно иметь очень большую таблицу или аппроксимировать между отсчетами полиномом. Но тут речь о чистых и относительно медленных синусах и косинусах с динамическим диапазоном под ~190dB. Это была первая причина.
Вторая причина заключается в том, откуда Вы будете выполнять код и откуда брать табличные значения. Если предполагается брать отсчеты sin/cos и флешки, может оказаться что скорость ее работы ограничена WaitState-ами, а пара умножений и пара сложений для значений, которые не выходят за пределы регистров ядра, выполняются несколько быстрее чем чтение из Flash, особенно если у Вас есть возможность работать с SIMD. Речь о том, что в ряде случаев крутить NCO может быть выгоднее чем читать память. Все зависит от памяти конечно....
Кроме того, даже если Вы разместили табличку синуса в SRAM а не Flash, ту же память зачастую использует DMA, из той же памяти Вы берете те самые отсчеты, которые читаете или которые записываете. Конфликты на шине могут снизить производительность значительно сильнее чем кажется на первый взгляд, особенно если у Вас нет ITCM и DTCM, как например в STM32H750.
Да, кстати, анализ динамики погрешности синтеза отсчетов в данной реализации следующий:
Т.е. если Вам требуется получить непрерывный синус с большим количеством отсчетов, лучше вызывать синтезатор несколько раз. Например 32 раза по 1024 отсчета для получения выборки из 32768 отсчетов. Это не сильно уменьшит производительность, но позволит минимизировать погрешность. Почему так? Потому что при вызове функции она каждый раз рассчитывает опорный вектор, на что уходит порядка 64 тактов на квадратуру. Т.е. 130-140 тактов в сумме. Сам же поворот вектора быстрый, там всего два умножения и два сложения. Однако, он вносит погрешность. А затраты на расчет начального значения размазываются по расчету последующих отсчетов, и после выборки ~128 отсчетов получается порядка одного дополнительного такта на отсчет.
Собственно, за счет этого и выигрыш.
Спасибо за подробный ликбез! Да, если речь идет о высокоточном сигнале, то простой табличкой на сотню синусов не обойдешься, приходится думать хитрые алгоритмы и обходить коварные подводные камни.
между табличными значениями можно считать линейные отрезки. Тогда в таблице будет не только значение функции, но и её производная в этой точке, для линейной экстаполяции
Да. Если хорошо подумать, то метод который использую для опорных значений и есть интерполяция. По сути это таблица из 8 точек, если смотреть на него под другим углом. Но каждый из отрезков приводится к интервалу 0..PI/4.
На этом отрезке, в целых числах можно было сэкономить 2 бита. Перешел на формат Q0.32 для вычисления ряда Тейлора. Но можно было использовать другие полиномы, например кубическую интерполяцию между отсчетами с минимаксными коэффициентами.
Текущий код опороного sinFixed/cosFixed выглядит приблизительно так, если оформить его на C:
#include <stddef.h>
#include <stdint.h>
#include "FskModem.h"
#define PHASE_INDEX_0PI4_TO_1PI4 (0x000)
#define PHASE_INDEX_1PI4_TO_2PI4 (0x001)
#define PHASE_INDEX_2PI4_TO_3PI4 (0x002)
#define PHASE_INDEX_3PI4_TO_4PI4 (0x003)
#define PHASE_INDEX_4PI4_TO_5PI4 (0x004)
#define PHASE_INDEX_5PI4_TO_6PI4 (0x005)
#define PHASE_INDEX_6PI4_TO_7PI4 (0x006)
#define PHASE_INDEX_7PI4_TO_8PI4 (0x007)
#define PHASE_OFFSET_0PI4_TO_1PI4 (0x00000000)
#define PHASE_OFFSET_1PI4_TO_2PI4 (0x40000000)
#define PHASE_OFFSET_2PI4_TO_3PI4 (0x40000000)
#define PHASE_OFFSET_3PI4_TO_4PI4 (0x80000000)
#define PHASE_OFFSET_4PI4_TO_5PI4 (0x80000000)
#define PHASE_OFFSET_5PI4_TO_6PI4 (0xC0000000)
#define PHASE_OFFSET_6PI4_TO_7PI4 (0xC0000000)
#define PHASE_OFFSET_7PI4_TO_8PI4 (0x00000000)
#define PHASE_SCALE_PI4 (0xC90FDAA2)
inline uint64_t cortexm_umull( uint32_t a, uint32_t b ) {
register uint64_t umull_res = ((uint64_t)(a)) * ((uint32_t)(b));
return umull_res;
}
inline int64_t cortexm_smull( int32_t a, int32_t b ) {
register int64_t smull_res = ((int64_t)(a)) * ((int32_t)(b));
return smull_res;
}
inline int32_t cortexm_smull_q31( int32_t a, int32_t b ) {
register int64_t smull_q31_res = ((int64_t)(a)) * ((int32_t)(b));
register int32_t smull_q31_hi = smull_q31_res >> 32;
register int32_t smull_q31_low = smull_q31_res;
register int32_t res = ( smull_q31_hi << 1 ) | ( smull_q31_hi >> (32-1) );
return res;
}
inline int16_t cortexm_smull_q15( int32_t a, int32_t b ) {
register int64_t smull_q15_res = ((int64_t)(a)) * ((int32_t)(b));
register int32_t smull_15_hi = smull_q15_res >> 32;
register int32_t res = smull_15_hi >> 15;
return res;
}
inline uint32_t scale_phase_to_pi4( uint32_t phase ) {
register uint32_t angleScale_0_to_PI4 = PHASE_SCALE_PI4;
register uint64_t mul_res = cortexm_umull( phase, angleScale_0_to_PI4);
register uint32_t mul_hi = mul_res >> 32;
register uint32_t mul_low = mul_res;
register uint32_t res = ( mul_hi << 3 ) | ( mul_low >> (32-3) );
return res;
}
uint32_t sinFixedPI4( uint32_t phase ) {
register uint32_t angleScale_0_to_PI4 = PHASE_SCALE_PI4;
register uint32_t phaseScaled = scale_phase_to_pi4( phase );
register uint32_t phaseQuad = cortexm_umull( phaseScaled, phaseScaled ) >> 32;
register uint32_t sum = phaseScaled;
register uint32_t x = cortexm_umull( phaseQuad, phaseScaled ) >> 32;
sum -= cortexm_umull( x, 0x2AAAAAAA ) >> 32;
x = cortexm_umull( x, phaseQuad ) >> 32;
sum += cortexm_umull( x, 0x02222222 ) >> 32;
x = cortexm_umull( x, phaseQuad ) >> 32;
sum -= cortexm_umull( x, 0x000D00D0 ) >> 32;
x = cortexm_umull( x, phaseQuad ) >> 32;
sum += cortexm_umull( x, 0x00002E3B ) >> 32;
x = cortexm_umull( x, phaseQuad ) >> 32;
sum -= cortexm_umull( x, 0x0000006B ) >> 32;
return sum;
}
uint32_t cosFixedPI4( uint32_t phase ) {
register uint32_t angleScale_0_to_PI4 = 0xc90fdaa2;
register uint32_t phaseScaled = scale_phase_to_pi4( phase );
register uint32_t phaseQuad = cortexm_umull( phaseScaled, phaseScaled ) >> 32;
register uint32_t sum = 0 - cortexm_umull( phaseQuad, 0x80000000 );
register uint32_t x = cortexm_umull( phaseQuad, phaseQuad ) >> 32;
sum += cortexm_umull( x, 0x0AAAAAAA ) >> 32;
x = cortexm_umull( x, phaseQuad ) >> 32;
sum -= cortexm_umull( x, 0x005B05B0 ) >> 32;
x = cortexm_umull( x, phaseQuad ) >> 32;
sum += cortexm_umull( x, 0x0001A01A ) >> 32;
x = cortexm_umull( x, phaseQuad ) >> 32;
sum -= cortexm_umull( x, 0x0000049F ) >> 32;
x = cortexm_umull( x, phaseQuad ) >> 32;
sum -= cortexm_umull( x, 0x00000008 ) >> 32;
sum -= 1;
return sum;
}
int32_t sinFixedQ31( uint32_t phase ) {
register uint32_t phaseSelector = phase >> 29;
switch ( phaseSelector ) {
case PHASE_INDEX_0PI4_TO_1PI4: return ( (int32_t)(sinFixedPI4( phase - PHASE_OFFSET_0PI4_TO_1PI4 ) / 2 ) );
case PHASE_INDEX_1PI4_TO_2PI4: return ( (int32_t)(cosFixedPI4( PHASE_OFFSET_1PI4_TO_2PI4 - phase ) / 2 ) );
case PHASE_INDEX_2PI4_TO_3PI4: return ( (int32_t)(cosFixedPI4( phase - PHASE_OFFSET_2PI4_TO_3PI4 ) / 2 ) );
case PHASE_INDEX_3PI4_TO_4PI4: return ( (int32_t)(sinFixedPI4( PHASE_OFFSET_3PI4_TO_4PI4 - phase ) / 2 ) );
case PHASE_INDEX_4PI4_TO_5PI4: return -( (int32_t)(sinFixedPI4( phase - PHASE_OFFSET_4PI4_TO_5PI4) / 2 ) );
case PHASE_INDEX_5PI4_TO_6PI4: return -( (int32_t)(cosFixedPI4( PHASE_OFFSET_5PI4_TO_6PI4 - phase ) / 2 ) );
case PHASE_INDEX_6PI4_TO_7PI4: return -( (int32_t)(cosFixedPI4( phase - PHASE_OFFSET_6PI4_TO_7PI4 ) / 2 ) );
case PHASE_INDEX_7PI4_TO_8PI4: return -( (int32_t)(sinFixedPI4( PHASE_OFFSET_7PI4_TO_8PI4 - phase ) / 2 ) );
}
return 0;
}
int32_t cosFixedQ31( uint32_t phase ) {
return sinFixedQ31( phase - PHASE_OFFSET_2PI4_TO_3PI4 );
}
Но если компилировать его как есть, то компилятор добавляет некоторое количеств инструкций умножения.
Причëм все производные синуса в таблице уже есть)
Для того чтобы получить подобные характеристики табличным методом, нужно иметь очень большую таблицу или аппроксимировать между отсчетами полиномом.
Знакомы с таким методом? Можно найти компромис между точностью, размером таблиц и количеством вычислений. На cortex-M3, где есть DSP инструкции, вычисляет довольно шустро, ссылка на примеры там присутствует.
хочу такой хабор, как эта статья ✨
Что только не придумают, лишь бы не возиться с плис. (шутка)
Как же давно я ничего не писал на Ассемблере и не разу не писал на нём ничего для STM32!
А ведь были времена писал программы в МАШИННЫХ КОДАХ под PDP11 !
Теперь только С и в большинстве проектов + RTOS.
В 2008 решал подобную задачу на LPC ещё. Но она была попроще - надо было чистый синус получать но в широком спектре частот. Пришлось поизвращаться но таки удалось обойти без ассемблера
Кстати до сих пор в ходу в промышленной автоматизации HART модемы с FSK модуляцией и стоит эта экзотика бешенных денег!
Считать точнее, чем разрешение ЦАП/АЦП смысла не вижу. Тем более, что частота, особенно на приëме не постоянная, и инкремент фазы всë время меняется.
Поэтому я всë-таки обычно использую таблицу, если надо - с интерполяцией линейной или квадратичной. Производные-то в таблице синуса автоматически есть.
Если Вы синтезируете сигнал для DAC или свертки с ADC, да особого смысла нет. Но в некоторых случаях смысл есть.
1. Представьте что Вам требуется принять сигнал, который утоплен в шум настолько глубоко, что влияет только на вероятностные характеристики младшего разряда ADC. В этом случае тепловой шум подталкивает младший разряд туда-сюда равновероятно. Однако за счет того что Вы умножаете сигнал на масштабируемый sin(x)/cos(x), фактически, вы работаете с масштабированием вероятности. В этом случае Вы можете собрать сигнал, уровень которого ниже чем разрядная сетка ADC. Однако, это возможно только в случае стохастических процессов. В этом случае динамически диапазон не будет лишним.
2. Представьте что Вы обрабатываете аудио, и в процессе обработки что-то теряете на округлении. Или наоборот сглаживаете сигнал, раскладывая его по определенному базису, который позволяет сделать его более правильным в последствии.
В следующем сообщении продолжу...
О чем я говорю (продолжение ответа).
Представьте что у Вас в наличии есть компакт-диск.
Это 44100 и 16 бит. И тут возникает вопрос: А можно ли получить качество аудио-сигнала лучше чем исходное?
В общем и целом, кажется что нельзя, но в некоторых случаях, наоборот, можно.
Как это сделать?
На самом деле если Вы используете компакт-диск музыкой записанной на настоящих музыкальных инструментах, то нет ничего невозможного. Вы понимаете что композиция состоит из:
1. Транзиентов (это момент, когда бьют по ударному инструменту).
2. Типовой схемы работы музыкальных инструментов (атака, стабильная фаза, затухание).
3. Гармоник музыкальных инструментов.
4. Голоса, у которого существует модель голосовых связок.
Для анализа музыки, требуется сделать несколько процедур:
1. Оценка мощности (дисперсия).
2. Оценка спектра (FFT).
3. Оценка транзиентов (удар по ударному инструменту с последующей затухающей по экспоненте синусойде).
4. Синтез звука для последующей свертки.
Если Вы сможете качественно разбить звук на частотные составляющие, и на затухающие квадратурные составляющие, то можете его отфильтровать. Причем это будет не FIR/IIR фильтр, который оставляет за собой хвост или формирует пред-эхо перед ударом барабана. Это реальное одностоннее преобразование (свертка с затухающий функцией).
Для этих целей, требуется синус более гладкий чем исходный, поскольку с ним Вы будете коррелировать и его будете восстанавливать. А дальше, после синтеза, остается только вычесть тот шум который остался.
И вот Вам либо улучшайзер музыки, который если правильно сделать будет восстанавливать то что потерялось, либо дополнительная стенография. Т.е. если Вы после такой обработки добавите в музыку на компакт-диске шум математическое ожидание которого стремится к нулю, а дисперсия занимает 1-2 младших разряда, то в этом шуме можете кодировать свою информацию, но или просто работать с ним как с шумом. Это может быть совместимый формат старых CD, которые на одной аппаратуре просто проигрываются, а на другой содержат текст караоке например, и при этом остаются совместимыми, и не теряют качество. А может и улучшают его.
Да, в том применении где я показал, больше чем 16 отсчетов вроде-бы не требуется, но в принципе, 24-32 бита в целых числах это хорошо. Но тогда аккумулятор тоже нужно брать не 32-бита, а 64 для того чтобы копить дробную часть.
Модуляция CPFSK и не только в квадратурах на микропроцессорах STM32