О пользе округления
Анатолий Сергиенко
Все основные операции в вычислителях с фиксированной запятой выполняются точно. Исключением является усечение n+m -разрядного операнда до n старших разрядов. Мы знаем,что финансовые расчеты выполняются точно, но при начислении процентов m разрядов, означающих доли копейки, при усечении отбрасываются и результат оказывается приближенным.
Правда, в добросовестном учреждении эти доли копейки собираются вместе и распределяются между вкладчиками по случайному закону или к сумме вкладчика прибавляется копейка, если отбрасываемое число означает больше, чем полкопейки.Так осуществляется округление.
Повсеместное использование микропроцессоров отучило разработчиков средств вычислительной техники анализировать ошибки вычислений, их зависимость от выбранной разрядности арифметики, так как выбирать-то особенно не приходится.
Микропроцессоры предоставляют любую разрядность данных из 8, 16 или 32, иногда 24 или 64. Так что хорошо загруженный 16-разрядный сигнальный микропроцессор даст такой динамический диапазон обрабатываемого сигнала, на какой он способен, то есть 60-90 дб, но не более. А если требуется поднять динамический диапазон, то придется перейти на обработку слов двойной длины с неизбежным резким падением пропускной способности.
Характерной чертой ПЛИС является то, что в них можно задавать произвольную разрядность операндов на каждом из этапов их обработки. Минимизировав эту разрядность до допустимых пределов можно как уменьшить аппаратурные затраты, так и повысить быстродействие. Выбором удачного способа округления можно дополнительно минимизировать разрядность или повысить точность вычислений.
Языки описания вычислительных схем VHDL и Verilog дают возможность моделировать эти схемы с варьируемыми параметрами разрядности. Таким образом, после серии экспериментов с моделями вычислителя можно выбрать нужную разрядность. Однако, если алгоритм обработки сложный, то таких экспериментов должно быть очень много.
Можно пойти другим путем и заранее аналитически оценить необходимую разрядность данных на всех этапах обработки. Наиболее эффективной по точности и трудоемкости является статистическая методика определения ошибки вычислений [1]. Эта методика основана на следующих упрощениях и положениях.
При усечении и округлении n+m -разрядного операнда до n старших разрядов возникает случайная ошибка с некоторым матожиданием M и cреднеквадратическим отклонением q от матожидания. То есть средняя ошибка равна M±q. Чаще всего ошибка имеет равномерную вероятность распределения, тогда максимальная ошибка равна M±SQRT(3)*q.
У всех способов округления практически одинаковая среднеквадратическая ошибка q, которая несколько зависит от m, и которую можно оценить как sqrt(1/12) от веса младшего разряда числа, полученного после округления. Квадрат среднеквадратической ошибки назывется дисперсией. Дисперсии операций округления можно масштабировать и складывать в соответствии с графом алгоритма вычислений, получая в конце дисперсию ошибки результата.
Однако, с оценкой матожидания результата вычислений дело сложнее, так как оно зависит от характера вычислений, исходных данных и, не в последнюю очередь, от способа округления.
Поэтому в данной заметке идет речь об оценке матожидания при различных способах округления и выборе лучшего способа округления.
Способы округления будут рассматриваться на примере округления 12-разрядного целого числа a со знаком в дополнительном коде до 8-разрядного числа b.
Округлить — значит подобрать такое число b , которое с минимальной ошибкой соответствовало бы числу а , то есть b приблизительно равно K*а, где K — масштабный коэффициент.В данном случае K=1/16.
Способы округления описываются на языке VHDL с применением библиотеки IEEE.std_logic_arith.
Округление с недостатком или округление к минус бесконечности , судя по названию, сводится к выбору ближайшего меньшего числа после отбрасывания младших разрядов. Это наиболее простой способ и реализуется как:
TO_NEGATIVE: b<=a(11 downto 4);
По сути это просто усечение.
Округление с избытком или округление к плюс бесконечности сводится к выбору ближайшего большего числа и выполняется при добавлении 1 к младшему разряду числа после усечения:
TO_POSITIVE: b<=SIGNED(a(11 downto 4))+1;
Округление к нулю означает выбор после усечения числа, ближайшего к нулю. Для положительных чисел это просто усечение, а для отрицательных — необходимо добавление единицы к младшему разряду:
TO_ZERO:process(a) begin if a(11)='1' then b<=SIGNED(a(11 downto 4))+1; else b<=a(11 downto 4); end if; end process;
Округление от нуля означает выбор после усечения числа, по модулю ближайшего к плюс бесконечности. Для отрицательных чисел это просто усечение, а для положительных - необходимо добавление единицы к младшему разряду:
FROM_ZERO:process(a)
begin
if a(11)='1' then
b<=a(11 downto 4);
else
b<=SIGNED(a(11 downto 4))+1;
end if;
end process;
Округление к ближайшему целому означает выбор после усечения ближайшего целого числа, если идет речь об округлении целых чисел. Это тот самый способ округления, который знаком всем со школы.
Так осуществляется округление.
Для двоичных чисел это означает прибавление 1 к старшему отбрасываемому разряду. Это то же самое, что прибавить 1 к младшему разряду после усечения, если старшая отбрасываемая цифра равна 1:
TO_INTEGER:process(a) variable one:std_logic; begin if a(3)='1' then b<=SIGNED(a(11 downto 4))+1; else b<=a(11 downto 4); end if; end process;
Уточненное округление к ближайшему целому. Вышеприведенное округление к ближайшему целому не всегда правильно выполняется для отрицательных чисел. Например, код 111111101000, равный -24/16= -1.5, после округления должен дать ближайшее целое -2, или 11111110 а в действительности получится 11111111, или -1.
Такая ситуация возникает всегда, когда старшая отбрасываемая цифра равна 1, а все цифры после нее равны 0. В этом случае для правильности округления 1 не прибавляется. Уточненное округление к ближайшему целому реализуется по следующему алгоритму:
TO_INTEGER_PRECISE:process(a)
variable one:std_logic;
begin
one:= a(2) or a(1) or a(0);
if ( a(11)='0' and a(3)='1' )
or ( a(11)='1' and a(3)='1' and one='1')
then
b6<=SIGNED(a(11 downto 4))+1;
else
b6<=a(11 downto 4);
end if;
end process;
Если округление выполняется с отбрасыванием большого числа разрядов, то ситуация "неправильного" округления встречается редко и тогда уточненное округление можно не делать.
Упрощенное округление к ближайшему целому. По сравнению с просто усечением для прибавления 1 при округлении обычно требуется дополнительный сумматор. Иногда применяют упрощенное округление, которое заключается в установке 1 в младшей цифре результата, если старшая цифра отбрасываемого кода равна 1. Тогда требуется всего лишь одна двухвходовая схема ИЛИ:
TO_INTEGER_SIMPLE: b<=a(11 downto 5)&(a(4) or a(3));
Ошибка упрощенного округления больше, чем ошибка округления к ближайшему целому, но меньше, чем округление к плюс бесконечности.
Вероятностное округление. При вычислении рекурсивных фильтров часто возникает такое явление, когда фильтр начинает генерировать заметный низкочастотный периодический шум. Такое явление возникает из-за нелинейности операции округления результатов фильтрации. Одним из действенных способов устранения такой генерации является вероятностное округление.
Оно заключается в случайном прибавлении единицы округления и может быть реализовано, например, по такому алгоритму:
PLUS_RANDOM: process(rst, clk,a) variable r: std_logic_vector(10 downto 0); variable one:std_logic; begin if rst='1' then r:=(others=>'1'); elsif rising_edge(clk) then r:=r(9 downto 0)&(r(10)xor r(8)); end if; if r(10)='1' then b8<=SIGNED(a(11 downto 4))+1; else b8<=a(11 downto 4); end if; end process;
Здесь переменная r генерируется с помощью генератора псевдослучайных чисел с равномерным законом распределения и управляет прибавлением 1.
Генератор псевдослучайных чисел построен на регистре сдвига с отводами и схеме Исключающее ИЛИ.
Как строить такие генераторы на ПЛИС, описано, например, в [2]. Данный генератор, построенный на 11 - разрядном регистре, выдает псевдослучайную последовательность длиной 211-1 с вероятностью появления единицы, равной 0.5.
Таким образом, к результату прибавляется, в среднем, число 0.5, что равно необходимому смещению для матожидания, которое в этом случае стремится к нулю.
Часто вероятностное округление применяют при квантовании данных в сигма-дельта АЦП при квантовании исходных данных во время вычисления корреляционных функций.
В этом случае заметно повышается точность результатов за счет фильтрации шума округления при дальнейшем накоплении суммы взвешенных данных.
Немного о реализации округления. Все вышеприведенные алгоритмы округления описаны синтезируемым стилем и требуют небольшого объема оборудования для своей реализации.
Так, для прибавления 1 требуется лишь схема распространения переноса. Можно было бы описать проще округление к ближайшему целому:
TO_INTEGER: b<=SIGNED(a(11 downto 4))+а(3);
однако компилятор-синтезатор распознает здесь схему полного сумматора, которая вдвое сложнее, чем схема распространения переноса.
Округление с прибавлением 1 все еще редко применяется в специализированных устройствах из-за мнения, что для такого округления требуется дополнительное оборудование или время.
В следующем примере показано, как выполнять округление к ближайшему целому, не добавляя оборудования. Пусть необходимо синтезировать накапливающий сумматор (интегратор), выполняющий сложение отсчетов сигнала A с накоплением суммы B, которая по сигналу RND округляется до результата C. Такое устройство можно описать следующей моделью:
entity ACC is
port( CLK: in std_logic;
RST: in std_logic;
RND: in std_logic;
A: in std_logic_vector(11 downto 0);
C: out std_logic_vector(7 downto 0));
end entity;
architecture ROUND of ACC is
signal B: std_logic_vector(11 downto 0);
begin
ACCUM:process(CLK,RST, RND,A,В)
variable ta,tb,tc: std_logic_vector(12 downto 0);
begin
ta(12 downto 5):=A(11 downto 4);
tb(12 downto 5):=B(11 downto 4);
ta(3 downto 0):=A(3 downto 0);
tb(3 downto 0):=B(3 downto 0);
ta(4):=RND;
tb(4):='1';
tc:= SIGNED(ta)+SIGNED(tb);
C<= tc(12 downto 5);
if RST='1' then -- register of the accumulator
B<=(others=>'0');
elsif rising_edge(CLK)
then B<= tc(12 downto 5)&tc(3 downto 0);
end if;
end process;
end ROUND;
В этой схеме четвертый разряд переменных ta,tb,tc введен искусственно для того, чтобы синтезатор распознал, что к четвертому разряду результата сложения нужно дополнительно прибавлять единицу округления.
Для проверки и оценки работы различных способов округления была составлена экспериментальная программа, приведенная в приложении. В ней процесс RANDOM генерирует последовательность случайных 12-разрядных данных а. При этом используется процедура UNIFORM из пакета IEEE.math_real для генерации случайных чисел с равномерным распределением в диапазоне 0,...,1.0.
Последовательность чисел а округляется восемью способами с получением последовательностей 8 - разрядных чисел b1,...,b8. Всего в каждой последовательности генерируется N=1 млн. чисел.
Среднеарифметическое значение достаточно большого количества знакопеременных чисел равно матожиданию этих чисел [1]. В программе числа из последовательностей bi складываются тремя способами в процессе ACCUM. Результаты сложения, будучи поделенными на N, дают оценки матожидания М+, М и М-.
Здесь М+ получено при сложении только положительных отсчетов с последующим вычитанием суммы этих же отсчетов, но не округленных, и является оценкой матожидания для вычислений с положительными числами. Оценка матожидания М получена при сложении знакопеременных чисел и является просто матожиданием.
Значение М- получено при суммировании четных отсчетов со знаком +, а нечетных - со знаком -. Во многих алгоритмах цифровой обработки сигналов, например, в быстрых ортогональных преобразованиях, выполняется парная операция: C=A+B; D=A-B. И если числа А,В - случайные, то результаты усреднения С и D должны быть приблизительно одинаковыми. Однако, из-за ошибки округления С и D могут отличаться на удвоенную величину матожидания.
Таким образом, М равно матожиданию результата С, а М- матожиданию результата D. Поэтому для быстрых ортогональных преобразований и других алгоритмов цифровой обработки сигналов необходимо выбирать такой способ округления, чтобы М и М- были близкими по величине.
Так, если быстрое преобразование Фурье реализовано с плохим округлением, то по этой причине при обработке сигнала ошибки округления отображаются в результирующем графике спектра как заметные ступеньки.
С другой стороны, поскольку М- получено при знакопеременном сложении, то в результате матожидания чисел взаимно скомпенсировались и М- без знака можно принять как некоторую оценку среднеквадратического отклонения.
Полученные значения оценки матожидания М+, М и М-, выраженные в доле веса младшего разряда числа после округления, приведены в следующей таблице:
Анализ таблицы показывает следующее. Округление может дать выигрыш в 2-4 двоичных разряда по сравнению с просто усечением. Для вычислений со знакопеременными числами желательно применять все способы округления, кроме округления с недостатком и округления с избытком. Если в алгоритме есть обработка положительных чисел, то также плохо подходят округление к нулю и округление от нуля. Упрощенное округление к ближайшему целому дает на порядок худшие результаты, чем, например, округление к ближайшему целому, но все же его результаты имеют вдвое меньшую ошибку, чем округление с недостатком, то есть просто усечение.
Наибольшую точность имеет уточненное округление к ближайшему целому и оно наиболее подходит для алгоритмов с операциями типа C=A+B; D=A-B и к любым другим, в которых необходимо накопление суммы большого числа операндов. Например, в интеграторах и чувствительных рекурсивных фильтрах. Интересным фактом является то, что вероятностное округление по точности приближается к округлению к ближайшему целому.
Литература Приложение
Способ округления
M+
M
M-
Округление с недостатком
-0,469
-0,448
-0,049
Округление с избытком
0,531
0,553
-0,049
Округление к нулю
-0,469
0,052
-0,048
Округление от нуля
0,531
0,053
-0,049
Округление к ближайшему целому
0,032
0,053
-0,049
Уточненное округление к ближайшему целому
0,032
0,022
-0,049
Упрощенное округление к ближайшему целому
-0,217
-0,197
-0,048
Вероятностное округление
0,032
0,053
-0,049
1. Соренков Э.И., Телига А.И., Шаталов А.С. Точность вычислительных устройств и алгоритмов. М.: -Машиностроение.-1976.-200с.
2. http:/www.xilinx.com/xapp/xapp052.pdf
library IEEE;
use IEEE.STD_logic_1164 .all;
use IEEE.STD_logic_arith.all;
use IEEE.math_real.all;
entity RNDER is
end RNDER;
architecture RNDER of RNDER is
signal a: std_logic_vector(11 downto 0);
signal b1,b2,b3,b4,b5,b6,b7,b8: std_logic_vector(7 downto 0);
signal n,acc1p,acc1m,acc2p,acc2m,acc3p,acc3m,acc4p,acc4m:integer;
signal acc5p,acc5m,acc6p,acc6m,acc7p,acc7m,acc8p,acc8m:integer;
signal accs,acc1s,acc2s,acc3s,acc4s,acc5s,acc6s,acc7s,acc8s:integer;
signal acc1,acc2,acc3,acc4,acc5,acc6,acc7,acc8:integer;
signal clk:std_logic;
signal rst,s,c: std_logic;
begin
RESET:process
begin
rst<='1';
wait for 11 ns;
rst<='0';
wait;
end process;
RANDOM:process(rst, clk)
variable r: std_logic_vector(11 downto 0);
variable x: real;
variable seed1,seed2: integer;
begin
UNIFORM(seed1,seed2, x); --random number generator
-- in range (0...+1.0)
if rst='1' then r:=(others=>'1');
elsif RISING_EDGE(clk) then
r:= CONV_STD_LOGIC_VECTOR(INTEGER((x-0.5)*4000.0),12);
end if;
a<=r; --random number in range (-0.99...+0.99)*2048
s<=r(11); -- sign of number a
end process;
TO_NEGATIVE: b1<=a(11 downto 4); --truncation
TO_POSITIVE: b2<=SIGNED(a(11 downto 4))+1;
TO_ZERO:process(a)
begin
if a(11)='1' then
b3<=SIGNED(a(11 downto 4))+1;
else
b3<=a(11downto 4);
end if;
end process;
FROM_ZERO:process(a)
begin
if a(11)='1' then
b4<=a(11 downto 4);
else b4<=SIGNED(a(11 downto 4))+1;
end if;
end process;
TO_INTEGER:process(a)
variable one:std_logic;
begin
if a(3)='1' then
b5<=SIGNED(a(11 downto 4))+1;
else
b5<=a(11 downto 4);
end if;
end process;
TO_INTEGER_PRECISE:process(a)
variable one:std_logic;
begin
one:= a(2) or a(1) or a(0);
if ( a(11)='0' and a(3)='1' )
or ( a(11)='1' and a(3)='1' and one='1')
then
b6<=SIGNED(a(11 downto 4))+1;
else
b6<=a(11 downto 4);
end if;
end process;
TO_INTEGER_SIMPLE: b7<=a(11 downto 5)&(a(4) or a(3));
PLUS_RANDOM: process(rst, clk,a)
variable r: std_logic_vector(10 downto 0);
variable one:std_logic;
begin
if rst='1'
then r:=(others=>'1');
elsif rising_edge(clk)
then r:=r(9 downto 0)&(r(10)xor r(8));
end if;
if r(10)='1'
then b8<=SIGNED(a(11 downto 4))+1;
else b8<=a(11 downto 4);
end if;
end process;
ACCUM:process (rst, clk)
variable t:bit;
begin
if rst='1'
then
accs<=0;
acc1p<=0; acc1m<=0; acc1s<=0;
acc2p<=0; acc2m<=0; acc2s<=0;
acc3p<=0; acc3m<=0; acc3s<=0;
acc4p<=0; acc4m<=0; acc4s<=0;
acc5p<=0; acc5m<=0; acc5s<=0;
acc6p<=0; acc6m<=0; acc6s<=0;
acc7p<=0; acc7m<=0; acc7s<=0;
acc8p<=0; acc8m<=0; acc8s<=0;
n<=0; t:='0';
elsif rising_edge(clk)
then -- Calculating M
acc1p<=acc1p+CONV_INTEGER(SIGNED(b1));
acc2p<=acc2p+CONV_INTEGER(SIGNED(b2));
acc3p<=acc3p+CONV_INTEGER(SIGNED(b3));
acc4p<=acc4p+CONV_INTEGER(SIGNED(b4));
acc5p<=acc5p+CONV_INTEGER(SIGNED(b5));
acc6p<=acc6p+CONV_INTEGER(SIGNED(b6));
acc7p<=acc7p+CONV_INTEGER(SIGNED(b7));
acc8p<=acc8p+CONV_INTEGER(SIGNED(b8));
if t='0' then -- Calculating M-
acc1m<=acc1m +CONV_INTEGER(SIGNED(b1));
acc2m<=acc2m +CONV_INTEGER(SIGNED(b2));
acc3m<=acc3m +CONV_INTEGER(SIGNED(b3));
acc4m<=acc4m +CONV_INTEGER(SIGNED(b4));
acc5m<=acc5m +CONV_INTEGER(SIGNED(b5));
acc6m<=acc6m +CONV_INTEGER(SIGNED(b6));
acc7m<=acc7m +CONV_INTEGER(SIGNED(b7));
acc8m<=acc8m +CONV_INTEGER(SIGNED(b8));
else
acc1m<=acc1m -CONV_INTEGER(SIGNED(b1));
acc2m<=acc2m -CONV_INTEGER(SIGNED(b2));
acc3m<=acc3m -CONV_INTEGER(SIGNED(b3));
acc4m<=acc4m -CONV_INTEGER(SIGNED(b4));
acc5m<=acc5m -CONV_INTEGER(SIGNED(b5));
acc6m<=acc6m -CONV_INTEGER(SIGNED(b6));
acc7m<=acc7m -CONV_INTEGER(SIGNED(b7));
acc8m<=acc8m -CONV_INTEGER(SIGNED(b8));
end if;
if s='0' then -- Calculating M+
accs<=accs+CONV_INTEGER(SIGNED(a)); --Exact sum of dates
acc1s<=acc1s +CONV_INTEGER(SIGNED(b1));
acc2s<=acc2s +CONV_INTEGER(SIGNED(b2));
acc3s<=acc3s +CONV_INTEGER(SIGNED(b3));
acc4s<=acc4s +CONV_INTEGER(SIGNED(b4));
acc5s<=acc5s +CONV_INTEGER(SIGNED(b5));
acc6s<=acc6s +CONV_INTEGER(SIGNED(b6));
acc7s<=acc7s +CONV_INTEGER(SIGNED(b7));
acc8s<=acc8s +CONV_INTEGER(SIGNED(b8));
end if;
n<=n +1;
t:= not t;
end if;
end process;
CLK_GEN: process --clock generator which stops
-- when end of simulation
begin
if n=1000000 then -- Calculating M+
acc1<=acc1s -accs/16;
acc2<=acc2s -accs/16;
acc3<=acc3s -accs/16;
acc4<=acc4s -accs/16;
acc5<=acc5s -accs/16;
acc6<=acc6s -accs/16;
acc7<=acc7s -accs/16;
acc8<=acc8s -accs/16;
wait for 100ns;
assert false
report "***END OF SIMULATION***"
severity FAILURE;
wait;
else
CLK<= '1';
wait for 10ns;
CLK<= '0';
wait for 10ns;
end if;
end process;
end rnder;