Autor Tema: [RESUELTO] XC8 - Mejorar precisión de funciones trigonométricas  (Leído 16623 veces)

0 Usuarios y 1 Visitante están viendo este tema.

Desconectado migsantiago

  • Colaborador
  • DsPIC33
  • *****
  • Mensajes: 8257
    • Sitio de MigSantiago
Hola

Estoy usando un PIC18F26K20 y el XC8 v1.30. Incluí la librería math.h y puse el linker para trabajar con float/double de 32 bits. Estoy llamando funciones sin/cos/acos en mi código y los resultados son aproximadamente correctos. Aún asú, la precisión no me es suficiente, es menor a la de una PC. Uso lo anterior en esta funcón de cálculo de distancias GPS:

Código: [Seleccionar]
float SIM_GPS_Distance(float lat1, float lon1, float lat2, float lon2)
{
   float theta;
   float dist;

   theta = lon1 - lon2;

   lat1 = SIM_Degree_to_Radian(lat1);
   lat2 = SIM_Degree_to_Radian(lat2);
   theta = SIM_Degree_to_Radian(theta);

   dist = (sin(lat1) * sin(lat2)) + (cos(lat1) * cos(lat2) * cos(theta));
   dist = acos(dist);

// dist = dist * 180.0 / SIM_PI_VALUE;
// dist = dist * 60.0 * 1.1515;
// /* Convert to km */
// dist = dist * 1.609344;

   dist *= 6370.693486F;

   return (dist);
}

Cuando depuro con MPLAB X, obtengo estos datos:

Código: [Seleccionar]
dist = (sin(lat1) * sin(lat2)) + (cos(lat1) * cos(lat2) * cos(theta)); -----> 0.99999994
dist = acos(dist); -----> 3.452301E-4
dist *= 6370.693486F; -----> 2.1993551

El último valor me indica 2.2km, pero debería ser 0.

Ahora, el mismo código en una PC regresa:

Código: [Seleccionar]
dist = (sin(lat1) * sin(lat2)) + (cos(lat1) * cos(lat2) * cos(theta)); -----> 1
dist = acos(dist); -----> 0
dist *= 6370.693486F; -----> 0

0km, que son correctos.

Sé que las funciones trigonométricas son aproximaciones con tablas y otros algoritmos que hacen su cálculo eficiente, pero no me molestaría que ocuparan mucha ROM ni que fueran lentas con tal de obtener más precisión.

¿XC8 tiene forma de lograr lo anterior? ¿Podría usar alguna librería externa para mejorarlas?

¿Conviene usar un dsPIC on punto fijo mejor?

Gracias  :mrgreen:
« Última modificación: 13 de Mayo de 2014, 21:53:41 por migsantiago »

Desconectado elgarbe

  • Moderadores
  • PIC24H
  • *****
  • Mensajes: 2178
Re: XC8 - Mejorar precisión de funciones trigonométricas
« Respuesta #1 en: 11 de Mayo de 2014, 20:45:34 »
Tiro una idea al aire, no sé si lo que voy a decir es correcto o no, pero aqui va:

yo provaría calcular la distancia entre dos latitutdes/longitudes distintas y compararía con lo obtenido en la PC. Puede que el problema sea en la funcion ACOS cuando el argumento se hacerca a 1 y cuando no está tan cerca quizá la funcion funcione bien. En cuyo caso quizá podrías poner un if despues del ACOS, si es menor a un valor determinado le asignas 0 directamente. lo que quiero decir es que no creo que sea problema de presicion de las funciones en todo el rango, sino que las funciones no deben funcionar bien cuando el argumento se hacerca a 1 o 0...

saludos!
-
Leonardo Garberoglio

Desconectado migsantiago

  • Colaborador
  • DsPIC33
  • *****
  • Mensajes: 8257
    • Sitio de MigSantiago
Re: XC8 - Mejorar precisión de funciones trigonométricas
« Respuesta #2 en: 11 de Mayo de 2014, 21:03:43 »
Hola Elgarbe

También lo había pensado, el descartar valores muy pequeños... el detalles es que esos valores tan pequeños podrían equivaler a unos 20m por ejemplo. El objeto que estamos rastreando podría perderse en una zona así y pues no me gustaría redondear tan bruscamente el dato.

 :police:

Desconectado elgarbe

  • Moderadores
  • PIC24H
  • *****
  • Mensajes: 2178
Re: XC8 - Mejorar precisión de funciones trigonométricas
« Respuesta #3 en: 11 de Mayo de 2014, 21:25:49 »
Entiendo... igual chequearía que realmente sea un problema de presicion del float...

otra cosa, los modulos GPS puede que no te den resolucion menor a 10-5 mts. tendrías que tener muchos stélites linkeados para que tengas buena resolucion...

Saludos!
-
Leonardo Garberoglio

Desconectado MGLSOFT

  • Moderadores
  • DsPIC33
  • *****
  • Mensajes: 7918
Re: XC8 - Mejorar precisión de funciones trigonométricas
« Respuesta #4 en: 11 de Mayo de 2014, 21:55:27 »
Es extraño, uso ese mismo PIC y calculo varias cosas con el (uso CCS) y nunca me dio problemas, es mas lo comparamos en los cálculos durante el desarrollo, haciendo los mismos cálculos en PC y la única diferencia la obtuvimos en el sexto decimal.
No habrá un bug en la librería del XC8 ?
Todos los dias aprendo algo nuevo, el ultimo día de mi vida aprenderé a morir....
Mi Abuelo.

Desconectado willynovi

  • Colaborador
  • PIC24F
  • *****
  • Mensajes: 546
Re: XC8 - Mejorar precisión de funciones trigonométricas
« Respuesta #5 en: 11 de Mayo de 2014, 22:21:46 »
Los FLoat siempre van a tener un error, y como dice mglsoft, en el sexto decimal ha comprobado eso, pero si tu multiplacas en el orden de los 1000, como lo haces en la última multiplicación, ese error lo tienes ahora en el segundo decimal.

Otro punto a tener en cuenta es que si multiplicas un float por otro float ese error se hace casa vez mayor.

Creo que no puedes comparar los datos con los que te da la PC, digo, no pretendas tener en el microcontrolador el mismo error que en la PC, los micros usados para PC hoy tienen funciones trigonometricas propias de él, y tu usas una libreria que implementa por software los calculos.

Estuve viendo que hay un algoritmo para estos cálculos, y que se usa justamente en microcontroladores que carecen de instrucciones especificas para calculos de trigonometrias, se llama CORDIC.

No se como será el algoritmo de la libreria math, quizás sea ese mismo.
Intento enseñarte a pescar, si solo quieres pescados, espera que un pescador te regale los suyos.

Desconectado willynovi

  • Colaborador
  • PIC24F
  • *****
  • Mensajes: 546
Re: XC8 - Mejorar precisión de funciones trigonométricas
« Respuesta #6 en: 11 de Mayo de 2014, 22:52:56 »
La AN660 habla sobre los errores que podes encontrarte cuando usas los float, Floating Point Math Functions

Intento enseñarte a pescar, si solo quieres pescados, espera que un pescador te regale los suyos.


Desconectado Picuino

  • Moderadores
  • DsPIC33
  • *****
  • Mensajes: 5892
    • Picuino
Re: XC8 - Mejorar precisión de funciones trigonométricas
« Respuesta #8 en: 12 de Mayo de 2014, 11:06:16 »
El problema está en los bits que utilizas.
Con 32 bits, tienes 8 de exponente y 24 de mantisa.

El error que vas a tener puedes calcularlo como:

1 - 2^-24 = 0,99999994039

Que es lo que has conseguido.

Para mejorar los resultados, puedes:
   1.- Pasarte a un PIC24 que creo que tiene soporte para coma flotante de 64 bits
   2.- Intentar calcular de otra forma el número, utilizando las reglas trigonométricas:
        http://people.math.sfu.ca/~cbm/aands/page_72.htm
       Teniendo en cuenta que la resolución se pierde en las sumas y restas
       (no se pierde resolución en las multiplicaciones y divisiones) y que la función
       coseno tiene una resta (pierde resolución cerca del valor cero, es mejor convertirla en ese caso a una función seno).

   3.- Buscar rutinas para coma flotante de 64 bits para PIC18 (creo que Suky había publicado unas, ahora no sé donde)
   4.- Hacer tu mismo las rutinas para 64 bits (poco recomendable)

La mejor resolución que vas a obtener para flotantes de 32 bits es de 1 en 16 millones.

Saludos.
« Última modificación: 12 de Mayo de 2014, 11:08:21 por Picuino »

Desconectado Picuino

  • Moderadores
  • DsPIC33
  • *****
  • Mensajes: 5892
    • Picuino
Re: XC8 - Mejorar precisión de funciones trigonométricas
« Respuesta #9 en: 12 de Mayo de 2014, 11:16:35 »
Algo debe estar mal en el código, porque el error máximo con precisión de coma flotante de 32bits, debería ser de unos 3 metros.

¿Qué valores de latitud y longitud estás poniendo?

Saludos.

Desconectado Picuino

  • Moderadores
  • DsPIC33
  • *****
  • Mensajes: 5892
    • Picuino
Re: XC8 - Mejorar precisión de funciones trigonométricas
« Respuesta #10 en: 12 de Mayo de 2014, 11:37:15 »
Fórmula:

http://es.wikipedia.org/wiki/F%C3%B3rmula_del_Haversine

d_lat = lat1 - lat2

d_lon = lon1 - lon2

cos(dist/R) = 1 - cos(d_lat) + cos(lat1) * cos(lat2) * (1-cos(d_lon))


dist = R * acos( 1 - cos(d_lat) + cos(lat1) * cos(lat2) * (1 - cos(d_lon)) )

R = radio de la tierra:

   Ecuatorial     6378100 m
   Polar            6356800 m 
   Medio           6371000 m

Saludos.

Desconectado migsantiago

  • Colaborador
  • DsPIC33
  • *****
  • Mensajes: 8257
    • Sitio de MigSantiago
Re: XC8 - Mejorar precisión de funciones trigonométricas
« Respuesta #11 en: 12 de Mayo de 2014, 22:26:56 »
Vaya que me dieron que leer. Gracias!

Tendré qué evaluar por qué camino irme. Por acá también me recomendaron el Haversine.

Empezaré a ver qué onda.

Saludos!

Desconectado migsantiago

  • Colaborador
  • DsPIC33
  • *****
  • Mensajes: 8257
    • Sitio de MigSantiago
Re: XC8 - Mejorar precisión de funciones trigonométricas
« Respuesta #12 en: 12 de Mayo de 2014, 23:10:31 »
otra cosa, los modulos GPS puede que no te den resolucion menor a 10-5 mts. tendrías que tener muchos stélites linkeados para que tengas buena resolucion...

Sí, ése detalle es muy importante. Por el momento me estoy enfocando en la precisión de las trigonométricas.

Es extraño, uso ese mismo PIC y calculo varias cosas con el (uso CCS) y nunca me dio problemas, es mas lo comparamos en los cálculos durante el desarrollo, haciendo los mismos cálculos en PC y la única diferencia la obtuvimos en el sexto decimal.
No habrá un bug en la librería del XC8 ?

De hecho las aproximaciones trigonométricas me están dando un error en el octavo decimal, lo cual es bueno, muy bueno, pero no suficientemente bueno para mi cálculo.

El PIC no sabe nada de trigonometría, así que todo viene dado por el compilador y cómo aproxima los cálculos trigonométricos para ahorrar ROM. Sería interesante correr el mismo ejemplo pero con CCS para ver qué tanta precisión ellos aplican en sus funciones trigonométricas.

Creo que no puedes comparar los datos con los que te da la PC, digo, no pretendas tener en el microcontrolador el mismo error que en la PC, los micros usados para PC hoy tienen funciones trigonometricas propias de él, y tu usas una libreria que implementa por software los calculos.

Estuve viendo que hay un algoritmo para estos cálculos, y que se usa justamente en microcontroladores que carecen de instrucciones especificas para calculos de trigonometrias, se llama CORDIC.

No se como será el algoritmo de la libreria math, quizás sea ese mismo.


De acuerdo Willynovi. Tomé como referencia los de la PC porque no sabía si el algoritmo como tal en el PIC estaba funcionando. Sí funciona, pero su margen de error viene dado por las aproximaciones trigonométricas usadas en XC8, no creo que sea por la precisión de 32 bits.

Voy a poner en segundo lugar al CORDIC en mis pruebas. Por ahora el Haversine tiene 2 votos a favor jeje.

En la que hay funciones trigonométricas:
Circular trigonometric functions

Gracias Manolo... el implementar las funciones trigonométricas desde cero será mi tercera opción. Tengo un aproximado de 25kB de ROM para quemar... a ver si caben.

El problema está en los bits que utilizas.
Con 32 bits, tienes 8 de exponente y 24 de mantisa.

El error que vas a tener puedes calcularlo como:

1 - 2^-24 = 0,99999994039

Se me hace extraño que el tamaño de 32 bits en los cálculos sea el que determine la calidad de los resultados. Afecta... sí, pero creo que va más por el lado de cómo XC8 implementa los cálculos trigonométricos. Una PC tiene su coprocesador que es súper preciso, pero XC8 implementa tablas de aproximaciones.

Algo debe estar mal en el código, porque el error máximo con precisión de coma flotante de 32bits, debería ser de unos 3 metros.

¿Qué valores de latitud y longitud estás poniendo?

Saludos.

El código es idéntico en el PIC y en la PC. En la PC también uso 32 bits. La diferencia es que la PC usa el coprocesador para trigonometría. Esto es lo que uso:

Código: [Seleccionar]
int main(int argc, char *argv[])
{
  float lat1 = 20.574570;
  float lon1 = -100.382116;
  float lat2 = 20.574570;
  float lon2 = -100.382116;
  
  float result;
  
  result = SIM_GPS_Distance(lat1, lon1, lat2, lon2);
  
  printf("%f\n", result);
  
  system("PAUSE");
  return 0;
}

Ambos datos son idénticos. No he hecho pruebas con metros de distancia, pero el cero de distancia es el que me preocupó.

Voy a hacer pruebas amigos. Gracias de nuevo por sus tips. Qué padre es el foro Todopic  :mrgreen:

Desconectado migsantiago

  • Colaborador
  • DsPIC33
  • *****
  • Mensajes: 8257
    • Sitio de MigSantiago
Re: XC8 - Mejorar precisión de funciones trigonométricas
« Respuesta #13 en: 13 de Mayo de 2014, 21:51:50 »
Hola

Acabo de probar el algoritmo de Haversine y me gustó mucho su rendimiento.

Les dejo el código. Hice pruebas con mapas de Google y calculadoras de distancia en línea. Hice los cálculos en la PC y en el PIC con 32 bits y ambos fueron muy aproximados entre ellos. Con double mejorarían, pero no se puede en el PIC.

El resultado sale en metros, pero lo pueden pasar a km modificando SIM_EARTH_RADIUS.

Saludos y muchas gracias  :mrgreen:

Código: [Seleccionar]
/**
 * Convert degrees to radians
 * @param x - degrees to convert
 * @return radians
 */
static inline float SIM_Degree_to_Radian(float x)
{
   return (x * 0.017453292F);
}

#define SIM_EARTH_RADIUS                     (6372797.560856F)

float SIM_GPS_Distance(float lat1, float lon1, float lat2, float lon2)
{
   float delta_lat;
   float delta_lon;
   float a;
   float c;
   float dist;

   lat1 = SIM_Degree_to_Radian(lat1);
   lat2 = SIM_Degree_to_Radian(lat2);

   delta_lat = SIM_Degree_to_Radian(lat2 - lat1);
   delta_lon = SIM_Degree_to_Radian(lon2 - lon1);

   a = sin(delta_lat / 2.0F) * sin(delta_lat / 2.0F) +
         cos(lat1) * cos(lat2) *
         sin(delta_lon / 2.0F) * sin(delta_lon / 2.0F);
   c = 2.0F * atan2(sqrt(a), sqrt(1.0F - a));
   dist = SIM_EARTH_RADIUS * c;

   return dist;
}

Desconectado Picuino

  • Moderadores
  • DsPIC33
  • *****
  • Mensajes: 5892
    • Picuino
Re: XC8 - Mejorar precisión de funciones trigonométricas
« Respuesta #14 en: 14 de Mayo de 2014, 05:19:47 »
Con double mejorarían, pero no se puede en el PIC.

Se puede hacer que el pic calcule con más precisión a costa de más esfuerzo programando tus librerías para manejar punto flotante en el pic18.

Librerías de ejemplo:
gnu glibc:      http://www.gnu.org/software/libc/libc.html                  http://ftp.gnu.org/gnu/glibc
softfloat:       http://www.jhauser.us/arithmetic/SoftFloat.html          http://svnweb.freebsd.org/base/head/lib/libc/softfloat
Fdlibm:         http://www.netlib.org/fdlibm
Crlibm:         http://lipforge.ens-lyon.fr/www/crlibm/

En realidad las funciones trascendentes (seno, coseno, tangente, logaritmos, etc) no se implementan con tablas, sino con polinomios.
Hay unos polinomios especiales, denominados polinomios de Chebyshev, que tienen un error absoluto mínimo en la aproximación de funciones.
Esos son los que se utilizan para simular las diferentes funciones

Un ejemplo de implementación de la función seno con la librería fdlibm:
Código: [Seleccionar]
* ====================================================
 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
 *
 * Developed at SunSoft, a Sun Microsystems, Inc. business.
 * Permission to use, copy, modify, and distribute this
 * software is freely granted, provided that this notice
 * is preserved.
 * ====================================================
 */

/* __kernel_sin( x, y, iy)
 * kernel sin function on [-pi/4, pi/4], pi/4 ~ 0.7854
 * Input x is assumed to be bounded by ~pi/4 in magnitude.
 * Input y is the tail of x.
 * Input iy indicates whether y is 0. (if iy=0, y assume to be 0).
 *
 * Algorithm
 * 1. Since sin(-x) = -sin(x), we need only to consider positive x.
 * 2. if x < 2^-27 (hx<0x3e400000 0), return x with inexact if x!=0.
 * 3. sin(x) is approximated by a polynomial of degree 13 on
 *    [0,pi/4]
 *            3            13
 *     sin(x) ~ x + S1*x + ... + S6*x
 *    where
 *
 * |sin(x)         2     4     6     8     10     12  |     -58
 * |----- - (1+S1*x +S2*x +S3*x +S4*x +S5*x  +S6*x   )| <= 2
 * |  x            |
 *
 * 4. sin(x+y) = sin(x) + sin'(x')*y
 *     ~ sin(x) + (1-x*x/2)*y
 *    For better accuracy, let
 *      3      2      2      2      2
 * r = x *(S2+x *(S3+x *(S4+x *(S5+x *S6))))
 *    then                   3    2
 * sin(x) = x + (S1*x + (x *(r-y/2)+y))
 */

#include "fdlibm.h"

#ifdef __STDC__
static const double
#else
static double
#endif
half =  5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */
S1  = -1.66666666666666324348e-01, /* 0xBFC55555, 0x55555549 */
S2  =  8.33333333332248946124e-03, /* 0x3F811111, 0x1110F8A6 */
S3  = -1.98412698298579493134e-04, /* 0xBF2A01A0, 0x19C161D5 */
S4  =  2.75573137070700676789e-06, /* 0x3EC71DE3, 0x57B1FE7D */
S5  = -2.50507602534068634195e-08, /* 0xBE5AE5E6, 0x8A2B9CEB */
S6  =  1.58969099521155010221e-10; /* 0x3DE5D93A, 0x5ACFD57C */

#ifdef __STDC__
double __kernel_sin(double x, double y, int iy)
#else
double __kernel_sin(x, y, iy)
double x,y; int iy; /* iy=0 if y is zero */
#endif
{
double z,r,v;
int ix;
ix = __HI(x)&0x7fffffff; /* high word of x */
if(ix<0x3e400000) /* |x| < 2**-27 */
   {if((int)x==0) return x;} /* generate inexact */
z =  x*x;
v =  z*x;
r =  S2+z*(S3+z*(S4+z*(S5+z*S6)));
if(iy==0) return x+v*(S1+z*r);
else      return x-((z*(half*y-v*r)-y)-v*S1);
}


Saludos.