eugenk 29 сен 2019 в 08:50

Один способ вычисления логарифма по основанию 2

5 мин

27K

Recovery Mode

Вычисление логарифмов довольно распространённая операция в цифровой обработке сигналов. Чаще пожалуй приходится считать только свёртки (умножение с накоплением) и амплитуды с фазами. Как правило для вычисления логарифмов на FPGA применяется алгоритм CORDIC в гиперболическом варианте, требующий только таблицы и простых арифметических операций. Однако это не всегда бывает удобно, особенно если проект большой, кристалл маленький и начинаются танцы с оптимизацией. Именно с такой ситуацией и пришлось мне однажды столкнуться. Оба порта блока RAM (Cyclone IV) уже плотненько были в работе, не оставляя свободных окон. Использовать ещё один блок под гиперболический CORDIC не хотелось. Зато был умножитель, для которого во временной диаграмме получалось приличное свободное окно. Денёк подумав, я сочинил следующий алгоритм, в котором не используется таблиц, но есть умножение, точнее возведение в квадрат. И поскольку схемотехнически возведение в квадрат проще общего случая умножения, возможно этот алгоритм представляет интерес для специализированных микросхем, хотя для FPGA разницы конечно нет. Подробнее под катом.

Объяснять что к чему, тут проще для действительных чисел. С них и начнём. К целочисленной реализации перейдём потом.

Пусть есть число X. Требуется найти число Y такое, чтобы $inline$ .
Положим так же что X лежит в интервале от 1 до 2. Это не слишком ограничивает общность, поскольку X всегда можно перевести в этот интервал умножением или делением на степень двойки. Для Y это будет означать добавление или вычитание целого числа, что делается легко. Итак X лежит в интервале от 1 до 2. Тогда Y будет лежать в интервале от 0 до 1. Запишем Y как бесконечную двоичную дробь:

$Y={b_0}2^0 + {b_1}2^{-1}+...+{b_n}2^{-n}+...$

Коэффициенты ${b_i}$ в этой записи есть ни что иное, как просто биты двоичного представления числа Y. Причём поскольку Y меньше 1, очевидно что ${b_0}$ =0.

Возведём наше первое уравнение в квадрат: $X^2=2^{2Y}$ и как и ранее, запишем двоичное представление числа 2Y. Очевидно, что

$2Y={b_1}2^0+{b_2}2^{-1}+...+{b_n}2^{-(n-1)}+...$

Т.е. биты ${b_i}$ остались теми же, просто сдвинулись степени двойки. Бита ${b_0}$ в представлении нет, потому что он равен нулю.

Возможны два случая:

1) $inline$ , 2Y > 1, ${b_1} = 1$

2) $inline$ , 2Y < 1, ${b_1} = 0$

В первом случае в качестве нового значения X примем $inline$ , во втором — $inline$ .

В итоге задача свелась к прежней. Новое 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%c\n",s,w,err,'%');
	return 0;
}

Мы посчитали логарифм с точностью до 16 бит и сравнили с тем, что даёт математическая библиотека. Программа вывела:

res=0.485413, log=0.485427, err=0.002931%

Результат совпал с библиотекой с точностью 0.003%, что показывает работоспособность нашего алгоритма.

Перейдём к целочисленной реализации.

Пусть N-разрядные двоичные беззнаковые числа, отображают интервал [0, 1]. Для удобства будем считать единицей число $inline$ , а не $inline$ , и соответственно двойкой число $2^{N+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("%X\n", (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=%d\n",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. Верилог выдал то же самое что С. Временная диаграмма тут довольно тривиальна и особого интереса не представляет. Поэтому её не привожу.

Сравним промежуточный вывод симулятора (acc) и программы на С (s):

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

Убеждаемся, что они совпадают бит в бит. Итого, реализация на верилоге повторяет с точностью до бита модель на С. Это и есть тот результат, которого следует достигать, реализуя алгоритмы в железе.

На этом пожалуй всё. Надеюсь кому-то этот мой опыт окажется полезен.

Теги:

Хабы:

Один способ вычисления логарифма по основанию 2

Публикации

Истории

Ближайшие события