Вычисление логарифма

Анатолий Сергиенко

Тот, кто еще три десятилетия назад изучал теоретические основы электротехники, должен был в совершенстве владеть логарифмической линейкой. Потому, что надо было много умножать, а в ту пору цена калькулятора равнялась цене сотни экземпляров 7-значной таблицы логарифмов. Логарифмическая линейка демонстрировала фокус, как простое сложение отрезков на логарифмических шкалах дает функции произведения, деления, извлечения корня, преобразования координат — что называется, ловкость рук.

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

Сейчас функция логарифма часто встречается в специализированных вычислителях. Например, при цифровой обработке сигналов логарифмическая шкала используется для сравнения уровней сигналов в децибелах, в простых алгоритмах сжатия сигналов, некоторых алгоритмах нелинейной обработки сигналов, таких как кепстральная обработка. Поэтому в данной заметке рассмотрены некоторые способы реализации этой функции.

В спецвычислителях выгодно вычислять функцию логарифма по основанию 2.

Во-первых, если аргумент x — целое число, то логарифм приближенно равен числу значащих цифр p. То есть,

log2(x) = log2( 2^p * xm ) = p + log2(xm),

где 1 &gt xm &gt 0.5 — мантисса числа.

Во-вторых, специалисту по цифровой обработке сигналов привычно оценивать уровень сигнала в децибелах как приблизительно равным 6 * p, если x — это амплитуда сигнала или приблизительно равным 3 * p, если x — мощность сигнала.

На вышекуказанном свойстве основаны следующие алгоритмы вычисления логарифма. Алгоритмы описаны на VHDL с применением библиотек IEEE.std_logic_1164 и IEEE.STD_logic_arith.

Модель PARITY предназначена для приближенного вычисления логарифма, то есть для вычисления числа значащих цифр p.

Она определяет номер старшего значащего разряда n — значного двоичного целого числа x. По сути это схема приоритетного шифратора, поэтому она называется PARITY.

entity LOG2 is
     generic(n: integer:=15; m: integer:= 4);
	port (X: in STD_LOGIC_VECTOR(n-1 downto 0);
		Y: out STD_LOGIC_VECTOR(m-1 downto 0));
end LOG2;

architecture PARITY of LOG2 is
begin 
    process(X)
    variable p:STD_LOGIC_VECTOR(m-1 downto 0);
    begin            
        p:= CONV_STD_LOGIC_VECTOR(n,m);
    for i in 1 to n loop
        exit when x(n-i)='1';
        p:= UNSIGNED(p)-1;       
    end loop;
      Y<=p; 
   end process;
end PARITY;

Вначале операнд-вектор p принимает значение числа n разрядов данного X. Для этого используется функция CONV_STD_LOGIC_VECTOR(n,m), преобразующая целое число разрядов n в m - разрядный вектор. Затем в цикле разряды вектора X последовательно проверяются, начиная со старшего, пока не будет найдена цифра 1. Одновременно из p вычитается 1. В результате, p становится равным номеру старшей значащей 1 в числе x, то есть приближенному значению логарифма.

Несмотря на сравнительную простоту описания алгоритма (попробуйте эту задачу решить с помощью логических функций), он дает приличные результаты при синтезе логической схемы. Правда, не все компиляторы-синтезаторы распознают правильно идиому таких циклов, поэтому необходимы эксперименты с каждым конкретным синтезатором. В таблице 1 приведены некоторые результаты синтеза архитектуры PARITY синтезатором SYNPLIFY для ПЛИС серии XC4000XLA-0.7.

Таблица 1.

Разрядность , n	7	15	24	31	48
Аппаратные затраты, LUT	4	10	24	27	55
Период вычислений, нс	10	13.2	16.6	13.9	17

Для более точного вычисления логарифма его можно представить как:

log2(x) = log2( 2^p * xm ) = p-1 + log2(xl),

где xl = 2 * xm; 2 > xl ≥ 1.

Тогда старшие m разрядов результата будут составлять разряды числа p-1, а младшие разряды - разряды положительной функции log2(xl). Функцию log2(xl) можно реализовать таблично, например, следующим процессом:

LOGTABL: process(xl)
  begin                  
   case xl is
	when "0000" => log<="0000";
	when "0001" => log<="0010";
	when "0010" => log<="0011";
	when "0011" => log<="0100";
	when "0100" => log<="0101";
	when "0101" => log<="0111";
	when "0110" => log<="1000";
	when "0111" => log<="1001";
	when "1000" => log<="1010";
	when "1001" => log<="1011";
	when "1010" => log<="1011";
	when "1011" => log<="1100";
	when "1100" => log<="1101";
	when "1101" => log<="1110";
	when "1110" => log<="1111";
	when "1111" => log<="1111";
	when others => 
  end case;                
end process;

где xl - 4 старшие значащие разряды исходного числа, log приблизительно равен log2(xl).

Для выделения разрядов xl необходимо построить сдвигатель, на входы которого подаются входное число X и номер старшей значащей единицы p.

Алгоритм работы сдвигателя можно совместить с алгоритмом вычисления числа p следующим образом:

 signal xl:STD_LOGIC_VECTOR (3 downto 0);      
 . . .
 SHIFT:process(X)
    variable xs:STD_LOGIC_VECTOR (n-1 downto 0);
    variable p:STD_LOGIC_VECTOR(m-1 downto 0);
  begin            
	xs:=X;
	p:= CONV_STD_LOGIC_VECTOR(n-1,m);
  for i in 1 to n-1 loop
	exit when x(n-i)='1';         
         xs:=xs(n-2 downto 0)&'0';--сдвиг влево
	  p:=unsigned(p)-1;       
		end loop;                                  
	  xl<=xs(n-2 downto n-5); 
	end process;

Полученная в результате архитектура DSPLOG приведена в приложении 1.

В таблице 2 показаны некоторые результаты синтеза архитектуры DSPLOG синтезатором SYNPLIFY для ПЛИС серии XC4000XLA-0.7.

Таблица 2.

Разрядность , n	8	16	24	32	48
Аппаратные затраты, LUT	26	66	109	147	233
Период вычислений, нс	19.8	22.4	23.9	25	27

Синтезированная структура показана на рис. 1. На входе и выходе устройства включены регистры для обеспечения высокой пропускной способности и повторяемости при реализации в ПЛИС. Устройство имеет 3 уровня логики. Первый уровень составляют схемы И сравнения соседних разрядов, второй уровень - два приоритетных коммутатора для определения веса p старшей единицы и сдвига четырех старших значащих разрядов, третий уровень - таблица логарифмов LOGTABL. Точности результатов этой схемы (7-9 двоичных разрядов) достаточно для многих приложений цифровой обработки сигналов.

Рис1.Структура синтезированного вычислителя логарифма 8 - разрядных чисел.

Если на вход архитектуры PARITY подавать квадрат числа х, то удвоится количество входных разрядов, зато разрядность результата увеличится на 1. Это означает, что удвоится точность результата. При q-кратном последовательном возведении в квадрат входного данного, количество точных разрядов результата увеличивается на q. Для упрощения таких вычислений возводить в квадрат можно только несколько старших значащих разрядов нормализованного операнда, выполняя слежение за масштабом результата.

На этом свойстве основан следующий алгоритм вычисления логарифма, который является модернизацией алгоритма, приведенного в [1].

entity LOG2 is
	generic(n: integer:=16;  --input width
		m: integer := 4;  -- m+q =output width
		q:integer :=12);
	port (X: in STD_LOGIC_VECTOR(n-1 downto 0);
		Y: out STD_LOGIC_VECTOR(m+q-1 downto 0));
end LOG2;

architecture LOGPRECISE of LOG2 is
	begin 
	LOGQUAD:process(X)
	     variable xs:STD_LOGIC_VECTOR(n-1 downto 0);
	     variable xl:STD_LOGIC_VECTOR(q downto 0);
	     variable xq:STD_LOGIC_VECTOR(2*q+1 downto 0);
	     variable  p:STD_LOGIC_VECTOR(m-1 downto 0);  
	     variable yl:STD_LOGIC_VECTOR(q-1 downto 0);
	begin 
		xs:=X;
		p:=CONV_STD_LOGIC_VECTOR(n-1, m);
		for i in 1 to n loop
			exit when xs(n-1)='1';      
			xs:=xs(n-2 downto 0)&'0';
			p:=unsigned(p)-1;
		end loop;                   
		xl:=xs(n-1 downto n-q-1);
		for i in 1 to q loop
			xq:=UNSIGNED(xl)*UNSIGNED(xl);
			if   xq(2*q+1)='1' then
				yl(q-i):='1';
				xl:=xq(2*q+1 downto q+1);
			else yl(q-i):='0';
				xl:=xq(2*q downto q);
     		        end if;
		end loop;
		Y<=p&yl;
	end process;
end LOGPRECISE;

Этот алгоритм может быть странслирован в логическую схему. Но для добавления q точных разрядов будут синтезированы q последовательно включенных умножителей.

Так, схема для вычисления 12 - разрядного логарифма от 16 - разрядного целого имеет затраты 1100 LUT и задержку 200 нс. Если алгоритм преобразовать так, чтоб он выполнялся на одном регистре сдвига, счетчике и умножителе за n+q итераций, то аппаратурные затраты уменьшатся в q раз.

В отличие от других алгоритмов вычисления логарифма, данный алгоритм не содержит никаких констант, кроме начального состояния переменной p.

Кроме того, количество точных разрядов результата m+q < n ничем не ограничено, так что алгоритм может быть использован при вычислениях со сверхвысокой точностью.

В приложении 2 представлен этот же алгоритм, реализованный как функция LOG2 для целых аргументов. Эта функция оформлена в пакете, который может быть помещен в рабочую библиотеку. Так что пользователь может в любом месте использовать вызов этой функции в своих программах, не помещая тело функции в описание архитектуры. Функция LOG2 дает 31 точный двоичный разряд результата, не считая ошибки в младшем разряде, что по крайней мере точнее, чем при вычислениях с плавающей запятой с одинарной точностью.

Конечно, проще всего использовать пакет MATH_REAL из библиотеки IEEE и сослаться на его функцию:

function LOG2(X: in REAL) return REAL;

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

Литература
1. Благовещенский Ю.В., Теслер Г.С. Вычисление элементарных функций на ЭВМ. Киев: Техніка. -1977. - 207с.

Приложение 1.

--^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
--  FILE: LOG2.VHD
--  PROJECT:     Logarithmic Function Unit
--  AUTHOR:       Anatoli Sergyienko
--  Email:        aser@comsys.kpi.ua
--
--  FUNCTION: - synthesable model for
--   calculating the binary logarithm function 
--  	from the positive integer data.
--  CONSTRAINTS: the input data X is in range:
--   1 to 2**n-1, 
--   generics 2**(m-1) log<="0000";
			when "0001" => log<="0010";
			when "0010" => log<="0011";
			when "0011" => log<="0100";
			when "0100" => log<="0101";
			when "0101" => log<="0111";
			when "0110" => log<="1000";
			when "0111" => log<="1001";
			when "1000" => log<="1010";
		        when "1001" => log<="1011";
			when "1010" => log<="1011";
			when "1011" => log<="1100";
			when "1100" => log<="1101";
			when "1101" => log<="1110";
			when "1110" => log<="1111";
			when "1111" => log<="1111";
			when others => 
   	         end case;
	end process;

RGS:process(CLK,log,yh)
		begin
		if CLK='1' and CLK'event then
		 Y<=yh & log;
		 xi<=X;
		end if;
	end process;
end DSPLOG;

Приложение 2.

--^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
--  FILE: LOG_FUNC.VHD
--  PROJECT:     Logarithmic function
--  AUTHOR:       Anatoli Sergyienko
--  Email:        aser@comsys.kpi.ua
--
--  FUNCTION: - calculating the binary logarithmic
--        function from the positive integer data.
--  CONSTRAINTS: the input data X is in range:
--   1 to 2**31-1, 
--   the output data is equal to  2**25*log2(X),
--   when X= 0 then the result is 0,  
--   and error report is given.  
--   Maximum error is equal to 1 l.s.b.
--^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
library IEEE;
use IEEE.std_logic_1164.all;
use IEEE.std_logic_arith.all;
package LOG_FUNC is
    function LOG2(X: INTEGER) return INTEGER;
end LOG_FUNC;

library IEEE;
use IEEE.std_logic_1164.all;
use IEEE.std_logic_arith.all;
package body LOG_FUNC is
  	function LOG2(x: INTEGER) return INTEGER is
	  variable xs: STD_LOGIC_VECTOR (30 downto 0);
	  variable xp: STD_LOGIC_VECTOR (26 downto 0);
	  variable xq: STD_LOGIC_VECTOR (53 downto 0);   
	  variable p: STD_LOGIC_VECTOR (4 downto 0);
	  variable yp: STD_LOGIC_VECTOR (25 downto 0);
	begin 
	   	assert  not(x=0)
		report "operand has value 0"
		severity ERROR;
		xs:=CONV_STD_LOGIC_VECTOR(x, 31); 
         	yt:="11110";
          for i in 1 to 31 loop
			exit when xs(30)='1';
			xs:=xs(29 downto 0)&'0';
			p:=unsigned(p)-1;
		end loop;                   
		xp:=xs(30 downto 4);
		for i in 1 to 26 loop
			xq:=UNSIGNED(xp) * UNSIGNED(xp);
			if   xq(53)='1' then
				yp(26-i):='1';
				xp:=xq(53 downto 27); 
			else    
				yp(26-i):='0';
				xp:=xq(52 downto 26); 
			end if;
		end loop;   
        if x=0 then return 0;
        else
		return CONV_INTEGER(UNSIGNED(p & yp));
        end if;
	end;      
    end package body;