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