Splitting NAV frames into fields

Pulling out the individual values out of the NAV frames, then converting them into the correct range of values is quite a struggle. Here is my solution!

First off, here are a few worker routines that pull a range of bits out of the NAV subfields (where the parity bits are in the lowest 6 bits of each value):

Return the lowest ‘n_bits’ bits of the value ‘u’:

unsigned int mask(unsigned u, int n_bits) {
  return u & ((1<<n_bits)-1);
}

Bit ‘len’ is the sign, and if set it needs to be extended into all the more significant bits – note this returns an ‘int’ type:

int sign_extend(unsigned u,int len) {
 if(len < 32 && u >0)
 if(u>>(len-1)&1) u |= 0xFFFFFFFF << len;
 return (int)u;
}

This is the key function to extract bits. It returns ‘len’ bits of the value ‘val’,starting with the bit at ‘offset’:

unsigned int bits(int val, int offset, int len) {
 return mask(val >> offset,len);
}

This extracts two sets of bits, from two different values, and then concatenates them together to build an unsigned value

unsigned join_bits_u(int val1, int offset1, int len1,
                     int val2, int offset2, int len2) {
  return (bits(val1, offset1, len1) << len2) | bits(val2, offset2, len2);
}

This also extracts two sets of bits and concatenates them together, however it also sign-extends the leftmost bit, making a signed value:

signed join_bits_s(int val1, int offset1, int len1, 
                   int val2, int offset2, int len2) {
  return sign_extend(join_bits_u(val1, offset1, len1, val2, offset2, len2),len1+len2);
}

With this groundwork in place, extracting the parameters becomes pretty easy, and by always calling “join_bits” even if only one set of bits is needed makes it quite clean:

/******** From FRAME 2 *******/
 iode1     = join_bits_u(           0, 0, 0, n->frame2[2], 22, 8);
 Crs       = join_bits_s(           0, 0, 0, n->frame2[2],  6,16) * pow(2.0,-5.0);
 M0        = join_bits_s(n->frame2[3], 6, 8, n->frame2[4],  6,24) * pow(2.0,-31.0) * PI;
 delta_n   = join_bits_s(           0, 0, 0, n->frame2[3], 14,16) * pow(2.0,-43.0) * PI;
 Cuc       = join_bits_s(           0, 0, 0, n->frame2[5],  6,16) * pow(2.0,-29.0);
 e         = join_bits_s(n->frame2[5], 6, 8, n->frame2[6],  6,24) * pow(2.0,-31.0);
 Cus       = join_bits_s(           0, 0, 0, n->frame2[7], 14,16) * pow(2.0,-29.0);
 sqrt_A    = join_bits_u(n->frame2[7], 6, 8, n->frame2[8],  6,24) * pow(2.0,-19.0);
 Toe       = join_bits_u(           0, 0, 0, n->frame2[9], 14,16);
 fit_flag  = join_bits_u(           0, 0, 0, n->frame2[9], 13, 1);
 aodo      = join_bits_u(           0, 0, 0, n->frame2[9],  8, 6);

/******** From FRAME 3 *******/
 Cic       = join_bits_s(           0, 0, 0, n->frame3[2], 14,16) * pow(2.0,-29.0);
 omega_0   = join_bits_s(n->frame3[2], 6, 8, n->frame3[3],  6,24) * pow(2.0,-31.0) * PI;
 Cis       = join_bits_s(           0, 0, 0, n->frame3[4], 14,16) * pow(2.0,-29.0);
 i_0       = join_bits_s(n->frame3[4], 6, 8, n->frame3[5],  6,24) * pow(2.0,-31.0) * PI;
 Crc       = join_bits_s(           0, 0, 0, n->frame3[6], 14,16) * pow(2.0,-5.0);
 w         = join_bits_s(n->frame3[6], 6, 8, n->frame3[7],  6,24) * pow(2.0,-31.0) * PI;
 omega_dot = join_bits_u(           0, 0, 0, n->frame3[8],  6,24) * pow(2.0,-43.0) * PI;
 iode2     = join_bits_u(           0, 0, 0, n->frame3[9], 22, 8);
 idot      = join_bits_s(           0, 0, 0, n->frame3[9],  8,14) * pow(2.0,-43.0) * PI;

 

Advertisements

Calculating orbits

Once you have acquired a GPS Space Vehicle, you are than faced with the daunting task of calculating where it is at any given instant.

NAV frames 2 and 3 provide information about the orbit, and NAV frame 1 provides information about the time.

Using the following untested code it is possible to print out the Earth Centered, Earth Fixed location of the Space Vehicle at any given time.

I haven’t verified it yet in a working design, but plotting the output makes it look correct. I’ll revisit this code later after my design is working.

/************************************************
* orbit.c - calculate a Space Vehicle's position 
* Last changed at 2016-09-24 
************************************************/
#include <stdio.h>
#include <math.h>

#define PI 3.1415926535898

double eccAnom(double ec, double mean_anomaly, double dp) {
 int iterations = 30;
 double delta=pow(10,-dp);
 double estimate, correction;

/* First estimate */
 estimate = (ec<0.8) ? mean_anomaly : PI;
 correction = estimate - ec*sin(mean_anomaly) - mean_anomaly;
 
 /* Solve iteratively */
 while ((abs(correction)>delta) && iterations > 0) {
 estimate -= correction/(1.0-ec*cos(estimate));
 correction = estimate - ec*sin(estimate) - mean_anomaly;
 iterations--;
 }
 return estimate;
}

int orbit_calc_position(
  double t_now,
  /* Orbital parameters */
  double sqrt_A,
  double e,
  double w,
  /* Position at ephemeris */
  double time_of_ephemeris,
  double mean_motion_at_ephemeris,
  double omega_at_ephemeris,
  double inclination_at_ephemeris,
 
 /* Corrections */
  double delta_mean_motion,
  double Cus, double Cuc,
  double Cis, double Cic,
  double Crs, double Crc,
  double IDOT,
  double omegaDot,

  /* Outputs */
  double *ecef_x, double *ecef_y, double *ecef_z)
{
 const double mu = 3.986005e14; /* Earth's universal gravitation parameter */
 const double omegaDot_e = 7.2921151467e-5; /* Earth's rotation rate */


 double semi_major_axis;
 double computed_mean_motion;
 double time_from_ephemeris;
 double correct_mean_motion;
 double mean_anomaly;
 double ek;
 double true_anomaly;
 double argument_of_latitude;
 double argument_of_latitude_correction;
 double radius_correction;
 double correction_of_inclination;
 double corrected_argument_of_latitude;
 double corrected_radius;
 double corrected_inclination;
 double pos_in_orbital_plane_x;
 double pos_in_orbital_plane_y;
 double corrected_angle_of_ascending_node;

/***********************
 * Calculate orbit
 ***********************/
 semi_major_axis = sqrt_A * sqrt_A;
 computed_mean_motion = sqrt(mu / pow(semi_major_axis,3.0));
 time_from_ephemeris = t_now - time_of_ephemeris;
 correct_mean_motion = computed_mean_motion + delta_mean_motion; 
 /* Add on how much we have moved through the orbit since ephemeris */
 mean_anomaly = mean_motion_at_ephemeris + correct_mean_motion * time_from_ephemeris;

 /* solve for eccentric anomaly */
 ek = eccAnom(e,mean_anomaly,10);
 /* Now calculate the first approximation of the latitude */
 argument_of_latitude = w + ek;

/*****************************************
 * Second Harmonic Perbturbations 
 *****************************************/
 argument_of_latitude_correction = Cus * sin(2*argument_of_latitude) + Cuc * cos(2*argument_of_latitude);
 corrected_argument_of_latitude = argument_of_latitude+argument_of_latitude_correction;

 radius_correction = Crc * cos(2*argument_of_latitude) + Crs * sin(2*argument_of_latitude);
 corrected_radius = semi_major_axis * (1- e * cos(ek)) + radius_correction;
 
 correction_of_inclination = Cic * cos(2*argument_of_latitude) + Cic * sin(2*argument_of_latitude);
 corrected_inclination = inclination_at_ephemeris + correction_of_inclination + IDOT*time_from_ephemeris;

 pos_in_orbital_plane_x = corrected_radius * cos(corrected_argument_of_latitude);
 pos_in_orbital_plane_y = corrected_radius * sin(corrected_argument_of_latitude);


corrected_angle_of_ascending_node = omega_at_ephemeris + (omegaDot - omegaDot_e)*time_from_ephemeris - omegaDot_e * time_of_ephemeris;

/******************************************************
 * Project into Earth Centered, Earth Fixed coordinates
 ******************************************************/
 *ecef_x = pos_in_orbital_plane_x * cos(corrected_angle_of_ascending_node)
         - pos_in_orbital_plane_y *cos(corrected_inclination) * sin(corrected_angle_of_ascending_node);
 *ecef_y = pos_in_orbital_plane_x * sin(corrected_angle_of_ascending_node) 
         + pos_in_orbital_plane_y *cos(corrected_inclination) * cos(corrected_angle_of_ascending_node);
 *ecef_z = pos_in_orbital_plane_y * sin(corrected_inclination);

#if 0
 printf("Semi major axis                   %20.10f\n", semi_major_axis);
 printf("Computed mean motion              %20.10f\n", computed_mean_motion);
 printf("Time from ephemeris               %20.10f\n", time_from_ephemeris);
 printf("Corrected mean_motion             %20.10f\n", correct_mean_motion);
 printf("mean_anomaly                      %20.10f\n", mean_anomaly);
 printf("ek                                %20.10f\n", ek);
 printf("true anomaly                      %20.10f\n", true_anomaly); 
 printf("argument of latitude              %20.10f\n", argument_of_latitude); 
 printf("corrected_argument_of_latitude    %20.10f\n", corrected_argument_of_latitude);
 printf("corrected_radius                  %20.10f\n", corrected_radius);
 printf("correction_of_inclination         %20.10f\n", correction_of_inclination);
 printf("pos_in_orbital_plane_x            %20.10f\n", pos_in_orbital_plane_x);
 printf("pos_in_orbital_plane_y            %20.10f\n", pos_in_orbital_plane_y);
 printf("corrected_angle_of_ascending_node %20.10f\n",corrected_angle_of_ascending_node);
#endif
#if 0
 /* Checks that nothing is amiss */
 printf("%12.2f ", sqrt(pos_in_orbital_plane_x*pos_in_orbital_plane_x+pos_in_orbital_plane_y*pos_in_orbital_plane_y));
 printf("%12.2f\n", sqrt((*ecef_x)*(*ecef_x)+(*ecef_y)*(*ecef_y)+(*ecef_z)*(*ecef_z)));
#endif
}

int main(int c, char *v[]) {
 double t_now = 7241; /* zero time */

/* Orbital parameters, transmitted by satellite */
 double sqrt_A = 5153.800608;
 double time_of_ephemeris = 7241;
 double delta_n = 4.279821129e-009;
 double M0 = -2.600551;
 double e = 0.030228948;
 double w = -1.048966714;
 double Cuc = -0.000020975247; 
 double Cus = 0.000004233792;
 double Cic = -0.000000014901;
 double Cis = 0.000000130385;
 double Crc = 305.562500;
 double Crs = 3.812500;
 double i0 = 0.978481871;
 double IDOT = 1.689356e-010;
 double omega_0 = -0.494186143;
 double omegaDot = 0.000005984;

/* Outputs */
 double ecef_x, ecef_y, ecef_z;

 /* Print out position 24 minutes apart */
 for(t_now == 7241; t_now < 7241+60*60*24*3; t_now += 60*24) {
   orbit_calc_position(
     t_now,
     /* orbital parameters */
     sqrt_A, e, w, 
     /* Time and position at ephemeris */
     time_of_ephemeris, M0, omega_0, i0,
     /* correction factors */
     delta_n, Cus, Cuc, Cis, Cic, Crs, Crc,
     IDOT, omegaDot,
     /* Outputs */ 
     &ecef_x, &ecef_y, &ecef_z);

   printf("ECEF at %10.1f = ( %12.2f, %12.2f, %12.2f )\n", t_now, ecef_x, ecef_y, ecef_z);
 }
return 0;
}

Three ways for verifying NAV subframe parity

The GPS specs are pretty average on the how the parity calculation for each subframe is performed, with things that I found confusing.

The most confusing thing for me is that the bit numbering within the fields is most significant bit first, so if the Time of Week is from bits 1 to 17 in the subframe bit 1 holds the  2^16 bit, and bit 17 holds the 2^0 bit.

To add to the fun, the bits are numbered from 1 to 30, rather than the more programmer-friendly 0 to 29.

And even worse, it is documented to be viewed as a hardware operation with lots of XOR gates, and that is not to be efficiently be implemented in software.

The routine below takes the last 32 bits to be received in ‘d’ and returns 1 if the parity is correct, or zero if false.

It does this three ways – a either a mass of bit XOR operations (which is great for hardware), or a set of test-bit-then-XOR operations, either from a constants or from a lookup table (which is great for software) – you can select which one to use by changing the “#if” values.

int nav_test_parity(unsigned int d)
{
 /* 'd' holds the 32 bits received, and this function will
 * return true if the low 30 bits is a valid subframe - 
 * the most recent bit is held in bit zero */
 
 /* If the last bit of the last subframe is set
 * then the data in this frame is flipped */
 if(d & 0x40000000)
   d ^= 0x3FFFFFC0;

#if 1
  static const unsigned char parity[32] = {
    0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x13, 0x25,
    0x0B, 0x16, 0x2C, 0x19, 0x32, 0x26, 0x0E, 0x1F,
    0x3E, 0x3D, 0x38, 0x31, 0x23, 0x07, 0x0D, 0x1A,
    0x37, 0x2F, 0x1C, 0x3B, 0x34, 0x2A, 0x16, 0x29}; 
 for(i = 6; i < 32; i++) {
   if(d&(1<<i)) d ^= parity[i]; 
 }  
 if((d & 0x3F) != 0) return 0; 
#endif 
#if 0
 /* Bits 3 through 0 */  
  if(d & 0x00000001) d ^= 0x00; // Does nothing  
  if(d & 0x00000002) d ^= 0x00; // Does nothing  
  if(d & 0x00000004) d ^= 0x00; // Does nothing
  if(d & 0x00000008) d ^= 0x00; // Does nothing
  /* Bits 7 through 6 */
  if(d & 0x00000010) d ^= 0x00; // Does nothing
  if(d & 0x00000020) d ^= 0x00; // Does nothing
  if(d & 0x00000040) d ^= 0x13;  
  if(d & 0x00000080) d ^= 0x25;  
  /* Bits 11 through 8 */  
  if(d & 0x00000100) d ^= 0x0B;  
  if(d & 0x00000200) d ^= 0x16;  
  if(d & 0x00000400) d ^= 0x2C;  
  if(d & 0x00000800) d ^= 0x31;  
  /* Bits 15 through 12 */  
  if(d & 0x00001000) d ^= 0x32;  
  if(d & 0x00002000) d ^= 0x26;  
  if(d & 0x00004000) d ^= 0x0E;  
  if(d & 0x00008000) d ^= 0x1F;  
  /* Bits 19 through 16 */  
  if(d & 0x00010000) d ^= 0x3E;  
  if(d & 0x00020000) d ^= 0x3D;   
  if(d & 0x00040000) d ^= 0x38;  
  if(d & 0x00080000) d ^= 0x31;  
  /* Bits 23 through 20 */  
  if(d & 0x00100000) d ^= 0x23;  
  if(d & 0x00200000) d ^= 0x07;   
  if(d & 0x00400000) d ^= 0x0D;   
  if(d & 0x00800000) d ^= 0x1A;   
  /* Bits 27 through 24 */  
  if(d & 0x01000000) d ^= 0x37;   
  if(d & 0x02000000) d ^= 0x2F;   
  if(d & 0x04000000) d ^= 0x1C;
  if(d & 0x08000000) d ^= 0x3B;
  /* Bits 29 through 28 */
  if(d & 0x10000000) d ^= 0x34;
  if(d & 0x20000000) d ^= 0x2A;
  /* Bits 1 through 0 of the previous subframe */ 
  if(d & 0x40000000) d ^= 0x16;
  if(d & 0x80000000) d ^= 0x29;
  /* Test all 6 parity bits (5 downto 0) at once */
  if((d & 0x3F) != 0) return 0; #
endif 
#if 0  char D24, D25, D26, D27, D28, D29;  /* A mass of XOR operations */  D24 = (l>> 1) ^ (d>>29) ^ (d>>28) ^ (d>>27) ^ (d>>25)
     ^ (d>>24) ^ (d>>20) ^ (d>>19) ^ (d>>18) ^ (d>>17) 
     ^ (d>>16) ^ (d>>13) ^ (d>>12) ^ (d>>10) ^ (d>> 7); 
 D25 = (l>> 0) ^ (d>>28) ^ (d>>27) ^ (d>>26) ^ (d>>24)
     ^ (d>>23) ^ (d>>19) ^ (d>>18) ^ (d>>17) ^ (d>>16)
     ^ (d>>15) ^ (d>>12) ^ (d>>11) ^ (d>> 9) ^ (d>> 6);
 D26 = (l>> 1) ^ (d>>29) ^ (d>>27) ^ (d>>26) ^ (d>>25)
     ^ (d>>23) ^ (d>>22) ^ (d>>18) ^ (d>>17) ^ (d>>16) 
     ^ (d>>15) ^ (d>>14) ^ (d>>11) ^ (d>>10) ^ (d>> 8);
 D27 = (l>> 0) ^ (d>>28) ^ (d>>26) ^ (d>>25) ^ (d>>24) 
     ^ (d>>22) ^ (d>>21) ^ (d>>17) ^ (d>>16) ^ (d>>15)
     ^ (d>>14) ^ (d>>13) ^ (d>>10) ^ (d>> 9) ^ (d>> 7);
 D28 = (l>> 0) ^ (d>>29) ^ (d>>27) ^ (d>>25) ^ (d>>24) 
     ^ (d>>23) ^ (d>>21) ^ (d>>20) ^ (d>>16) ^ (d>>15) 
     ^ (d>>14) ^ (d>>13) ^ (d>>12) ^ (d>> 9) ^ (d>> 8) ^ (d>> 6);
 D29 = (l>> 1) ^ (d>>27) ^ (d>>25) ^ (d>>24) ^ (d>>22)
     ^ (d>>21) ^ (d>>20) ^ (d>>19) ^ (d>>17) ^ (d>>15)
     ^ (d>>11) ^ (d>> 8) ^ (d>> 7) ^ (d>> 6);
 /* mask the parity values to a single bit */
 D24 &=1; D25 &=1; D26 &=1; D27 &=1; D28 &=1; D29 &=1;
 /* Test and reject */
 if(((d>>5) & 1) != D24) { return 0; }
 if(((d>>4) & 1) != D25) { return 0; }
 if(((d>>3) & 1) != D26) { return 0; } 
 if(((d>>2) & 1) != D27) { return 0; }
 if(((d>>1) & 1) != D28) { return 0; }
 if(((d>>0) & 1) != D29) { return 0; }
#endif
 return 1; /* Success! */
}

 

Correlation time! Come on! Lets correlate

(Apologies to K. C. and the Sunshine band).

The key maths behind receiving GPS is cross-correlation – the testing of two signals to see how alike they are, in frequency and phase. For a GPS receiver this technique is used to see if the pseudo-random bit pattern is being transmitted on a given frequency, and when required, measure the phase  difference between the local sample clock and the down-converted signal.

The process can be quite easily described:

  1. Stretch the 1023-bit Gold code of the Space Vehicle to match the number of samples per millisecond of the incoming data stream.
  2. Multiply the sample stream with the Gold code you are hunting for, treating the 0s in the stretched Gold code as -1s. If the Gold code is in the data stream, and they are aligned in phase, then this should remove the Gold code leaving just the intermediate frequency signal and any other noise.
  3. Multiply the resulting stream with the sin() and cos() functions of the intermediate frequency that you are probing, to give you an ‘I’ (in phase) and ‘Q’ (quadrature) stream.
  4. Sum up the I and Q streams for a sensible period of time (usually a multiple of length of time it takes to transmit a Gold code. These two numbers will indicate how much signal is in phase with the sin() and cos() functions used in step three. This can be treated as a 2D vector to measure phase
  5. Finally work out the magnitude of the vector – i^2+q^2 will give you the power level for that Gold code, with that alignment, at that frequency

This all sees like a lot of work – for each sample, there is a sin(), a cos(), three multiplies and two additions – tracking 5 space vehicles using a 5456kHz sample rate will require at least a hundred of megaflops.

Work avoidance scheme

A simple trick can save the day. For one-bit samples (where ‘1’ indicates a positive value, and ‘0’ indicates a negative value), a binary XOR is equivalent to multiplication by either 1 or -1 as it flips the bit.

This reduces the process to:

  1. Stretch the Gold code from 1023 bits to match the sample rate (in bits per millisecond)
  2. XOR the raw bitstream with the stretched Gold code.
  3. XOR the bitstream from step 2 with a binary string, which is ‘1’ where sine of the frequency being tested is positive (and ‘0’ where it is negative) to give the I stream. Also XOR the bitstream from step 2 with a binary string, which is ‘1’ where cosine is positive (and ‘0’ where it is negative) to give the Q stream.
  4. Count up the ‘1’s and ‘0’s in the I and Q streams for a sensible time period, then subtract the ‘1’s count from the ‘0’s count.
  5. Finally work out the magnitude of the vector – i^2+q^2 will give you the power of that Gold code, with that alignment, at that frequency

It should also be pretty self-evident step 1 (stretching the Gold code) only needs to be performed once, and that steps 2 and 3 can be performed in any convenient order. A less obvious optimization is that in step 4, you only need to count the ‘1’s, as you can always calculate how many zeros there must have been.

So, with all that explanation out of the way, here is my code that teases out GPS signals from the bit stream from a single bit ADC:

 cos_total = 0;
 sin_total = 0;

 for(i = 0; i < MAX_WINDOW_USED; i++) {
   unsigned int gc_bit = sv->gold_code_stretched[i];
   if(data[i])
     gc_bit ^= 1;

   if(gc_bit) {
     sin_total += intermediate_freq_sin_table[i];
     cos_total += intermediate_freq_cos_table[i];
   }
 }
 total = sin_total*sin_total+cos_total*cos_total;

So that is it – the key to detecting GPS signals in 14 lines of C!

 

What are Gold codes?

The key to how GPS works are “Gold codes” – a 1023-bit binary number that each Space Vehicle transmits once every millisecond to allow the presence and phase of the signal to be accurately detected.

Here is the Gold code for Space Vehicle 1:

110010000011100101001001111001010001001111101010110100010001010101
011001000111101001111110110111001101111100101010100001000000001110
101001000100110111100000111101011100110011110110000000101111001111
101010011000101101110001101111010100010101100000100000000100000011
000111011000000111000110111111111010011101001011011000010101011000
100111001011011101100011101110111100001101100001100100100100000110
110100101101111000101110000001010010011111100000101010111001111101
011111001100110001110001101101010101101100011011101110000000000010
110011011001110110100000101010111010111010010100011100111000100101
000101001011010000101011011010110110001110011110110010000111111001
011010001000011111010101110011001001001001011111111110000111110111
100011011100101100001110010101000010100101011111100011110110100111
011001111110111110100011000111110000000100101000101101000100010011
011000000111011010001101000100100011100010110011001001111001101111
110011001010011010011010111100110110101001110111100011010100010000
100010010011100001110010100010000

Testing the correlation of two binary numbers is pretty simple – XOR the two values together, and then count the numbers of 1s and 0s in the resulting string.

int correlate(unsigned char *a, unsigned char *b; int length) {
  for(i = 0; i < length; i++) {
    if(a[i] ^ b[i] == 0)
      match++
    else
      match--
  }
  return (match < 0) ? -match : match;
}

To prove how little they match, here are graphs of Space Vehicles 1’s Gold code correlated against itself and that of the Gold code for Space Vehicle 2:

correlation

The single peak can clearly be seen standing out from the noise, SV1’s code XORed with SV1’s code, with an phase offset of zero.

This special property of Gold codes allows a GPS receiver to check if a Space Vehicle’s signal is being received, and the relative phase to with sub-microsecond microsecond accuracy. However, this comes at a price. You can’t just ‘tune in’ to a Space Vehicle, you have to do a lot of work, deliberately looking for the correct single, at the correct frequency and with the correct phase.

Here is the C program used to produce the data for the graphs. You might find it useful:

/*****************************************************************************
* correlate.c - Testing Gold code correlation
*****************************************************************************/

unsigned char gold_sv1[1023] =
 "110010000011100101001001111001010001001111101010110100010001010101"
 "011001000111101001111110110111001101111100101010100001000000001110"
 "101001000100110111100000111101011100110011110110000000101111001111"
 "101010011000101101110001101111010100010101100000100000000100000011"
 "000111011000000111000110111111111010011101001011011000010101011000"
 "100111001011011101100011101110111100001101100001100100100100000110"
 "110100101101111000101110000001010010011111100000101010111001111101"
 "011111001100110001110001101101010101101100011011101110000000000010"
 "110011011001110110100000101010111010111010010100011100111000100101"
 "000101001011010000101011011010110110001110011110110010000111111001"
 "011010001000011111010101110011001001001001011111111110000111110111"
 "100011011100101100001110010101000010100101011111100011110110100111"
 "011001111110111110100011000111110000000100101000101101000100010011"
 "011000000111011010001101000100100011100010110011001001111001101111"
 "110011001010011010011010111100110110101001110111100011010100010000"
 "100010010011100001110010100010000";

unsigned char gold_sv2[1023] =
 "111001000011100000111110100110010110111111001011001011111111010010"
 "110000100010001011000111100011000001100111000010001110001000111010"
 "011110000010110100001101001011110000111000001111100011011100110001"
 "101011100000001111000110100010111001000110011001101101001100000101"
 "111000100100101010001100000010101010000010011110000010010111111111"
 "111011101010110101010000010101001010101010000100100011001000010100"
 "001001010011011100000011100101010100110111000001111011101100000010"
 "001000101101011100110101001001100100111011011111010110000010010111"
 "100111000010110000110100011010110000001001011101100010100101101110"
 "100101101100111110000010001000101011010010100011101001110110010010"
 "001111101011111011110110001110101110101001101111001101010010111010"
 "111100000111101111101101111010001111100110110111100010000101101011"
 "010000111100110010110011000001100001000011111101011001111101001100"
 "111010100010111101000011001111101000001001111010011001011111011100"
 "101101101001111011001110011110010110111001100110101111101111011010"
 "010100101101000000101000011001000";

/****************************************************************************/
int main(int argc, char *argv[]) {
 int offset;
 for(offset = 0; offset < 1023; offset++)
 {
  int i,match;

  match = 0;
  for(i = 0; i < 1023; i++) {
   if(gold_sv1[i%1023] == gold_sv2[(i+offset)%1023])
    match++;
   else
    match--;
  }
  printf("%4i: %i\n",offset,match);
 }
}
/*****************************************************************************
* End of file!
*****************************************************************************/

 

 

 

Why GPS Demystified?

GPS is now a key technology for many diverse fields, from controlling autonomous vehicles, location aware games on cellphones to controlling continental power grids. I want to understand it.

One unique feature of GPS is the ability for it to provide a highly accurate time reference, and it was thinking about how this feature was implemented that started me wondering how GPS actually works. How can a small electronic device determine an accurate time, based on the transmission of satellites that are 20,200 km away?

After investigating on Google, I’ve reached the conclusion that one way to understand everything about it is to make my own functional GPS receiver.

In a act of serendipity I am just about to receive the hardware for this project – the KiwiSDR cape for the Beagle Bone. It handles the messy RF front-end, delivers the output of a Skyworks single chip GPS front-end to an Xilinx Spartan 6 FPGA for decoding. You can find the full details on the project’s web site – http://kiwisdr.com/kiwisdr/ including the full set of design files.

For this blog I will focus on the most primitive, understandable techniques I can find to decode this raw ADC output into usable navigation information. Rather than efficiency the focus will be on being understandable and testable with the minimal advanced maths skills I have.