Możesz użyć logarytmów, aby pozbyć się podziału. Dla ( x , y) w pierwszej ćwiartce:
z= log2)( y) - log2)( x )atan2 ( y, x ) = atan ( y/ x)=atan( 2z)
Rycina 1. Wykres atan ( 2z)
Aby uzyskać wymaganą dokładność 1E-9, musisz zbliżyć atan ( 2z) w zakresie - 30 < z< 30 . Możesz skorzystać z symetrii atan ( 2- z) = π2)- atan ( 2z)lub alternatywnie upewnij się, że( x , y)jest w znanej liczbie oktanowej. Aby przybliżyćlog2)( ):
b = podłoga ( log2)( a ) )c = a2)blog2)( a ) = b + log2)( c )
b można obliczyć, znajdując lokalizację najbardziej znaczącego niezerowego bitu. do można obliczyć przez przesunięcie bitowe. Konieczne byłoby przybliżenielog2)( c ) w zakresie1 ≤ c < 2 .
Ryc. 2. Wykres log2)( c )
Dla twoich wymagań dokładności wystarczy interpolacja liniowa i jednorodne próbkowanie, 2)14+ 1 = 16385 próbek log2)( c ) i 30 × 212+ 1 = 122881 próbek atan ( 2z) dla 0 < z< 30 . Ten drugi stół jest dość duży. Dzięki niemu błąd spowodowany interpolacją zależy w dużej mierze od z :
Ryc. 3. największy błąd bezwzględny aproksymacji atan(2z) dla różnych zakresów z (oś pozioma) dla różnych liczb próbek (od 32 do 8192) na przedział jednostkowy z . Największy błąd bezwzględny dla 0≤z<1 (pominięty) jest nieco mniejszy niż dla floor(log2(z))=0 .
atan(2z) tabela może być podzielony na wiele subtables które odpowiadają 0≤z<1 i innym floor(log2(z)) z z≥1 , który jest łatwy do obliczenia. Długości tabeli można wybrać zgodnie z rys. 3. Indeks wewnątrz podtabeli można obliczyć za pomocą prostej manipulacji ciągiem bitów. Dla potrzeb dokładności atan(2z) będą zawierały łącznie 29217 próbek, jeśli rozszerzysz zakres z do 0≤z<32 dla uproszczenia.
Do późniejszego wykorzystania oto niezgrabny skrypt Pythona, którego użyłem do obliczenia błędów aproksymacji:
from numpy import *
from math import *
N = 10
M = 20
x = array(range(N + 1))/double(N) + 1
y = empty(N + 1, double)
for i in range(N + 1):
y[i] = log(x[i], 2)
maxErr = 0
for i in range(N):
for j in range(M):
a = y[i] + (y[i + 1] - y[i])*j/M
if N*M < 1000:
print str((i*M + j)/double(N*M) + 1) + ' ' + str(a)
b = log((i*M + j)/double(N*M) + 1, 2)
err = abs(a - b)
if err > maxErr:
maxErr = err
print maxErr
y2 = empty(N + 1, double)
for i in range(1, N):
y2[i] = -1.0/16.0*y[i-1] + 9.0/8.0*y[i] - 1.0/16.0*y[i+1]
y2[0] = -1.0/16.0*log(-1.0/N + 1, 2) + 9.0/8.0*y[0] - 1.0/16.0*y[1]
y2[N] = -1.0/16.0*y[N-1] + 9.0/8.0*y[N] - 1.0/16.0*log((N+1.0)/N + 1, 2)
maxErr = 0
for i in range(N):
for j in range(M):
a = y2[i] + (y2[i + 1] - y2[i])*j/M
b = log((i*M + j)/double(N*M) + 1, 2)
if N*M < 1000:
print a
err = abs(a - b)
if err > maxErr:
maxErr = err
print maxErr
y2[0] = 15.0/16.0*y[0] + 1.0/8.0*y[1] - 1.0/16.0*y[2]
y2[N] = -1.0/16.0*y[N - 2] + 1.0/8.0*y[N - 1] + 15.0/16.0*y[N]
maxErr = 0
for i in range(N):
for j in range(M):
a = y2[i] + (y2[i + 1] - y2[i])*j/M
b = log((i*M + j)/double(N*M) + 1, 2)
if N*M < 1000:
print str(a) + ' ' + str(b)
err = abs(a - b)
if err > maxErr:
maxErr = err
print maxErr
P = 32
NN = 13
M = 8
for k in range(NN):
N = 2**k
x = array(range(N*P + 1))/double(N)
y = empty((N*P + 1, NN), double)
maxErr = zeros(P)
for i in range(N*P + 1):
y[i] = atan(2**x[i])
for i in range(N*P):
for j in range(M):
a = y[i] + (y[i + 1] - y[i])*j/M
b = atan(2**((i*M + j)/double(N*M)))
err = abs(a - b)
if (i*M + j > 0 and err > maxErr[int(i/N)]):
maxErr[int(i/N)] = err
print N
for i in range(P):
print str(i) + " " + str(maxErr[i])
Lokalna maksymalna błędu z zbliżony funkcję f(x) przez interpolację liniową f ( x ) z próbek f ( x ) , wykonane z jednolitego próbek z próbek przedział Δ x może być w przybliżeniu analitycznie:f^(x)f(x)Δx
fˆ(x)−f(x)≈(Δx)2limΔx→0f(x)+f(x+Δx)2−f(x+Δx2)(Δx)2=(Δx)2f′′(x)8,
f′′(x)f(x)x
atanˆ(2z)−atan(2z)≈(Δz)22z(1−4z)ln(2)28(4z+1)2,log2ˆ(a)−log2(a)≈−(Δa)28a2ln(2).
Ponieważ funkcje są wklęsłe, a próbki pasują do funkcji, błąd jest zawsze w jednym kierunku. Lokalny maksymalny błąd bezwzględny można by zmniejszyć o połowę, gdyby znak błędu pojawiał się naprzemiennie tam iz powrotem raz na każdy interwał próbkowania. Dzięki interpolacji liniowej można uzyskać prawie optymalne wyniki, wstępnie filtrując każdą tabelę poprzez:
y[k]=⎧⎩⎨⎪⎪b2x[k−2]c1x[k−1]+b1x[k−1]b0x[k]+c0x[k]+b0x[k]+b1x[k+1]+c1x[k+1]+b2x[k+2]if k=0,if 0<k<N,if k=N,
xy0≤k≤Nc0=98,c1=−116,b0=1516,b1=18,b2=−116c0,c1N
(Δx)NlimΔx→0(c1f(x−Δx)+c0f(x)+c1f(x+Δx))(1−a)+(c1f(x)+c0f(x+Δx)+c1f(x+2Δx))a−f(x+aΔx)(Δx)N=⎧⎩⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪(c0+2c1−1)f(x)01+a−a2−c02(Δx)2f′′(x)if N=0,∣∣∣c1=1−c02if N=1,if N=2,∣∣∣c0=98
0≤a<1f(x)f(x)=exb0,b1,b2
(Δx)NlimΔx→0(b0f(x)+b1f(x+Δx)+b2f(x+2Δx))(1−a)+(c1f(x)+c0f(x+Δx)+c1f(x+2Δx))a−f(x+aΔx)(Δx)N=⎧⎩⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪(b0+b1+b2−1+a(1−b0−b1−b2))f(x)(a−1)(2b0+b1−2)Δxf′(x)(−12a2+(2316−b0)a+b0−1)(Δx)2f′′(x)if N=0,∣∣∣b2=1−b0−b1if N=1,∣∣∣b1=2−2b0if N=2,∣∣∣b0=1516
0≤a<1
log2(a)
W tym artykule prawdopodobnie przedstawiono bardzo podobny algorytm: R. Gutierrez, V. Torres i J. Valls, „ Implementacja FPGA atan (Y / X) w oparciu o transformację logarytmiczną i techniki oparte na LUT ”, Journal of Systems Architecture , tom . 56, 2010. Streszczenie mówi, że ich implementacja przewyższa poprzednie algorytmy oparte na CORDIC pod względem prędkości, a algorytmy oparte na LUT pod względem wielkości powierzchni.