Receiver Operation

The image below shows the computer screen when the receiver first starts up.  Note the channels are listed from 0 to 11, since the position and time were known it is attempting a warm start.  It has predicted 9 satellites in view and has determined their azimuth, elevation and doppler shift.  Note the state is 1 so the receiver is in the very beginning stages of acquisition.  The other pages in figures 2 through 4 can be reached by pushing the "P(p)" key.

Figure 1 

Default page 1 Computer Display

Page 2 shows the TLM, TOW and health of each satellite tracked.  The TOW is especially useful to compare to what the software thinks it is (in top left corner).

Figure 2 

Page 2 Computer Display

Page 3 shows the ionosphere and troposphere corrections being used for each satellite.

Figure 3 

Page 3 Computer Display

Page 4 shows the pseudorange and delta pseudorange from each satellite.

Figure 4

Page 4 Computer display

As seen in the diagram below 4 correlation counts are available from each channel.  The prompt and dither (either ahead or behind the prompt correlator by 1/2 chip) both have in-phase "I" and quadrature "Q" correlators.  By forming an RSS (root sum square) of the in-phase and quadrature values we can determine if we are lined up with the PRN code.  The in-phase and quadrature correlator values tell us what the phase of the signal is.  Although the source code includes comments I think it will be very helpful to explain in more detail how the receiver software works.

Figure 2 Channel Digital Processing Schematic

To track the signal we do not directly move the correlators, we speed up or slow down the carrier and code DCO's (digitally controlled oscillators) to keep the channel locked in code and carrier phase.

The correlator read timing of the receiver is set by taking over the PC's Int 0 which is normally used to keep the computer real time clock up to date.  The 8xxx normally interrupts the computer every xxx ms.  The software replaces this interrupt routine with its own (GPS_interrupt) and sets the interrupt time to about 500 us.  Whenever the interrupt occurs the channel data register is checked to see which channels have dumped correlation data.  The correlator data is read and a check is made to see if a 99.9999 ms "tic" has occurred.  If it has, the measurement data is stored.

This state diagram (figure 3) illustrates how the receiver works.

Figure 3 State Diagram

State 1 is the acquisition state.  It does the code and frequency (doppler) search to find a high correlation peak.  If a peak in either the prompt or delay correlator (RSS of in-phase and quadrature) is above threshold it switches to state 2.  When in warm or cold start the doppler is centered on the expected value.  The search is conducted with two loops, the inner loop is the code search and the outer is the frequency search.  The code search pattern is simply "slewing" the PRN generator by adding an extra chip of delay after each correlation sample.  Thus after 1023 samples we have searched every 1/2 chip in code space.  The frequency search pattern checks on either side of the expected frequency in integer multiples of the doppler search range i.e. ( 0, +1, -1, +2, -2, ......)*200 Hz.  The maximum range depends on what the receiver is doing.  The nominal value is for a warm start, it is larger if attempting a cold start and smaller if the receiver has a navigation fix.  Samples are taken every millisecond.  This allows each frequency bin to be searched in about 1 second (if it does not go into state 2).

State 2 is the confirmation state.  The confirmation state stops the search and dwells at the code and doppler where the high correlation peak was found in state 1 to confirm the presence of the signal in order to reduce the false alarm rate.   If n of m samples (e.g. 8 of 10) are above threshold it switches to state 3.

State 3 is the pull-in state.  The signal has been confirmed but the frequency may be up to ~500 Hz off.  The pull-in state attempts to start tracking the signal in order to pull the frequency in close enough that the carrier phase can be tracked.  This state is enabled for about 1500 ms, during the last 500 ms of this time the C/No and phase errors are measured and it attempts to synchronize on the edge of a data bit.  If it confirms carrier and code tracking it switches to state 4.  Since the frequency is likely to be very far off the receiver uses a combination of a frequency locked loop FLL and a phase locked loop PLL.  As the time progresses into the state the effect of the FLL is decreased so that at the end of the state the PLL is dominant.  The figures below are actual data from a pull-in attempt and illustrate how it works.  In Figure 4 it can be seen that the phase is changing rapidly at the start.  As time progresses to about 200 ms the phase is changing more slowly.  By the time we get to about 300 ms we are close to phase lock.  From this point on we can see the data bits as +90 and -90 degree phase transitions.  

Figure 4  Carrier Phase during Pull-In

In figure 5 we can see the noise in the carrier tracking loop.  It starts out high (due to the FLL) and transitions down in three steps every 500 ms.  At the end the effect of the FLL is negligible and the PLL dominates.

Figure 5 Carrier loop frequency command during Pull-in

In addition note that in order to reduce start-up transients the code loop is closed after 2 ms and the carrier loop is closed after 5 ms.

Pull-in Carrier Loop code

q_sum        = prompt_Q + dither_Q

i_sum         = prompt_I + dither_I

theta           = fix_atan2(q_sum,-i_sum)

theta_e        = theta - sign(theta)*25736

theta_dot    = theta - old_theta

wdot_gain   = (ms/499)4

dcarr           = pull_carr_k*(theta_dot*5)/(1 + wdot_gain)+theta_e

ddcar          = (dcarr - dcarr1)*pull_carr_d

carrier_freq = (dcarr + ddcar)/16384 + carrier_freq

dcarr1         = dcarr

old_theta     = theta


Pull-in Code loop code

prompt_mag  = RSS(prompt_I, prompt_Q)

dith_mag       = RSS(dither_I, dither_Q)

dfreq             = pull_code_k*(prompt_mag-dith_mag)*16384/(prompt_mag+dith_mag)

ddf                = (dfreq - dfreq1)*pull_code_d

code_freq     = (dfreq + ddf)*16384+code_freq

dfreq1           = dfreq

State 4 is the normal tracking state.  The tracking loops integrate over a data bit (20 ms) to track code and 1ms to track phase.  The need to track phase at 1 ms may have been due to excessive phase noise from the previous clock I had used.  I have not yet tried to track carrier at 20ms on the new receiver.  The data message is recorded and the time is synchronized to the TOW (time of week) of the data message.

Track Carrier Loop code

q_sum        = prompt_Q + dither_Q

i_sum         = prompt_I + dither_I

dcarr           = i_sum*16384*trk_carr_k*sign(q_sum)/RSS(q_sum,i_sum)

ddcar          = (dcarr - dcarr1)*trk_carr_d

carrier_freq = (dcarr + ddcar)/16384 + carrier_freq

dcarr1         = dcarr

Track Code loop code

prompt_mag  = RSS(prompt_I, prompt_Q)

dith_mag       = RSS(dither_I, dither_Q)

dfreq             = trk_code_k*(prompt_mag-dith_mag)

ddf                = (dfreq - dfreq1)*trk_code_d

code_freq     = (dfreq + ddf)/trk_div+code_freq

dfreq1           = dfreq

While in any of states 2 through 4 if a problem occurs the channel transitions back to state 1.

If the C/No drops below about 35 dBHz and the satellite is below the mask angle the channel is deallocated.

Once every minute the receiver checks to see if a new satellite has gone above the mask angle, assigns a channel to it, and starts it out in state 1.

Computing Position

Recall we are solving this equation:

(HtH)-1Hr = dx

This is a condensed version of the code used

for (i=1;i<=nsl;i++)  // nsl = number of satellites with pseudorange data
// Compute range in ECI at the time of arrival at the receiver
alpha = (dt[i] - t)*omegae;  // rotate to ECI at time of tic

                                      // omegae = rotation rate of earth in rad/sec

r = sqrt(

         pow(track_sat[i].x*cos(alpha) - track_sat[i].y*sin(alpha) - x,2)+
         pow(track_sat[i].y*cos(alpha) + track_sat[i].x*sin(alpha) - y,2)+
         pow(track_sat[i].z - z,2)); // compute range

bm[i] = r - (dt[i] - t)*c;  // bm = difference of range from pseudorange

ms[1][i] = (track_sat[i].x*cos(alpha) - track_sat[i].y*sin(alpha) - x)/r;
ms[2][i] = (track_sat[i].y*cos(alpha) + track_sat[i].x*sin(alpha) - y)/r;
ms[3][i] = (track_sat[i].z - z)/r;
ms[4][i] = 1.0;     // ms = row vector of direction cosines of H

                        // Note in this software version G is a I matrix and is not used

for (k=1;k<=nsl;k++)
a1 += ms[1][k]*ms[1][k];
b1 += ms[1][k]*ms[2][k];
c1 += ms[1][k]*ms[3][k];
d1 += ms[1][k]*ms[4][k];
// e1 += ms[2][k]*ms[1][k];
f1 += ms[2][k]*ms[2][k];                   // compute HtH  the result is arranged as:
g1 += ms[2][k]*ms[3][k];                  
h1 += ms[2][k]*ms[4][k];                  //  a1   b1  c1  d1 \\
// i1 += ms[3][k]*ms[1][k];               //   e1   f1   g1  h1   \\
// j1 += ms[3][k]*ms[2][k];               \\   i1    j1   k1  l1   //
k1 += ms[3][k]*ms[3][k];                  \\  m1  n1  o1  p1 //
l1 += ms[3][k]*ms[4][k];
// m1 += ms[1][k];
// n1 += ms[2][k];
// o1 += ms[3][k];
p1 += ms[4][k];
o1=l1;m1=d1;n1=h1;e1=b1;i1=c1;j1=g1;  // matrix symmetry


Not that it's a great idea to do it this way but since it is only a 4x4 matrix and we don't have to worry about using someone's math functions

denom = (k1*p1-l1*o1)*(a1*f1-b1*e1)+(l1*n1-j1*p1)*(a1*g1-c1*e1)+(j1*o1-k1*n1)*(a1*h1-d1*e1)
dd[1][1] = f1*(k1*p1-l1*o1)+g1*(l1*n1-j1*p1)+h1*(j1*o1-k1*n1);
dd[1][2] = e1*(l1*o1-k1*p1)+g1*(i1*p1-l1*m1)+h1*(k1*m1-i1*o1);
dd[1][3] = e1*(j1*p1-n1*l1)-i1*(f1*p1-n1*h1)+m1*(f1*l1-j1*h1);
dd[1][4] = e1*(n1*k1-j1*o1)+i1*(f1*o1-n1*g1)+m1*(j1*g1-f1*k1);
dd[2][1] = b1*(l1*o1-k1*p1)+j1*(c1*p1-d1*o1)+n1*(d1*k1-c1*l1);
dd[2][2] = a1*(k1*p1-l1*o1)+c1*(l1*m1-i1*p1)+d1*(i1*o1-k1*m1);
dd[2][3] = a1*(l1*n1-j1*p1)+i1*(b1*p1-n1*d1)+m1*(j1*d1-b1*l1);
dd[2][4] = a1*(j1*o1-n1*k1)-i1*(b1*o1-n1*c1)+m1*(b1*k1-c1*j1);
dd[3][1] = b1*(g1*p1-h1*o1)-f1*(c1*p1-o1*d1)+n1*(c1*h1-d1*g1);
dd[3][2] = a1*(o1*h1-g1*p1)+e1*(c1*p1-o1*d1)+m1*(d1*g1-c1*h1);
dd[3][3] = a1*(f1*p1-h1*n1)+b1*(h1*m1-e1*p1)+d1*(e1*n1-f1*m1);
dd[3][4] = a1*(n1*g1-f1*o1)+e1*(b1*o1-c1*n1)+m1*(c1*f1-b1*g1);
dd[4][1] = b1*(h1*k1-g1*l1)+f1*(c1*l1-d1*k1)+j1*(d1*g1-c1*h1);
dd[4][2] = a1*(g1*l1-h1*k1)-e1*(c1*l1-d1*k1)+i1*(c1*h1-d1*g1);
dd[4][3] = a1*(j1*h1-f1*l1)+e1*(b1*l1-d1*j1)+i1*(d1*f1-b1*h1);
dd[4][4] = a1*(f1*k1-g1*j1)+b1*(g1*i1-e1*k1)+c1*(e1*j1-f1*i1);
if ( denom <= 0.0 )
    result.x = 0.0; // something went wrong
    result.y = 0.0; // set solution to center of earth
    result.z = 0.0;
    result.dt = 0.0;
for (i = 1;i <= 4;i++)
     for (j = 1;j <= 4;j++) dd[i][j] = dd[i][j]/denom;  // dd is now the inverse
for (i = 1;i <= nsl;i++)
   for (j = 1;j <= 4;j++)
       pm[j][i] = 0.0;
       for (k = 1;k <= 4;k++) pm[j][i] += dd[j][k]*ms[k][i];  // pm = (HtH)-1H
for (i=1;i<=4;i++)
    br[i] = 0.0;
    for (j = 1;j <= nsl;j++) br[i] += bm[j]*pm[i][j];  // br = (HtH)-1H dr
x = x + br[1];  // add the correction vector to the position and time vector
y = y + br[2];
z = z + br[3];
t = t - br[4]/c;
correct_mag = sqrt(br[1]*br[1] + br[2]*br[2] + br[3]*br[3]);
} while ( correct_mag > 0.01 && correct_mag < 1.e8 && nits < 10);
result.x = x;
result.y = y;
result.z = z;
result.dt = t         // result is in ecef

I have found that it generally converges in 3-4 iterations if started at the center of the earth and 2 iterations if a recent solution is used for the starting point.

Computing Velocity

Recall we are solving this equation:

(HtH)-1Hdf = v

This is a condensed version of the code used

for (i=1;i<=nsl;i++)
   alpha = (dt[i] - t)*omegae; // rotate to ECI at time of tic
   r = sqrt(

       pow(track_sat[i].x*cos(alpha) -  track_sat[i].y*sin(alpha) - x,2)+
       pow(track_sat[i].y*cos(alpha) + track_sat[i].x*sin(alpha) - y,2)+
       pow(track_sat[i].z - z,2));

   bm[i] = ((track_sat[i].x*cos(alpha)-track_sat[i].y*sin(alpha)-x)*d_sat[i].x+
   (track_sat[i].z-z)*d_sat[i].z)/r-meas_dop[i]*lambda;// df in meters  lambda is L1

                                                                                  // wavelength of about 19 cm

for (i=1;i<=4;i++)
   for (j = 1;j <= nsl;j++) br[i] += bm[j]*pm[i][j];  //  br = (HtH)-1H df
result.xv = br[1] + y*omegae; // get rid of earth
result.yv = br[2] -  x*omegae; // rotation rate
result.zv = br[3];
result.df = br[4]/c*1000000.0;  // clock frequency offset in ppm (parts per million)
return(result);  // result is in ecef

Last Updated September 6, 2003