Wednesday, February 26, 2014

Removing the DC component on the fly

Correlator

For the GPS-frontend I need a correlator. The ADC-output from the LMS is 12 bits. That is a lot of bits and makes the correlator huge. Therefore I want to reduce the ADC-output to a single analogue value: -1 or +1.
This reduction is not too difficult to put in an FPGA. All you have to do is:
If sample > average then output := 1 else output := -1  (or perhaps 0 in the latter case?)

0/1 or -1/+1 ?

Multiply
+1 * +1 = +1
+1 * -1 = -1
-1 * +1 = -1
-1 * -1 = +1

This looks like an EXOR:
1 XOR 1 = 0
1 XOR 0 = 1
0 XOR 1 = 1
0 XOR 0 = 0

Average or DC-component

If you have a file with samples, calculating the average is very simple. Put all samples in an accumulator and divide the result by the number of samples. But waiting for a zillion of samples is not necessary. If you average the current 100 samples you have a fairly accurate estimation of the average. The average, by the way, is the DC-component. With an FFT, or DFT for that matter, X[0] is the DC-component, the average. But that is a little bit overdone. How many samples do you need for a reasonable estimate of the average? I have to divide the accumulated sum of the samples by the number of samples. If the number of samples is a power of 2 then division is simply removing some zeroes or shifting to the right!

So I need two simple parallel processes:

reducer: Process()
    If sample > average
    then output := 1
    else output := -1
end reducer

calc average: Process()
accumaverage := 0
counter := 0
    accumaverage := accumaverage + sample
    counter := counter + 1
    if counter == 128
    then
        average := accumaverage/128
        counter := 0
        accumaverage := 0
end calc average

This should be done for I and for Q.  
                                                                                                                                                                  What number system is used by the ADC? Two’s complement
How many samples: 128, more? less?

If the average is reasonable constant, this is a good estimation if you have ‘enough’ samples. If the average fluctuates a bit, you want as few samples as possible. How few? Furthermore you want the estimation of the average asap, so fewer samples is better.

Simulations in python

# binaryAverage01.py
# Feb-2014 Kees de Groot
#
# remove average (DC-component)
# by accumulating 128 samples and divide by 128
# perform some experiments/simulations
   
import pylab as pl
import numpy as np
import math
import sys

######################################################################
############# parameter section ######################################
######################################################################

path = '/Temp/'
filename = 'DDOUD.csv'
f = open(path + filename, 'rb')
print "filename = ", path + filename
skip = 0 # 10000 # there is a discontinuity in the file
            # so skip half of it
print "skip ", skip, "samples of input-data from datafile"

###########################################
#############  main loop ##################
###########################################

dataI = []
dataQ = []
averI = 0.0
averQ = 0.0

# read (I,Q) from bladeRF file
n = 1024
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
       
pl.figure(4)
pl.title("I/Q, DC-component removed")
pl.plot(dataI, dataQ, 'ro')
pl.show()

sys.exit('debug')

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 ================================
>>>
filename =  /Temp/DDOUD.csv
skip  0 samples of input-data from datafile
read  1024  samples from file
averI =  518907.0
averQ =  628875.0
averI/n =  506.745117188
averQ/n =  614.135742188


The center of the cluster datapoints is at about (500,600)


With the DC component removed, the center is at (0, 0).

Moving average


# binaryAverage02.py
# Feb-2014 Kees de Groot
#
# remove average (DC-component)
# by accumulating 128 samples and divide by 128
# perform some experiments/simulations
   
import pylab as pl
import numpy as np
import math
import sys

######################################################################
############# parameter section ######################################
######################################################################

path = '/Temp/'
filename = 'DDOUD.csv'
f = open(path + filename, 'rb')
print "filename = ", path + filename
skip = 0 # 10000 # there is a discontinuity in the file
            # so skip half of it
print "skip ", skip, "samples of input-data from datafile"

###########################################
#############  main loop ##################
###########################################

n = 1024 # number of samples to process
countermax = 128 # number of samples to average
counter = 0
dataI = []
dataQ = []
averI = 0.0
averQ = 0.0
accQ = 0.0
accI = 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])
    accQ += Q
    I = int(list[1])
    accI += I
    dataI.append(I - averI)
    dataQ.append(Q - averQ)
    i += 1
    if i >= n:
        break
    counter += 1
    if counter == countermax:
        averQ = accQ / countermax
        accQ = 0
        averI = accI / countermax
        accI = 0
        counter = 0
       
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
       
pl.figure(4)
pl.title("I/Q, DC-component removed")
pl.plot(dataI, dataQ, 'ro')
pl.show()

sys.exit('debug')


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 ================================
>>>
filename =  /Temp/DDOUD.csv
skip  0 samples of input-data from datafile
read  1024  samples from file
averI =  510
averQ =  621
averI/n =  0
averQ/n =  0



The first 128 points have the original DC-component and are topright in the figure.
The remaining sample have the DC-component removed, so are upperleft, around (0,0)

Plot of moving average


# binaryAverage03.py
# Feb-2014 Kees de Groot
#
# remove average (DC-component)
# by accumulating 128 samples and divide by 128
# perform some experiments/simulations
# plot average
   
import pylab as pl
import numpy as np
import math
import sys

######################################################################
############# parameter section ######################################
######################################################################

path = '/Temp/'
filename = 'DDOUD.csv'
f = open(path + filename, 'rb')
print "filename = ", path + filename
skip = 0 # 10000 # there is a discontinuity in the file
            # so skip half of it
print "skip ", skip, "samples of input-data from datafile"

###########################################
#############  main loop ##################
###########################################

n = 1024*256 # number of samples to process
print "n = (number of samples to process)", n
countermax = 32 # number of samples to average
counter = 0
dataI = []
dataQ = []
averI = 0.0
averQ = 0.0
accQ = 0.0
accI = 0.0
averageI = []

# 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])
    accQ += Q
    I = int(list[1])
    accI += I
    dataI.append(I - averI)
    dataQ.append(Q - averQ)
    i += 1
    if i >= n:
        break
    counter += 1
    if counter == countermax:
        averQ = accQ / countermax
        accQ = 0
        averI = accI / countermax
        averageI.append(averI)
        accI = 0
        counter = 0
       
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("average I")
pl.plot(averageI)
##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
       
pl.figure(4)
pl.title("I/Q, DC-component removed")
pl.plot(dataI, dataQ, 'ro')
pl.show()

sys.exit('debug')

With output

>>> ================================ RESTART ================================
>>>
filename =  /Temp/DDOUD.csv
skip  0 samples of input-data from datafile
n = (number of samples to process) 262144
read  262144  samples from file
averI =  487
averQ =  597
averI/n =  0
averQ/n =  0



Ah, interesting. Why peaks?
Need some more experiments. At the moment enough interesting things for my blog…

Summary

With a python program I simulated the behaviour of a piece of VHDL-code that I want to try out in the next instalment of this blog. This looks promising but at the end I found some strange peaks in the plot that I have to understand.

As always: if you find errors, if you have suggestions, if you see what I don't see immediately, please comment. 

No comments:

Post a Comment