Микроконтроллер + DAC + табличный синус

В этом документе пойдёт речь об использовании микроконтроллера для прямого цифрового синтеза сигнала (DDS, direct digital synthesis). Например, с помощью микроконтроллера STM32F100RB с тактовой частотой 24 МГц или аналогичного ему, можно синтезировать сигнал с произвольной частотой от 0 до 100 кГц, а при оптимизации метода по скорости и до 500 кГц, задавая требуемое значение частоты с точностью до долей герца и даже до тысячных далей (при наличии достаточно стабильного тактового генератора для микроконтроллера, конечно). При этом перестройка частоты осуществляется крайне просто, её можно сделать очень плавной, при необходимости - очень быстрой.

Используя специализированные микросхемы для DDS, можно получить сигналы с гораздо более высокими частотами (до сотен МГц, например: AD9858), зато микроконтроллер более гибок и обеспечивает синтез сигналов чисто программными средствами, т.е. "даром".

Общие вопросы прямого цифрового синтеза сигнала (DDS)

В статье "Микроконтроллер + DAC как синтезатор сетки частот" рассматривались некоторые теоретические вопросы прямого цифрового синтеза сигнала (DDS), в частности, был исследован спектр синтезированного синусоидального сигнала и было показано, на основе опорной частоты f₀ можно получить сигнал с частотой $$ f=f_0 \frac m n, \tag {*} $$ где m, n - любые целые числа при условии, что m<n/2 (на практике требуется хотя бы трёхкратное отличие). Для синтеза синусоиды с любой требуемой частотой, на цифро-аналоговый преобразователь (DAC) с частотой f₀ подаются значения вида \(a\sin\frac{2\pi mi}n,\) где i - номер точки, в которой выполняется цифро-аналоговое преобразование, a - амплитуда сигнала.

В соответствии с формулой (*), частотой синтезируемого сигнала можно управлять тремя способами: изменяя опорную частоту f₀; изменяя величину m, определяющую количество полных периодов в синтезируемом отрезке или изменяя общее количество точек n, приходящихся на синтезируемый отрезок сигнала. Лучше всего для управления частотой сигнала изменять величину m, что даёт возможность линейной перестройки частоты от 0 до максимального возможного значения при данной f₀.

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

Величина n плохо подходит для управления частотой: зависимость частоты от n нелинейная, кроме того, уменьшая n, мы уменьшаем разрешение синтезатора по частоте, так что выгодно выбрать максимально возможное значение n и тем самым получить максимальное разрешение. При условии, что используются 32-разрядные целые числа, таким значением будет n=2³², которому соответствует точность задания частоты сигнала 2^-32 или около 2*10^-10 от опорной частоты. Это более чем достаточное разрешение для большинства случаев, так как нестабильность даже хорошего генератора, стабилизированного кварцем, на несколько порядков выше.

Способы расчёта значений сигнала

Будем рассматривать синтез синусоидального колебания. Естественно, в случае необходимости, могут быть синтезированы и прочие сигналы, с учётом ограничений на частоту высших гармоник со стороны используемых аппаратных средств.

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

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

Табличный метод синтеза в своём крайнем варианте, когда значения заранее вычисляются для всех n точек синтезируемого сигнала, практически полностью разгружает процессор от вычислений. Особенно если задействовать DMA для пересылки данных в DAC. Однако, данный способ является затратным в отношении используемой памяти. Для n=2³² он вообще становится неприменимым, поскольку микроконтроллер с объёмом памяти в 4 Гслов пока ещё не самое рядовое явление.

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

Оптимизация табличного метода синтеза

Самое очевидное и неправильное решение задачи синтеза с использованием табличного метода состоит в использовании таблицы с размещением элементов в том порядке, в котором они загружаются в DAC, т.е. когда элемент с индексом i содержит значение синуса от аргумента \(\frac{2\pi m}n i\). Такая таблица требует пересчёта при каждом изменении m, а это никуда не годится - любое изменение частоты генерируемого сигнала будет сопровождаться чудовищными задержками на пересчёт таблицы. Совсем другое дело, если мы составим таблицу из n элементов, в которой элемент k содержит значение синуса в точке \(\frac{2\pi}n k\), т.е. таблицу значений синуса для разных фаз в пределах только одного периода. Таблица ставит в соответствие целочисленному индексу фазы k значение синуса для соответствующей этому индексу фазы \(\phi=\frac{2\pi}n k\). Для того, чтобы найти нужное нам значение \(\sin\frac{2\pi m}n i\), просто берём из таблицы значение элемента с индексом m*i, точнее с индексом, равным остатку от деления m*i на n, так как само произведение может быть больше n, а размер таблицы равен n.

В самом деле, пусть при делении числа m*i нацело на n, получается частное - целое число p и остаток - целое число q, меньшее n. Это можно записать в следующем виде
m*i=p*n+q, 0<=q<n.
Тогда интересующее нас значение сигнала в точке i, с учётом периодичности синуса, равно $$ \sin\frac{2\pi}n(mi)=\sin\left(\frac{2\pi}n pn+\frac{2\pi}n q\right)=\\ \sin\left(2\pi p+\frac{2\pi}n q\right)=\sin\frac{2\pi}n q, $$ т.е. равно значению элемента таблицы с индексом q=(m*i)%n.

Причём вычислять остаток от деления на практике не требуется, для перехода от одного элемента таблицы к следующему, достаточно прибавить к текущему индексу число m, если результат больше или равен n, вычитаем из него n и получаем следующий индекс. В случае 32-битовых регистров и n=2³² алгоритм ещё проще: на каждом шаге добавляем к текущему индексу величину m, в случае переполнения бит переноса будет автоматически отброшен и получен результат по модулю 2³², т.е. не превосходящий n.

Теперь пересчитывать таблицу при изменении частоты не требуется, нужно лишь изменить шаг приращения индекса, равный m. Сразу же начнётся генерация сигнала с новой частотой f=f₀*m/n.

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

Перейдём к уменьшению таблицы значений до разумных размеров. Легко заметить, что при большом количестве n элементов в таблице, разница между фазами для соседних элементов очень мала, а значения функции, соответственно, очень близки. С учётом ограниченной разрядности используемого DAC, множеству соседних элементов в таблице будет соответствовать одно и то же загружаемое в цифро-аналоговый преобразователь значение. А значит таблица избыточна и может быть "прорежена" для уменьшения количества одинаковых элементов. Легко показать, что в абсолютных величинах изменение функции sin не превышает вызвавшего его изменения аргумента: $$ \Delta y=|\sin(\phi+\Delta\phi)-\sin\phi| \le |\Delta\phi|. $$ Предположим, что мы синтезируем синусоиду, используя DAC с однополярным выходом, напряжение на котором может быть в пределах от 0 до V_REF. В таком случае необходимо наличие постоянной составляющей в сигнале, чтобы сигнал не уходил в область отрицательных значений: $$ u=a_0+a_1\sin\phi, $$ где \(\phi\) - фаза. Размах сигнала не может превышать диапазона выходных напряжений, максимальная амплитуда вдвое меньше размаха, \(a_1 \le V_{REF}/2\). Ошибка, обусловленная погрешностью табличного определения sin из сокращённого варианта таблицы $$ \Delta u=a_1 \Delta y \le a_1 \Delta\phi \le \frac{V_{REF}}2 \Delta\phi, $$ наибольшее значение ошибки тогда: $$ \Delta u_{max} = \frac{V_{REF}}2 \Delta\phi. $$ Допустим, требуется, чтобы ошибка вычислений не превышала половины значения младшего бита DAC (1/2 LSB), то $$ \Delta u_{max} \le \frac 1 2 \frac{V_{REF}}{2^N}, $$ тогда $$ \Delta \phi \le \frac 1 {2^N}. $$ Для того, чтобы заполнить таблицу значениями от 0 до 2π с найденным нами шагом, потребуется L элементов: $$ L=\frac{2\pi}{\Delta \phi}=2\pi 2^N \approx 2^{N+3}. $$ Итак, если используется DAC с разрядностью N, нет смысла делать таблицу с количеством элементов, превосходящим 2^N+3..2^N+4. Размер выбирается равным степени двойки, в этом случае индексация таблицы оказывается наиболее простой - для получения индекса в таблице достаточно индекс фазы сдвинуть вправо на (32-Nt), где Nt определяет размер таблицы, L=2^Nt. При N=12 размер таблицы составит 32768 или 65536 элементов, что уже вполне реализуемо с использованием обычных микроконтроллеров.

Таблицу легко уменьшить ещё в 4 раза с учётом повторяемости значений синуса по четвертям: $$ \text{если } \pi \le x \lt 2\pi, \\ \sin x=-\sin(2\pi-x), \text{ причём, } 0 \lt 2\pi-x \le \pi, \\ \text{если } \pi/2 \le x \lt \pi, \\ \sin x=\sin(\pi-x), \text{ причём, } 0 \lt \pi-x \le \pi/2. \\ $$ Размер такой таблицы составит 8192+1 или 16384+1 элементов, но потребуются дополнительные вычисления, если получен индекс за пределами укороченного варианта таблицы (размер таблицы L/4+1, так как мы должны включить в таблицу элемент с индексом L/4, соответствующий фазе π/2).

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

Табличный метод и интерполяция

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

Здесь не будем приводить расчётов по оценке точности метода интерполяции, подробно эта задача рассматривается в документе "Табличное вычисление синуса с использованием линейной интерполяции", где показано, что при шаге аргумента в таблице, равном d, получаем точность вычисления синуса любого аргумента не хуже d²/8. Для получения точности вычисления значений сигнала не хуже 1/2 LSB, потребуется таблица размером не менее 2.22*2^N/2.

Так, использование 12-разрядного DAC потребует таблицы синусов размером 256 элементов, причём благодаря повторяемости значений синуса по четвертям, её можно уменьшить ещё в 4 раза (за счёт использования формул приведения), в данном примере это будет 64+1 элементов, что в 128..256 раз меньше размера, необходимого без применения интерполяции (как описано в предыдущем пункте)! Такая крошечная таблица уже по силам практически любому микроконтроллеру.

Использование усечённой таблицы синусов

Рассмотрим, каким образом определить значение элемента с индексом t в полной таблице синусов (размера L), имея доступ только к первым L/4+1 элементам таблицы. Это позволит нам использовать усечённый вариант таблицы, который почти в 4 раза меньше по размеру.

Индексу t=0..L-1 соответствует аргумент \(\phi=2\pi t/L\), например, индексу t=0 соответствует аргумент 0, t=L/4 - аргумент π/2, t=L/2 - аргумент π, t=3*L/4 - аргумент 3π/2. Изменение аргумента при переходе от элемента с индексом t таблицы к следующему элементу t+1, очевидно, равно \(d=2\pi/L\). Это означает, что добавить/вычесть из аргумента угол \(\Delta\phi\) эквивалентно перемещению по таблице на \(\Delta\phi/d=\frac{\Delta\phi}{2\pi}L\) в ту или иную сторону. Пользуясь перечисленными соотношением, мы можем записать формулы приведения в применении к индексам таблицы, а формулы приведения позволяют нам "понизить" значение аргумента.

Допустим, $$ L/2 \le t \lt L \text{ или } \pi \le \phi \lt 2\pi, \\ $$ тогда $$ \sin \phi=-\sin(2\pi-\phi), 0 \lt 2\pi-\phi \le \pi, $$ что в индексной форме можно записать как $$ \sin \frac{2\pi t}L=-\sin\left(2\pi-\frac{2\pi t}L\right)= -\sin\left(\frac{2\pi}L\left(L-t\right)\right), \\ \text{где } 0 \lt (L-t) \le L/2, $$ т.е. можно "отразить" значения из второй половины таблицы на первую половину, для чего следует взять из таблицы элемент с противоположным знаком элемент с индексом L-t. Проблема возникает только с точкой t=L/2, которая отражается сама на себя: L-L/2=L/2, этот случай следует учесть отдельно; табличное значение для t=L/2 равно 0: \(\sin\pi=0\), а значит элемент t=L/2 может быть отображён на элемент t=0.

Рассмотренное преобразование позволяет уменьшить размеры таблицы вдвое. Теперь уменьшим таблицу ещё в 2 раза. Для этого, если $$ L/4 \le t \lt L/2 \text{ или } \pi/2 \le \phi \lt \pi, $$ воспользуемся формулой приведения $$ \sin(\pi-\phi)=\sin\phi, \text{ или} \\ \sin\phi=\sin(\pi-\phi), \text{ где } 0 \lt (\pi-\phi) \le \pi/2. $$ Можно записать то же самое в индексной форме: $$ \sin \frac{2\pi t}L=\sin\left(\pi-\frac{2\pi t}L\right)= \sin\left(\frac{2\pi}L\left(\frac L 2 -t\right)\right), \\ 0 \lt (L/2-t) \le L/4. $$ То есть, чтобы получить значение элемента с индексом t из второй четверти таблицы, следует взять значение с индексом (L/2-t), который находится либо в первой четверти, либо является особым случаем, когда t=L/4 (аргумент равен π/2). В этом случае L/2-t=L/2-L/4=L/4. Эта ситуация требует отдельного учёта, либо следует использовать таблицу с индексами от 0 до L/4, т.е. имеющую всего L/4+1 элементов.

Для того, чтобы рассмотренные преобразования были возможны, размер полной таблицы L должен быть кратен 4. А если, кроме того, L является степенью двойки (L=2^p, p - некоторое целое число), то вычисления ещё более упрощаются. В этом случае индекс t в полной таблице может быть представлен в виде p-разрядного двоичного числа из диапазона 0..(2^p-1), в котором старшие 2 бита определяют расположение по четвертям аргумента, соответствующего этому индексу. Всего возможны 4 комбинации из двух битов, рассмотрим их все.

00
Индекс t имеет вид в двоичной форме 00'xx...x, может изменяться в пределах от 00'0...0 до 00'1...1 или от 0 до 2^p-2-1=L/4-1, всего L/4 возможных значения индекса. Этим индексам соответствуют аргументы \(0 \le \phi_t \lt \pi/2\). Для определения значения синуса по индексу с использованием усечённого варианта таблицы никаких преобразований не требуется - индекс не выходит за её границы.

01
Индекс t имеет вид в двоичной форме 01'xx...x, может изменяться в пределах от 01'0...0 до 01'1...1 или от 2^p-2=L/4 до 2^p-1-1=L/2-1, всего L/4 возможных значения индекса. Этим индексам соответствуют аргументы \(\pi/2 \le \phi_t \lt \pi\). Пользуемся преобразованием \(\sin\phi=\sin(\pi-\phi)\), т.е. значение элемента с индексом t равно значению элемента с индексом t'=L/2-t, причём \(0 \lt t' \le L/4\) - индекс t' находится в пределах усечённого варианта таблицы.

10
Индекс t имеет вид в двоичной форме 10'xx...x, может изменяться в пределах от 10'0...0 до 10'1...1 или от 2^p-1=L/2 до 2^p-1+2^p-2-1=L/2+L/4-1, всего L/4 возможных значения индекса. Этим индексам соответствуют аргументы \(\pi \le \phi_t \lt 3\pi/2\). Пользуемся преобразованием \(\sin \phi=-\sin(2\pi-\phi)\), т.е. значение элемента с индексом t равно взятому с противоположным знаком значению элемента с индексом t'=L-t, причём \(L/4 \lt t' \le L/2\). Видим, что индекс t' выходит за пределы усечённого варианта таблицы, но он отвечает требованиям предыдущего пункта (за исключением значения t=L/2). Применяя преобразование из предыдущего пункта, получаем, что значение элемента t' равно элементу с индексом t''=L/2-t', причём t'' находится в пределах усечённого варианта таблицы. Выполняем суперпозицию двух преобразований и получаем, что индексу t соответствует взятый с противоположным знаком элемент усечённой таблицы с индексом t''=L/2-(L-t)=t-L/2. Причём, оказывается, что это справедливо и в случае t=L/2, которому соответствует t''=0 и значение синуса 0.

11
Индекс t имеет вид в двоичной форме 11'xx...x, может изменяться в пределах от 11'0...0 до 11'1...1 или от 2^p-1+2^p-2=L/2+L/4 до 2^p-1=L-1, всего L/4 возможных значения индекса. Этим индексам соответствуют аргументы \(3\pi/2 \le \phi_t \lt 2\pi\). Пользуемся преобразованием \(\sin \phi=-\sin(2\pi-\phi)\), т.е. значение элемента с индексом t равно взятому с противоположным знаком значению элемента с индексом t'=L-t, причём \(0 \lt t' \le L/4\), индекс t' находится в пределах усечённого варианта таблицы.

Обобщим полученные результаты в виде таблицы.

Как видно из таблицы, старший (используемый) бит индекса [p-1] определяет знак результата (0 - '+', 1 - '-'). Предшествующий старшему бит [p-2] определяет преобразование индекса для получения нового индекса, который указывает на элемент в первой четверти таблицы с таким же по модулю значением синуса. Если бит равен 0, то преобразование либо не требуется (первая строчка таблицы), либо преобразование имеет вид t'=t-L/2, что эквивалентно сбросу старшего бита индекса (третья строчка таблицы). Бит можно сбросить, применяя к индексу побитовую операцию 'И' с маской, состоящей из (p-1) единичных битов, т.е. содержащей биты 1 во всех используемых разрядах индекса, кроме старшего, определяющего знак результата. Для единообразия, маску можно применять и для случая в первой строчке таблицы, это не изменит значения индекса, так старший бит и так равен 0. На C++ это можно оформить в таком виде:

// В данном случае предполагается размер полной таблицы равным 256,
// значит индекс будет содержать 8 бит.
const unsigned int p=8;

// Определяем маску из (p-1)=7 единичных младших битов.
const int mask=((1<<Nt)-1)>>1;

// В случаях, соответствующих строкам 1, 3 таблицы, выполняем
// преобразование индекса:
t&=mask;

Если бит, определяющий преобразование индекса равен 1, то выполняется преобразование t'=L/2-t или t'=L-t (вторая и четвёртая строки таблицы). Можно показать, что и в том, и в другом случае t'=(-t)&mask. В самом деле, так как \(0 \le t' \le L/4<L/2\), значит старший используемый в индексах бит в t' сброшен, следовательно
t'=t'&mask,
если t'=L/2-t, то
t'=t'&mask=(L/2-t)&mask=(((L/2)&mask)-(t&mask))&mask=(0-(t&mask))&mask=
=(-(t&mask))&mask=(-t)&mask.
Аналогично L-t=(-t)&mask.

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

const unsigned int p=8;
const int mask=((1<<Nt)-1)>>1;

int b0=(t>>(p-2))&1; // Выделяем бит преобразования индекса.
int b1=(t>>(p-1))&1; // Выделяем знаковый бит результата.

// Преобразование индекса, используется код без ветвлений:
// если b0==0, t'=(t^0)+0=t,
// если b0==1, t'=(t^(-1))+1=~t+1=-t.
t=(t^(-b0))+b0;
t&=mask;

// Извлекаем значение из таблицы. Здесь 0<=t<=L/4, поэтому
// может использоваться короткий вариант таблицы.
Type_of_result r=table[t];
// Теперь, в зависимости от используемого типа значений для
// хранения элементов таблицы Type_of_result, используем тот
// или иной способ присвоения результату r знака, определяемого
// битом b1. Например:
r*=1-2*b1;

Использование арифметики с фиксированной точкой

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

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

Умножение также выполняется как целочисленное, но точка в результате перемещается - количество битов в дробной части равно сумме количеств битов в дробных частях множителей. Поэтому результат целочисленного умножения требует коррекции - сдвига вправо до требуемого положения точки. Кроме того, следует избегать переполнения при умножении целых. Если используется "длинное" умножение, когда результат умножения имеет двойную длину, переполнения никогда не возникает, однако, даже если в процессоре аппаратно реализовано "длинное" умножение, оно не поддерживается языками высокого уровня, требует ассемблерных вставок в код и делает код аппаратно-зависимым, что нехорошо. При использовании обычного умножения, можно предложить проводить коррекцию в два этапа: сначала выполнить предкоррекцию одного или обоих множителей - тем самым теряем точность, но избегаем переполнения; затем выполнить посткоррекцию для сдвига точки в требуемую позицию.

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

Рассмотрим реализацию линейной интерполяции с использованием арифметики с фиксированной точкой. Размер полной таблицы синусов принимаем равным L=2^p, тогда индекс таблицы t=0..L-1 будет содержать p используемых бит. Например, если используется 12-разрядный DAC, то получим точность вычислений не хуже 1/2 LSB при выборе значений L=256 (p=8, для индексации таблицы достаточно 8 битов).

Индекс фазы k, имеющий размер 32 бита, может принимать одно из 2³² значений. Следовательно, на 2³²/L значений индекса фазы приходится одно табличное значение синуса. И полный диапазон значений индекса фазы, и полный диапазон индексов таблицы соответствуют одному и тому же диапазону значений фазы [0, 2π); $$ \phi=\frac{2\pi}L t, \\ \phi=\frac{2\pi}{2^{32}} k, \\ t=\frac L{2\pi}\phi=\frac L{2\pi} \frac{2\pi}{2^{32}} k=\frac L {2^{32}} k, \\ L=2^p, \\ t=\frac{2^p}{2^{32}}k=\frac 1{2^{32-p}}k. $$ Видим, для того, чтобы получить индекс в таблице, необходимо разделить k на 2^32-p или сдвинуть k на (32-p) битов вправо. В нашем примере с 12-разрядным DAC, (32-p)=24: сдвиг производится на 24 бита вправо.

Точное целочисленное деление возможно только в том случае, если все (32-p) младших битов в индексе фазы k являются нулевыми. Только в этом случае фаза, соответствующая индексу фазы k точно равна фазе, соответствующей строчке таблицы t. В любом другом случае, фаза, соответствующая индексу фазы k находится "между строчками" таблицы с индексами t=k>>(32-p) и t+1 (если t+1==L, то в качестве t+1 используем индекс 0). Шаг изменения аргумента при переходе от одного элемента таблицы к следующему в единицах индекса фазы равен d=1<<(32-p). Разность между текущим индексом фазы и индексом фазы, соответствующим элементу таблицы с индексом t равен k&mask, где mask=(1<<(32-p))-1 - маска, состоящая из (32-p) единичных битов в младших разрядах. С учётом этого, формула для линейной интерполяции примет вид:
s(k)=T[t]+(T[t+1]-T[t])*(k&mask)/d,
где под T подразумевается таблица синусов. Выражение можно преобразовать следующим образом: $$ s(k)=T[t]\left(1-\frac{k\&mask}d\right)+T[t+1]\frac{k\&mask}d= \\ =T[t]\left(1-\frac{k\&mask}{2^{32-p}}\right)+T[t+1]\frac{k\&mask}{2^{32-p}}, $$ $$ s(k)=\frac{T[t](2^{32-p}-(k\&mask))+T[t+1](k\&mask)}{2^{32-p}}. \tag 1 $$ Здесь числа (k&mask) и (2^32-p-(k&mask)) являются целыми. Умножение числа с фиксированной точкой на целое выполняется как перемножение целых чисел, результат будет являться числом с фиксированной точкой с таким же положением точки, как у исходного множителя с фиксированной точкой. Поэтому, в полученном выражении, если значения таблицы T представлены числами с фиксированной точкой, то все операции при вычислении числителя выполняются в целых числах. Для получения окончательного результата, следует найденный числитель разделить на целое число 2^32-p; деление числа с фиксированной точкой выполняется как деление целых, полученный результат будет числом с фиксированной точкой, положение которой такое же, как у делимого. В данном случае деление может быть выполнено как сдвиг вправо на (32-p) битов.

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

Например, пусть таблица синусов хранится в формате с фиксированной точкой с 16 битами в дробной части. Максимальное по модулю значение синуса равно 1, в используемом нами формате с фиксированной точкой оно будет представлено целым числом 2¹⁶. В соответствии с формулой (1), нам потребуется умножать значение синуса из таблицы на числа (2^32-p-(k&mask)) и (k&mask). Если полная таблица синусов имеет размер 256 элементов, то p=8 и указанные числа не превышают 2²⁴.

Для того, чтобы не происходило переполнения 32-битовых целых со знаком, нам потребуется корректирующий множитель K, такой что
$$ K \cdot 2^{16} \cdot 2^{24} \le \text{0x7FFFFFFF} $$ или $$ K \cdot 2^{40} \lt \text{0x80000000}=2^{31}, \\ K \lt 2^{-9}, \text{ например,} \\ K=2^{-10}, m1=10. $$ Для предкоррекции потребуется сдвинуть вправо множители (2^32-p-(k&mask)) и (k&mask) на 10 битов вправо, при этом их эффективная разрядность уменьшиться с (32-p)=24 битов до 14 битов. Это допустимо, так как новая разрядность сопоставима с разрядностью значений в таблице синусов (16 разрядов после точки) и превышает разрядность DAC. Можно распределить коррекцию между двумя множителями: табличное значение синуса сдвигать на 1 бит, а фазовые множители - на 9 битов, тогда эффективные разрядности обоих сомножителей будут равны по 15 бит.

Предкоррекция множителей должна быть учтена при выполнении деления в соответствии (1) на коэффициент 2^32-p. То есть, вычисленное по формуле значение числителя (с выполнением предкоррекции путём сдвига вправо на m1), надо будет сдвигать вправо уже не на (32-p) бита вправо, а на (32-p-m1), в нашем примере - на 14.

Пример практической реализации

Скачать исходный код примера

В этом примере реализован синтез синусоидального сигнала с использованием табличного метода вычисления синуса и линейной интерполяции. Код отлаживался на микроконтроллере STM32F100RB, тактовая частота 24 МГц. Аналоговый сигнал формируется на выходе первого канала DAC, вывод PA4 микроконтроллера. Для увеличения быстродействия, вычисляемые значения сигнала помещаются в буфер и передаются в DAC с использованием механизма DMA.

В режиме отладки, без оптимизации, быстродействие микроконтроллера позволило достичь значений опорной частоты F0 до 80 кГц (при этом можно синтезировать сигнал с частотой до 20 кГц).
С оптимизацией по скорости (даже при выполнении из RAM) легко достижимы значения F0 до 400 кГц (частота синтезируемого сигнала до 100..150 кГц). При выполнении программы из ROM достижимо и более высокое быстродействие.

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

Смотрите пояснения к коду после листинга.

// Генерация синусоиды методом DDS; демонстрируется
// табличное вычисление синуса с линейной интерполяцией.
// Отлажено на STM32F100RB (24МГц).
// Используется: DAC канал 1 (PA4), DMA1 канал 3, TIM6.

// F0 до 80кГц при компиляции без оптимизации;
// до 400кГц при компиляции с оптимизацией по скорости.

#include <stm32f10x.h>

#include <math.h>

// Максимальный выходной сигнал, постоянное смещение и амплитуда
// синтезируемой синусоиды в единицах Vref/65536.
const int out_max=0xFFFF;
int out_a0=out_max/2;
int out_a1=out_max/3;

// Тактовая частота сигнала на входе TIM6.
#define FSYSCLK 24000000
// Задаём желаемую опорную частоту (максимальное значение
// зависит от быстродействия микроконтроллера).
// Возможны только значения вида FSYSCLK/K, где
// K-целое из диапазона 2..65536.
// Можно задать требуемый делитель непосредственно в коде,
// отвечающем за конфигурирование таймера.
#define F0 80000

// Текущее значение индекса фазы.
unsigned int ph_ind=0;

// Приращение индекса фазы в каждом цикле синтеза;
// частота синтезируемого сигнала f=F0*ph_m/(2.0**32),
// где F0 - опорная частота.
unsigned int ph_m=0;

// Количество бит в индексе таблицы синусов.
#define Nt 8
// Размер полной таблицы.
const unsigned int L=1<<Nt;
// Размер используемой усечённой таблицы.
const unsigned int tsize=L/4+1;

// Вспомогательные константы.
const unsigned int Nt_1=Nt-1;
const unsigned int Nt_2=Nt-2;
const unsigned int m=32-Nt;

// Маска, выделяющая используемую часть индекса для
// доступа к усечённой таблице: (Nt-1) единиц в младших битах.
const unsigned int ind_mask=((1<<Nt)-1)>>1;

// Таблица значений sin для индексов 0..L/4 (аргументы 0..M_PI/2);
// формат с фиксированной точкой, 15 бит в дробной части.
unsigned short int table[tsize];    //fix15 format

// Функция инициализации таблицы синусов.
void sin_table_init()
{
    for(unsigned int i=0; i<tsize; i++)
        table[i]=0x8000*sinf(2*(float)M_PI*i/L);
}

// Функция, используя усечённую таблицу, возвращает значение sin для
// заданного индекса в полной таблице: t=0..L-1
// (допустимы значения t>=L, данная функция периодична с периодом L).
inline int select(unsigned int t)
{
    unsigned int res_sign=(t>>Nt_1)&1;
    unsigned int arg_sign=(t>>Nt_2)&1;
    t=(t^(-arg_sign))+arg_sign;
    t&=ind_mask;
    int r=table[t];
    return (r^(-res_sign))+res_sign;
}

// Функция, используя табличные значения sin, выполняет линейную
// интерполяцию и возвращает sin для заданного индекса фазы
// (индекс фазы - любое целое 32-битное число без знака).
inline int interp(unsigned int k)
{
    const unsigned int mask2=(1<<m)-1;
    const unsigned int m1=m-15;
    const int c=1<<15;

    unsigned int t1=k>>m;
    int y1=select(t1);
    int y2=select(t1+1);  // Случай t1==L-1 => t1+1==L допустим.

    int dk=(k&mask2)>>m1;
    int r=(y1*(c-dk)+y2*dk)>>14;
    return r;
}

// Буфер значений, загружаемых в DAC с использованием DMA.
const unsigned int dma_buf_length=64;
short int dma_buf[dma_buf_length];

// Функции для обновления буфера.
// Используют глобальные переменные ph_ind,  ph_m.

// Обновить заданный фрагмент буфера.
inline void dma_buf_fragment_init(short int *begin, const short int *end)
{
    int a1=out_a1>>1; // Предкоррекция во избежание переполнения.
    while(begin!=end)
    {
        *begin=out_a0+
            (a1*interp(ph_ind)>>15); // Умножение и посткоррекция.
        ph_ind+=ph_m;
        begin++;
    }
}

// Вычислить значения для первой половины буфера.
inline void dma_buf_first_half_init()
{
    dma_buf_fragment_init(dma_buf, dma_buf+dma_buf_length/2);
}

// Вычислить значения для второй половины буфера.
inline void dma_buf_second_half_init()
{
    dma_buf_fragment_init(dma_buf+dma_buf_length/2, dma_buf+dma_buf_length);
}

// Обработчик прерываний от DMA (канал 3).
extern "C" void DMA1_Channel3_IRQHandler()
{
    if(DMA1->ISR&DMA_ISR_HTIF3)
        dma_buf_first_half_init();
    else if(DMA1->ISR&DMA_ISR_TCIF3)
        dma_buf_second_half_init();
    DMA1->IFCR|=                // Сброс соотв. флагов записью 1.
            DMA_IFCR_CTEIF3|
            DMA_IFCR_CHTIF3|
            DMA_IFCR_CTCIF3|
            DMA_IFCR_CGIF3;
}

int main()
{
    // Настраиваем приоритет о разрешаем обработку прерываний DMA.
    NVIC_SetPriority(DMA1_Channel3_IRQn, 0);
    NVIC_EnableIRQ(DMA1_Channel3_IRQn);

    // Включаем тактовые сигналы используемой периферии:
    // DMA, DAC, TIM6, GPIOA.
    RCC->AHBENR|=RCC_AHBENR_DMA1EN;
    RCC->APB1ENR|=RCC_APB1ENR_DACEN|RCC_APB1ENR_TIM6EN;
    RCC->APB2ENR|=RCC_APB2ENR_IOPAEN;

    // PA4 - аналоговый вход в соответствии с требованиями DAC.
    GPIOA->CRL&=~(0xf<<4*4);

    // Инициализация таблицы синусов.
    sin_table_init();

    // Настраиваем и включаем используемый для пересылки данных
    // в DAC канал DMA.
    DMA1_Channel3->CPAR=(uint32_t)&DAC->DHR12L1; // ПУ
    DMA1_Channel3->CMAR=(uint32_t)dma_buf;       // Адрес буфера в памяти.
    DMA1_Channel3->CNDTR=dma_buf_length;         // Кол-во элементов.

    uint32_t tmp=DMA1_Channel3->CCR&~0x7FFF|
            DMA_CCR3_TCIE|      // Прерывание после завершения передачи.
            DMA_CCR3_HTIE|      // Прерывание после полузавершения.
            DMA_CCR3_DIR|       // Направление: память->периферия.
            DMA_CCR3_CIRC|      // Включить циклический режим.
            // DMA_CCR3_PINC=0  // Не инкрементировать адрес периферии.
            DMA_CCR3_MINC|      // Инкрементировать адрес в памяти.
            DMA_CCR3_PSIZE_1|   // Размер слова периферии=10: 32 бита.
            DMA_CCR3_MSIZE_0|   // Размер слова в памяти=01: 16 бит.
            DMA_CCR3_PL;        // 11 - наибольший приоритет.
    DMA1_Channel3->CCR=tmp;
    DMA1_Channel3->CCR|=DMA_CCR3_EN;

    // Настраиваем и включаем DAC.
    DAC->CR|=
        DAC_CR_DMAEN1| // Включить генерацию запросов DMA.
        DAC_CR_TEN1;   // Включить запуск DAC преобразования по внешнему сигналу.
        // TSEL[2:0]=000: запуск сигналом Timer 6 TRGO event.
    DAC->CR|=DAC_CR_EN1;

    // Настраиваем и включаем источник опорного сигнала -
    // таймер TIM6. После включения таймера начнётся генерация
    // сигнала с частотой f=F0*ph_m/(2.0**32) =>
    // ph_m=(uint32_t)((2.0**32)*f/F0);

    // Примеры:
    ph_m=2684354; // 50 Гц
    //ph_m=2689723; // 50.1 Гц
    //ph_m=268435456; // 5000 Гц
    //ph_m=268435456*2; // 10000 Гц
    //ph_m=662803111; // 12345.67Hz

    // Начальная инициализация буфера данных для загрузки в DAC,
    // в дальнейшем буфер будет обновляться по прерываниям от DMA.
    dma_buf_first_half_init();
    dma_buf_second_half_init();

    TIM6->ARR=FSYSCLK/F0-1;   // Частота переполнений: FSYSCLK/(ARR+1)=F0.
    TIM6->CR2&=~TIM_CR2_MMS;
    TIM6->CR2|=TIM_CR2_MMS_1; //MMS[2:0]=010: как TRGO используем UEV.
    TIM6->CR1|=TIM_CR1_CEN;   // Включаем таймер.

    while(true);
}

Пояснения к коду.

В используемом микроконтроллере имеется 12-разрядный цифро-аналоговый преобразователь, однако, он может работать с данными, выровненными влево, т.е. позволяет загружать 16-разрядные числа. В нашем случае этот вариант удобен, он позволяет непосредственно загружать в DAC числа с фиксированной точкой с 16 битами в дробной части. Значению 0 соответствует выходной сигнал 0, а значению 0xFFFF - максимальное выходное напряжение V_REF(-1 LSB). При этом, оставаясь по-прежнему 12-разрядным, DAC на самом деле отбрасывает последние 4 бита в 16-битовым значении, так что от нас требуется обеспечить лишь такую точность, чтобы старшие 12 бит в 16-битном значении были верными.

В данном коде имеется возможность в процессе выполнения изменять амплитуду синтезируемого сигнала и величину постоянной составляющей, изменяя значения переменных out_a1, out_a0. Значения сигнала вычисляются по формуле u(ph)=a0+a1*sin(ph), ph - значение фазы; значение синуса находится табличным способом.

Для хранения текущего индекса фазы и шага приращения индекса используются глобальные переменные ph_ind, ph_m. Это 32-битные целые без знака. Изменяя ph_m, можем задать желаемую частоту синтезируемого сигнала: f=F0*ph_m/(2.0**32), f<F0/(3..4). Да, использовать подобным образом глобальные переменные - плохой стиль, однако в демонстрационных целях - допустимо.

Для вычисления значений синуса используется табличный метод с линейной интерполяцией. С учётом разрешающей способности DAC (12-разрядный DAC), выбираем размер полной таблицы синусов в L=256 элементов, размер усечённого варианта составит L/4+1=65 элементов. Для индексации полной таблицы требуется Nt=8 битов (2⁸=256). Используем 16-битовый целый тип для хранения элементов таблицы, элементы таблицы интерпретируются как значения с фиксированной точкой, 15 битов в дробной части. Выбор формата, с одной стороны, позволяет получить максимально возможную точность, с другой стороны, имеет 1 бит в целой части, что даёт возможность хранить значение 1 (sin π/2=1).

Определение значения синуса по индексу в полной таблице с использованием только значений из усечённой подробно описано ранее. Здесь это реализовано в функции select. Линейная интерполяция значений реализована в функции interp, аргументом которой является целое без знака 32-битовое число - индекс фазы, значениям которого от 0 до 0xFFFFFFFF соответствует фаза [0, 2π).

Для достижения максимального быстродействия используется механизм DMA, который обеспечивает загрузку значений в DAC из буфера без участия процессора. Буфер обновляется функциями dma_buf_first_half_init(), dma_buf_second_half_init(), которые соответственно обновляют первую и вторую половину буфера. Функции используют значения глобальных переменных ph_ind, ph_m и обновляют значение ph_ind. В первый раз буфер инициализируется в начале выполнения программы, а затем - по прерываниям от устройства DMA. Прерывание возникает после пересылки половины буфера, что служит сигналом для обновления первой половины буфера и после пересылки в DMA всех значений из буфера, что служит сигналом к обновлению второй половины буфера (к этому моменту первая половина буфера будет уже обновлена в результате обработки предыдущего прерывания). Устройство DMA настроено для работы в циклическом режиме, так что после завершения передачи всех данных из буфера, передача начинается сначала.

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