Вычисление логарифмов довольно распространённая операция в цифровой обработке сигналов. Чаще пожалуй приходится считать только свёртки (умножение с накоплением) и амплитуды с фазами. Как правило для вычисления логарифмов на FPGA применяется алгоритм CORDIC в гиперболическом варианте, требующий только таблицы и простых арифметических операций. Однако это не всегда бывает удобно, особенно если проект большой, кристалл маленький и начинаются танцы с оптимизацией. Именно с такой ситуацией и пришлось мне однажды столкнуться. Оба порта блока RAM (Cyclone IV) уже плотненько были в работе, не оставляя свободных окон. Использовать ещё один блок под гиперболический CORDIC не хотелось. Зато был умножитель, для которого во временной диаграмме получалось приличное свободное окно. Денёк подумав, я сочинил следующий алгоритм, в котором не используется таблиц, но есть умножение, точнее возведение в квадрат. И поскольку схемотехнически возведение в квадрат проще общего случая умножения, возможно этот алгоритм представляет интерес для специализированных микросхем, хотя для FPGA разницы конечно нет. Подробнее под катом.
Объяснять что к чему, тут проще для действительных чисел. С них и начнём. К целочисленной реализации перейдём потом.
Пусть есть число X. Требуется найти число Y такое, чтобы 
.
Положим так же что X лежит в интервале от 1 до 2. Это не слишком ограничивает общность, поскольку X всегда можно перевести в этот интервал умножением или делением на степень двойки. Для Y это будет означать добавление или вычитание целого числа, что делается легко. Итак X лежит в интервале от 1 до 2. Тогда Y будет лежать в интервале от 0 до 1. Запишем Y как бесконечную двоичную дробь:
Коэффициенты 
в этой записи есть ни что иное, как просто биты двоичного представления числа Y. Причём поскольку Y меньше 1, очевидно что 
=0.
Возведём наше первое уравнение в квадрат: 
и как и ранее, запишем двоичное представление числа 2Y. Очевидно, что
Т.е. биты 
остались теми же, просто сдвинулись степени двойки. Бита 
в представлении нет, потому что он равен нулю.
Возможны два случая:
1) 
, Y > 1, 
2) 
, Y < 1, 
В первом случае в качестве нового значения X примем 
, во втором — 
.
В итоге задача свелась к прежней. Новое X опять лежит в интервале от 1 до 2, новое Y от 0 до 1. Но мы узнали один бит результата. Делая такие же шаги в дальнейшем, мы можем получить сколько угодно битов Y.
Посмотрим как это работает в программе на С:
#include <stdio.h>
#include <math.h>
int main()
{
double w=1.4;
double s=0.0;
double a=0.5;
double u=w;
for(int i=0; i<16; i++)
{
u=u*u;
if(u>2)
{
u=u/2;
s+=a;
}
a*=0.5;
}
w=log2(w);
double err=100*abs(2*(s-w)/(s+w));
printf("res=%f, log=%f, err=%f%cn",s,w,err,'%');
return 0;
}
Мы посчитали логарифм с точностью до 16 бит и сравнили с тем, что даёт математическая библиотека. Программа вывела:
res=0.485413, log=0.485427, err=0.002931%
Результат совпал с библиотекой с точностью 0.003%, что показывает работоспособность нашего алгоритма.
Перейдём к целочисленной реализации.
Пусть N-разрядные двоичные беззнаковые числа, отображают интервал [0, 1]. Для удобства будем считать единицей число 
, а не 
, и соответственно двойкой число 
. Напишем программу по образу и подобию предыдущей, но работающую с целыми числами:
#include <stdio.h>
#include <math.h>
#define DIG 18 //разрядность чисел
#define N_BITS 16 //число бит которое считаем
unsigned ONE=1<<(DIG-1); //единица
unsigned TWO=ONE<<1; //двойка
unsigned SCALE=1<<(N_BITS+1); //масштаб логарифмов
unsigned myLog(unsigned w)
{
unsigned s=0;
unsigned long long u=w;
for(int i=0; i<N_BITS+1; i++)
{
s<<=1;
u=(u*u)>>(DIG-1); //берём только старшие разряды !
if(u&TWO) //тут достаточно одного бита для сравнения
{
u>>=1;
s+=1;
}
printf("%Xn", (int)u);
}
return s;
}
int main()
{
double w=1.2345678;
unsigned iw=(unsigned)(ONE*w);
double dlog=log2(w);
unsigned ilog=myLog(iw);
unsigned test=(unsigned)(SCALE*dlog);
int err=abs((int)(ilog-test));
printf("val=0x%X, res=0x%X, log=0x%X, err=%dn",iw,ilog,test,err);
return 0;
}
Поиграв в программе с разными разрядностью чисел (DIG), точностью вычислений (N_BITS), и аргументами логарифма (w), видим, что всё вычисляется правильно. В частности с теми параметрами, которые заданны в этом исходнике, программа выдаёт:
val=0x27819, res=0x9BA5, log=0x9BA6, err=1
Теперь всё готово для того, чтобы реализовать на верилоге железку, делающую точно то же самое, что и функция myLog на С. Переменные s и u в нашей функции можно распечатать в цикле и сравнивать с тем, что выдает симулятор верилога. Соответствие этих переменных железной реализации весьма прозрачно и понятно. u это рабочий регистр, принимающий во время итераций новые значения X. s это сдвиговый регистр, в котором накапливается результат. Интерфейс нашего модуля будет выглядеть вот так:
module logarithm(
input clk, //клок
input wr, //сигнал записи данных
input[17:0] din, //входная шина данных
output[nbits-1:0] dout, //выходная шина данных
output rdy //сигнал готовности
);
parameter nbits=16; //точность вычисления
Входная шина принята 18-разрядной, соответственно разрядности умножителей в Cyclone IV. Числа на наш модуль должны поступать нормализованные. Т.е. со старшим битом, равным единице. В моём проекте это выполнялось автоматически. Но в случае чего реализовать нормализатор, думаю никому труда не составит. Точность вычислений задаётся параметром nbits, по умолчанию равным 16. Модуль считает по одному биту за такт, и за 16 тактов вычислит логарифм с точностью до 16 битов. Если надо быстрее с той же точностью или точнее с той же скоростью, надеюсь никому не составит особого труда разделить модуль на несколько устройств и конвейеризовать.
//--------------------- logarithm.v ------------------------------//
module logarithm(
input clk, //клок
input wr, //сигнал записи данных
input[17:0] din, //входная шина данных
output[nbits-1:0] dout, //выходная шина данных
output rdy //сигнал готовности
);
parameter nbits=16; //точность вычисления
reg[4:0] cnt; //счетчик
reg[17:0] acc; //рабочий регистр-аккумулятор
reg[nbits-1:0] res; //результат
always @(posedge clk)
if(wr) cnt<=nbits+1;
else
if(cnt != 0)
cnt<=cnt-1;
wire[35:0] square=acc*acc; //квадрат аккумулятора
wire bit=square[35]; //бит результата
wire[17:0] next = bit ? square[35:18] : square[34:17]; //новый аккумулятор
always @(posedge clk)
if(wr) acc<=din;
else
if(cnt != 0) begin
acc<=next;
#10 $display("%X", acc);
end
always @(posedge clk)
if(wr) res<=0;
else
if(cnt != 0) begin
res[nbits-1:1]<=res[nbits-2:0];
res[0]<=bit;
end
assign dout=res;
assign rdy=(cnt==0);
endmodule
//======================== testbench.v =====================//
module testbench();
reg clk; //клок
always
#100 clk=~clk;
reg wr; //запись
reg[17:0] din; //вход
wire rdy; //готовность
wire[15:0] dout; //выход
logarithm log2(
.clk (clk),
.wr (wr),
.din (din),
.dout (dout),
.rdy (rdy)
);
// ждёт n клоков и немного ещё
task skipClk(integer n);
integer i;
begin
for(i=0; i<n; i=i+1)
@(posedge clk);
#10 ;
end
endtask
initial
begin //тест
$dumpfile("testbench.vcd");
$dumpvars(0, testbench);
clk=0;
wr=0;
din=18'h27819;
skipClk(3);
wr=1;
skipClk(1);
wr=0;
@(rdy);
skipClk(3);
$display("value=%X, result=%X", din, dout);
$display("Done !");
$finish;
end
endmodule
Запустим тест вот таким скриптом:
#!/bin/sh
rm -f *.vvp
rm -f *.vcd
iverilog -o testbench.vvp logarithm.v testbench.v
vvp testbench.vvp
gtkwave testbench.vcd testbench.gtkw
Запустив тест, видим окончательный вывод симулятора — value=27819, result=9ba5. Верилог выдал то же самое что С. Временная диаграмма тут довольно тривиальна и особого интереса не представляет. Поэтому её не привожу.
Verilog C
30c5d 30C5D
252b1 252B1
2b2bc 2B2BC
3a3dc 3A3DC
35002 35002
2be43 2BE43
3c339 3C339
38a0d 38A0D
321b0 321B0
273a3 273A3
30163 30163
24214 24214
28caf 28CAF
34005 34005
2a408 2A408
37c9d 37C9D
30a15 30A15
Убеждаемся, что они совпадают бит в бит. Итого, реализация на верилоге повторяет с точностью до бита модель на С. Это и есть тот результат, которого следует достигать, реализуя алгоритмы в железе.
На этом пожалуй всё. Надеюсь кому-то этот мой опыт окажется полезен.
Автор: Калягин Евгений Игоревич
