méthodes de calcul du point fixe atan2 sur FPGA

12

J'ai besoin de calculer atan2(x,y)sur un FPGA avec un flux continu de données d'entrée / sortie. J'ai réussi à l'implémenter à l'aide de noyaux CORDIC déroulés et pipelinés, mais pour obtenir la précision dont j'avais besoin, j'ai dû effectuer 32 itérations. Cela a conduit à consacrer une assez grande quantité de LUT à cette seule tâche. J'ai essayé de changer le flux pour utiliser des noyaux CORDIC partiellement déroulés, mais j'avais ensuite besoin d'une fréquence d'horloge multipliée pour exécuter des boucles répétées tout en maintenant un flux d'entrée / sortie continu. Avec cela, je ne pouvais pas respecter le calendrier.

Alors maintenant, je cherche des moyens de calcul alternatifs atan2(x,y).

J'ai pensé à utiliser des tables de recherche de bloc-RAM avec interpolation, mais comme il y a 2 variables, j'aurais besoin de 2 dimensions de tables de recherche, ce qui nécessite beaucoup de ressources en termes d'utilisation de bloc-RAM.

J'ai alors pensé à utiliser le fait atan2(x,y)lié à l' atan(x/y)ajustement quadrant. Le problème avec cela est qu'il a x/ybesoin d'une vraie division car ce yn'est pas une constante, et les divisions sur les FPGA sont très gourmandes en ressources.

Existe-t-il d'autres façons nouvelles d'implémenter atan2(x,y)sur un FPGA qui entraîneraient une utilisation moindre des LUT, mais fourniraient toujours une bonne précision?

user2913869
la source
2
Quelle est votre fréquence d'horloge de traitement et votre taux de données d'entrée?
Jim Clay
Quelle est la précision dont vous avez besoin? Je suppose également que vous utilisez le calcul à virgule fixe. Quelle profondeur de bits utilisez-vous? Une approximation polynomiale (ou LUT) avec ajustement en quadrant est une méthode courante à mettre en œuvre atan2. Je ne sais pas si vous pouvez vous en tirer sans division, cependant.
Jason R
L'horloge d'entrée est de 150 MHz, le débit de données d'entrée est de 150 MSamps / sec. En gros, j'obtiens une nouvelle entrée à chaque cycle d'horloge. La latence, c'est bien, mais je dois aussi produire une sortie à 150 MSamps / sec.
user2913869
Mes simulations montrent que je peux vivre avec environ 1 * 10 ^ -9. Je ne suis pas sûr du nombre absolu de bits à virgule fixe, mais j'ai simulé avec un format à virgule fixe
Q10.32
Cet article explique une implémentation à virgule fixe de atan2. Mais vous aurez toujours besoin d'une division.
Matt L.

Réponses:

20

Vous pouvez utiliser des logarithmes pour vous débarrasser de la division. Pour (x,y) dans le premier quadrant:

z=log2(y)log2(x)atan2(y,x)=atan(y/x)=atan(2z)

atan (2 ^ z)

Figure 1. Diagramme d' atan(2z)

Vous auriez besoin d'approximativement un atan(2z) dans la plage 30<z<30 pour obtenir la précision requise de 1E-9. Vous pouvez profiter de la symétrie atan(2z)=π2atan(2z)ou bien assurez-vous que(x,y)est dans un octant connu. Pour approximer lelog2(a):

b=floor(log2(a))c=a2blog2(a)=b+log2(c)

b peut être calculé en trouvant l'emplacement du bit non nul le plus significatif. c peut être calculé par un décalage de bits. Vous devez approximerlog2(c) dans la plage1c<2 .

log2 (c)

Figure 2. Diagramme du log2(c)

214+1=16385log2(c)30×212+1=122881atan(2z)0<z<30z

Erreur d'approximation atan (2 ^ z)

atan(2z)zz0z<1floor(log2(z))=0

atan(2z)0z<1floor(log2(z))z1atan(2z)z0z<32

Pour référence ultérieure, voici le script Python maladroit que j'ai utilisé pour calculer les erreurs d'approximation:

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

f(x)f^(x)f(x)Δx

f^(x)f(x)(Δx)2limΔx0f(x)+f(x+Δx)2f(x+Δx2)(Δx)2=(Δx)2f(x)8,

où est la dérivée seconde de et est à un maximum local de l'erreur absolue. Avec ce qui précède, nous obtenons les approximations:f(x)f(x)x

atan^(2z)atan(2z)(Δz)22z(14z)ln(2)28(4z+1)2,log2^(a)log2(a)(Δa)28a2ln(2).

Étant donné que les fonctions sont concaves et que les échantillons correspondent à la fonction, l'erreur se produit toujours dans une direction. L'erreur absolue maximale locale pourrait être réduite de moitié si le signe de l'erreur était alterné d'avant en arrière une fois à chaque intervalle d'échantillonnage. Avec une interpolation linéaire, des résultats proches de l'optimisation peuvent être obtenus en préfiltrant chaque table en:

y[k]={b0x[k]+b1x[k+1]+b2x[k+2]if k=0,c1x[k1]+c0x[k]+c1x[k+1]if 0<k<N,b2x[k2]+b1x[k1]+b0x[k]if k=N,

où et sont l'original et la table filtrée s'étendant sur et les poids sont . Le conditionnement final (première et dernière ligne de l'équation ci-dessus) réduit l'erreur aux extrémités du tableau par rapport à l'utilisation d'échantillons de la fonction en dehors du tableau, car le premier et le dernier échantillon n'ont pas besoin d'être ajustés pour réduire l'erreur d'interpolation entre elle et un échantillon juste à l'extérieur de la table. Les sous-tables avec différents intervalles d'échantillonnage doivent être préfiltrées séparément. Les valeurs des poids ont été trouvées en minimisant séquentiellement pour augmenter l'exposantxy0kNc0=98,c1=116,b0=1516,b1=18,b2=116c0,c1N la valeur absolue maximale de l'erreur approximative:

(Δx)NlimΔx0(c1f(xΔx)+c0f(x)+c1f(x+Δx))(1a)+(c1f(x)+c0f(x+Δx)+c1f(x+2Δx))af(x+aΔx)(Δx)N={(c0+2c11)f(x)if N=0,|c1=1c020if N=1,1+aa2c02(Δx)2f(x)if N=2,|c0=98

pour les positions d'interpolation entre échantillons , avec une fonction concave ou convexe (par exemple ). Une fois ces poids résolus, les valeurs des poids de conditionnement d'extrémité ont été trouvées en minimisant de manière similaire la valeur absolue maximale de:0a<1f(x)f(x)=exb0,b1,b2

(Δx)NlimΔx0(b0f(x)+b1f(x+Δx)+b2f(x+2Δx))(1a)+(c1f(x)+c0f(x+Δx)+c1f(x+2Δx))af(x+aΔx)(Δx)N={(b0+b1+b21+a(1b0b1b2))f(x)if N=0,|b2=1b0b1(a1)(2b0+b12)Δxf(x)if N=1,|b1=22b0(12a2+(2316b0)a+b01)(Δx)2f(x)if N=2,|b0=1516

pour . L'utilisation du préfiltre réduit de moitié l'erreur d'approximation et est plus facile à faire que l'optimisation complète des tableaux.0a<1

Erreur d'approximation avec et sans préfiltre et conditionnement d'extrémité

Figure 4. Erreur d'approximation de partir de 11 échantillons, avec et sans préfiltre et avec et sans conditionnement d'extrémité. Sans conditionnement final, le préfiltre a accès aux valeurs de la fonction juste à l'extérieur du tableau.log2(a)

Cet article présente probablement un algorithme très similaire: R. Gutierrez, V. Torres et J. Valls, « FPGA-implementation of atan (Y / X) based on logarithmic transformation and LUT-based techniques », Journal of Systems Architecture , vol . 56, 2010. Le résumé indique que leur mise en œuvre bat les précédents algorithmes basés sur CORDIC en vitesse et les algorithmes basés sur LUT en taille d'encombrement.

Olli Niemitalo
la source
3
Matthew Gambrell et moi avons procédé à une ingénierie inverse de la puce audio Yamaha YM3812 1985 (par microscopie) et y avons trouvé des tables similaires de mémoire morte (ROM) log / exp. Yamaha avait utilisé une astuce supplémentaire pour remplacer chaque entrée sur deux dans chaque tableau par une différence avec l'entrée précédente. Pour des fonctions fluides, la différence prend moins de bits et de surface de puce à représenter que la fonction. Ils avaient déjà un additionneur sur la puce qu'ils ont pu utiliser pour ajouter la différence à l'entrée précédente.
Olli Niemitalo
3
Merci beaucoup! J'adore ce genre d'exploits de propriétés mathématiques. Je vais certainement développer quelques sims MATLAB de cela, et si tout va bien, passer au HDL. Je ferai rapport de mes économies LUT lorsque tout sera fait.
user2913869
J'ai utilisé votre description comme guide et je suis heureux de rester que j'ai réduit les LUT de près de 60%. J'avais un besoin de réduire les BRAM, alors j'ai compris que je pouvais obtenir une erreur maximale cohérente sur ma table ATAN en faisant un échantillonnage non uniforme: j'avais plusieurs BRAM LUT (tous le même nombre de bits d'adresse), le plus proche de zéro, plus l'échantillonnage est rapide. J'ai choisi mes plages de table pour être des puissances de 2 afin que je puisse facilement détecter dans quelle plage j'étais également et faire l'indexation automatique des tables via la manipulation de bits. J'ai également appliqué une symétrie, donc je n'ai enregistré que la moitié de la forme d'onde.
user2913869
De plus, j'ai peut-être manqué certaines de vos modifications, mais j'ai réussi à implémenter 2 ^ z en le divisant en 2 ^ {if} = 2 ^ i * 2 ^ {0.f}, où i est la partie entière et f est la partie fractionnaire. 2 ^ i est simple, juste une manipulation de bits, et 2 ^ {0.f} avait une plage limitée, donc il se prêtait bien à LUT avec interpolation. J'ai également traité le cas négatif: 2 ^ {- if} = 2 ^ {- i} * 1 / (2 ^ {0.f}. Donc, un autre tableau pour 1/2 / ^ ^ {0.f}. Ma prochaine étape pourrait être d'appliquer la puissance de 2 échantillonnages de portée / non uniforme sur les LUT log2 (y), car cela semble être la forme d'onde candidate parfaite pour ce genre de chose. A
bientôt
1
Lol yup j'ai totalement raté cette étape. Je vais essayer ça maintenant. Cela devrait me faire économiser encore plus de LUT et encore plus de BRAM
user2913869