Costas loop

[hamster_nz]

Hi,

I'm prototyping/building my own GPS receiver in S/W and have hit a bit of trouble.

- I can acquire the GPS space vehicles (work out rough frequency/phase of the L1 signal)
- I can follow the Gold Codes as the move in phase relative to the local osc.
- I can almost recover the BPSK NAV messages

But I am having trouble with recovering the exact carrier frequency - so my I/Q plot looks like a circle (approx 300 units in radius), not the  nice clean

The standard technique is a Costas loop (https://en.wikipedia.org/wiki/Costas_loop), which I've tried implementing many different ways without success. Does anybody have any hints or tips? My questions are mainly about the bandwidth of the arm/loop filters.

[uncle_bob]

Hi

How stable is your local oscillator?

What are you doing to estimate doppler change?

There is nothing crazy about the loop estimation process. You feed the phase in and it effectively runs a PLL that locks up your software tuning process. The bandwidth is adjusted based on where you are in the acquisition process. Wide band is in the 100 Hz range in this case. Locked / narrow band can easily be < 10 Hz. In precision work with a proper doppler estimation and a good local oscillator, you might get below 0.1 Hz on the control loop bandwidth. Your plots look like either the bandwidth is way to high or there is a glitch in your phase estimation.

Bob

[hamster_nz]

Hi Bob

How stable is your local oscillator?

Unsure - I am using a 1-bit ADC sample file off of the interweb while I wait for the real hardware to arrive. However others have successfully decoded it.

What are you doing to estimate doppler change?

At the moment the initial estimate comes from a 10,912 point DTF (2ms of data) around the 4092 MHz IF, giving bins that are 500 Hz wide. Then by looking at the relative power levels in the bins either side it
narrows it down to under 200 Hz. Given that the C/A codes are every 1ms, this should give (and does seem to give) an initial guess that is within 2pi/5 radians (72 degrees) per sample of the ideal lock.

Given the 'openness' of the I/Q plot (except for the BPSK 180 degree changes) it looks like there is plenty of signal, just that I am doing something dumb.

There is nothing crazy about the loop estimation process. You feed the phase in and it effectively runs a PLL that locks up your software tuning process. The bandwidth is adjusted based on where you are in the acquisition process. Wide band is in the 100 Hz range in this case. Locked / narrow band can easily be < 10 Hz. In precision work with a proper doppler estimation and a good local oscillator, you might get below 0.1 Hz on the control loop bandwidth. Your plots look like either the bandwidth is way to high or there is a glitch in your phase estimation.

I know my phase estimation is wrong - at the moment I'm multiplying the I & Q values to find out which quadrant I'm in, and then using that to tweak the 32-bit phase accumulator, but I can't work out what to do in the loop filter. - I'm guessing if filtered(I*Q) > 0 then I need to reduce the phase accumulator's step, if filtered(I*Q) < 0 I adjust it the other way. But the behavior of the filter is what is escaping me - how much bandwidth, how much gain... IIR,.. FIR...

For testing I've taken a very long DFT, and used that to narrow the initial Phase Accumulator step to be as close as possible (< 2Hz off, less than 0.01 radians per sample), so it is start where the Costas loop where it should be in lock, and watched what is happening. It just drifts away...

After a couple of nights playing with it I am thinking I am just doing it wrong - either the gain or bandwidth must be wrong.

[uncle_bob]

Hi

Ok, let's back off a bit:

What you are building is a phase locked loop. That's all it is. Your I/Q and loop process forms a phase detector. That detector has a sensitivity. Simply put, X radians of phase change gives you Y bits of output. The I/Q needs to be demodulated into phase to make this happen. The Arc Tangent function is your friend in this case.

Once you have phase to bits worked out, you have a linear change in bits out as the phase shifts. That gets translated into a "signal" that changes your NCO in order to line it up with the incoming data. You process this bit stream through a very conventional PID (or PI) loop. The bandwidth of this loop can be less ( = a good thing) if you do forward correction on the doppler (it's shifting XXX Hz per sample).

The control loop that drives the NCO is like any other. To much gain and it goes nuts. To little gain and it does not track. Exactly how that is all accomplished depends on your implementation.

Bob

[hamster_nz]

Hi

Ok, let's back off a bit:

What you are building is a phase locked loop. That's all it is. Your I/Q and loop process forms a phase detector. That detector has a sensitivity. Simply put, X radians of phase change gives you Y bits of output. The I/Q needs to be demodulated into phase to make this happen. The Arc Tangent function is your friend in this case.

Once you have phase to bits worked out, you have a linear change in bits out as the phase shifts. That gets translated into a "signal" that changes your NCO in order to line it up with the incoming data. You process this bit stream through a very conventional PID (or PI) loop. The bandwidth of this loop can be less ( = a good thing) if you do forward correction on the doppler (it's shifting XXX Hz per sample).

The control loop that drives the NCO is like any other. To much gain and it goes nuts. To little gain and it does not track. Exactly how that is all accomplished depends on your implementation.

Bob

Thanks for the very clear thoughts!

So to problems I am working over in my head are:

The BSPK gives a flip to the other side of the I/Q graph, which can look to be a very large phase shift if I am naively tracking phase. As GPS NAV is 50 bps, a bit is 20ms long, so these can occur every 20 samples. This throws my phase detector nuts - how do I keep these out of the loop? Something like "if the shift is > pi/2, subtract pi, else if < pi/2 add pi" (which I will now have to try)

Even if I got a stable frequency lock, the key is to get the points pretty much lying on the I axis, so that the I channel encodes the BPSK data. Achieving this requires some way to feed the absolute phase into the loop some how, so it locks down the absolute phase as well.

Also, I'm trying to avoid trig functions so I can implement in programmable logic...

Humm... plenty of thinking (and hopefully learning) going on here...

[uncle_bob]

Hi

Ummm ..... errrrr ....

The heart of a Costas loop is the squaring function followed by a bandpass filter. Once you double the frequency (by squaring), all the phase shift goop goes away ....

Bob

[6thimage]

It might be worth you having a look at 'understanding gps: principles and applications' by Kaplan and Hegarty (you can find a pdf of it online quite easily). In it they list 4 costas loops to use, in order of complexity - the first is IxQ, second is sign(I)xQ, third is Q/I and fourth is atan(Q/I). The fourth is by far the best, but the second is incredibly easy to get working in hardware.

One thing to mention is that a lot of GPS receivers use software for the error calculation and loop filters, and so don't care too much on the time it takes to compute them (or can use large tables). It is possible to do the error calculation and loop filtering in an FPGA, as long as you keep it to fixed points.

Now for some links:

GPS data is a pain to find, but the SETI project have two very strong GPS signals - they only contain single satellites but can be very useful to test your code against http://setiquest.org/wiki/index.php/SetiQuest_Data_Links
This website - http://www.acasper.org/2011/11/07/gps-signal-analysis/ - has some information on how someone else has tracked them

There is some short baseband data available at http://gfix.dk/matlab-gnss-sdr-book/gnss-signal-records/

The only fully open source design I have ever found is Andrew Holme's (it uses a raspberry pi and an FPGA) - http://www.aholme.co.uk/GPS/Main.htm

There is also the piksi - where the software is open sourced (which runs on an arm cortex m4) but the FPGA is closed source - https://github.com/swift-nav/piksi_firmware
Swift-nav (who make the piksi) also have an open sourced software receiver written in python and c (using their own c library, which is also open sourced), but it is quite difficult to follow what part does what - https://github.com/swift-nav/peregrine
In fact, swift-nav's library is worth its own link - https://github.com/swift-nav/libswiftnav

[edit]
I have collected some baseband data from a spirent simulator, but the files are incredibly big (roughly 32 GB per hour, and most of the recordings were for 30 minutes to 2 hours). I've uploaded one of the smaller files I have to mega (never used it before, but they have a large upload size on their free tier) so hopefully it might be useful

data file - https://mega.nz/#!athwHBQB!WTLbYw8jwf0sdxUQ3TQqH5wmHpV5JVPzVHgtyeZcxQE
simulator logs - https://mega.nz/#!6kIUjZpb!GznR4jTWg9Loe36Q2qx3o6e18rSg6-kr0uJCSmA4GGc

The data file is for a test I was performing, so it does not have any antenna attenuation (which is a little unrealistic, as patch antenna's receive very little from the horizon). It is recorded at 16.368 MHz with a centre frequency of 4.092 MHz. The first 512 bytes of the data file are blank (due to the way I designed the recorder - SD card connected to an FPGA), the data is 3 bits (without any complex data) and is stored in the top of each nibble. Below is a c function I use to load the data (my programs, separately to this function, open the file and fseek 512 bytes):

Code: [Select]

/* nibble 3 data format
 * two samples, stored in nibbles of byte, with no imaginary data
 *
 *     76543210
 *     smm*SMM*
 *
 * s - sign, m - magnitude, * - unused, cases indicate sample
 *
 * 7-5 first sample
 * 3-1 second sample
 */
static double n3_sample_table_double[]={ 1., /* 000 */  3., /* 001 */  5., /* 010 */  7., /* 011 */
                                        -1., /* 100 */ -3., /* 101 */ -5., /* 110 */ -7.  /* 111 */
                                       };
int32_t data_n3_read(FILE *data_file, uint32_t num_samples, uint32_t pad, double complex **data,
                     data_malloc dmalloc)
{
    size_t bytes_to_read=num_samples/2, data_counter=0;

    /* malloc the data */
    *data=(double complex*)dmalloc(sizeof(double complex)*(num_samples+pad));
    if(!*data)
        return 0;

    /* read data */
    while(bytes_to_read)
    {
        size_t i, block_size=4096, bytes_read;
        uint8_t data_block[4096];

        if(bytes_to_read<block_size)
            block_size=bytes_to_read;

        /* perform read */
        bytes_read=fread(data_block, 1, block_size, data_file);
        if(bytes_read<1)
            return -(int32_t)data_counter;

        /* store data */
        for(i=0; i<bytes_read; ++i, data_counter+=2)
        {
            (*data)[data_counter]=(1. + I)*n3_sample_table_double[(data_block[i]>>5)&0x7];
            (*data)[data_counter+1]=(1. + I)*n3_sample_table_double[(data_block[i]>>1)&0x7];
        }

        /* update counter */
        bytes_to_read-=bytes_read;
    }

    /* pad data */
    for(; data_counter<(num_samples+pad); ++data_counter)
        (*data)[data_counter]=0.;

    return num_samples+pad;
}

Ignore the pad and dmalloc - they are there so that I could use it with the fftw library (which uses its own aligned version of malloc), with the padding being so I could use a 2^n length FFT.

[hamster_nz]

It might be worth you having a look at 'understanding gps: principles and applications' by Kaplan and Hegarty (you can find a pdf of it online quite easily). In it they list 4 costas loops to use, in order of complexity - the first is IxQ, second is sign(I)xQ, third is Q/I and fourth is atan(Q/I). The fourth is by far the best, but the second is incredibly easy to get working in hardware.

Great stuff - where is the "Thanks a million" button?

I've been experimenting with the data file from http://www.jks.com/gps/gps.html, which is about 80 seconds of single-bit data, with an 5.456MHz sample rate, with a LO of 4.092MHz (yes it aliases down to a LO of 1,364 MHz

I think sign(I)*Q is what I've been looking for....

Being digital I can also adjust the LO phase retrospectively,

[hamster_nz]

I got somewhere eventually, using a Frequency Locked Loop.

For me, the key was to reduce the I/Q into the -3 to +3 range (just by left shifting, then use the following lookup table to convert to an angle (scaled to 0-255)

Code: [Select]

unsigned char angle_lookup[64] =
         
          { // i = -4  -3  -2  -1   0   1   2   3   
                    0,  0,  0,  0,  0,  0,  0,  0,  // q = -4
                    0,160,168,178,192,206,216,224,  //     -3
                    0,152,160,173,192,211,224,232,  //     -2
                    0,142,147,160,192,224,237,242,  //     -1
                    0,128,128,128,  0,  0,  0,  0,  //     -0
                    0,114,109, 96, 64, 32, 19, 14,  //      1
                    0,104, 96, 83, 64, 45, 32, 24,  //      2
                    0, 96, 88, 78, 64, 50, 40, 32}; //      3

I'm still trying to get the Costas filter tracking to eliminate the approx 2Hz drift, but I can at least strip out the NAV messages.

[hamster_nz]

Might not have got my Costas loop to work, but am starting to extract orbital and time data from the signal:

Code: [Select]

=========== SV 29 ==============

Orbit parameters:
iode    31
Cfb     -0.000001522
delta_n 12572
M0         61.616402
Cuc      0.000000000
e        0.002519569
Cus      0.000009390
sqrt(A)  5217.634533 m^1/2
A       27223710.119208 m
Toe     7242
Fit     0
aodo    0F
=========== SV 30 ==============

Orbit parameters:
iode    3A
Cfb     -0.000004211
delta_n 13465
M0         16.222033
Cuc      0.000000000
e        0.013084712
Cus      0.000006450
sqrt(A)  5219.346394 m^1/2
A       27241576.776231 m
Toe     7242
Fit     0
aodo    0C
=========== SV 31 ==============

Orbit parameters:
iode    13
Cfb      0.000000227
delta_n 11983
M0        -52.977987
Cuc      0.000000000
e        0.007557237
Cus      0.000004234
sqrt(A)  5153.800608 m^1/2
A       26561660.703736 m
Toe     7241
Fit     0
aodo    0E

I guess I now need to learn a bit of Orbital Mechanics...

[hamster_nz]

It has taken a month of Sundays, but I learned a bit of Orbital Mechanics and can now go from a file of raw 1-bit ADC samples to a reasonable position fix. Considering that the Space Vehicles orbit at 27,200km I am dead chuffed with myself.

Code is is about 1600 lines of C, using only standard libraries. Now to tweak, tidy, add maybe add ionosphere corrections, then convert the heavy bits to hardware.

But man, the precision and accuracy makes it a real pain to debug this stuff. One little error and you are stuffed. I spent best part of a week tracking down where I had extracted the wrong number of bits for a clock correction constant...

[uncle_bob]

Hi

You are only just at the start

How about SBAS?

Bob

[hamster_nz]

Hi

You are only just at the start

How about SBAS?

Bob

Being in the darkest reaches of the South Pacific I don't think that there is any SBAS service out here to make it useful!

[uncle_bob]

Hi

You are only just at the start

How about SBAS?

Bob

Being in the darkest reaches of the South Pacific I don't think that there is any SBAS service out here to make it useful!

Hi

Making getting it working even *more* complicated

Bob