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

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!