Звезда не активнаЗвезда не активнаЗвезда не активнаЗвезда не активнаЗвезда не активна
 

  (заметка будет дописываться)

Цель: аппроксимация функции asin целочисленными линейными полиномами совместно с табличным методом с целью увеличения скорости вычислений и уменьшения расхода памяти при вычислениях на процессоре с фиксированной точкой PIC18F. Нотация представления 1.0 float eq 16384 int.

1. Разобьем функцию в диапазоне аргументов x=  0.. 1.0 на более мелкие участки, длина которых обратно пропорциональна log2N. Для этого используем MATLAB.

Весь участок выглядит так, как показано на рисунке 1.

Арксинус

Рис.1 Арксинус (в радианах) на участке  x=  0.. 1.0. arcsin(x). 

M-файл:

clear;
indx=0;
for argum= 0.00:1./16384.:16384./16384.
%нумерация для массивов, которые будем использовать для аппроксимации
insideidx0 = mod(indx,8192)+1;
insideidx1 = mod(indx,4096)+1;
insideidx2 = mod(indx,2048)+1;
insideidx3 = mod(indx,1024)+1;
insideidx4 = mod(indx,512)+1;
insideidx5 = mod(indx,256)+1;
insideidx6 = mod(indx,128)+1;
insideidx7 = mod(indx,64)+1;
insideidx8 = mod(indx,32)+1;
 %сложность аппроксикации возникает в конце, где возрастает нелинейность, переворачиваем нумерацию для построения условий
indx_c = 16384 - indx ;
if(indx_c > 8192)
  asinout_0(insideidx0) = asin(argum);
  xcoord_0(insideidx0) = argum;
else
  if(indx_c > 4096)
  asinout_1(insideidx1) = asin(argum);
  xcoord_1(insideidx1) = argum;
  else
  if(indx_c > 2048)
  asinout_2(insideidx2) = asin(argum);
  xcoord_2(insideidx2) = argum;
  else
  if (indx_c > 1024)
  asinout_3(insideidx3) = asin(argum);
  xcoord_3(insideidx3) = argum;
  else
  if (indx_c > 512)
  asinout_4(insideidx4) = asin(argum);
  xcoord_4(insideidx4) = argum;
  else
if (indx_c > 256)

  asinout_5(insideidx5) = asin(argum);

  xcoord_5(insideidx5) = argum;

 

  else

 

  if (indx_c > 128)

  asinout_6(insideidx6) = asin(argum);

  xcoord_6(insideidx6) = argum;

 

  else

 

  if (indx_c > 64)

  asinout_7(insideidx7) = asin(argum);

  xcoord_7(insideidx7) = argum;

 

else

 

  if (indx_c > 32)

  asinout_8(insideidx8) = asin(argum);

  xcoord_8(insideidx8) = argum;

font-size: 10.0pt; font-family: 'Courier New'; color: blue; background: white; mso-highlight: white; mso-ansi-language: EN-US;

 

 end

 end

 end

 end

 end

 end

 end

end

 

end

 

indx = indx+1;

full_asin(indx)= asin(argum);

xcoord_full(indx) = argum;

 

end

Теперь мы имеем 8 массивов с длинами, кратными 2n : 8192, 4096, 2048, 1024, 512, 256, 128, 64, 32. Внутри этих участков функция ведёт себя близко к линейной.

2. Используем Curve Fitting Tool для аппроксимации полиномов (1-й степени):
Полином первой степени имеет вид:

f(x) = p1*x + p2

>> cftool

 

Рис.2 Настройки аппроксимации

 

Результат аппроксимации предпоследнего из участков на рис.3:

 

Линейная аппроксимация

Рис.3 Результат линейной аппроксимации предпоследнего из участков 

SSE (сумма квадратов ошибок) = 1.35*10-5

SSE: 1.35e-005
R-square: 0.9981
Adjusted R-square: 0.9981
RMSE: 0.0004667

Для участка аргументов в целочисленной нотации  16256... 16384

 Для остальных участков получили:

Участки целочисленного параметра float приближение коэффициентов полинома
0..8191 p1 = 1.041 (1.041, 1.042)
p2 = -0.004679 (-0.004799, -0.00456)
(в целочисленной нотации 17056 и -77)
Goodness of fit:
SSE: 0.06196
R-square: 0.9997
Adjusted R-square: 0.9997
RMSE: 0.002751
8192...12287 p1 = 1.291 (1.29, 1.292)
p2 = -0.1282 (-0.129, -0.1273)
 (в целочисленной нотации 21152 и -2100)  
Goodness of fit:
SSE: 0.04157
R-square: 0.9988
Adjusted R-square: 0.9988
RMSE: 0.003187
12288... 14335 p1 = 1.721 (1.718, 1.724)
p2 = -0.4479 (-0.4504, -0.4453)
(в целочисленной нотации 28197 и -7338)
Goodness of fit:
SSE: 0.01376
R-square: 0.9983
Adjusted R-square: 0.9983
RMSE: 0.002593
  14336... 15359 p1 = 0.04305 (0.04294, 0.04316)
p2 = 1.136 (1.136, 1.136)
(в целочисленной нотации 705 и 18612)
Goodness of fit:
SSE: 0.003348
R-square: 0.9982
Adjusted R-square: 0.9982
RMSE: 0.00181
                 15360... 15871 p1 = 0.03009 (0.02997, 0.0302)
p2 = 1.265 (1.265, 1.265)
(в целочисленной нотации 493 и 20726)
Goodness of fit:
SSE: 0.0008592
R-square: 0.9981
Adjusted R-square: 0.9981
RMSE: 0.001298
                15872... 16127 p1 = 0.02116 (0.02105, 0.02128)
p2 = 1.355 (1.355, 1.355)
(в целочисленной нотации 347 и 22200)
Goodness of fit:
SSE: 0.000217
R-square: 0.9981
Adjusted R-square: 0.9981
RMSE: 0.0009244
              16128... 16255 p1 = 0.01494 (0.01482, 0.01506)
p2 = 1.418 (1.418, 1.418)
(в целочисленной нотации 245 и 23233)
Goodness of fit:
SSE: 5.434e-005
R-square: 0.9981
Adjusted R-square: 0.9981
RMSE: 0.0006567
16256...16319

 p1 = 0.01056 (0.01044, 0.01069)
 p2 = 1.463 (1.463, 1.463)
(в целочисленной нотации 173 и 23970)
SSE: 1.35e-005
R-square: 0.9981
Adjusted R-square: 0.9981
 RMSE: 0.0004667

 16320...16351

p1 =     0.00751  (0.007388, 0.007631)
p2 =       1.494  (1.494, 1.494)
(в целочисленной нотации 123 и 24478)
Goodness of fit:
  SSE: 3.313e-006
  R-square: 0.9981
  Adjusted R-square: 0.998
  RMSE: 0.0003323

3. Чтобы не разбивать далее последовательность на участки с ещё меньшей длиной "хвост" (последние 33 отсчётов (не 32, т.к. границы включены (нуль и единица) )) аппроксимируем с помощью табличных значений.

16352 => 1.5083
16353 => 1.5093
16354 => 1.5103
16355 => 1.5113
16356 => 1.5123
16357 => 1.5134
16358 => 1.5145
16359 => 1.5155
16360 => 1.5167
16361 => 1.5178
               1.5190 
               1.5202
               1.5214
               1.5226
               1.5239
               1.5252
               1.5266
               1.5280
               1.5295
               1.5310
               1.5325
               1.5342
               1.5359
               1.5376
               1.5395
               1.5416
               1.5437
               1.5461
               1.5487
               1.5517
               1.5552
               1.5597
16384 => 1.5708 

 В целочисленном виде, начиная с конца (для удобства имплементации на языке Си (смотрите ниже)):

Добавим переменную в начало программы:
insideidx9 =1;

Модифицируем участок программы:
if (indx_c > 32)
asinout_8(insideidx8) = asin(argum);
xcoord_8(insideidx8) = argum;

заменим на
if (indx_c > 32)
asinout_8(insideidx8) = asin(argum);
xcoord_8(insideidx8) = argum;
else

FixedOutput(insideidx9)= round(asin(argum)*16384.0);

 

 Получим
25736,     25555,     25480,     25422,     25374,     25331,     25293,     25257,     25224,     25193,     25163,     25136,     25109,     25083,     25059,     25035,     25012,     24990,     24968,     24947,     24926,     24906,     24887,     24868,     24849,     24831,     24813,     24795,     24778,     24761,     24744,     24728,     24712

Си имплементация может в этом случае выглядеть примерно так:
---------------------------------------------------------------------------
static union pol
{
  int16_t halfs_inp[2];
  int32_t inp;
}
result;

int16_t fixed_asin(int16_t x){
  if(x&16384)return get_table(x);
  if(!(x&8192))return get_poly(x, 0);
  if(!(x&4096))return get_poly(x, 1);
  if(!(x&2048))return get_poly(x, 2);
  if(!(x&1024))return get_poly(x, 3);
  if(!(x&512))return get_poly(x, 4);
  if(!(x&256))return get_poly(x, 5);
  if(!(x&128))return get_poly(x, 6);
  if(!(x&64))return get_poly(x, 7);
  if(!(x&32))return get_poly(x, 8);
return get_table(x);
}

int16_t get_poly(int16_t x, int8_t num){
const int16_t p1[9] = {17056,21152,28197,705, 493, 347, 245, 173, 123};
const int16_t p2[9] = {-77, -2100,-7338,18612,20726, 22200,23233,23970,24478};
result.inp = x;
result.inp*=p1[num];
result.inp <<= 2;
result.halfs_inp[1]+=p2[num];
 return result.halfs_inp[1];
}

int16_t get_table(int16_t x){
const int16_t x_num[33]={25736,25555,25480,25422,25374,
25331,25293,25257,25224,25193,
25163,25136,25109,25083,25059,
25035,25012,24990,24968,24947,
24926,24906,24887,24868,24849,
24831,24813,24795,24778,24761,
24744,24728,24712};
return x_num[16384-x];
}

 Однако, погрешность при такой кусочно-линейной аппроксимации значительная.

---------------------------------------------------------------------------

Ещё аппроксимации из сети (внимание, не на все

х участках сходимость может быть хорошая)...

tan(x) = (x - x*x*x/9) / (1 - 4*x*x/9 + x*x*x*x/63);

float facos( float x ) {
return 0.93 * (1.68902830838 + x * (x * x - 1.07526881720));
}

float fasin( float x ) {
return 0.25f * x * (4.0 + x * x);
}

 

-----------

VHDL:

function arcsin(X:double):double;
Type ArrND=array[1..30] of double;
var I:integer;
P,Q:double;
A, B:ArrND;
begin
A[1]:= 1.95387311727746E-0005;
A[2]:=-1.52852221394629E+0000;
A[3]:=-2.01776736316873E-0003;
A[4]:= 7.76296426337929E-0001;
B[1]:=-1.52652894854093E+0000;
B[2]:=-1.74609738502257E-0003;
P:=A[4];
for i:=3 downto 1 do P:=x*P+A[I];
Q:=X+B[2];
for I:=1 downto 1 do Q:=x*Q+B[I];
Result:=P/Q;
end;

nbsp;

х участках сходимость может быть хорошая)...nbsp; font-size: 10.0pt; font-family: 'Courier New'; color: blue; background: white; mso-highlight: white; mso-ansi-language: EN-US; RMSE: 0.0009244