Wednesday, November 27, 2013

Without bladeRF: correlation with FFT and problem: python class returns by value

My first styled page

Python problem: class returns by reference, not by value

I made some changes to my python-class Gold.

The problem is that in python an attribute of a class is returned by reference. This might be nasty in some situations.

#import old class Gold goldutil myGold1 = Gold(13) prn = myGold1. prnRegAnalogue prn[0] = 100 prn2 = myGold1. prnRegAnalogue

if I modify prn I modify the internal list prnRegAnalogue because they are identical: they have the same address!
So now prn[0] and prn2[0] both have the same value 100.

Therefore I changed my class. Now I hide prnRegAnalogue by prefixing the name with a double-underscore: __ prnRegAnalogue and I defined a method prnReg() to return a copy of prnRegAnalogue. The effect is that now the list prnRegAnalogue is called by value.

# import modified class Gold goldutil myGold1 = Gold(13) prn = myGold1. prnRegAnalogue prn[0] = 100 prn2 = myGold1. prnRegAnalogue

Now prn2[0] has not been changed.

The new class now reads:

class Gold: # Gold-sequence goldutil.py # Nov-2013 Kees de Groot # call: # myprn1 = Gold(3, 7) for code phase selection = (3,7) # or # myprn2 = Gold(4) for the 4th from the list # = satellite ID number = GPS PRN signal number # start with (0,0) to start numbering with 1 # note that PRN 34 and 37 are identical prnlist = [ (0,0), (2,6), (3,7), (4,8), (5,9), (1,9), (2,10), (1,8), (2,9), (3,10), (2,3), (3,4), (5,6), (6,7), (7,8), (8,9), (9,10), (1,4), (2,5), (3,6), (4,7), (5,8), (6,9), (1,3), (4,6), (5,7), (6,8), (7,9), (8,10), (1,6), (2,7), (3,8), (4,9), (5,10), (4,10), (1,7), (2,8), (4,10)] def __init__(self, a, b = 0): self.reg1 = [1]*11 # don't use first bit-position [0] self.reg2 = [1]*11 # so count bit-positions from 1 to 10 self.__prnRegBinary = [] # contains 0 and 1 self.__prnRegAnalogue = [] # contains -1 and +1 if b == 0: # a = GPS PRN signal number or satellite ID (self.prn1, self.prn2) = self.prnlist[a] else: # direct code phase selection (self.prn1, self.prn2) = (a, b) for i in range(1024): G1 = self.reg1[10] G2 = self.reg2[self.prn1] ^ self.reg2[self.prn2] out = G1 ^ G2 self.__prnRegBinary.append(out) # values 0 and 1 self.__prnRegAnalogue.append(2*out-1) # values -1 and +1 val1 = self.reg1[3] ^ self.reg1[10] val2 = (self.reg2[2] ^ self.reg2[3] ^ self.reg2[6] ^ self.reg2[8] ^ self.reg2[9] ^ self.reg2[10]) # shift one position to the right for j in range(9): # j = 0,1,2,3,4,5,6,7,8 k = 10 - j # k = 10,9,8,7,6,5,4,3,2 self.reg1[k] = self.reg1[k-1] self.reg2[k] = self.reg2[k-1] self.reg1[1] = val1 self.reg2[1] = val2 # return list by value, not by reference def prnReg(self): alist = [] for val in self.__prnRegAnalogue: alist.append(val) return alist

Correlation with fft


A very interesting property of the Fourier transform is that convolution in the time-domain can be done by multiplication in the frequency-domain.
That means that to calculate a @ b, in which @ is meant to stand for convolution, you can do that calculation as follows:
ifft( fft(a) * fft(b) )
which means do an fft of a, do an fft of b, multiply and inverse-fft the result.
But the result of fft(a) has to be complex conjugated first, which means that all imaginary values has to change sign.
So, the procedure is: ifft(conjugate(fft(a)) * fft(b))
It sound strange perhaps, but this complex (sic) procedure might be very efficient to calculate because in practice the part conjugate(fft(a)) has to be calculated only once during initialization. The remaining calculation boils down to
c = conjugate(fft(a)

ifft(c * fft(b))

To calculate the correlation between a and b.

I made a small python program to experiment with this procedure:

import goldutil as gl # home-made Gold-sequence generator import pylab as pl # module for plotting import numpy as np # module for fft # CorrelFFTdemo2s.py # Simulation of message constructed using Gold-sequences # - Define an 8-bit message: eg 0 0 1 1 0 1 0 0 1 = -1, 1, 1, -1, 1, -1, -1, 1 # - Construct an array with 9*1024 samples of the same Gold-sequence, # eg. satellite id 12 # - Multiply the first 1024 samples with 1, the 2n set with -1, # the 3rd set with -1, etc to code our message # - Now try to decode our message using correlation implemented with FFT # nov 2013 Kees de Groot messageBin = (0, 0, 1, 1, 0, 1, 0, 0, 1) print "messageBin = ", messageBin messageAna = [] for val in messageBin: messageAna.append(2*val-1) print "messageAna = ", messageAna sat1 = 12 # this is the sat we are looking for myGold1 = gl.Gold(sat1) pprn1 = myGold1.prnReg() # get list with 1024 bits of Goldsequence # code message into codedmsg by multiplying codedmsg = [] for mesbit in messageAna: for prnbit in pprn1: codedmsg.append(mesbit * prnbit) codedmsg = [0]*200 + codedmsg # add some zeroes to improve readability of plot # decode codedmsg by correlation with origional pprn1 Gold-sequence # cross-correlation with the right Gold-sequence # correlate codedmsg with pprn1 and put result in ac # once-only calculate complex conjugate of prn-sequence pprn1 # problem is that I want the fft to have the same length as codedmsg # so, just add zeros at the end of pprn1 numzeros = len(codedmsg) - len(pprn1) pprn1 = pprn1 + [0]*numzeros fftprn = np.fft.rfft(pprn1) # complex conjugate (there must be a better method) for i in range(len(fftprn)): fftprn[i] = fftprn[i].conjugate() fftcodedmsg = np.fft.rfft(codedmsg) # convolution (correlation) == multiplying with complex conjugate # after fft and inverse-fft the result # multiply fftcodedmsg with fftprn for i in range(len(fftcodedmsg)): fftcodedmsg[i] *= fftprn[i] # inverse-fft ac = np.fft.irfft(fftcodedmsg) pl.figure(1) pl.title("correlation, calculated with fft-method") pl.ylim(-1500, 1500) pl.stem(ac) pl.show()
This program might look daunting, but there is really a lot of comment and not so much code.

The output is:

Python 2.7.3 (default, Apr 10 2012, 23:31:26) [MSC v.1500 32 bit (Intel)] on win32 Type "copyright", "credits" or "license()" for more information. >>> ================================ RESTART ================================ >>> messageBin = (0, 0, 1, 1, 0, 1, 0, 0, 1) messageAna = [-1, -1, 1, 1, -1, 1, -1, -1, 1]
Graphics


It is obvious that this is exactly the same result as in the previous blog-message where I programmed plain convolution.

Correlation with fft, block-mode

I expect to handle my real-time samples in blocks of 1024 chips. A chip is one bit of a Gold-sequence. In this example I simulated this procedure. In a loop I correlate 1024 chips from my input and decode the message bit by bit.

import goldutil as gl # home-made Gold-sequence generator import pylab as pl # module for plotting import numpy as np # module for fft # CorrelFFTblock.py # Simulation of message constructed using Gold-sequences # - Define an 8-bit message: eg 0 0 1 1 0 1 0 0 1 = -1, 1, 1, -1, 1, -1, -1, 1 # - Construct an array with 9*1024 samples of the same Gold-sequence, # eg. satellite id 12 # - Multiply the first 1024 samples with 1, the 2n set with -1, # the 3rd set with -1, etc to code our message # - Now try to decode our message using correlation implemented with FFT # do this piece-wise in blocks of 1024 bits # nov 2013 Kees de Groot def minmax(ac): minval = min(ac) maxval = max(ac) if abs(minval) > abs(maxval): result = 0 else: result = 1 return result messageBin = (0, 0, 1, 1, 0, 1, 0, 0, 1) print "messageBin = ", messageBin messageAna = [] for val in messageBin: messageAna.append(2*val-1) print "messageAna = ", messageAna sat1 = 12 # this is the sat we are looking for myGold1 = gl.Gold(sat1) pprn1 = myGold1.prnReg() # get list with 1024 bits of Goldsequence # code message into codedmsg by multiplying codedmsg = [] for mesbit in messageAna: for prnbit in pprn1: codedmsg.append(mesbit * prnbit) codedmsg = [0]*200 + codedmsg + [0]*824 # add some zeroes # decode codedmsg by correlation with origional pprn1 Gold-sequence # cross-correlation with the right Gold-sequence # correlate codedmsg with pprn1 and put result in ac # once-only calculate complex conjugate of prn-sequence pprn1 fftprn = np.fft.rfft(pprn1) # complex conjugate (there must be a better method) for i in range(len(fftprn)): fftprn[i] = fftprn[i].conjugate() # the above calculation has to be done once during initialization # now take 1024 bit blocks and correlate using fft N = len(codedmsg) / len(pprn1) print "N = ", N for n in range(N): fftcodedmsg = np.fft.rfft(codedmsg[n*1024:(n+1)*1024]) for i in range(len(fftcodedmsg)): fftcodedmsg[i] *= fftprn[i] ac = np.fft.irfft(fftcodedmsg) print "block ", n, ", message-bit = ", minmax(ac) pl.figure(n+1) pl.title("correlation, calculated with fft-method " + "block " + str(n)) pl.ylim(-1500, 1500) pl.stem(ac) pl.show()
with output:

Python 2.7.3 (default, Apr 10 2012, 23:31:26) [MSC v.1500 32 bit (Intel)] on win32 Type "copyright", "credits" or "license()" for more information. >>> ================================ RESTART ================================ >>> messageBin = (0, 0, 1, 1, 0, 1, 0, 0, 1) messageAna = [-1, -1, 1, 1, -1, 1, -1, -1, 1] N = 10 block 0 , message-bit = 0 block 1 , message-bit = 0 block 2 , message-bit = 1 block 3 , message-bit = 1 block 4 , message-bit = 0 block 5 , message-bit = 1 block 6 , message-bit = 0 block 7 , message-bit = 0 block 8 , message-bit = 1 block 9 , message-bit = 1
and a lot of nice graphics:












Conclusion

To find a message coded with a Gold-sequence one uses correlation. One method is to use convolution. It is a simple process but might cost a lot of time.

Another method is to make use of a property of the Fourier-transform:
Convolution in the time-domain is the same as multiplication in the frequency-domain
I experimented with this method and showed it to be possible.
Now I am confident to use this method with real samples from the bladeRF.

Keep tuned!


Thursday, November 21, 2013

For GPS on a bladeRF you have to start with Gold-sequences to find signals buried in noise: it is all about correlation

My first styled page

Gold-sequence

GPS-signals is all about correlation. The signals are completely buried in the noise. Only by searching (or knowing) the code-sequences and the Doppler frequencies is it possible to decode the signals. And they are all transmitted at that same frequency!

I did some experiments in python to get aquainted (again) with Gold-sequences and correlation. So, I programmed a class Gold and wrote some test-programs. Very nice experiments!

I had a lot of fun with it. Python is a very simple language. Plotting is a breeze. And with a class that you can import you get nice short programs. The IDLE editor is sufficient to get experiments running.


A very nice article is: “Understanding Spread Spectrum for Communications” at http://www.ni.com/white-paper/4450/en/

A python class "Gold"


In python I programmed the following class:

class Gold: # Gold-sequence goldutil.py # Nov-2013 Kees de Groot # call: # myprn1 = Gold(3, 7) for code phase selection = (3,7) # or # myprn2 = Gold(4) for the 4th from the list # = satellite ID number = GPS PRN signal number # start with (0,0) to start numbering with 1 # note that PRN 34 and 37 are identical prnlist = [ (0,0), (2,6), (3,7), (4,8), (5,9), (1,9), (2,10), (1,8), (2,9), (3,10), (2,3), (3,4), (5,6), (6,7), (7,8), (8,9), (9,10), (1,4), (2,5), (3,6), (4,7), (5,8), (6,9), (1,3), (4,6), (5,7), (6,8), (7,9), (8,10), (1,6), (2,7), (3,8), (4,9), (5,10), (4,10), (1,7), (2,8), (4,10)] def __init__(self, a, b = 0): self.shreg1 = [1]*11 # don't use first bit-position [0] self.shreg2 = [1]*11 # so count bit-positions from 1 to 10 self.prnRegBinary = [] # contains 0 and 1 self.prnRegAnalogue = [] # contains -1 and +1 if b == 0: # a = GPS PRN signal number or satellite ID (self.prn1, self.prn2) = self.prnlist[a] else: # direct code phase selection (self.prn1, self.prn2) = (a, b) for i in range(1024): G1 = self.shreg1[10] G2 = self.shreg2[self.prn1] ^ self.shreg2[self.prn2] out = G1 ^ G2 self.prnRegBinary.append(out) self.prnRegAnalogue.append(2*out-1) val1 = self.shreg1[3] ^ self.shreg1[10] val2 = (self.shreg2[2] ^ self.shreg2[3] ^ self.shreg2[6] ^ self.shreg2[8] ^ self.shreg2[9] ^ self.shreg2[10]) # shift one position to the left for j in range(9): k = 10 - j self.shreg1[k] = self.shreg1[k-1] self.shreg2[k] = self.shreg2[k-1] self.shreg1[1] = val1 self.shreg2[1] = val2
A first test-program:

import goldutil as gl myGold37 = gl.Gold(3, 7) prn1 = myGold37.prnRegBinary prn2 = myGold37.prnRegAnalogue print "the first 20 points of the binary version of Gold37:" for i in range(20): print prn1[i], print "\n..." print "\nthe first 20 points of the analogue version of Gold37:" for i in range(20): print prn2[i], print "\n..."

with output:

Python 2.7.3 (default, Apr 10 2012, 23:31:26) [MSC v.1500 32 bit (Intel)] on win32 Type "copyright", "credits" or "license()" for more information. >>> ================================ RESTART ================================ >>> the first 20 points of the binary version of Gold37: 1 1 1 0 0 1 0 0 0 0 1 1 1 0 0 0 0 0 1 1 ... the first 20 points of the analogue version of Gold37: 1 1 1 -1 -1 1 -1 -1 -1 -1 1 1 1 -1 -1 -1 -1 -1 1 1 ... >>>
With this framework I can write simple programs by importing the class Gold with the command
import goldutil as gl

Autocorrelation

As a first test I calculate autocorrelation.

# Gold02b.py # autocorrelation # nov-2013 Kees de Groot import goldutil as gl import pylab as pl myGold37 = gl.Gold(3, 7) prn2 = myGold37.prnRegAnalogue # autocorrelation ac = [0]*1024 for m in range(1024): for l in range(1024): ac[m] += prn2[l] * prn2[(l + m)%1024] print "\nthe first 20 points of the autocorrelation of myGold37:" for i in range(20): print ac[i], print "\n..." pl.figure(1) pl.title("autocorrelation, note big peak of 1024 at x = 0") pl.vlines(range(1024),[0]*1024,ac) pl.show()
with output:

Python 2.7.3 (default, Apr 10 2012, 23:31:26) [MSC v.1500 32 bit (Intel)] on win32 Type "copyright", "credits" or "license()" for more information. >>> ================================ RESTART ================================ >>> the first 20 points of the autocorrelation of myGold37: 1024 0 0 0 -68 -4 0 0 60 -8 4 4 -4 0 -4 60 64 60 -8 4 ...


There is nice peak at x=0, indicating positive correlation, the codes overlaps perfectly

Now correlate two different Gold-sequences.

# Gold02c.py # correlation between different Gold-sequences # nov-2013 Kees de Groot import goldutil as gl import pylab as pl myGold37 = gl.Gold(3, 7) # code (3,7) prn1 = myGold37.prnRegAnalogue myGold15 = gl.Gold(15) # satellite 15 prn2 = myGold15.prnRegAnalogue # correlation ac = [0]*1024 for m in range(1024): for l in range(1024): ac[m] += prn1[l] * prn2[(l + m)%1024] print "\nthe first 20 points of the correlation:" for i in range(20): print ac[i], print "\n..." pl.figure(1) pl.title("correlation between two different Gold-codes") pl.vlines(range(1024),[0]*1024,ac) pl.show()
with output:

Python 2.7.3 (default, Apr 10 2012, 23:31:26) [MSC v.1500 32 bit (Intel)] on win32 Type "copyright", "credits" or "license()" for more information. >>> ================================ RESTART ================================ >>> the first 20 points of the correlation: 0 0 0 -64 0 0 0 0 0 -4 64 0 4 -4 -64 0 -4 -68 0 64 ...



There are some peaks, but not a single strong one, so there is no correlation.

The real thing


Now it is time for the real thing.
  • Take a Gold-sequence of sat 18
  • Shift it to the right 200 chips
  • Add two more Gold-sequences: sat 4 and sat 17
  • Now I have a shifted contaminated satellite signal.
First correlate with a Gold-sequence sat=20. Of course there is not any correlation, only noise
Then, the wonder: correlate with a Gold-sequence sat=18. We see an obvious peak at x=200.
So we know we found the satellite and we know the delay of the signal!

import goldutil1 as gl import pylab as pl # Gold04.py # Generate a Gold-sequence, shift it 200 chips # and add two more different Gold-sequences # then, correlate it with an arbitrary Gold-sequence # and finally correlate it with the same Gold-sequence # nov 2013 Kees de Groot sat1 = 18 # this is the sat we are looking for myGold1 = gl.Gold(sat1) pprn1 = myGold1.prnRegAnalogue sat2 = 20 # this is our guess1 myGold2a = gl.Gold(sat2) pprn2a = myGold2a.prnRegAnalogue sat2 = 18 # this is our guess2 myGold2b = gl.Gold(sat2) pprn2b = myGold2b.prnRegAnalogue # shift pprn1 200 bits to the right k = 200 for i in range(1024): temp = pprn1[(i + k)%1024] pprn1[(i + k)%1024] = pprn1[i] pprn1[i] = temp # add two different Gold-sequences to pprn1 myGold3 = gl.Gold(4) myGold4 = gl.Gold(17) for i in range(1024): pprn1[i] = (pprn1[i] + myGold3.prnRegAnalogue[i] + myGold4.prnRegAnalogue[i]) # now pprn1 is contaminated with two other sequences # cross-correlation with wrong Gold-sequence ac = [0]*1024 for m in range(1024): for l in range(1024): ac[m] += pprn1[l] * pprn2a[(l + m)%1024] print "\n 20 points around the peak at 200:" for i in range(190, 210): print ac[i], print "\n..." pl.figure(1) pl.title("correlation with different sequence") pl.vlines(range(1024),[0]*1024,ac) # cross-correlation with the right Gold-sequence ac = [0]*1024 for m in range(1024): for l in range(1024): ac[m] += pprn1[l] * pprn2b[(l + m)%1024] print "\n 20 points around the peak at 200:" for i in range(190, 210): print ac[i], print "\n..." pl.figure(2) pl.title("correlation with the right sequence") pl.vlines(range(1024),[0]*1024,ac) pl.show()
With output:

Python 2.7.3 (default, Apr 10 2012, 23:31:26) [MSC v.1500 32 bit (Intel)] on win32 Type "copyright", "credits" or "license()" for more information. >>> ================================ RESTART ================================ >>> 20 points around the peak at 200: 12 116 112 28 4 -4 -8 12 0 -48 -48 -36 -8 128 8 -4 -88 -20 -24 -28 ... 20 points around the peak at 200: 48 104 124 184 72 -12 -32 0 60 8 708 -40 -20 -32 -32 80 12 -24 -80 -48 ...
Note there is no peak in the first run, but in the second run there is a peak with a value of 708, indicating strong correlation around delay = 200 chips.



Only noisy correlation, not one obvious peak



And there it is: a peak around x=200. A clear indication of the existence of a signal in the noise.

Conclusion:

I implemented the Gold-sequence as a class Gold in python.
With that class I first demonstrated a working implementation.
Then I did some experiments with (auto)correlation and finally demonstrated that you can find a certain code-sequence in a noisy contaminated signal and even calculate the delay of that signal.


Thursday, November 7, 2013

BladeRF data around GPS frequency of 1575.42 MHz

My first styled page

GPS revisited and some confessions

The book

An interesting book is:
 "A software defined GPS and Galileo receiver" , a single frequency approach, by Kai Borre et al.
In the back of the book is a CDROM with a dataset. That dataset contains the if-output of a real SDR tuned to a GPS-signal. I copied the test-file to a convenient place and renamed it to testdata.dat.
That started my struggle with GPS. In that timeperiod I did a course on java and wanted to do something ‘useful’ with that knowledge. I decided to build a GPS-receiver and program it in java. So, I bought some books, read a lot about it on the internet and really got a deep understanding of the system and the signals used in a GNSS-system.
I decided to first develop the software and later concentrate on the hardware. The software was difficult enough for me…
To make a really long story short I finally managed to read almanac-data from the GPS-signals. However I got stuck in the tracking loop. I never got the loop stable. That last problem never made it to my blog at http://www.qsl.net/pa1kdg/ by the way.
So, I also did not succeed in calculating the position from the data in that huge data-file. But, what the heck, I learned soo much of that whole project…

Next step

My next step now is restart that whole project with a bladeRF-board. I want to use python now as a programming language. Just because I can ;>))

TANSTAAFL There aint no such thing as a free lunch

Being a retired engineer, the Mathworks, the developers of MATLAB, is no longer interested in me. When I worked at Wageningen University and Research I got an enormous lot of software free. I was addicted to MATLAB (and LABVIEW by the way), with that software you can solve all your scientific problems. After my retirement I soon discovered that the love of the Mathworks I felt was quite different than real life love I perceive on a daily basis from the people around me. Just before retirement I got an SDR-board with a huge Xilinx FPGA. That would be my big and only project after retirement! Nope, Xilinx too discovered I am not longer interesting for them. Real-life for a CEO is about earning money, not about supporting retired engineers. I have not any problem to pay for hardware or software but 2000 bucks for my hobby is in-explainable to my wife (ok, this whole hobby is inexplainable…) OK, I am back on earth now with two feet on the ground, with my bladeRF and open source software. This is my real life now!!

Life is complicated though

On my way along my GPS-project I discovered a lot about linux, VMWare, github, python-eggs, C, VHDL, makefiles, Quartus, SDR-Radio, gnuradio, Octavo and finally plotting in python (and kind of python-MATLAB too). And I found many sidepaths, DesignSpark Mechanical e.g. and soldering SMD’s.

What now

Back to GPS. I want to record a data-file with real-bladeRF-data and try to process that file. First discover what satellites are visible, later read the almanac data. Finally try to get past that data-tracking-loop problem and solve the position-equations. In the meantime I want to try to get some calculations done in the FPGA. It should be possible to do the correlation-part in the FPGA. Quite a challenge!

Some facts

I concentrate on GPS, not GLONASS. GPS transmits on one frequency. All satellites transmit on the same frequency, that is to say, because of Doppler effects there is a slight deviation of that frequency for every satellite. Each satellite has a different code. If you know the code and the Doppler-frequency you can decode the signal of that very satellite. The signal is simply almanac-data at 50 baud. The satellites tell you constantly where they are. Because you know exactly which bit of the code-sequence you receive and because you ‘see’ 4 satelites or more you can calculate where your GPS-antenna is. For 50 $ you can buy a module that does this trick in a very good way, so…

Practical now

Ok, I tune bladeRF to 1575.42 MHz. With a bandwidth of 2 MHz I might expect to get a weak noisy signal, barely a signal, -160dBW. So, only noise.
There are 32 different code sequences. Just try all those sequences one by one.
There are about 60 different frequencies to consider because of the unknown Doppler shift. So try all frequencies.
In a 3D-display one can find some peaks. Then you know some codes and some Doppler shifts. From then on you can decode the satellite signals and get the almanac data. Then you know te distance to each found satellite. In 3D, if you know the distance to 3 or more satellites, you can calculate where you are.  Because you know where the satellites are you know where you are. Simple eh? That is GPS

GLONASS, the Russian system, uses different frequencies. No problem, bladeRF can handle that too. And, with an FPGA you can put the most time-consuming calculations on that chip.

That is the theory.

Tune to 1575.42 MHz. Bandwidth 2 MHz, so samplerate 4 MHz. VGA-gains at 30 dB. Get 1M samples into a file. Perform an FFT. Hopefully something interesting will show up.
Filename A1575.csv without GPSantenna connected and B1575.csv with the GPSantenna connected and pearing out of my window. In the meantime I have a TomTom and a Garmin device locked at the satellites in the sky so I know some data already.

Got my first two files with bladeRF-cli. GPS-systems say: satellites in the sky: 2,4,5,8,9,10,13,16 of which 2,4, 9 are strongest. So I might see satellites 2, 4 and 9 in my datafiles.


[INFO] Using libusb version 1.0.16.10774 [INFO] Found a bladeRF [INFO] Claimed all inferfaces successfully [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... [INFO] Change to alternate interface 3 [INFO] Change to alternate interface 1 [INFO] Changed into RF link mode: LIBUSB_SUCCESS / LIBUSB_TRANSFER_COMPLETED [INFO] Setting integer sample rate: 1000000 [INFO] Found r value of: 4 [INFO] MSx a + b/c: 316 + 4/5 [INFO] MSx a + b/c: 316 + 4/5 [INFO] MSx P1: 0x00009c66 (40038) P2: 0x00000002 (2) P3: 0x00000005 (5) [INFO] Calculated samplerate: 1000000 + 0/1 [INFO] Set actual integer sample rate: 1000000 [INFO] Setting integer sample rate: 1000000 [INFO] Found r value of: 4 [INFO] MSx a + b/c: 316 + 4/5 [INFO] MSx a + b/c: 316 + 4/5 [INFO] MSx P1: 0x00009c66 (40038) P2: 0x00000002 (2) P3: 0x00000005 (5) [INFO] Calculated samplerate: 1000000 + 0/1 [INFO] Set actual integer sample rate: 1000000 Done. bladeRF> print rxvga1 RXVGA1 Gain: 33 bladeRF> print rxvga2 RXVGA2 Gain: 3dB bladeRF> set rxvga2 33 [WARNING] Setting gain above 30dB? You crazy!! bladeRF> set rxvga2 30 bladeRF> print rxvga2 RXVGA2 Gain: 30dB bladeRF> set bandwidth 2000000 Set RX bandwidth - req: 2000000Hz actual: 2500000Hz Set TX bandwidth - req: 2000000Hz actual: 2500000Hz bladeRF> set frequency 1575420000 Set RX frequency: 1575420000Hz Set TX frequency: 1575420000Hz bladeRF> set bandwidth 2000000 Set RX bandwidth - req: 2000000Hz actual: 2500000Hz Set TX bandwidth - req: 2000000Hz actual: 2500000Hz bladeRF> set samplerate 4000000 [INFO] Setting integer sample rate: 4000000 [INFO] Found r value of: 1 [INFO] MSx a + b/c: 316 + 4/5 [INFO] MSx a + b/c: 316 + 4/5 [INFO] MSx P1: 0x00009c66 (40038) P2: 0x00000002 (2) P3: 0x00000005 (5) [INFO] Calculated samplerate: 4000000 + 0/1 [INFO] Set actual integer sample rate: 4000000 Setting RX sample rate - req: 4000000Hz, actual: 4000000Hz [INFO] Setting integer sample rate: 4000000 [INFO] Found r value of: 1 [INFO] MSx a + b/c: 316 + 4/5 [INFO] MSx a + b/c: 316 + 4/5 [INFO] MSx P1: 0x00009c66 (40038) P2: 0x00000002 (2) P3: 0x00000005 (5) [INFO] Calculated samplerate: 4000000 + 0/1 [INFO] Set actual integer sample rate: 4000000 Setting TX sample rate - req: 4000000Hz, actual: 4000000Hz bladeRF> set bandwidth 2000000 Set RX bandwidth - req: 2000000Hz actual: 2500000Hz Set TX bandwidth - req: 2000000Hz actual: 2500000Hz bladeRF> rx config format=csv n=1000000 file=/temp/A1575.csv bladeRF> rx start bladeRF> rx start bladeRF> rx start bladeRF> rx config format=csv n=1000000 file=/temp/B1575.csv bladeRF> rx start bladeRF> rx start bladeRF> rx start bladeRF>

Note the warning about gain higher than 30 dB ;>))
# FFTSDRdata03.py # Nov-2013 Kees de Groot # # show and analyse snapshot file from bladeRF # remove DC-component # show data in different ways; I/Q, realtimeplot, FFT import pylab as pl import numpy as np import math filename = 'B1575.csv' fsample = 4E6 timestep = 1 / fsample f = open('/Temp/' + filename, 'r') print f print "fsample = ", fsample dataI = [] dataQ = [] n = 0 sumI = 0 sumQ = 0 # read I, Q - values into memory for line in f: list = line.split(',') # two values with comma inbetween Q = int(list[0]) I = int(list[1]) # print 1st 10 lines of data-file if n == 0: print '1st 10 lines of I/Q-values' if n < 10: print I, Q dataI.append(I) # save data in memory dataQ.append(Q) sumI += I sumQ += Q n += 1 averI = sumI / n # calculate average averQ = sumQ / n #averI = 0 # debug #averQ = 0 # debug print "n = ", n # remove DC-component and construct complex dataC dataC = [] for i in range(n): I = dataI[i] - averI Q = dataQ[i] - averQ dataI[i] = I dataQ[i] = Q dataC.append(complex(I,Q)) # this one has complex data pl.figure(1) pl.title("Q versus I of " + filename) pl.xlabel('=======> I') pl.ylabel('=======> Q') pl.plot(dataI, dataQ, 'ro') pl.figure(2) pl.subplot(211) pl.title("I-values of " + filename) pl.ylabel('=======> I') pl.plot(dataI) pl.subplot(212) pl.title("Q-values of " + filename) pl.ylabel('=======> Q') pl.plot(dataQ) pl.figure(3) sp = np.fft.fft(dataC) # use fft, not rfft, since dataC is complex freq = np.fft.fftfreq(n, d=timestep) # for convenient x-axis-ticks pl.title('fft of ' + filename) pl.xlabel('frequency') pl.ylabel('magnitude') pl.plot(freq, abs(sp)/n) # abs converts to magnitude; /n normalizes pl.show()
Run the program for the file without antenna connected: I would expect simply noise but there are three peaks at -583kHz, 0 kHz and 1022 kHz. Perhaps the gain was too big



I would expect simply noise but there are three peaks at
-583kHz, 0 kHz and 1022 kHz.

Now connect GPS-antenna





Perhaps the gain was too big



bladeRF> set rxvga2 27 bladeRF> set rxvga1 27 bladeRF> rx config format=csv n=1000000 file=/temp/C1575.csv bladeRF> rx start bladeRF> rx start bladeRF> rx start bladeRF>







I am still puzzled about the fact that I see a signal when no antenna is connected. I remove the SMB cable from RX and set gains to 0 dB. Then create file D1575.csv:





Biggest peak is at 0 Hz exactly, the other is at 1.01882E6, so around 1 MHz

Again without antenna but a 50 ohm dummyload directly connected to bladeRF-RX-input.
Something went wrong with the board. Error because of physical touching it with too big dummyload directly connected?


[INFO] Using libusb version 1.0.16.10774 [INFO] Found a bladeRF [INFO] Claimed all inferfaces successfully [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... [INFO] Change to alternate interface 3 [INFO] Change to alternate interface 1 [INFO] Changed into RF link mode: LIBUSB_SUCCESS / LIBUSB_TRANSFER_COMPLETED [INFO] Setting integer sample rate: 1000000 [INFO] Found r value of: 4 [INFO] MSx a + b/c: 316 + 4/5 [INFO] MSx a + b/c: 316 + 4/5 [INFO] MSx P1: 0x00009c66 (40038) P2: 0x00000002 (2) P3: 0x00000005 (5) [INFO] Calculated samplerate: 1000000 + 0/1 [INFO] Set actual integer sample rate: 1000000 [INFO] Setting integer sample rate: 1000000 [INFO] Found r value of: 4 [INFO] MSx a + b/c: 316 + 4/5 [INFO] MSx a + b/c: 316 + 4/5 [INFO] MSx P1: 0x00009c66 (40038) P2: 0x00000002 (2) P3: 0x00000005 (5) [INFO] Calculated samplerate: 1000000 + 0/1 [INFO] Set actual integer sample rate: 1000000 Done. bladeRF> print rxvga1 RXVGA1 Gain: 33 bladeRF> set rxvga1 0 bladeRF> set rxvga2 0 bladeRF> set samplerate 1000000 [INFO] Setting integer sample rate: 1000000 [INFO] Found r value of: 4 [INFO] MSx a + b/c: 316 + 4/5 [INFO] MSx a + b/c: 316 + 4/5 [INFO] MSx P1: 0x00009c66 (40038) P2: 0x00000002 (2) P3: 0x00000005 (5) [INFO] Calculated samplerate: 1000000 + 0/1 [INFO] Set actual integer sample rate: 1000000 Setting RX sample rate - req: 1000000Hz, actual: 1000000Hz [INFO] Setting integer sample rate: 1000000 [INFO] Found r value of: 4 [INFO] MSx a + b/c: 316 + 4/5 [INFO] MSx a + b/c: 316 + 4/5 [INFO] MSx P1: 0x00009c66 (40038) P2: 0x00000002 (2) P3: 0x00000005 (5) [INFO] Calculated samplerate: 1000000 + 0/1 [INFO] Set actual integer sample rate: 1000000 Setting TX sample rate - req: 1000000Hz, actual: 1000000Hz bladeRF> set frequency 1575420000 Set RX frequency: 1575420000Hz Set TX frequency: 1575420000Hz bladeRF> set bandwidth 2000000 Set RX bandwidth - req: 2000000Hz actual: 2500000Hz Set TX bandwidth - req: 2000000Hz actual: 2500000Hz bladeRF> rx config format=csv n=1000000 file=/temp/E1575.csv bladeRF> start Unrecognized command: start bladeRF> rx start bladeRF> rx start Error: Operation invalid in current state bladeRF> rx start bladeRF> rx start bladeRF> rx start bladeRF>
Beware: sample rate = 1 Mhz


One peak at zero Hz, the other is at 19.354 Hz.
Now samplerate = 2 MHz





Now same with samplerate = 4MHz






Still, no antenna connected, but dummyload directly at RX

Conclusion

I have a problem. I expected a flat noise floor without any antenna connected. But I see some spurious signals. I might be able to explain the signal at 0 Hz. But the others are strange.
I have to check again without any antenna connected and vary the parameters that I can set:
n, the number of samples, it might be a hickup of the interface
samplerate, there is a clock, so there are some signals
frequency, mixer-products?
rxvga1 the more gain, the more spurious signals
rxvga2 same
bandwidth, dont think so

Well I had a nice evening trying to get a GPS-satellite-signal but discovered my bladeRF also is very creative?




Monday, November 4, 2013

Data analysis of a file of bladeRF samples with python

My first styled page

Finally, some python-programming

With bladeRF-CLI, the bladeRF-control-program, one can collect received data into a file. That file contains I- and Q-samples. With a very simple python-program I try to get some insight in the received data.

First python-program for bladeRF-datafile

Create a data-file with bladeRF-cli:
bladeRF> version [WARNING] FPGA currently does not have a version number. Serial #: b436de8c8212b9aeaaeba852246866e7 VCTCXO DAC calibration: 0x9a5d FPGA size: 115 KLE Firmware version: 1.5 FPGA version: 0.0 bladeRF> print bandwidth RX Bandwidth: 28000000Hz TX Bandwidth: 28000000Hz bladeRF> print samplerate [INFO] Calculated samplerate: 1000000 + 0/1 [INFO] Calculated samplerate: 1000000 + 0/1 RX sample rate: 1000000 TX sample rate: 1000000 bladeRF> set frequency 800000000 Set RX frequency: 800000000Hz Set TX frequency: 800000000Hz bladeRF> rx config format=csv n=1024 file=C:/temp/A800MHzCW.csv bladeRF> rx start bladeRF>
And now analyse the datafile with a python-program

f = open('/Temp/A800MHzCW.csv', 'r') n = 0 averageI = 0 averageQ = 0 for line in f: if n < 10: print line, list = line.split(',') I = int(list[0]) Q = int(list[1]) #print "I =",I, "Q =",Q averageI = averageI + I averageQ = averageQ + Q n = n + 1 averageI = averageI / n averageQ = averageQ / n print f print "averageI = ", averageI, " averageQ = ", averageQ
and the output is:
Python 2.7.3 (default, Apr 10 2012, 23:31:26) [MSC v.1500 32 bit (Intel)] on win32
Type "copyright", "credits" or "license()" for more information.
>>> ================================ RESTART ================================
>>> 
110, 100
113, 100
112, 103
114, 104
112, 105
110, 101
114, 100
112, 100
115, 97
113, 97

averageI =  110  averageQ =  98
>>>

This may seem a silly simple program but:

  • I can read a file
  • I can interpret the content and get I and Q samples
  • I can print values
  • I can calculate some data
So, this is a sort of 'Hello world' program.

Plotting in python

How?
I need to install matplotlib for python 2.7
Download from http://sourceforge.net/projects/matplotlib/postdownload?source=dlp
And execute matplotlib-1.3.1.win-amd64-py2.7.exe


Helloo, I _have_ python 2.7 installed!!!

I have python 2.7 running on my machine, obviously not registered properly?
I’’ try matplotlib-1.3.1.win32-py2.7.exe instead


Seems to be ok.

But...


So...


From http://www.lfd.uci.edu/~gohlke/pythonlibs/
python-dateutil-1.5.win32-py2.7.exe

>>> import matplotlib

Traceback (most recent call last):
  File "<pyshell#4>", line 1, in <module>
    import matplotlib
  File "C:\Python27\lib\site-packages\matplotlib\__init__.py", line 125, in <module>
    raise ImportError("matplotlib requires pyparsing")
ImportError: matplotlib requires pyparsing
>>>

From http://www.lfd.uci.edu/~gohlke/pythonlibs/
pyparsing-2.0.1.win32-py2.7.exe

>>> import matplotlib
>>>

Quite encouraging!

All examples from internet about plotting with matplotlib function very well!

Nice playing by the way!!

Sometimes it is difficult to go persistently only into the direction of your goal. Before my retirement I had to focus on result, that was very important. But there are so many distracting sidepaths... Nowadays I play a lot, go into an enormous lot of sidepaths, and that is great fun! So, progress is not the issue anymore...

Some data-files

summary:
freq = 1575.450.000 Hz
bandwidth = 28.000.000 Hz
samplerate = 1.000.000 Hz
file=/temp/AA.csv
n = 5000 samples

# FFTSDRdata01.py # Oct-2013 Kees de Groot # # show data collected from bladeRF import pylab as pl import numpy as np import math filename = 'AA.csv' f = open('/Temp/' + filename, 'r') print f n = 0 dataI = [] dataQ = [] dataC = [] for line in f: # print 1st 10 lines of data-file if n < 10: print line, list = line.split(',') I = int(list[0]) Q = int(list[1]) dataI.append(I) dataQ.append(Q) dataC.append(complex(I,Q)) n += 1 print "n = ", n pl.figure(1) pl.title("Q versus I of " + filename) pl.xlabel('=======> I') pl.ylabel('=======> Q') pl.plot(dataI, dataQ, 'ro') pl.figure(2) pl.subplot(211) pl.title(filename) pl.ylabel('=======> I') pl.plot(dataI) pl.title("I-values of " + filename) pl.subplot(212) pl.title(filename) pl.ylabel('=======> Q') pl.plot(dataQ) pl.figure(3) fsample = 1E6 timestep = 1 / fsample sp = np.fft.fft(dataC) freq = np.fft.fftfreq(n, d=timestep) pl.title(filename) pl.plot(freq, abs(sp)) pl.show()


Python 2.7.3 (default, Apr 10 2012, 23:31:26) [MSC v.1500 32 bit (Intel)] on win32
Type "copyright", "credits" or "license()" for more information.
>>> ================================ RESTART ================================
>>>
<open file '/Temp/AA.csv', mode 'r' at 0x0343D390>
80, 80
82, 77
81, 77
77, 79
80, 77
80, 81
78, 79
81, 76
84, 80
80, 76
n =  5000




No rocket-science, but quite encouraging.
I see some interruptions/setting-changes/discontinuities, but what-the-heck: nice pictures!
(I may be a strange person, but I really became quite excited when I first got these pictures!)

Settable gains inside the bladeRF-board

From the log of SDR-Radio console:

16:33:51> bladeRF
16:33:51>     --Starting sending
16:33:51>     Set RX VGA gain 2: value 30dB, status = Success
16:33:51>     Set RX LPF mode: Normal, status = Success
16:33:51>     Set RX VGA gain 1: value 30dB, status = Success
16:33:51>     Set LNA gain: Maximum, status = Success
16:33:51> Input Reader (Part 1)
16:33:51>     Starting listener, Radio defn 04A686B0, buffer size 65536

There is a third settable gain: lna-gain.


[INFO] Using libusb version 1.0.16.10774 [INFO] Found a bladeRF [INFO] Claimed all inferfaces successfully [INFO] Change to alternate interface 1 [INFO] Changed into RF link mode: LIBUSB_SUCCESS / LIBUSB_TRANSFER_COMPLETED [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> print bandwidth RX Bandwidth: 28000000Hz TX Bandwidth: 28000000Hz bladeRF> print frequency RX Frequency: 1575450000Hz TX Frequency: 1575450000Hz bladeRF> print samplerate [INFO] Calculated samplerate: 38400000 + 0/1 [INFO] Calculated samplerate: 1000000 + 0/1 RX sample rate: 38400000 TX sample rate: 1000000 bladeRF> print vga1gain print: Invalid parameter (vga1gain) bladeRF> print gainvga1 print: Invalid parameter (gainvga1) bladeRF> help set set The set command takes a parameter and an arbitrary number of arguments for that particular command. The parameter is one of: bandwidth Bandwidth settings config Overview of everything frequency Frequency settings lmsregs LMS6002D register dump loopback Loopback settings mimo MIMO settings pa PA settings pps PPS settings refclk Reference clock settings rxvga1 Gain setting of RXVGA1 in dB (range: ) rxvga2 Gain setting of RXVGA2 in dB (range: ) samplerate Samplerate settings trimdac VCTCXO Trim DAC settings txvga1 Gain setting of TXVGA1 in dB (range: ) txvga2 Gain setting of TXVGA2 in dB (range: ) bladeRF> print rxvga1 RXVGA1 Gain: 30 bladeRF> print rxvga2 RXVGA2 Gain: 30dB bladeRF>
nope, no other gain-setting.
What is the maximum data-rate / samplerate for USB2?

A lot of different experiments lead to the following python-program:

A.csv is a file without any connection from bladeRF to the signal generator, so, pure noise

The python-program:

# FFTSDRdata03.py # Nov-2013 Kees de Groot # # show and analyse snapshot file from bladeRF # remove DC-component # show data in different ways; I/Q, realtimeplot, FFT import pylab as pl import numpy as np import math filename = 'A.csv' fsample = 1E6 timestep = 1 / fsample f = open('/Temp/' + filename, 'r') print f print "fsample = ", fsample dataI = [] dataQ = [] n = 0 sumI = 0 sumQ = 0 # read I, Q - values into memory for line in f: list = line.split(',') # two values with comma inbetween Q = int(list[0]) I = int(list[1]) # print 1st 10 lines of data-file if n == 0: print '1st 10 lines of I/Q-values' if n < 10: print I, Q dataI.append(I) # save data in memory dataQ.append(Q) sumI += I sumQ += Q n += 1 averI = sumI / n # calculate average averQ = sumQ / n print "n = ", n # remove DC-component and construct complex dataC dataC = [] for i in range(n): I = dataI[i] - averI Q = dataQ[i] - averQ dataI[i] = I dataQ[i] = Q dataC.append(complex(I,Q)) # this one has complex data pl.figure(1) pl.title("Q versus I of " + filename) pl.xlabel('=======> I') pl.ylabel('=======> Q') pl.plot(dataI, dataQ, 'ro') pl.figure(2) pl.subplot(211) pl.title("I-values of " + filename) pl.ylabel('=======> I') pl.plot(dataI) pl.subplot(212) pl.title("Q-values of " + filename) pl.ylabel('=======> Q') pl.plot(dataQ) pl.figure(3) sp = np.fft.fft(dataC) # use fft, not rfft, since dataC is complex freq = np.fft.fftfreq(n, d=timestep) # for convenient x-axis-ticks pl.title('fft of ' + filename) pl.xlabel('frequency') pl.ylabel('magnitude') pl.plot(freq, abs(sp)/n) # abs converts to magnitude; /n normalizes pl.show()
With output
A.csv is a file without any connection from bladeRF to the signal generator, so, pure noise:





There is an enormous peak at 200kHz which was not there before on other pictures. I have to find out what is happening there.

I did a lot of more measurements:

BladeRF is setup:
RX Bandwidth:  28.000.000Hz
RX sample rate: 1.000.000
RX frequency:  646.000.000 Hz
RXVGA1 Gain:  33
RXVGA2 Gain:   3dB
A.csv is a file without any connection from bladeRF to the signal generator, so, pure noise
Now the signal generator, a HP8656B is directly connected to RX of bladeRF
B.csv is 646 MHz, AM 1 kHz 75% modulation, 1 µV amplitude
C.csv is 3,98 µV
D.csv is 7.94 µV
E.csv is 31.6 µV
F.csv is 63,1 µV

Back to 1 µV and now FM 1 kHz sweep 75 kHz wide

AF.csv is 1 µV
BF.csv is 3.98 µV
CF.csv is 7.94 µV
DF.csv is 31,6 µV

This last signal, DF.csv, is very nice:




Conclusion

I managed to get a file with data (I and Q-samples) from bladeRF-board. With a program written in python I managed to read these files, interpret and show the data. Nice results, and a big question: what is that signal at 200kHz? So, every experiment gives you new things to investigate...

Wife: what are you doing?
Me: Well, you know, my project..
Wife: What are you trying to achieve?
Me: Well, at the end this board will tell me where I am..
Wife: Do you know where you are at the end?
Me: Yes..
Wife: So, ???

This whole project, that keeps me out of all evil places of the world, that keeps me home, ok, that keeps me behind my laptop the whole evening and the first part of the night, that keeps me active,
h o w  c a n  I  e v e r  e x p l a i n  t o  a n  a l i e n  w h a t  I  a m  d o i n g  h e r e?
My wife stops asking when I speak the magic spell: "it is just my hobby!"