métodos de cálculo de punto fijo atan2 en FPGA

12

Necesito computar atan2(x,y)en un FPGA con un flujo continuo de datos de entrada / salida. Logré implementarlo usando núcleos CORDIC enrollados y canalizados, pero para obtener la precisión que necesitaba, tuve que realizar 32 iteraciones. Esto llevó a que se dedicara una gran cantidad de LUT a esta tarea. Intenté cambiar el flujo para usar núcleos CORDIC parcialmente desenrollados, pero luego necesitaba una frecuencia de reloj multiplicada para ejecutar bucles repetidos mientras mantenía un flujo continuo de entrada / salida. Con esto, no pude cumplir con el tiempo.

Así que ahora estoy buscando formas alternativas de computación atan2(x,y).

Pensé en usar tablas de búsqueda de RAM en bloque con interpolación, pero dado que hay 2 variables, necesitaría 2 dimensiones de las tablas de búsqueda, y esto requiere muchos recursos en términos de uso de RAM en bloque.

Luego pensé en usar el hecho que atan2(x,y)está relacionado atan(x/y)con el ajuste del cuadrante. El problema con esto es que x/ynecesita una verdadera división ya yque no es una constante, y las divisiones en FPGAs requieren muchos recursos.

¿Hay más formas novedosas de implementar atan2(x,y)en un FPGA que resulten en un menor uso de LUT, pero que aún así proporcionen una buena precisión?

usuario2913869
fuente
2
¿Cuál es su velocidad de reloj de procesamiento y su velocidad de datos de entrada?
Jim Clay
¿Cuál es la precisión que necesita? Supongo que también estás usando computación de punto fijo. ¿Qué profundidad de bits estás usando? Una aproximación polinómica (o LUT) con ajuste de cuadrante es un método común para implementar atan2. Sin embargo, no estoy seguro de poder sobrevivir sin una división.
Jason R
El reloj de entrada es de 150MHz, la velocidad de datos de entrada es de 150 MSamps / seg. Básicamente obtengo una nueva entrada en cada ciclo de reloj. Tener latencia está bien, pero también debo producir una salida a 150 MSamps / seg.
user2913869
Mis simulaciones muestran que puedo vivir con aproximadamente 1 * 10 ^ -9. No estoy seguro de los bits de punto fijo mínimo absoluto, pero he estado simulando con un formato de punto fijo
Q10.32
Este artículo explica una implementación de punto fijo de atan2. Sin embargo, aún necesitarás una división.
Matt L.

Respuestas:

20

Puedes usar logaritmos para deshacerte de la división. Para (x,y) en el primer cuadrante:

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

atan (2 ^ z)

Figura 1. Parcela de atan(2z)

atan(2z) aproximar atan ( 2 z ) en el rango 30<z<30 para obtener la precisión requerida de 1E-9. Puede aprovechar la simetría atan(2z)=π2atan(2z)o alternativamente asegúrese de que(x,y)esté en un octante conocido. Para aproximar ellog2(a):

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

b se puede calcular encontrando la ubicación del bit distinto de cero más significativo. c puede calcularse mediante un cambio de bit. Debería aproximar ellog2(c) en el rango1c<2 .

log2 (c)

Figura 2. Gráfico del log2(c)

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

Error de aproximación atan (2 ^ z)

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

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

Para referencia posterior, aquí está el script torpe de Python que usé para calcular los errores de aproximación:

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,

donde es la segunda derivada de y está en un máximo local del error absoluto. Con lo anterior obtenemos las aproximaciones: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).

Debido a que las funciones son cóncavas y las muestras coinciden con la función, el error siempre es en una dirección. El error absoluto máximo local podría reducirse a la mitad si el signo del error se alternara una y otra vez cada intervalo de muestreo. Con la interpolación lineal, se pueden lograr resultados casi óptimos filtrando previamente cada tabla de la siguiente manera:

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,

donde e son la tabla original y la filtrada que abarcan y los pesos son . El acondicionamiento final (primera y última fila en la ecuación anterior) reduce el error en los extremos de la tabla en comparación con el uso de muestras de la función fuera de la tabla, porque la primera y la última muestra no necesitan ajustarse para reducir el error de la interpolación entre él y una muestra justo afuera de la mesa. Las subtablas con diferentes intervalos de muestreo deben prefiltrarse por separado. Los valores de los pesos se encontraron minimizando secuencialmente para aumentar el exponentexy0kNc0=98,c1=116,b0=1516,b1=18,b2=116c0,c1N El valor absoluto máximo del error aproximado:

(Δ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

para las posiciones de interpolación entre muestras , con una función cóncava o convexa (por ejemplo ). Con esos pesos resueltos, los valores de los pesos de acondicionamiento final se encontraron minimizando de manera similar el valor absoluto máximo 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

para . El uso del prefiltro acerca de la mitad del error de aproximación y es más fácil de hacer que la optimización completa de las tablas.0a<1

Error de aproximación con y sin prefiltro y acondicionamiento final

Figura 4. Error de aproximación de de 11 muestras, con y sin prefiltro y con y sin acondicionamiento final. Sin el acondicionamiento final, el prefiltro tiene acceso a los valores de la función justo fuera de la tabla.log2(a)

Este artículo probablemente presenta un algoritmo muy similar: R. Gutiérrez, V. Torres y J. Valls, " Implementación FPGA de atan (Y / X) basada en transformación logarítmica y técnicas basadas en LUT " , Journal of Systems Architecture , vol . 56, 2010. El resumen dice que su implementación supera los algoritmos anteriores basados ​​en CORDIC en velocidad y los algoritmos basados ​​en LUT en tamaño de huella.

Olli Niemitalo
fuente
3
Matthew Gambrell y yo hemos realizado una ingeniería inversa del chip de sonido Yamaha YM3812 de 1985 (por microscopía) y encontramos en él tablas de memoria de solo lectura (ROM) de registro / exp similares. Yamaha había utilizado un truco adicional para reemplazar cada segunda entrada en cada tabla por una diferencia con respecto a la entrada anterior. Para funciones suaves, la diferencia requiere menos bits y área de chip para representar que la función. Ya tenían un sumador en el chip que podían usar para agregar la diferencia a la entrada anterior.
Olli Niemitalo
3
¡Muchas gracias! Me encantan este tipo de hazañas de propiedades matemáticas. Definitivamente desarrollaré algunos simulacros de MATLAB de esto, y si todo se ve bien, pasaré a HDL. Informaré de mis ahorros de LUT cuando todo esté listo.
user2913869
Utilicé su descripción como guía y estoy feliz de quedarme porque reduje los LUT en casi un 60%. Tenía la necesidad de reducir los BRAM, así que descubrí que podía obtener un error máximo constante en mi tabla ATAN haciendo un muestreo no uniforme: tenía múltiples BRAM LUT (todos el mismo número de bits de dirección), cuanto más cerca estaba cero, cuanto más rápido sea el muestreo. Elegí mis rangos de tabla para que fueran potencias de 2 para poder detectar fácilmente en qué rango estaba también y hacer una indexación automática de tablas a través de la manipulación de bits. También apliqué simetría atan, así que solo almacené la mitad de la forma de onda.
user2913869
Además, es posible que me haya perdido algunas de sus ediciones, pero logré implementar 2 ^ z dividiéndolo en 2 ^ {if} = 2 ^ i * 2 ^ {0.f}, donde i es la parte entera y f es La porción fraccional. 2 ^ i es simple, solo manipulación de bits, y 2 ^ {0.f} tenía un rango limitado, por lo que se prestó bien a LUT con interpolación. También manejé el caso negativo: 2 ^ {- if} = 2 ^ {- i} * 1 / (2 ^ {0.f}. Entonces, una tabla más para 1/2 ^ {0.f}. Mi siguiente paso podría ser la aplicación de la potencia de 2 que van de muestreo / no uniforme en las tablas de búsqueda (y) log2, ya que parece como que sería la forma de onda candidato perfecto para ese tipo de cosas Saludos.!
user2913869
1
Jaja, sí, me perdí totalmente ese paso. Voy a intentar eso ahora. Debería ahorrarme aún más LUT e incluso más BRAM
user2913869