**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)^{-1}Hr = 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

}

a1=0.;b1=0.;c1=0.;d1=0.;

e1=0.;f1=0.;g1=0.;h1=0.;

i1=0.;j1=0.;k1=0.;l1=0.;

m1=0.;n1=0.;o1=0.;p1=0.;

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

/*

SOLVE FOR THE MATRIX INVERSE

*/

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)

+(l1*m1-i1*p1)*(c1*f1-b1*g1)+(i1*o1-k1*m1)*(d1*f1-b1*h1)+(i1*n1-j1*m1)*(c1*h1-d1*g1);

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;

}

else

{

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)^{-1}H

}

}

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)^{-1}H dr

}

nits++;

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)^{-1}Hdf = 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].y*cos(alpha)+track_sat[i].x*sin(alpha)-y)*d_sat[i].y+

(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++)

{

br[i]=0.0;

for (j = 1;j <= nsl;j++) br[i] += bm[j]*pm[i][j]; //
br = (HtH)^{-1}H 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