El Templete de la melancolía

Hoy es un día apropiado para volver a la poesía con Li Po (o Li Bai, Li Bo, Li Tai-Pei, Li Tai Po, vaya con la transliteración)

El Templete de la Melancolía

Será el lugar más triste de la tierra
el pabellón donde nos despedimos.
¡Oh viento de primavera!
Qué bien conoces tú
la tristeza de la separación:
Dejas sin reverdecer
los sauces en derredor.

Y como no, aprovechar para acompañar del mismo "Poesía clásica china", Cátedra, 2001, edición de Guojian Chen, otra versión del

Sentado, solo, en la montaña de Jingting

Los pájaros han vuelto a sus nidos en bandadas.

Perezosa, la última nube se aleja.

Oh montaña, eres mi única compañera.

Ni a ti ni a mí el mirarnos nos cansa.

Octave y el PRN

Y ahora los detalles escabrosos. En Octave se puede generar la secuencia PRN del código C/A de los satélites GPS con la siguiente función

function [xca,g1,g2] = xcacode(s1, s2, lon)
 xca=zeros(1,lon);
 g1=zeros(1,lon);
 g2=zeros(1,lon);
 r1=ones(1,10);
 r2=ones(1,10);
 for ind = [1:lon]
  g2i=bitxor(r2(s1),r2(s2));
  g1(ind)=r1(10);
  g2(ind)=r2(10);
  xca(ind)=bitxor(g1(ind),g2i);
  new1=bitxor(r1(3),r1(10));
  new2=bitxor(r2(2),r2(3));
  new2=bitxor(new2,r2(6));
  new2=bitxor(new2,r2(8));
  new2=bitxor(new2,r2(9));
  new2=bitxor(new2,r2(10));
  r1=[new1 r1(1:9)];
  r2=[new2 r2(1:9)];
 endfor
endfunction

La secuencia de interés se obtiene en el primer elemento de retorno, xca.
r1 y r2 son los dos registros de desplazamiento de bits, con diez bits cada uno. Los parámetros de entrada s1 y s2 determinan la secuencia concreta que se devuelve, indicando lon la longitud. Para los satélites 05, 19, y 25 sus PRN se obtienen ( Interface Specification IS-GPS-200, pg. 7 y 8) así

xca05=xcacode(1,9,1023);
xca19=xcacode(3,6,1023);
xca25=xcacode(5,7,1023);
...
xca19_ = de01a1m1(xca19);
...

Para pasar de la versión "lógica" de 0 y 1 a la versión de señal -1 y +1 se usa esta función bien corta (gracias a la sugerencia de Misra&Enge):

function v1m1 = de01a1m1(v01)
 v1m1 = (v01==0)*(1) + (v01==1)*(-1);
endfunction

Las siguientes funciones desplazan (shift) circularmente (circ) a la izquierda o a la derecha una secuencia ini un número de pasos n

function  res = shiftcircizq(ini,n)
 siz = max(size(ini));
 if (abs(n)>= siz )
  res=ini;
 elseif (n>=0)
  res = ini([n+1:siz 1:n]);
 else
  n=-n;
  res = ini([siz-n+1:siz 1:siz-n]);
 endif
endfunction

function  res = shiftcircder(ini,n)
 res = shiftcircizq(ini,-n);
endfunction

Por ejemplo

> shiftcircizq([0 0 1 0 1 1 1 0],2)
ans =
1   0   1   1   1   0   0   0
> shiftcircder([0 0 1 0 1 1 1 0],2)
ans =
1   0   0   0   1   0   1   1

La función de correlación circular de dos vectores x e y se calcula con la siguiente función

function res = correlacir(x,y)
 siz=max(size(x));
 res=zeros(1,siz);
 for i = [1:siz]
  res(i) = sum(x.*shiftcircizq(y,i-1));
 endfor
endfunction

que corresponde a esta definición matemática

La primera línea es la definición, y la siguiente describe propiedades de simetría al intercambiar el orden de las funciones a correlacionar. Según esta definición para la correlación con n=10 de x con y, multiplico término a término los elementos de x con una copia de y desplazada a la izquierda 10 posiciones. x(1) se multiplica con y(1+10), etc.

No es por tanto lo mismo la correlación de x con y que la de y con x, para un mismo valor de n. Pero sí es lo mismo la correlación de x con y para n=10, que la de y con x para n=-10.

Hay que andarse con cuidado a la hora de poner en el orden adecuado los vectores en la llamada a la función de correlación. En la anterior entrada del blog, la gráfica de la correlación con el pico en el paso 350 se hacía mostrando

cc19compruido=correlacir(xca19_, xcomp+ruido);

xca19_representa al PRN19 generado en el receptor(con valores +1 y -1), xcomp+ruido es la suma de las señales de los tres satélites y el ruido aleatorio, la señal total recibida. A la vista del orden, lo que se hace es coger el vector PRN19, luego el xcomp+ruido al que se desplaza a la izquierda entre 0 y 1022 pasos, obteniendo la función de autocorrelación para n=0 hasta n=1022. Pero resulta que eso es lo mismo que dejar xcomp+ruido quieto, y desplazar a la derecha xca19_.

Sin embargo si se calcula la correlación con correlacir(xcomp+ruido, xca19_), el pico del máximo no sale en el paso 350, sino en el 673 (es decir, 1023-350). ESto significa lo mismo pues por la periodicidad en 1023 de la función de correlación el máximo en 673 también está en -350, y desplazar a la izquierda -350 pasos el segundo parámetro, xca19_ , es justamente desplazarlo +350 a la derecha. Es lo mismo, pero si se interpreta bien.

En cuanto al ruido, se "genera" usando una función de octave ruido=4*randn(1,1023);, y hay que recordar que cada vez que se invoca la función se obtienen diferentes valores.

Para hacer las gráficas de un vector v se puede usar en octave el plot(v) sin más, y se pueden usar símbolos discretos para cada valor, como se hizo en la figura de autocorrelación de las PRN 05, 19 y 25. Si en vez de símbolos se usa una línea tendremos que esta consiste en segmentos rectilineos que unen los diferentes valores. Eso no sirve directamente para sacar un vector de valores +1 o -1 como pulsos rectangulares. Hay un tipo de gráfica en octave, el bar(...) que tampoco es lo que necesitaba, de modo que lo mejor es construir a partir de un vector de valores, otro con los puntos adecuados que unidos den los pulsos rectangulares deseados, una función como

#{
Una secuencia de valores en un vector v se quiere representar en forma de pulsos rectangulares,
que se inicien en x0 y de anchura xdelta. Por cada valor de v se emplean dos puntos desde el inicio al final
del intervalo correspondiente, a la altura del valor v. El primer punto se asume que es (x0,0)
#}
function [x,y] = graficabin(v, x0, xdelta)
 n=max(size(v));
 lon=2*n+1;
 x=zeros(1,lon);
 y=zeros(1,lon);
 x(1)=x0;
 y(1)=0;
 for i = [2:lon]
  y(i)=v(floor(i/2));
  x(i)=x0+xdelta*floor((i-1)/2);
 endfor 
endfunction

function plotbin(v,x0, xdelta)
 [x,y]=graficabin(v, x0, xdelta);
 margenx=(max(x)-min(x))/100;
 margeny=(max(y)-min(y))/50;
 plot(x,y);set (gca (), "xlim", [min(x)-margenx, max(x)+margenx]);set (gca (), "ylim", [min(y)-margeny, max(y)+margeny]);
endfunction

La segunda función es útil como sustituto rápido del plot.

Octave sigue siendo una herramienta inmejorable en relación calidad/precio para cacharrear con números, y hacer deberes. Aunque a veces hasta que consigues hacer algo pasas un buen (o mal) rato, como la simple tarea de comprobar los 10 primeros dígitos de las secuencias PRN del código C/A en octal, o los 16 primeros en hexadecimal. Pero al final es algo tan sencillo como
convertir de número a texto, quitar los espacios, pasar de binario a decimal y de decimal a octal o hexadecimal

> xca19(1:16)
ans =
1 1 1 0 0 1 1 0 1 1 0 1 0 1 1 0
> dec2base(bin2dec(strrep(num2str(xca19(1:10))," ","")),8)
ans = 1633
> dec2hex(bin2dec(strrep(num2str(xca19(1:16))," ","")))
ans = E6D6


La magia del PRN

La magia del PRN consigue que un cacharro "gepeese" sea capaz de escuchar a los satélites que allá arriba le mandan sus señales.

Un satélite GPS gira en una órbita casi circular a unos 20.180 km de altura sobre la superficie terrestre. En enviar su señal no consume demasiado, unos 27 vatios, como una bombillita de poca luz. Emplea una buena antena, que enfocada hacia la Tierra distribuye algo más de energía hacía las zonas más alejadas, respecto de las que le quedan justo debajo. Esto compensa en parte el hecho de que el debajo "sólo" está a 20.180km, y la parte rasante a unos 25.700km. Con todo, dada la lejanía del satélite respecto de un receptor típico, a este le llegan aproximadamente 10^-16 vatios. La señal es tan rematadamente débil que está inmersa en medio del ruido radioeléctrico existente debido a todo tipo de fuentes humanas presentes en la propia superficie.

Es como estar al fondo de una larga sala repleta de gente y esperar que nuestras palabras sean entendidas por alguién en el otro extremo, rodeado de personas en animada conversación.

Obviando muchos detalles, el caso es que podemos considerar que además de otra información, un satélite GPS repite sin cesar, cada milisegundo, una secuencia binaria (de ceros y unos) de 1023 dígitos. Cada satélite emplea una secuencia única, diferente a la empleada por el resto. El orden de los ceros y unos en la secuencia parece a primera vista aleatorio, sin una pauta clara, por lo que se les conoce como secuencias PRN.

PRN son las siglas inglesas de Pseudo Random Noise, ruido pseudo aleatorio, aunque puede usarse Number en vez Noise, entendiendo que no es pseudo aleatorio el número sino la secuencia. La señal formada por la repetición de esas secuencias de 1023 dígitos es lo que se llama el código C/A, pensado para uso civil, y encargado de ayudar a que el receptor pueda Adquirir (la A de C/A) la señal del satélite, además de hacer con él una estimación burda (la C en C/A es de Coarse) de la pseudodistancia. Hay otros códigos con secuencias PRN mucho más largas, como el P (de Preciso, en contraste a burdo) pensado para uso militar, o su versión encriptada Y (encrYpted).

El "Pseudo" indica que realmente la secuencia de ceros y unos no es al azar sino que se determina exactamente mediante ciertas reglas. Sin embargo tiene unas propiedades muy similares a las secuencias genuinamente aleatorias, y eso es lo que importa para su aplicación. Claro que el código C/A nos puede dar el pego de aleatorio como mucho 1 milisegundo, pues luego la secuencia se repite machaconamente.

¿Cuales son las reglas de creación de cada secuencia única? Se pueden explicar en términos de registros de desplazamiento de bits, y de operaciones lógicas. El caso es que esas reglas que usa el satélite, son conocidas y empleadas por los receptores. Y por cualquiera interesado si se lee por ejemplo la Interface Specification IS-GPS-200.

Yo he preferido seguir un excelente libro de texto sobre el GPS, el "Global Positioning System. Signals, Measurements, and Performance", segunda edición, de Pratap Misra y Per Enge, y hacer los "deberes" del final del capítulo 2, "GPS in 2005: An overview". Allí se pide calcular las secuencias PRN de 3 satélites distintos, el 05, el 19 y el 25, según las reglas explicadas. Las reglas son comunes a todos los satélites, cada satélite se distingue por usar unos parámetros con valores específicos. Una etiqueta como por ejemplo PRN05 designa tanto una secuencia PRN concreta, como el único satélite GPS que la emite.

Usando de nuevo octave, las gráficas de las secuencias de esos 3 satélites tienen esta pinta:



La secuencia de 1023 ceros y unos lógicos se muestra como una señal digital de pulsos, usando un valor de la señal de +1 para representar los ceros, y una señal con valor -1 para representar los unos. Cada dígito se convierte así en lo que se llama un "chip", un pulso rectangular de amplitud positiva o negativa. Hay 1023 chips en cada secuencia repetida. Visto como en la imagen anterior parece que se trata de un código de barras, pero acercandonos para ver sólo los 200 primeros chips de cada satélite, se aprecia mejor el sube y baja de la señal, sobre todo si se abre la imagen en una ventana para ella sola.



En estas figuras en el eje horizontal se representa el tiempo de emisión del chip. Si el primer chip se emite en t=0, el segundo se emite en t=1, ..., y el último chip de la secuencia, el 1023, se emite en t=1022, todo expresado en unidades de chip (de duración de chip), donde 1 chip equivale a 1/1023 milisegundos.

Y ahora llega la propiedad mágica crucial de las secuencias PRN. La correlación entre dos secuencias PRN de dos satélites cualesquiera es casi nula. La autocorrelación de cada secuencia PRN consigo misma desplazada es también casi nula, salvo para el desplazamiento nulo.

La correlación entre secuencias cualesquiera de ceros y unos, de la misma longitud, es algo muy sencillo. Poniendolas en paralelo se van comparando los elementos que ocupan la misma posición. Si sus valores coinciden, sumo 1. Si difieren, resto 1. Lo que me dé la suma al final es la correlación, si bien puedo normalizar diviendo entre la longitud de las secuencias. Es como hacer una operación de lógica booleana, el OR exclusivo o XOR.

Otra forma de cálcular el mismo resultado, empleando los niveles del "chip", +1 para los ceros y -1 para los unos, es la multiplicación. Cuando hay coincidencia tanto 1x1 como (-1)x(-1) dan 1. Cuando hay discrepancia (-1)x1 = 1x(-1)= -1. Agregando los resultados de tales multiplicaciones llegamos a lo mismo que antes, aunque así tenemos una definición más general, aplicable a cualquier señal con amplitudes tanto negativas como positivas, lo que nos será útil más adelante.

Si comparamos una secuencia con una copia de sí misma desplazada un número variable de pasos, se habla de autocorrelación. Si la copia no se desplaza (o se "desplaza" 0 pasos), siempre se suma 1 por cada dígito o chip, nunca se resta. Para una secuencia genuinamente aleatoria de 1023 dígitos la autocorrelación (sin normalizar) será de 1023. Por otro lado si la comparo con una copia exacta pero desplazada a derecha o izquierda, por la misma aleatoriedad de la situación de ceros y unos, casi tantas veces tendré coincidencias, sumando 1, como discrepancias, restando 1, y la suma total estará muy cerca de 0.

Este tipo de correlación es "circular". Al comparar una secuencia con otra desplazada a la derecha 10 pasos (cada paso corresponde a un dígito o a un chip), entendemos que el valor de la posición 1 pasa a la posición 11, el de la 2 a la 12, ..., el de la 1013 a la 1023, el de la 1014 a la 1, ..., y el de la 1023 a la 10. Si el desplazamiento es de 10 pasos a la izquierda, la posición 1 pasa a la 1014, ..., la 10 a la 1023, la 11 a la 1, ..., y la 1023 a la 1013.

La siguiente imagen refleja los valores de la autocorrelación para las PRN 05, 19 y 25, cada una con una copia desplazada un número de pasos (eje horizontal) entre 0 y 50. El aspecto para los pasos del 51 al 1022 es parecido al de los pasos 1 a 50. Un desplazamiento de 1023 pasos lleva a la configuración inicial.



Por cierto, esos pasos ¿son hacia la derecha o hacia la izquierda? En el caso de la autocorrelación da igual, por la simetría. Si comparo A con una copia C, de A desplazada a la izquierda 10 pasos, es igual que comparar C con una copia suya desplazada a la derecha 10 pasos, que resulta ser A.

Hay que fijarse que para desplazamiento cero siempre se obtiene la máxima correlación de 1023, faltaría más. Las funciones de autocorrelación en cada uno de los tres casos no coinciden exactamente, pero todas (esas tres y las del resto de satélites) comparten algo en común: los únicos valores de autocorrelación (sin normalizar) que se dan para cualquier paso entre 1 y 1022 son los tres siguientes: 63, -1 y -65. La clave es que esos valores son mucho más pequeños que el máximo, 1023.

Además la correlación entre la PRN de un satélite y la de cualquier otro satélite distinto, para cualquier paso (incluido el cero) también tiene siempre uno de esos tres valores: 63, -1, -65.

Las PRN de los satélites GPS para el código C/A son un tipo de secuencias matemáticas estudiadas por el Dr. Robert Gold, por lo que se conocen como códigos de Gold. Las propiedades enunciadas arriba se pueden probar rigurosamente.

El caso es que cada satélite emite una especie de huella dactilar. Y los receptores GPS tienen un "detector de huellas", un correlacionador. Es un circuito electrónico que genera la secuencia PRN de chips de un satélite dado, o su versión desplazada cierto número de pasos, y durante un milisegundo calcula la correlación de dicha secuencia con la señal que le llega.

Imaginemos que el correlacionador va a ver si recibe algo del satélite 19 en la señal que, captada por la antena del receptor, es acondicionada previamente por otros componentes del receptor. Para ello calculará la correlación entre la PRN19 generada localmente en el receptor y la señal. Por cada valor de la función de correlación para un cierto número de pasos de desplazamiento, hay que generar durante un milisegundo la secuencia PRN19 desplazada ese número de pasos, multiplicarla chip a chip por la señal recibida e ir sumando para tener al final el valor buscado. Este cálculo hay que repetirlo para los pasos del cero al 1022, hasta obtener los 1023 puntos de la función correlación.

Imaginemos también (problema 2.2(f) de Misra&Enge) que la señal que le llega es una suma de una señal completamente aleatoria, con distribución normal, media cero y desviación estándar de 4 unidades, junto con tres señales de los satélites 05, 19 y 25, de amplitud 1 (inferior pues a la del ruido aleatorio).

Como la distancia de cada satélite hasta el receptor será diferente, el origen de la secuencia de cada satélite (emitido cada milisegundo del reloj del satélite), llegará al receptor con un retraso variable. Supongamos que respecto a la marca de milisegundo del reloj del receptor, cuando comienza este a generar su PRN particular, el origen de la PRN05 recibida del satélite va retrasada 75 chips, el origen de la PRN19 recibida va retrasada 350 chips, y el origen de la PRN25 recibida va retrasada 905 chips. (Aunque es equivalente decir, por ejemplo, que el origen de la PRN05 desde su satélite va adelantada 1023-75 = 948 chips).

Bueno, pues la mezcla de las señales indicadas tendrá el siguiente aspecto, que dificilmente da idea de satélite alguno.




Teniendo que calcular 1023 valores, empleando un milisegundo por cada uno, en poco más de un segundo el correlacionador obtiene ... ¡algo así!:



Tatacháaaaaaan. Un maravilloso pico que destaca primoroso ¡justo en el paso 350! ¡Magia!

Ese pico por un lado descubre el código PRN del satélite 19 entre el aparente ruido informe, y por otro no menos importante, mide la pseudo distancia entre el satélite 19 y el receptor, salvo por un número entero de milisegundos.

Con la información del correlacionador, el receptor no sabe si la señal que le llega salió del satélite 19 hace 66 milisegundos y pico, o hace 67 milisegundos y pico, o hace 68 milisegundos y pico, o... pero sí sabe que el "pico" es de 350/1023 milisegundos. Cada chip dura 1/1023 milisegundos, tiempo en que la luz recorre "sólo" unos 293 metros. Los circuitos electrónicos pueden afinar y precisar el paso fraccionario en que se da el máximo de correlación, llegando a las milésimas, por ejemplo 350,276 en vez de 350. En una milésima de chip la luz recorre casi 30 centímetros, y ese es el orden de error en la estimación de la pseudo distancia basada en el código C/A (sin olvidar que es pseudo, no la distancia real, por el tema de la falta de sincronía entre el reloj del receptor y del satélite, y que hay otras fuentes de errores).

Vamos a ver esto del "pico". La función de correlación de la figura se obtiene al comparar la señal recibida con una copia de la secuencia PRN19 desplazada n pasos a la derecha. En términos temporales "desplazar a la derecha" significa "emitir con retraso": cada valor en la copia se emite n chips más tarde que en el original. El máximo encontrado en el paso 350 se da al utilizar una copia de PRN19 retrasada 350 chips. Eso indica que hay una coincidencia con el PRN19 original del satélite cuando esperamos 350 chips tras la marca de milisegundo del receptor. Esos 350 chips son el intervalo de espera entre una marca de milisegundo del reloj del receptor y la siguiente marca de milisegundo del reloj del satélite.

Tras captar el máximo de correlación y adquirir la señal del satélite 19, consideremos el momento en que por el reloj del receptor llega la marca de milisegundo X. El receptor sabe que llegará una marca de milisegundo del reloj del satélite dentro de 350/1023 milisegundos, es decir, cuando sea el tiempo de recepción X+(350/1023)ms. Pero no sabe qué marca de milisegundo es la que llega, cuándo fué emitida esa marca por el satélite. Necesita aún que alguién le "sople" que esa marca del satélite es la del milisegundo Y. Entonces el receptor ya podrá estimar la pseudodistancia como c por el intervalo temporal entre emisión y recepción: 0,001.c.(X+(350/1023)-Y).

¿Quién le sopla al receptor el tiempo exacto de emisión, según el reloj del satélite, de la marca de milisegundo que le llega? Pues la propia señal del satélite, que además del código C/A con su PRN, incluye bits de datos con todo lo necesario. Aunque como se transmiten a un ritmo muy lento de 50 bits por segundo, hace falta un poco de paciencia para poder seguir la señal (tracking). Pero antes de poder seguirla, y de poder acceder a esos bits de datos, hay que adquirir (adquisition) la señal, que es donde interviene el correlacionador.

Pero nunca es tan facil. Hay tantos detalles... Hemos asumido que la antena capta la señal, algún sistema electrónico la "acondiciona", y ya puede el correlacionador calcular punto a punto la correlación. Ya es bastante trabajo que esa función tenga 1023 puntos, y que haya que probar las PRN de todos los satelites, porque si el receptor de entrada no sabe dónde está, no sabe qué satélites pueden estar a la vista. Es que además la señal que el satélite emite en una frecuencia muy determinada se modifica por el efecto Doppler, debido al movimiento del satélite y el posible del receptor. Este tiene que buscar en varios intervalos de frecuencias, sin saber exactamente en cúal estará la señal de cada satélite, y por cada intervalo elegido hay que repetir los cálculos del correlacionador.

Es como si éste fuese un labriego que busca un anillo perdido en un surco, recorriendo con cuidado el mismo a ver si descubre su brillo, y recibe una aclaración de su amada, y propietaria del anillo: "... pero no estoy segura de que fuese en ese surco, salté por todo el terreno, y perdí no uno sino todos mis anillos, 4 ... o 10, no estoy segura".

Menos mal que estos cacharricos electrónicos son muy rápidos, y no tienen venas.

La magia del GPS

Parece que Arthur C. Clarke, además de escribir libros de Ciencia Ficción y de idear el satélite de comunicaciones geoestacionario, dijo eso de que

Cualquier tecnología suficientemente avanzada es indistinguible de la magia.

A los efectos prácticos de cualquier usuario de navegador parlanchín que nos manda girar a izquierda o derecha, que ese cacharro sepa dónde estamos y por dónde ir es asimilable a la magia.

¿Qué "magia" usa el navegador o el móvil o cualquier dispositivo con GPS para saber dónde está?

Aunque llamamos a menudo "el gepeese" a los cacharritos que nos indican dónde estamos, realmente GPS hace referencia al Sistema de Posicionamiento Global, propiedad del ejército USA. Los soviéticos dieron la réplica con su GLONASS, y la Unión Europea anda tras una versión propia y civil, Galileo, al que citaba en noticias del espacio. También China está poniendo en marcha su propio sistema Compass.

Todos esos sistemas se engloban bajo unas siglas, GNSS, Sistemas de Navegación Global por Satélite. Permiten determinar (para ayudar en la Navegación) la posición y velocidad de un receptor, en cualquier parte de la Tierra (por eso lo de Global), con la ayuda de una constelación de Satélites. Para conseguirlo se hace uso no de una, sino de muchas, muchas magias acumuladas por los hechiceros científicos a lo largo de los siglos.

La magia más esencial y poderosa, que sustenta al resto de magias, es la Matemática. Por ejemplo, la geometría de un espacio euclídeo tridimensional nos permite resolver el problema siguiente: Dados tres puntos distintos S₁, S₂ y S₃ no colineales, determinar el punto R tal que su distancia a S₁ es D₁, a S₂ es D₂ y a S₃ es D₃. Basta considerar las esferas centradas en cada punto S_i, de radios dados por las distancias D_i, y en su intersección común, de existir, se encontrará R. Por la simetría de las 3 esferas consideradas respecto al plano determinado por sus 3 centros, si hay una solución a un lado del plano, habrá otra "reflejada", y si la solución resulta estar en ese plano, será única.

Nuestros modernos cacharros electrónicos no saben hacer construcciones geométricas, pero se les da de maravilla operar con números. Un gran paso en la historia de las Matemáticas fue el desarrollo de la geometría analítica. Las construcciones geométricas se convierten en ecuaciones algebraicas, y la solución concreta de encontrar las 3 coordenadas cartesianas de R dadas las 9 de S₁, S₂ y S₃, y las tres distancias D₁, D₂, D₃, implica resolver 3 ecuaciones (de segundo grado) con 3 incógnitas.

(R_x-S_ix)²+(R_y-S_iy)²+(R_z-S_iz)² = D_i²   (i=1,2,3)

Ya podemos usar los rápidos calculistas electrónicos, cada vez más chicos y poderosos gracias a las magias cuánticas y lógicas, para resolver las ecuaciones.

Pero antes de resolver las ecuaciones y despejar las incógnitas, hay que obtener los datos, ¡y esa es la parte más difícil! Esos 3 puntos S₁, S₂ y S₃ no colineales van a ser 3 satélites, no 3 puntos fijos. La razón para usar una constelación de satélites esta en la G de GPS, que allá donde vayamos sobre la Tierra y sus inmediaciones aéreas, encontremos los puntos de referencia necesarios, sobre nuestras cabezas. Resulta que con unos 24 a 30 satélites bien repartidos hacemos el apaño.

Mediante las magias correspondientes, largas de enumerar, se construyen y lanzan los satélites, y se les "sigue" estrechamente para saber no sólo dónde están sino por dónde van a estar, decírselo a los propios satélites, y que estos a su vez se lo digan a todo receptor que quiera y pueda entenderlos. Un satélite de un sistema GNSS en esencia va emitiendo continuamente una información como "soy satelifulanito de tal, ahora que son las tantas estoy justo en tales coordenadas".

Bueno, si soy el punto R del problema geométrico y a la vez un Receptor que entiende las señales del satélite S₁, del S₂ y del S₃, podría en un instante dado tener los datos de las 9 coordenadas de los 3 puntos que ocupaban los satélites cuando salieron los 3 mensajes que me llegan a la vez. Solo me hace falta obtener la distancia a cada uno de esos 3 puntos.

Aquí la geometría recibe una ayuda de la física, que le dice "esa información desde el satélite llega al receptor por una onda electromagnética que viaja a la velocidad de la luz, c=299.792.458 m/s". Como el espacio recorrido es igual a la velocidad por el intervalo de tiempo empleado, y en el mensaje desde el satélite me dan "la hora" en que estaba en "esas coordenadas", miro el momento en que me ha llegado el mensaje (más tarde que cuando salió del satélite), hago la diferencia entre tiempo de salida y tiempo de llegada, multiplico por c, y ya tengo (repitiendo el proceso para cada satélite) D₁, D₂, D₃, y a resolver las ecuaciones.

Llegados a este punto nos damos de bruces con un problemilla. Dar con la distancia entre S₁ y R implica restar tiempos, medidos por diferentes relojes. S₁ y R cada uno tiene su reloj. Si ambos estuviesen marcando exactamente la hora en perfecta sincronía, nada que objetar. Los relojes que van en los satélites son de lo mejorcito, relojes atómicos de Cesio o de Rubidio, varios en cada satélite por redundancia. Pero el reloj de un receptor, sin ánimo de ofender, no les llega a las suelas de los zapatos.

¿Entonces? Imaginemos un caso aproximadamente real. El satélite S₁ está sobre la cabeza de nuestro receptor R a 20.100km de altura, y pasa por allí a las 12h 0minutos 0segundos 0 microsegundos mandando su mensaje, que le llega a R cuando por el reloj de R han pasado 66 milisegundos de las 12h en punto. Redondeando c=300.000km/s, la distancia estimada por R será 0,066.c = 19.800km. Por llevar el reloj adelantado 1 milisegundo (marca 66 cuando son 67), R comete un error en su estimación de la distancia a S₁ de 300km. Incluso un adelanto tan pequeño como 10 microsegundos en el reloj de R implicaría un error en D₁ de 3 km.

La conclusión es que R no puede estimar directamente la distancia verdadera a un satélite, D, pero sí una pseudodistancia pseudoD, que encierra una discrepancia c. deltaR respecto a la verdadera por culpa del adelanto (deltaR > 0) o atraso (deltaR < 0) de su reloj.

La nueva incógnita deltaR es desconocida pero la misma para cualquier distancia a diferentes satélites. Una nueva incógnita requiere, para resolver el sistema de ecuaciones, una nueva ecuación, que se obtiene fácilmente si además de los satélites S₁, S₂ y S₃ interviene un cuarto satélite S₄ favorablemente situado. El nuevo sistema de ecuaciones a resolver es

(R_x-S_ix)²+(R_y-S_iy)²+(R_z-S_iz)² = (c . deltaR)² + pseudoD_i²   (i=1,2,3,4)

En el ejemplo de antes, uno de los datos obtenidos por R sería que pseudoD₁=19.800km. Al resolver el sistema con las otras 3 pseudodistancias y las posiciones de los 4 satélites, debería obtenerse la posición de R y que deltaR es aproximadamente 1 milisegundo de adelanto sobre el tiempo verdadero.

En esencia así se obtiene la posición de R. Aunque nunca es tan fácil. En la vida real todo es fuente de error. Las posiciones de los satélites no se conocen con total exactitud. Las ondas electromagnéticas que llevan información del satélite al receptor pasan por la ionosfera y la estratosfera, donde su velocidad varía según las condiciones ambientales. Las características de esa señal se ven influidas además por el movimiento tanto del satélite como del receptor, con el efecto Doppler cambiando frecuencias. También hay que contar con los efectos relativistas. Y con el ruido y demoras en los circuitos electrónicos. Incluso los avanzados relojes atómicos de los satélites tienen sus errores.

Ante tanta "fuente de error", la Física y la Matemática contraatacan modelizando esas fuentes y sus "errores". En general se sigue una estrategia como la anterior cuando se ampliaban las ecuaciones de 3 a 4, con una nueva incógnita. Cada fuente de variabilidad adecuadamente modelizada conlleva una serie de nuevos parámetros con valores desconocidos pero que se pueden estimar incorporando nuevas ecuaciones. Cuanto más datos y más ecuaciones, más se mejoran las estimaciones de los parámetros de los modelos, y más se reducen los errores.

La resolución de ecuaciones y estimación de parámetros es un asunto denso y un tanto aburrido, no toda la magia es igual de atrayente. Pero hay números de magia especialmente sorprendentes, como el que hace cada Satelifulanito para decir quién es y, el más difícil todavía, con el que un Receptor avispado le escucha e identifica incluso en medio de un ensordecedor ruido. Es ... la magia del PRN.