Tuesday, January 7, 2014

Finally I can see GPS-satellites with bladeRF!

My first styled page

Summary

Happy New Year to all of you, or better perhaps, I wish you all a Better New Year!

In my struggle to decode a bladeRF-datafile with GPS-data I rewrote a working java-program into python. When I was able to decode a previous old data-file and got the same results as with my previous old java-program I felt confident to go back to my dear bladeRF again.

Luck was with me. To escape the Dutch newyear-firework me and my wife and our dog Boris went for a long fortnight in our caravan somewhere in the north of the Netherlands in the middle of nowhere.  It was January 3rd when I put my active GPS-antenna in front of the window together with my Garmin GPS. My Garmin told me to see satellites 5,7,13,16 and 23. I made a recording into a few data-files and managed to have my python-program really consistently find those satellites! My bladeRF got DC-current from a separate adaptor, an old fashioned one with a transformer happily delivering current at 5 V. 

I got some problems that I will describe here together with my python program.


# findPRNs05.py # Jan-2014 Kees de Groot # # find satellites in bladeRf-file # search for prn's in a big loop over all possible frequencies: # loop over all prn's: # calculate list with prn[i] # fft prn[i] # multiply with fft'd conjugated signal # inverse-fft # search for maximum # print i, maximum import pylab as pl import numpy as np import math import sys import goldutil as gl # home-made Gold-sequence generator ###################################################################### ############# parameter section ###################################### ###################################################################### path = '/Temp/' filename = 'DKIEL.csv' centerfreq = 0 fsample = 4E6 timestep = 1.0 / fsample f = open(path + filename, 'rb') print "filename = ", path + filename print "fsample = ", fsample print "centerfreq = ", centerfreq # get 1 msec data n = int(1.0E-3 / timestep) print "1 msec data is ", n, " samples" skip = 10000 # there is a discontinuity in the file # so skip half of it print "skip ", skip, "samples of input-data from datafile" ###################################################################### ############# routines/functions section ############################ ###################################################################### def getPrnTimed(pprn, t): "return prn-chip at time t" # there are 1023 chips in 1 msec # normalize t while t >= 1E-3: t -= 1E-3 index = int((t / 1E-3) * 1023.0) return pprn[index] def getFFTPrn(prn): myGold = gl.Gold(prn) pprn = myGold.prnReg() samprn = [] t = 0 for i in range(n): # I have n data-samples samprn.append(getPrnTimed(pprn,t)) t += timestep fftprn = np.fft.fft(samprn) return fftprn def findPeak(timeresult): index = 0 # find peak in inverse-fft max = abs(timeresult[0]) for i in range(len(timeresult)): if abs(timeresult[i]) > max: max = int(abs(timeresult[i])) index = i res = [] res.append(int(max)) res.append(index) return res # [max, index] def findMax(peaklist): # find maximum in peaklist i = 0 max = peaklist[0][0] index = 0 prn = 0 for list in peaklist: i += 1 if list[0] > max: max = list[0] index = list[1] prn = i res = [] res.append(max) res.append(index) res.append(prn) return res # [max, index, prn] def getFFTdata(freq, dataI, dataQ, timestep): # demodulate by multiplying with freq # use complex mixer t = 0 w = 2 * math.pi* freq # omega dataCC = [] for i in range(n): valI = dataI[i] valQ = dataQ[i] wt = w * t factc = math.cos(wt) facts = math.sin(wt) valI = valI * factc - valQ * facts valQ = valI * facts + valQ * factc dataCC.append(complex(valI, valQ)) t += timestep sp = np.fft.fft(dataCC) # conjugate sp for i in range(len(sp)): sp[i] = sp[i].conjugate() return sp ########################################### ############# main loop ################## ########################################### dataI = [] dataQ = [] #dataC = [] averI = 0.0 averQ = 0.0 # read (I,Q) from bladeRF file i = 0 for line in f: skip = skip - 1 if skip > 0: continue list = line.split(',') Q = int(list[0]) averQ += Q I = int(list[1]) averI += I dataI.append(I) dataQ.append(Q) i += 1 if i >= n: break print "read ", n, " samples from file" print "averI = ", averI print "averQ = ", averQ averI = averI/n averQ = averQ/n print "averI/n = ", averI print "averQ/n = ", averQ ##pl.figure(1) ##pl.title("raw I/Q") ##pl.plot(dataI, dataQ, 'ro') ##pl.figure(2) ##pl.title("dataI") ##pl.plot(dataI) ##pl.figure(3) ##pl.title("dataQ") ##pl.plot(dataQ) ##pl.show() ##sys.exit('debug') # remove DC-component dataC = [] for i in range(len(dataI)): dataI[i] -= averI dataQ[i] -= averQ # dataC.append(complex(dataI[i], dataQ[i])) ##pl.figure(1) ##pl.plot(dataI, dataQ, 'ro') ##pl.show() # sys.exit('debug') #test = np.fft.fft(dataI) maxmaxmax = [] for freq in range(centerfreq - 30*500, centerfreq + 30*500, 500): peaklist = [] # contains peak-values [max, codeindex] # get fft'd data in sp sp = getFFTdata(freq, dataI, dataQ, timestep) # conjugated fft-data for prn in range(1,33): # prn = 1 .. 32 fftprn = getFFTPrn(prn) # get array with fft'd prn # get convolution in time-domain by multiplying in freq domain freqresult = [] for i in range(len(fftprn)): freqresult.append(fftprn[i]*sp[i]) timeresult = np.fft.ifft(freqresult) # find peak in timeresult # findPeak returns (max, codeindex) maxcode = findPeak(timeresult) peaklist.append(maxcode) # peaklist[i] = (max, codeindex) at prn = i+1 elem = findMax(peaklist) # find prn with largest peak # elem = [max, codeindex, prn] print "freq, [peak, index, prn] = ", freq, elem elem.append(freq) maxmaxmax.append(elem) # list of [peak, index, prn, freq] maxmaxmax.sort() maxmaxmax.reverse() # max at the top of the list print print "sorted list, max at top of the list" print "[peak, index, prn, freq]" aver = 0 num = 0 for elem in maxmaxmax: print elem aver += elem[0] num += 1 aver = aver / num print print " bigger than 1.5 * average: " print "[peak, index, prn, freq[" for elem in maxmaxmax: if elem[0] > 1.5 * aver: print elem #sys.exit('debug')

Description of the program


The program has an outer loop over the frequencies from -30*500 Hz to +30*500Hz.
Why? Depending on the speed of the GPS-receiver (e.g. in a car) and the speed of the GPS-satellites there is a certain Doppler-shift of the received frequency. So I try every possible frequency by mixing the received data in a complex mixer. The result is a complex signal that has been shifted in frequency.
That signal has to be convoluted with every possible prn-sequence (about 32 ot them). Convolution in the time domain can be achieved by multiplying in the frequency domain.
So, I fft the complex mixed signal, conjugate it and multiply with the fft-ed prn-sequence. After inverse-fft I look for a peak. The peak indicates correlation and at the same time indicate the phase of the prn.
Complicated? Yes, it took me some time to get used to that idea.

DC component

The signal in the bladeRF-datafile has a huge DC-component. I presume that the DA-converters have an output-range from 0 V to, perhaps, +3.3 V. I did not check that in the schematics yet.
The DC-component results in a big peak at 0 Hz, ruining my scaling for the plots. So I write a simple loop to calculate the average for I and Q and in a second loop I subtract the average from the data. This can also be done in one loop but I prefer to make programs simple at first. Later I might change the code if necessary.

Discontinuity

The first 500 samples or so, sometimes show a discontinuity. I solve that problem by simply skipping the first 5000 samples. Later I might go into that problem but for now it is solved.

Commands at bladerf-cli:



[WARNING] extract_field: Field checksum mismatch [WARNING] Could not extract VCTCXO trim value [WARNING] extract_field: Field checksum mismatch [WARNING] Could not extract FPGA size bladeRF> load fpga hostedx115.rbf Loading fpga from hostedx115.rbf... Done. bladeRF> print frequency RX Frequency: 1000000000Hz TX Frequency: 1000000000Hz bladeRF> set frequency 1575420000 Set RX frequency: 1575420000Hz Set TX frequency: 1575420000Hz bladeRF> set samplerate 4000000 Setting RX sample rate - req: 4000000Hz, actual: 4000000Hz Setting TX sample rate - req: 4000000Hz, actual: 4000000Hz bladeRF> set bandwidth 2.5E6 Set RX bandwidth - req: 2500000Hz actual: 2500000Hz Set TX bandwidth - req: 2500000Hz actual: 2500000Hz bladeRF> set rxvga1 30 bladeRF> set rxvga2 30 bladeRF> rx config format=csv n=40000 file=/temp/AKIEL.csv bladeRF> rx start bladeRF> rx config format=csv n=40000 file=/temp/BKIEL.csv bladeRF> rx start bladeRF> rx config format=csv n=40000 file=/temp/CKIEL.csv bladeRF> rx start bladeRF> rx start bladeRF> rx start bladeRF> rx config format=csv n=40000 file=/temp/DKIEL.csv bladeRF> rx start bladeRF> rx start bladeRF> rx start bladeRF> rx config format=csv n=80000 file=/temp/CKIEL.csv bladeRF> rx config format=csv n=80000 file=/temp/DKIEL.csv bladeRF> rx start
bladeRF> rx start bladeRF>

As you can see I set the frequency at 1575420000 Hz. This is the main frequency for GPS. I choose a samplerate of 4 MHz, USB2 can handle this. I choose 2.5 MHz for the bandwidth. The gains rxvga1 and rxvga2 are put to the maximum: 30 dB.

Try to set a higher value, bladeRf software will protest with a funny remark!

I type start a few time to get rid of the discontinuities, but this might be stupid as well and cause problems. Again something I will have to investigate.

The output of my python-program



>>> ================================ RESTART ================================ >>> filename = /Temp/DKIEL.csv fsample = 4000000.0 centerfreq = 0 1 msec data is 4000 samples skip 10000 samples of input-data from datafile read 4000 samples from file averI = 1944131.0 averQ = 2448936.0 averI/n = 486.03275 averQ/n = 612.234 freq, [peak, index, prn] = -15000 [13745, 2616, 18] freq, [peak, index, prn] = -14500 [15293, 1835, 21] freq, [peak, index, prn] = -14000 [17838, 64, 23] freq, [peak, index, prn] = -13500 [15970, 2033, 22] freq, [peak, index, prn] = -13000 [14932, 0, 0] freq, [peak, index, prn] = -12500 [15607, 1528, 32] freq, [peak, index, prn] = -12000 [18891, 64, 23] freq, [peak, index, prn] = -11500 [15460, 329, 9] freq, [peak, index, prn] = -11000 [15352, 803, 11] freq, [peak, index, prn] = -10500 [14976, 1098, 27] freq, [peak, index, prn] = -10000 [14131, 1312, 14] freq, [peak, index, prn] = -9500 [15606, 2364, 21] freq, [peak, index, prn] = -9000 [18250, 64, 23] freq, [peak, index, prn] = -8500 [15687, 3959, 28] freq, [peak, index, prn] = -8000 [14267, 2166, 26] freq, [peak, index, prn] = -7500 [16446, 318, 31] freq, [peak, index, prn] = -7000 [18228, 64, 23] freq, [peak, index, prn] = -6500 [14471, 3802, 31] freq, [peak, index, prn] = -6000 [16286, 64, 23] freq, [peak, index, prn] = -5500 [14694, 1, 11] freq, [peak, index, prn] = -5000 [14824, 662, 10] freq, [peak, index, prn] = -4500 [14487, 3566, 11] freq, [peak, index, prn] = -4000 [15549, 1679, 17] freq, [peak, index, prn] = -3500 [16070, 2864, 21] freq, [peak, index, prn] = -3000 [16865, 144, 12] freq, [peak, index, prn] = -2500 [14733, 896, 15] freq, [peak, index, prn] = -2000 [19962, 1668, 10] freq, [peak, index, prn] = -1500 [17504, 64, 23] freq, [peak, index, prn] = -1000 [17899, 64, 23] freq, [peak, index, prn] = -500 [17439, 2403, 6] freq, [peak, index, prn] = 0 [32753, 64, 23] freq, [peak, index, prn] = 500 [26798, 64, 23] freq, [peak, index, prn] = 1000 [28387, 782, 13] freq, [peak, index, prn] = 1500 [18584, 63, 23] freq, [peak, index, prn] = 2000 [18123, 64, 23] freq, [peak, index, prn] = 2500 [38468, 1668, 10] freq, [peak, index, prn] = 3000 [15833, 64, 23] freq, [peak, index, prn] = 3500 [13939, 970, 20] freq, [peak, index, prn] = 4000 [21279, 388, 7] freq, [peak, index, prn] = 4500 [24322, 387, 7] freq, [peak, index, prn] = 5000 [16517, 64, 23] freq, [peak, index, prn] = 5500 [16828, 55, 19] freq, [peak, index, prn] = 6000 [28200, 1154, 5] freq, [peak, index, prn] = 6500 [21470, 1153, 5] freq, [peak, index, prn] = 7000 [16290, 1466, 28] freq, [peak, index, prn] = 7500 [14327, 1471, 14] freq, [peak, index, prn] = 8000 [15135, 2887, 15] freq, [peak, index, prn] = 8500 [14257, 1340, 24] freq, [peak, index, prn] = 9000 [15520, 2217, 4] freq, [peak, index, prn] = 9500 [14830, 1473, 2] freq, [peak, index, prn] = 10000 [17611, 64, 23] freq, [peak, index, prn] = 10500 [14564, 1794, 22] freq, [peak, index, prn] = 11000 [15766, 1328, 13] freq, [peak, index, prn] = 11500 [15376, 2485, 28] freq, [peak, index, prn] = 12000 [17781, 64, 23] freq, [peak, index, prn] = 12500 [15759, 808, 10] freq, [peak, index, prn] = 13000 [14780, 124, 26] freq, [peak, index, prn] = 13500 [15648, 64, 23] freq, [peak, index, prn] = 14000 [18595, 63, 23] freq, [peak, index, prn] = 14500 [14543, 1736, 15] sorted list, max at top of the list [peak, index, prn, freq] [38468, 1668, 10, 2500] [32753, 64, 23, 0] [28387, 782, 13, 1000] [28200, 1154, 5, 6000] [26798, 64, 23, 500] [24322, 387, 7, 4500] [21470, 1153, 5, 6500] [21279, 388, 7, 4000] [19962, 1668, 10, -2000] [18891, 64, 23, -12000] [18595, 63, 23, 14000] [18584, 63, 23, 1500] [18250, 64, 23, -9000] [18228, 64, 23, -7000] [18123, 64, 23, 2000] [17899, 64, 23, -1000] [17838, 64, 23, -14000] [17781, 64, 23, 12000] [17611, 64, 23, 10000] [17504, 64, 23, -1500] [17439, 2403, 6, -500] [16865, 144, 12, -3000] [16828, 55, 19, 5500] [16517, 64, 23, 5000] [16446, 318, 31, -7500] [16290, 1466, 28, 7000] [16286, 64, 23, -6000] [16070, 2864, 21, -3500] [15970, 2033, 22, -13500] [15833, 64, 23, 3000] [15766, 1328, 13, 11000] [15759, 808, 10, 12500] [15687, 3959, 28, -8500] [15648, 64, 23, 13500] [15607, 1528, 32, -12500] [15606, 2364, 21, -9500] [15549, 1679, 17, -4000] [15520, 2217, 4, 9000] [15460, 329, 9, -11500] [15376, 2485, 28, 11500] [15352, 803, 11, -11000] [15293, 1835, 21, -14500] [15135, 2887, 15, 8000] [14976, 1098, 27, -10500] [14932, 0, 0, -13000] [14830, 1473, 2, 9500] [14824, 662, 10, -5000] [14780, 124, 26, 13000] [14733, 896, 15, -2500] [14694, 1, 11, -5500] [14564, 1794, 22, 10500] [14543, 1736, 15, 14500] [14487, 3566, 11, -4500] [14471, 3802, 31, -6500] [14327, 1471, 14, 7500] [14267, 2166, 26, -8000] [14257, 1340, 24, 8500] [14131, 1312, 14, -10000] [13939, 970, 20, 3500] [13745, 2616, 18, -15000] bigger than 1.5 * average: [peak, index, prn, freq[ [38468, 1668, 10, 2500] [32753, 64, 23, 0] [28387, 782, 13, 1000] [28200, 1154, 5, 6000] [26798, 64, 23, 500] >>>

Results

Those results are encouraging. It is clear that prn 10, 23, 13, 5 and 23 are prominently present.

From my notes:

According to my Garmin, there are
PRN = 5,7,13,16,23
Smaller amplitudes: 2,8, 10 and 29
All satellites 2,4,5,7,8,9,10,13,16,23,29 
26 is invisible
Current location N 52 51 25.6 E 6 44 53.3
Time 17:05 Amsterdam Time (CET?)

 OK, now you have the GPS-coordinates of that beautiful minicamping somewhere in the Netherlands!

Epilogue

From now on I might go further with this python-program. I might generate a huge data-file, a recording of 20 seconds and try to decode the almanac-data from the satellites. But I want to do something else, something more hardware like. I think I am gonna reprogram the FPGA. First some study, later some experiments as a preparation for the real thing: a hardware multichannel GPS-correlator. But that is far in the future!

8 comments:

  1. Hi Kees de Groot ,

    This is a nice project. Can you the the received from your bladerf, published on the webseite for download?

    greeting Marco

    ReplyDelete
    Replies
    1. Marco,
      dunno what you want exactly. You can simply copy/paste the python-programs and run them in your own python-environment. The *.csv-files you can make yourself, just type the bladerf-cli-commands. I can send you the csv files or the python-files if you want. Just ask again what you want exactly and I'll try to help!

      Delete
    2. Hi Kees de Groot,
      Oh sorry, i have forget to write. I haven't a bladeRF. I need only the csv-files.
      greeting Marco

      Delete
  2. The file DKIEL.CSV is at:
    http://jmp.sh/fg3YypW

    ReplyDelete
  3. Hello, master.
    I am graduate student and I am very interested in GPS security. I will study the gps record&replay, spoofing, countermeasure. Firstly, I am implementing the gps signal record&replay. So, I buy the BladeRFx115, bias-tee, LNA, gps active antenna like unicorn team in DEF CON 23 (https://www.youtube.com/watch?v=jwJKMti_aw0)
    And I connected a tools like this to record the GPS signal.

    But, I could't received the gps signal enough to positioning. My rx is average I/Q= 130, 122 using your program. So, I think that the active antenna isn't operate well. Because your average I/Q = 486,612 in DKEIL.csv. Therefore,I connected the USB-DC cable to Bias-Tee to energize a DC to the antenna. And I get the I/Q = 270, 224. I couldn't increase the I/Q. If I couldn't record the gps signal enough to positioning, I can't implement the tx test.
    This is result of your program. I just edit the sentence about print.

    1) 1.1 <= average < 1.2:
    [529615, 3010, 23, -1500]
    [520345, 2915, 29, -500]
    [515706, 167, 4, 14500]

    2) 1.2 <= average < 1.3:

    3) 1.3 <= average < 1.4:
    [635206, 1278, 25, 500]

    4) 1.4 <= average < 1.5:

    5) 1.5 <= average :

    So, how can I solve this problem and I want to know about your tips in this test.

    ReplyDelete
    Replies
    1. Hi,

      I have run the python program on the Raspberry Pi 3. I have the same result. The runtime is a round 5 min.

      Delete
  4. Hi, could you please provide me an e-mail address to contact with you about some issues in my GPS project.

    Thanks.
    Baris Can.

    ReplyDelete