The GPS PRN (Gold Codes)

First read the previous post about the basics of GPS encoding.

PRN (“Psuedo Random Noise”) sequences are used as part of the CDMA scheme in GPS. But they don’t just use any random number. In fact, it’s not random, it’s a kind of pseudo-random number. One of the reasons you want a spreading code that looks random is to ensure you have good autocorrelation properties. It also helps the spreading properties.

A signal that’s statistically random doesn’t repeat itself at any point in its sequence and so it has very good autocorrelation properties. We also want the different codes for each satellite that are nearly orthogonal so the cross correlation of one code will be very high even when many satellite codes are mixed in. This is the scheme used to create pseudo-random sequences with all these properties for GPS:

Shift Registers

The Gold Codes are generated using a pair of shift registers with feedback. Let’s make a software emulation of a shift register in python! A shift register is just a list of bits that has an input end and an output end. The GPS registers are length 10. We’ll start them out with all 1’s, calculate the feedback, and then shift all the numbers over to the right. Taking it step by step:

# Make a shift register of length 10
# initilize with all 1's
G1 = [1 for i in xrange(10)]
print G1

[1, 1, 1, 1, 1, 1, 1, 1, 1, 1]

Feedback

Now we take feedback from a couple of positions in the register and modulo 2 add them together. Specifically for GPS, the first shift register (G1) is defined to have feedback from positions 3 and 10 (note the 1 indexing):

fb3 = G1[3-1]           # position 3
fb10 = G1[10-1]         # position 10
fb = (fb3 + fb10) % 2   # modulo 2 add
print fb

0

Now we have the number to push back into the fist position. Let’s “clock” the shift register and put the feedback in:

# shift everything one place to the right
for i in reversed(range(len(G1[1:]))):
    G1[i+1] = G1[i]

# feedback
G1[0] = fb

print G1

[0, 1, 1, 1, 1, 1, 1, 1, 1, 1]

Neat! Now just keep going…

fb = (G1[3-1] + G1[10-1]) % 2
for i in reversed(range(len(G1[1:]))):
    G1[i+1] = G1[i]
G1[0] = fb
print G1

[0, 0, 1, 1, 1, 1, 1, 1, 1, 1]

Etc., etc., etc. Lets put this in a function and clean things up.

def shift(register, feedback, output):
    """GPS Shift Register
    
    :param list feedback: which positions to use as feedback (1 indexed)
    :param list output: which positions are output (1 indexed)
    :returns output of shift register:
    
    """
    
    # calculate output
    out = [register[i-1] for i in output]
    if len(out) > 1:
        out = sum(out) % 2
    else:
        out = out[0]
        
    # modulo 2 add feedback
    fb = sum([register[i-1] for i in feedback]) % 2
    
    # shift to the right
    for i in reversed(range(len(register[1:]))):
        register[i+1] = register[i]
        
    # put feedback in position 1
    register[0] = fb
    
    return out

# example:
print shift(G1, [3,10], [10])
print G1

1 [0, 0, 0, 1, 1, 1, 1, 1, 1, 1]

Output

Notice we have a place defined to have some output. The shift registers can be tapped at different locations to define an output. These spots are modulo 2 added together to make an output from the register. This function returns that value.

Since we defined this generically we can use it for the other shift register used for the C/A code, G2, which has feedback from positions 2,3,6,8,9, and 10:

G2 = [1 for i in range(10)]
shift(G2, [2,3,6,8,9,10], [10])
print G2

[0, 1, 1, 1, 1, 1, 1, 1, 1, 1]

The C/A Code

The “coarse acquisition” code, or C/A code, in GPS is made of two shift registers, G1 and G2. G1 has feedback from position 3 and 10, and G2 has feedback from 2,3,6,8,9, and 10. The output of each is modulo 2 added to produce the code. To make different satellite codes different spots on the G2 register are tapped.

To produce a bit for PRN #1 we tap position 2 and 6 as the output (according to the GPS spec). We modulo 2 add this to the output of G1 (always position 10).

Lets make the first 10 bits of PRN #1:

# init registers
G1 = [1 for i in range(10)]
G2 = [1 for i in range(10)]

# work out first 10 bits
ca = []
for i in range(10):
    g1 = shift(G1, [3,10], [10]) #feedback 3,10, output 10
    g2 = shift(G2, [2,3,6,8,9,10], [2,6]) #feedback 2,3,6,8,9,10, output 2,6 for sat 1
    ca.append((g1 + g2) % 2)

print ca

[1, 1, 0, 0, 1, 0, 0, 0, 0, 0]

Visualization

To make things a little clearer, we can draw a diagram of this shift register feedback and output system.

Feedback lines are in blue. Output (in red) from G1 is always from position 10. Different choices for the two G2 outputs are assigned to different satellites. Everything is modulo 2 added. In this example to current C/A output is a 1:

from IPython.display import SVG
SVG(filename='images/gps_prn.svg')

Satellite Codes

Each satellite is assigned a different 1023 bit long PRN (a Gold Code). This is decided by choosing different spots to tap the G2 shift register. In the above example we are looking at PRN #1, which is defined as being tapped from position 2 and 6. To make PRN #2 we would tap 3 and 7. Here is the list of PRN assignments from the GPS spec:

SV = {
   1: [2,6],
   2: [3,7],
   3: [4,8],
   4: [5,9],
   5: [1,9],
   6: [2,10],
   7: [1,8],
   8: [2,9],
   9: [3,10],
  10: [2,3],
  11: [3,4],
  12: [5,6],
  13: [6,7],
  14: [7,8],
  15: [8,9],
  16: [9,10],
  17: [1,4],
  18: [2,5],
  19: [3,6],
  20: [4,7],
  21: [5,8],
  22: [6,9],
  23: [1,3],
  24: [4,6],
  25: [5,7],
  26: [6,8],
  27: [7,9],
  28: [8,10],
  29: [1,6],
  30: [2,7],
  31: [3,8],
  32: [4,9],
}

Final PRN Code

So now we have all the pieces to write our PRN creator! Given the input of satellite number (SV) we can now compute the 1023 bit long chip for that satellite:

def PRN(sv):
    """Build the CA code (PRN) for a given satellite ID
    
    :param int sv: satellite code (1-32)
    :returns list: ca code for chosen satellite
    
    """
    
    # init registers
    G1 = [1 for i in range(10)]
    G2 = [1 for i in range(10)]

    ca = [] # stuff output in here
    
    # create sequence
    for i in xrange(1023):
        g1 = shift(G1, [3,10], [10])
        g2 = shift(G2, [2,3,6,8,9,10], SV[sv]) # <- sat chosen here from table
        
        # modulo 2 add and append to the code
        ca.append((g1 + g2) % 2)

    # return C/A code!
    return ca

And now test with, say, satellite #24:

%matplotlib inline
import chart
from matplotlib import rcParams
rcParams.update({'font.size': 12})

f_prn = 10.23e6 / 10  # chipping frequency

# find ca code for sat 24, and make 0 into -1 to use in BPSK
sat_24 = [-1 if x==0 else x for x in PRN(24)]
prn = lambda x: sat_24[int(x*f_prn)%1023]

chart.Plot('GPS Sat 24 PRN Code', tend=1023/f_prn, sr=1.0e-7).signal(prn, lw=0.7).show()

See previous post and the next post to continue learning about GPS sinals.