Methoden zur Berechnung von Fixpunkt atan2 auf FPGA

Ich muss atan2(x,y)auf einem FPGA mit einem kontinuierlichen Eingabe- / Ausgabestrom von Daten rechnen. Ich habe es geschafft, es mit nicht gerollten CORDIC-Kerneln in Pipelines zu implementieren, aber um die Genauigkeit zu erreichen, die ich benötige, musste ich 32 Iterationen durchführen. Dies führte dazu, dass eine große Anzahl von LUTs dieser einen Aufgabe gewidmet war. Ich habe versucht, den Fluss so zu ändern, dass teilweise entrollte CORDIC-Kernel verwendet werden, aber dann brauchte ich eine multiplizierte Taktfrequenz, um wiederholte Schleifen auszuführen und trotzdem einen kontinuierlichen Eingabe- / Ausgabe-Fluss aufrechtzuerhalten. Damit konnte ich das Timing nicht einhalten.

Jetzt suche ich nach alternativen Möglichkeiten für die Datenverarbeitung atan2(x,y).

Ich dachte darüber nach, Block-RAM-Nachschlagetabellen mit Interpolation zu verwenden, aber da es 2 Variablen gibt, würde ich 2 Dimensionen von Nachschlagetabellen benötigen, und dies ist im Hinblick auf die Block-RAM-Verwendung sehr ressourcenintensiv.

Ich habe dann darüber nachgedacht, die Tatsache zu nutzen, atan2(x,y)die atan(x/y)mit der Quadrantenanpassung zusammenhängt. Das Problem dabei ist, dass x/yeine echte Teilung erforderlich ist, da dies ykeine Konstante ist und die Teilung auf FPGAs sehr ressourcenintensiv ist.

Gibt es neuere Möglichkeiten zur Implementierung atan2(x,y)auf einem FPGA, die zu einer geringeren LUT-Nutzung führen, aber dennoch eine gute Genauigkeit bieten?

Sie können Logarithmen verwenden, um die Division zu entfernen. Für $(x, y)$ im ersten Quadranten:

$$z = \log_2(y)-\log_2(x)\\ \text{atan2}(y, x) = \text{atan}(y/x) = \text{atan}(2^z)$$

Abbildung 1. Diagramm von $\text{atan}(2^z)$

Sie müssten sich einem ( 2 z ) im Bereich von - 30 < z < 30 $\text{atan}(2^z)$ , um die erforderliche Genauigkeit von 1E-9 zu erhalten. Sie können die Symmetrie atan ( 2 - z ) = π ausnutzen$-30 < z < 30$$\text{atan}(2^{-z}) = \frac{\pi}{2}-\text{atan}(2^z)$oder stellen Sie alternativ sicher, dass$(x, y)$in einem bekannten Oktanten ist. $\log_2(a)$approximieren:

$$b = \text{floor}(\log_2(a))\\ c = \frac{a}{2^b}\\ \log_2(a) = b + \log_2(c)$$

$b$ kann berechnet werden, indem der Ort des höchstwertigen Nicht-Null-Bits ermittelt wird. $c$ kann durch eine Bitverschiebung berechnet werden. Sie müssten$\log_2(c)$ im Bereich$1 \le c < 2$ approximieren.

Abbildung 2. Diagramm von $\log_2(c)$

Für Ihre Genauigkeitsanforderungen, lineare Interpolation und gleichmäßige Abtastung, $2^{14} + 1 = 16385$ Abtastungen von $\log_2(c)$ und $30\times 2^{12} + 1 = 122881$ Abtastungen von$\text{atan}(2^z)$ für$0 < z < 30$ ausreichen. Der letztere Tisch ist ziemlich groß. Der Fehler durch Interpolation hängt dabei stark von$z$ :

Abbildung 3. Größter absoluter Fehler bei einer ( 2 z ) $\text{atan}(2^z)$ für verschiedene Bereiche von $z$ (horizontale Achse) für verschiedene Anzahlen von Stichproben (32 bis 8192) pro Einheitsintervall von $z$ . Der größte absolute Fehler für $0 \le z < 1$ (weggelassen) ist geringfügig kleiner als für $\text{floor}(\log_2(z)) = 0$ .

Die $\text{atan}(2^z)$ -Tabelle kann in mehrere Untertabellen aufgeteilt werden, die $0 \le z < 1$ und in verschiedene $\text{floor}(\log_2(z))$ mit $z \ge 1$ , was einfach zu berechnen ist. Die Tabellenlängen können gemäß Abb. 3 gewählt werden. Der Index innerhalb der Untertabelle kann durch einfache Manipulation von Bitfolgen berechnet werden. Für Ihre Genauigkeitsanforderungen haben die $\text{atan}(2^z)$ -Untertabellen insgesamt 29217 Samples, wenn Sie den Bereich von $z$ bis $0 \le z < 32$ zur Vereinfachung.

Zum späteren Nachschlagen hier das klobige Python-Skript, mit dem ich die Approximationsfehler berechnet habe:

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])    

Der lokale Maximalfehler von Approximieren einer Funktion $f(x)$ durch lineare Interpolation f$\hat{f}(x)$ aus Proben von$f(x)$ , mit Abtastintervall durch gleichförmige Abtastung gemacht$\Delta x$ kann durch analytisches angenähert werden:

$$\widehat{f}(x) - f(x) \approx (\Delta x)^2\lim_{\Delta x\rightarrow 0}\frac{\frac{f(x) + f(x + \Delta x)}{2} - f(x + \frac{\Delta x}{2})}{(\Delta x)^2} = \frac{(\Delta x)^2 f''(x)}{8},$$

wo $f''(x)$ die zweite Ableitung von$f(x)$ und$x$ ist ein lokales Maximum des absoluten Fehlers. Mit dem obigen erhalten wir die Annäherungen:

$$\widehat{\text{atan}}(2^z) - \text{atan}(2^z) \approx \frac{(\Delta z)^2 2^z(1 - 4^z)\ln(2)^2}{8(4^z + 1)^2},\\ \widehat{\log_2}(a) - \log_2(a) \approx \frac{-(\Delta a)^2}{8 a^2\ln(2)}.$$

Da die Funktionen konkav sind und die Abtastwerte mit der Funktion übereinstimmen, liegt der Fehler immer in einer Richtung. Der lokale maximale absolute Fehler könnte halbiert werden, wenn das Vorzeichen des Fehlers einmal in jedem Abtastintervall hin und her wechselt. Mit linearer Interpolation können nahezu optimale Ergebnisse erzielt werden, indem jede Tabelle vorgefiltert wird durch:

$$y[k] = \cases{\begin{array}{rrrrrl}&&b_0x[k]&\negthickspace\negthickspace\negthickspace+ b_1x[k+1]&\negthickspace\negthickspace\negthickspace+ b_2x[k+2]&\text{if } k = 0,\\ &c_1x[k-1]&\negthickspace\negthickspace\negthickspace+ c_0x[k]&\negthickspace\negthickspace\negthickspace+ c_1x[k+1]&&\text{if }0 < k < N,\\ b_2x[k-2]&\negthickspace\negthickspace\negthickspace+ b_1x[k-1]&\negthickspace\negthickspace\negthickspace+ b_0x[k]&&&\text{if } k = N, \end{array}}$$

Dabei sind $x$ und $y$ die ursprüngliche und die gefilterte Tabelle, die beide $0 \le k \le N$$c_0 = \frac{9}{8}, c_1 = -\frac{1}{16}, b_0 = \frac{15}{16}, b_1 = \frac{1}{8}, b_2 = -\frac{1}{16}$$c_0, c_1$$N$

$$(\Delta x)^N\lim_{\Delta x\rightarrow 0}\frac{\left(c_1f(x - \Delta x) + c_0f(x) + c_1f(x + \Delta x)\right)(1-a) + \left(c_1f(x) + c_0f(x + \Delta x) + c_1f(x + 2 \Delta x)\right)a - f(x + a\Delta x)}{(\Delta x)^N} =\left\{\begin{array}{ll}(c_0 + 2c_1 - 1)f(x) &\text{if } N = 0, \bigg| c_1 = \frac{1 - c_0}{2}\\ 0&\text{if }N = 1,\\ \frac{1+a-a^2-c_0}{2}(\Delta x)^2 f''(x)&\text{if }N=2, \bigg|c_0 = \frac{9}{8}\end{array}\right.$$

$0 \le a < 1$$f(x)$$f(x) = e^x$$b_0, b_1, b_2$

$$(\Delta x)^N\lim_{\Delta x\rightarrow 0}\frac{\left(b_0f(x) + b_1f(x + \Delta x) + b_2f(x + 2 \Delta x)\right)(1-a) + \left(c_1f(x) + c_0f(x + \Delta x) + c_1f(x + 2 \Delta x)\right)a - f(x + a\Delta x)}{(\Delta x)^N} =\left\{\begin{array}{ll}\left(b_0 + b_1 + b_2 - 1 + a(1 - b_0 - b_1 - b_2)\right)f(x) &\text{if } N = 0, \bigg| b_2 = 1 - b_0 - b_1\\ (a-1)(2b_0+b_1-2)\Delta x f'(x)&\text{if }N = 1,\bigg|b_1=2-2b_0\\ \left(-\frac{1}{2}a^2 + \left(\frac{23}{16} - b_0\right)a + b_0 - 1\right)(\Delta x)^2f''(x)&\text{if }N=2, \bigg|b_0 = \frac{15}{16}\end{array}\right.$$

for $0 \le a < 1$. Use of the prefilter about halves the approximation error and is easier to do than full optimization of the tables.

Figure 4. Approximation error of $\log_2(a)$ from 11 samples, with and without prefilter and with and without end conditioning. Without end conditioning the prefilter has access to values of the function just outside of the table.

This article probably presents a very similar algorithm: R. Gutierrez, V. Torres, and J. Valls, “FPGA-implementation of atan(Y/X) based on logarithmic transformation and LUT-based techniques,” Journal of Systems Architecture, vol. 56, 2010. The abstract says their implementation beats previous CORDIC-based algorithms in speed and LUT-based algorithms in footprint size.