(заметка будет дописываться)
Цель: аппроксимация функции 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
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) |
| 16320...16351 |
p1 = 0.00751 (0.007388, 0.007631) |
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 |
|
В целочисленном виде, начиная с конца (для удобства имплементации на языке Си (смотрите ниже)):
Добавим переменную в начало программы:
Модифицируем участок программы:
заменим на FixedOutput(insideidx9)= round(asin(argum)*16384.0);
|
|
Получим |
Си имплементация может в этом случае выглядеть примерно так:
---------------------------------------------------------------------------
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