We want to find a satellite in our test data. For each satellite in the GPS constellation we can "search" for the hidden GPS signal by trying a correlation technique. We know exactly what signal we want to find is (the PRN of a given sat) so we can "ask" if it's there by looking for a peak in a cross correlation graph. There is a trick though: it only works if you run it at the exact right frequency. We can't know this ahead of time though, because the recieved signal may (in fact, will!) be doppler shifted by some amount. So we try at each doppler shift in some range.
Lets start by reading in data:
In [1]:
from gps import data
from gps import aquire
from numpy import argmax, absolute, divide
# read in 20 milliseconds of data assuming the spain data format
init = data.read('../../data/2013_04_04_GNSS_SIGNAL_at_CTTC_SPAIN.dat', data.SPAIN, 0.020)
print "read in", len(init), "samples"
In [2]:
# It turns out that we can choose an optimal step size
test_step = int(data.SPAIN['sample-rate'] / float(len(init)))
test_range = range(-15000, 15000, test_step)
# Stuff tests in here:
frequencies = []
max_snr = 0
min_snr = 1000
for doppler in test_range:
correlation = aquire.test(data.SPAIN['sample-rate'], init, 23, doppler)
# Compute stats about this test
tot_pwr = 0
snr = [0 for i in range(1023)]
for i, pwr in enumerate(correlation):
# bin by code phase
phase = int(((i/4.0e6)*1.023e6) % 1023)
tot_pwr += pwr
snr[phase] += pwr
# normilize to total power
snr_norm = divide(snr, tot_pwr)
max_norm_snr = max(snr_norm)
min_norm_snr = min(snr_norm)
if max_snr < max_norm_snr:
max_snr = max_norm_snr
if min_snr > min_norm_snr:
min_snr = min_norm_snr
frequencies.append(snr_norm)
Now we can look a the results.
In [3]:
%matplotlib inline
import matplotlib.pyplot as plt
from matplotlib import rcParams
rcParams.update({'font.size': 12})
yticks = range(0, len(frequencies)+1, 100)
ytick_labels = ["%d kHz" % ((i-len(frequencies)/2)/20) for i in yticks]
fig = plt.figure(figsize=(11,11))
ax = fig.add_subplot(111)
img = ax.imshow(frequencies, interpolation='nearest', aspect=1.6, cmap='hot')
plt.yticks(yticks, ytick_labels)
plt.xticks(range(0, 1025, 128))
plt.ylabel("Doppler")
plt.xlabel("Code Phase [chips]")
plt.title("PRN 23 Search Space")
ax.tick_params(axis='x', colors='white', labelcolor='black', length=10)
ax.tick_params(axis='y', colors='white', labelcolor='black', length=10)
cbar = plt.colorbar(img, ticks=[min_snr, (min_snr+max_snr)/2.0, max_snr], fraction=0.043, pad=0.04)
cbar.ax.set_ylabel('Correlation Power')
cbar.ax.set_yticklabels(["Low", "Medium", "High"])
plt.annotate("Found it!", color='#eecccc', xy=[617, 510], xycoords='data', xytext=[-10,50],
textcoords='offset points', arrowprops=dict(color='#eecccc',arrowstyle="->",
connectionstyle="arc3,rad=%0.1f" % 0.2) )
plt.show()
Looking at just the area near our suspected correlation peak:
In [4]:
yticks = range(518-40, 518+41, 20)
ytick_labels = ["%d kHz" % ((i-len(frequencies)/2)/20) for i in yticks]
fig = plt.figure(figsize=(11,11))
ax = fig.add_subplot(111)
img = ax.imshow(frequencies, interpolation='nearest', aspect=1, cmap='hot')
plt.yticks(yticks, ytick_labels)
plt.xticks(range(617-40, 617+41, 16))
plt.ylabel("Doppler")
plt.xlabel("Code Phase [chips]")
plt.title("PRN 23 Search, Detail")
ax.set_xlim([617-40, 617+40])
ax.set_ylim([518-40, 518+40])
ax.tick_params(axis='x', colors='white', labelcolor='black', length=10)
ax.tick_params(axis='y', colors='white', labelcolor='black', length=10)
cbar = plt.colorbar(img, ticks=[min_snr, (min_snr+max_snr)/2.0, max_snr], fraction=0.045, pad=0.04)
cbar.ax.set_ylabel('Correlation Power')
cbar.ax.set_yticklabels(["Low", "Medium", "High"])
plt.show()
Neat.