Finally I can generate a fresh new blinking warm hostedx115.rbf from the bladerf hdl sources!

Brian was so kind to answer my last "mayday, mayday, mayday" message.

To keep you informed I just copy and paste the last info from my private-diary:

Get a stable starting-point from scratch

Get the latest bladerf software.

Nuand.com, support, alle relevant … and end at https://github.com/Nuand/bladeRF

At right of the page is download zip (yes, I know, I should clone, but don’t yet know how, sorry…)

bladeRF-master (3).zip is in my download-area.

Rename old bladeRF-master to bladeRF-master03

Extract to C:\bladeRF-master

May be wise to get the latest firmware too.

The latest image is pointed to by: http://www.nuand.com/fx3/latest.img

Download latest.img and copy to C:\Program Files (x86)\bladeRF for loading with bladerf-cli

To be sure download FPGA-firmware from http://nuand.com/fpga/ choose V0.0.3 hostedx115.rbf

Rename current hostedx115.rbf to hostedx115.rbf_01

Now copy hostedx115.rbf to directory C:\Program Files (x86)\bladeRF

Now connect bladerf-board using USB2-port, and flash latest.img, power remove/connect, flash with hostedx115.rbf

Connect power (I use external poweradaptor), connect to USB2

One green LED D1 shows a nice constant green light.

Startup bladerf-cli (I use windows by the way)

[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 fx3 latest.img

Flashing firmware from latest.img...

[INFO] Erasing 0x20000 bytes starting at address 0x00

[INFO] Erased sector at 0x00...

[INFO] Erased sector at 0x10000...

[INFO] Writing 0x20000 bytes starting at address 0x00

[INFO] Verifying 0x20000 bytes starting at address 0x00

Done. Cycle power on the device.

bladeRF>

Remove/reconnect power

bladeRF> probe

    Backend:        libusb

    Serial:         b436de8c8212b9aeaaeba852246866e7

    USB Bus:        2

    USB Address:    6

bladeRF> version

  bladeRF-cli version:        0.7.0-git

  libbladeRF version:         0.9.0-git

[ERROR] status < 0: LIBUSB_ERROR_IO

Error: File or device I/O failure

bladeRF> exit

oeps, error, disconnect/reconnect USB and restart bladeRF

[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>

YEP! LED2 blinks about 40 times in 10 seconds, LED1 is off and LED3 is on.

Now try to create hostedx115.rbf from scratch

From this stable situation I am trying to create a new hostedx115.rbf with Quartus.

Little bit confusion about the LEDs

In the file C:\bladeRF-master\hdl\fpga\platforms\bladerf\vhdl\ bladerf-hosted.vhd:

…

                count := 100_000_00 ;

                led(1) <= not led(1) ;

            end if ;

        end if ;

    end process ;

    led(2) <= tx_underflow_led ;

    led(3) <= rx_overflow_led ;

…

I would expect led(1) to blink, not LED2

But, I might presume that LED1 = led(1), LED2 = led(2) and LED3 = led(3) but is that true?

The file C:\bladeRF-master\hdl\fpga\platforms\bladerf\constraints\pins.tcl contains:

set_location_assignment PIN_AA7 -to led[1]

set_location_assignment PIN_AB7 -to led[2]

…

set_location_assignment PIN_AB10 -to led[3]

Now to the schematics to find:

Sheet 3 of 14 FPGA “LEFT” BANK

LED1 = D11 AA7

LED2 = D12 AB7

AB10 is not connected on the schematic

On the bladeRF-board I read

LED1 = D12 (left)

LED2 = D11 (middle)

LED3 = D13 (right)

So, the middle blinking LED must be led(1)

I guess that the left LED must be led(2)

I’ll try to change the blinking rate of led(1) and set led(2) to the on state ‘1’

But firstly try to get a home-made hostedx115.rbf without any changes.

In Windows, you should be able to open the Nios II Command Shell (look in your Start menu). You can then navigate to the hdl/quartus directory and use the build_bladerf.sh script to build your project.

You can open the resulting project using the GUI in the work directory that gets created by the script.

Also note there is a README in the HDL directory to help with compilation. We currently use Quartus II 13.1 for building. Please make sure you are up to date.

I keep an eye on your blog and love what you're doing. Please keep it up! :)

Posted by Brian to SDR with BladeRF at February 1, 2014 at 3:43 AM

Start, Altera 13.1.0.162 Web Edition, NIOS II EDS 13.1.0.162, NIOS II 13.1 Command Shell

------------------------------------------------

Altera Nios2 Command Shell [GCC 4]

Version 13.1, Build 162

------------------------------------------------

lonneke@lonneke-VAIO /cygdrive/c/altera/13.1

$ cd ..

lonneke@lonneke-VAIO /cygdrive/c/altera

$ dir

13.0sp1  13.1

lonneke@lonneke-VAIO /cygdrive/c/altera

$ cd ..

lonneke@lonneke-VAIO /cygdrive/c

$ dir

$Recycle.Bin              END         Music                  Quartus       System\ Volume\ Information  WirelessDiagLog.csv  default_xdb   magazines     version

BBHuur                    Elektor     POI                    RHDSetup.log  Temp                         Xilinx               eclipse       pagefile.sys  vmware

Backup\ docs\ iPad        Firefox     PerfLogs               SPLASH.000    Update                       _FS_SWRINFO          fftrlog.txt   prive

CampingWijzer             ISOimages   Photo                  SPLASH.SYS    Users                        altera               hiberfil.sys  splash.idx

Cypress                   Intel       Program\ Files         Sdr           Utrecht                      bladeRF-master       homepage      t

DesignSparkMech           Keil        Program\ Files\ (x86)  SiLabs        VAIO\ Sample\ Contents       bladeRF-master_01    java          test.xml

Documentation             LinuxBoard  ProgramData            Softw\ Radio  WinRadio                     bladeRF-master_02    javaprojects  usb_driver

Documents\ and\ Settings  MSOCache    Python27               SwSetup       Windows                      bladeRF-master_03    lv.log        user.js

lonneke@lonneke-VAIO /cygdrive/c

$ cd bladerf-master

lonneke@lonneke-VAIO /cygdrive/c/bladerf-master

$ dir

CMakeLists.txt  CONTRIBUTORS  COPYING  README.md  debian  firmware_common  fx3_firmware  hdl  host  legal

lonneke@lonneke-VAIO /cygdrive/c/bladerf-master

$ cd hdl

lonneke@lonneke-VAIO /cygdrive/c/bladerf-master/hdl

$ dir

README.md  fpga  quartus

lonneke@lonneke-VAIO /cygdrive/c/bladerf-master/hdl

$ cd quartus

lonneke@lonneke-VAIO /cygdrive/c/bladerf-master/hdl/quartus

$ dir

bladerf.tcl  build.tcl  build_bladerf.sh  constraints  ip.ipx  signaltap

lonneke@lonneke-VAIO /cygdrive/c/bladerf-master/hdl/quartus

$ bash build_bladerf.sh -h

bladeRF FPGA build script

Usage: build_bladerf.sh -r <rev> -s <size>

Options:

    -r <rev>       FPGA revision

    -s <size>      FPGA size

    -a <stp>       SignalTap STP file

    -h             Show this text

Supported revisions:

    hosted

    headless

    fsk_bridge

    qpsk_tx

Supported sizes (kLE)

    40

    115

lonneke@lonneke-VAIO /cygdrive/c/bladerf-master/hdl/quartus

$ bash build_bladerf.sh -r hosted -s 115

Press RETURN or ENTER and go out with the dog for a short walk or go and listen to your partner

An enormous amount of information appears in the window, after a few minutes the last few lines read:

    Info: Elapsed time: 00:00:14

    Info: Total CPU time (on all processors): 00:00:12

Info (23030): Evaluation of Tcl script ../build.tcl was successful

Info: Quartus II 32-bit Shell was successful. 0 errors, 89 warnings

    Info: Peak virtual memory: 195 megabytes

    Info: Processing ended: Sun Feb 02 23:03:26 2014

    Info: Elapsed time: 00:06:05

    Info: Total CPU time (on all processors): 00:00:02

/cygdrive/c/bladerf-master/hdl/quartus

##########################################################################

    Done! Image copied to: hostedx115.rbf

##########################################################################

lonneke@lonneke-VAIO /cygdrive/c/bladerf-master/hdl/quartus

$

There it is: hostedx115.rbf, I rename the file to hostedx115new.rbf and copy it to C:\Program Files (x86)\bladeRF

The proof of the pudding is in the eating:

Everything is still fine on my bladerf-board, the LEDs are still blinking. SDR-RADIO.COM SDR radio console works like a charm.

Now edit a file with the GUI of Quartus

You can open the resulting project using the GUI in the work directory that gets created by the script.

Double-click on C:\bladeRF-master\hdl\quartus\work\ bladerf.qpf

The GUI of Quartus opens!!!

I edit the file ../../fpga/platforms/bladerf-hosted.vhd to read:

    toggle_led1 : process(fx3_pclk)

        variable count : natural range 0 to 100_000_000 := 100_000_000 ;

    begin

        if( rising_edge(fx3_pclk) ) then

            count := count - 1 ;

            if( count = 0 ) then

--                count := 100_000_00 ;

                count := 100_000_000 ;                                -- KdG 2-Feb-2014                                        

                led(1) <= not led(1) ;

            end if ;

        end if ;

    end process ;

--    led(2) <= tx_underflow_led ;

    led(2) = '1';                             -- KdG 2-Feb-2014

    led(3) <= rx_overflow_led ;

--    toggle_led2 : process(rx_clock)

--        variable count : natural range 0 to 38_400_00 := 38_400_00 ;

--    begin

--        if( rising_edge(rx_clock) ) then

--            count := count - 1 ;

--            if( count = 0 ) then

--                count := 38_400_00 ;

--                led(2) <= not led(2) ;

--            end if ;

--        end if ;

--    end process ;

Now try the wholce circus again after having saved this file.

Just click on save all and exit

What a pity, at the end I get:

2014.02.02.23:31:37 Progress: Done reading input file

2014.02.02.23:31:38 Info: nios_system: Generating nios_system "nios_system" for QUARTUS_SYNTH

2014.02.02.23:31:38 Info: pipeline_bridge_swap_transform: After transform: 12 modules, 41 connections

2014.02.02.23:31:38 Info: No custom instruction connections, skipping transform

2014.02.02.23:31:38 Info: merlin_initial_interconnect_transform: After transform: 12 modules, 35 connections

2014.02.02.23:31:39 Info: merlin_translator_transform: After transform: 25 modules, 74 connections

2014.02.02.23:31:39 Info: merlin_domain_transform: After transform: 50 modules, 207 connections

2014.02.02.23:31:39 Info: merlin_router_transform: After transform: 63 modules, 246 connections

2014.02.02.23:31:39 Info: merlin_burst_transform: After transform: 64 modules, 249 connections

2014.02.02.23:31:40 Info: merlin_network_to_switch_transform: After transform: 89 modules, 301 connections

2014.02.02.23:31:40 Info: merlin_width_transform: After transform: 91 modules, 307 connections

2014.02.02.23:31:40 Info: merlin_clock_and_reset_bridge_transform: After transform: 93 modules, 399 connections

2014.02.02.23:31:41 Info: merlin_hierarchy_transform: After transform: 13 modules, 43 connections

2014.02.02.23:31:41 Info: merlin_mm_transform: After transform: 13 modules, 43 connections

2014.02.02.23:31:41 Info: merlin_interrupt_mapper_transform: After transform: 14 modules, 46 connections

2014.02.02.23:31:41 Info: reset_adaptation_transform: After transform: 16 modules, 60 connections

2014.02.02.23:31:43 Warning: nios_system: "No matching role found for uart_0:s1:dataavailable (dataavailable)"

2014.02.02.23:31:43 Warning: nios_system: "No matching role found for uart_0:s1:readyfordata (readyfordata)"

2014.02.02.23:31:44 Error: Couldn't write C:/bladeRF-master/hdl/fpga/ip/altera/nios_system/synthesis/nios_system.v

2014.02.02.23:31:44 Error: Generation stopped, 14 or more modules remaining

2014.02.02.23:31:44 Info: nios_system: Done "nios_system" with 15 modules, 1 files, 53301 bytes

lonneke@lonneke-VAIO /cygdrive/c/bladerf-master/hdl/quartus

$

Perhaps I should do some cleaning?

From C:\bladeRF-master\hdl\quartus delete directory work and try again

Same problem:

------------------------------------------------

Altera Nios2 Command Shell [GCC 4]

Version 13.1, Build 162

------------------------------------------------

lonneke@lonneke-VAIO /cygdrive/c/altera/13.1

$ cd ..

lonneke@lonneke-VAIO /cygdrive/c/altera

$ cd ..

lonneke@lonneke-VAIO /cygdrive/c

$ cd bladerrf-master/hdl/quartus

bash: cd: bladerrf-master/hdl/quartus: No such file or directory

lonneke@lonneke-VAIO /cygdrive/c

$ cd bladerf-master/hdl/quartus

lonneke@lonneke-VAIO /cygdrive/c/bladerf-master/hdl/quartus

$ dir

bladerf.tcl  build.tcl  build_bladerf.sh  constraints  ip.ipx  signaltap

lonneke@lonneke-VAIO /cygdrive/c/bladerf-master/hdl/quartus

$ bash build_bladerf.sh -r hosted -s 115

## Non-Linux OS Detected (Windows?)

## Forcing QUARTUS_BINDIR to C:/altera/13.1/quartus/bin

##########################################################################

    Generating NIOS II Qsys for bladeRF ...

##########################################################################

2014.02.02.23:38:28 Progress: Loading nios_system/nios_system.qsys

2014.02.02.23:38:28 Progress: Reading input file

2014.02.02.23:38:28 Progress: Adding clk_0 [clock_source 13.1]

2014.02.02.23:38:29 Progress: Parameterizing module clk_0

2014.02.02.23:38:29 Progress: Adding nios2_qsys_0 [altera_nios2_qsys 13.1]

2014.02.02.23:38:29 Progress: Parameterizing module nios2_qsys_0

2014.02.02.23:38:29 Progress: Adding spi_0 [altera_avalon_spi 13.1]

2014.02.02.23:38:29 Progress: Parameterizing module spi_0

2014.02.02.23:38:29 Progress: Adding jtag_uart_0 [altera_avalon_jtag_uart 13.1]

2014.02.02.23:38:29 Progress: Parameterizing module jtag_uart_0

2014.02.02.23:38:29 Progress: Adding onchip_memory2_0 [altera_avalon_onchip_memory2 13.1]

2014.02.02.23:38:29 Progress: Parameterizing module onchip_memory2_0

2014.02.02.23:38:29 Progress: Adding timer_0 [altera_avalon_timer 13.1]

2014.02.02.23:38:29 Progress: Parameterizing module timer_0

2014.02.02.23:38:29 Progress: Adding uart_0 [altera_avalon_uart 13.1]

2014.02.02.23:38:29 Progress: Parameterizing module uart_0

2014.02.02.23:38:29 Progress: Adding spi_1 [altera_avalon_spi 13.1]

2014.02.02.23:38:29 Progress: Parameterizing module spi_1

2014.02.02.23:38:29 Progress: Adding pio_0 [altera_avalon_pio 13.1]

2014.02.02.23:38:29 Progress: Parameterizing module pio_0

2014.02.02.23:38:29 Progress: Adding iq_corr_rx_phase_gain [altera_avalon_pio 13.1]

2014.02.02.23:38:29 Progress: Parameterizing module iq_corr_rx_phase_gain

2014.02.02.23:38:29 Progress: Adding iq_corr_tx_phase_gain [altera_avalon_pio 13.1]

2014.02.02.23:38:29 Progress: Parameterizing module iq_corr_tx_phase_gain

2014.02.02.23:38:29 Progress: Adding bladerf_oc_i2c_master_0 [bladerf_oc_i2c_master 1.0]

2014.02.02.23:38:30 Progress: Parameterizing module bladerf_oc_i2c_master_0

2014.02.02.23:38:30 Progress: Building connections

2014.02.02.23:38:30 Progress: Parameterizing connections

2014.02.02.23:38:30 Progress: Validating

2014.02.02.23:38:30 Progress: Done reading input file

2014.02.02.23:38:34 Progress: Loading nios_system/nios_system.qsys

2014.02.02.23:38:34 Progress: Reading input file

2014.02.02.23:38:34 Progress: Adding clk_0 [clock_source 13.1]

2014.02.02.23:38:34 Progress: Parameterizing module clk_0

2014.02.02.23:38:34 Progress: Adding nios2_qsys_0 [altera_nios2_qsys 13.1]

2014.02.02.23:38:35 Progress: Parameterizing module nios2_qsys_0

2014.02.02.23:38:35 Progress: Adding spi_0 [altera_avalon_spi 13.1]

2014.02.02.23:38:35 Progress: Parameterizing module spi_0

2014.02.02.23:38:35 Progress: Adding jtag_uart_0 [altera_avalon_jtag_uart 13.1]

2014.02.02.23:38:35 Progress: Parameterizing module jtag_uart_0

2014.02.02.23:38:35 Progress: Adding onchip_memory2_0 [altera_avalon_onchip_memory2 13.1]

2014.02.02.23:38:35 Progress: Parameterizing module onchip_memory2_0

2014.02.02.23:38:35 Progress: Adding timer_0 [altera_avalon_timer 13.1]

2014.02.02.23:38:35 Progress: Parameterizing module timer_0

2014.02.02.23:38:35 Progress: Adding uart_0 [altera_avalon_uart 13.1]

2014.02.02.23:38:35 Progress: Parameterizing module uart_0

2014.02.02.23:38:35 Progress: Adding spi_1 [altera_avalon_spi 13.1]

2014.02.02.23:38:35 Progress: Parameterizing module spi_1

2014.02.02.23:38:35 Progress: Adding pio_0 [altera_avalon_pio 13.1]

2014.02.02.23:38:35 Progress: Parameterizing module pio_0

2014.02.02.23:38:35 Progress: Adding iq_corr_rx_phase_gain [altera_avalon_pio 13.1]

2014.02.02.23:38:35 Progress: Parameterizing module iq_corr_rx_phase_gain

2014.02.02.23:38:35 Progress: Adding iq_corr_tx_phase_gain [altera_avalon_pio 13.1]

2014.02.02.23:38:35 Progress: Parameterizing module iq_corr_tx_phase_gain

2014.02.02.23:38:35 Progress: Adding bladerf_oc_i2c_master_0 [bladerf_oc_i2c_master 1.0]

2014.02.02.23:38:35 Progress: Parameterizing module bladerf_oc_i2c_master_0

2014.02.02.23:38:35 Progress: Building connections

2014.02.02.23:38:35 Progress: Parameterizing connections

2014.02.02.23:38:35 Progress: Validating

2014.02.02.23:38:36 Progress: Done reading input file

2014.02.02.23:38:37 Info: nios_system: Generating nios_system "nios_system" for QUARTUS_SYNTH

2014.02.02.23:38:37 Info: pipeline_bridge_swap_transform: After transform: 12 modules, 41 connections

2014.02.02.23:38:37 Info: No custom instruction connections, skipping transform

2014.02.02.23:38:37 Info: merlin_initial_interconnect_transform: After transform: 12 modules, 35 connections

2014.02.02.23:38:38 Info: merlin_translator_transform: After transform: 25 modules, 74 connections

2014.02.02.23:38:38 Info: merlin_domain_transform: After transform: 50 modules, 207 connections

2014.02.02.23:38:38 Info: merlin_router_transform: After transform: 63 modules, 246 connections

2014.02.02.23:38:39 Info: merlin_burst_transform: After transform: 64 modules, 249 connections

2014.02.02.23:38:39 Info: merlin_network_to_switch_transform: After transform: 89 modules, 301 connections

2014.02.02.23:38:39 Info: merlin_width_transform: After transform: 91 modules, 307 connections

2014.02.02.23:38:39 Info: merlin_clock_and_reset_bridge_transform: After transform: 93 modules, 399 connections

2014.02.02.23:38:40 Info: merlin_hierarchy_transform: After transform: 13 modules, 43 connections

2014.02.02.23:38:40 Info: merlin_mm_transform: After transform: 13 modules, 43 connections

2014.02.02.23:38:40 Info: merlin_interrupt_mapper_transform: After transform: 14 modules, 46 connections

2014.02.02.23:38:41 Info: reset_adaptation_transform: After transform: 16 modules, 60 connections

2014.02.02.23:38:42 Warning: nios_system: "No matching role found for uart_0:s1:dataavailable (dataavailable)"

2014.02.02.23:38:42 Warning: nios_system: "No matching role found for uart_0:s1:readyfordata (readyfordata)"

2014.02.02.23:38:43 Error: Couldn't write C:/bladeRF-master/hdl/fpga/ip/altera/nios_system/synthesis/nios_system.v

2014.02.02.23:38:43 Error: Generation stopped, 14 or more modules remaining

2014.02.02.23:38:43 Info: nios_system: Done "nios_system" with 15 modules, 1 files, 53301 bytes

lonneke@lonneke-VAIO /cygdrive/c/bladerf-master/hdl/quartus

$

 Summary

For me all this was a giant step, perhaps a small step for Brian.

I created a stable start-point.

I managed to generate hostedx115.rbf using the NIOS shell (note the use of the command bash).

I managed to create a Quartus project and could use the GUI.

I managed to change bladerf-hosted.vhdl.

But I did not manage to compile again.

I should do some cleaning perhaps, or is there a command from the GUI?

Anyway, I had a great great evening this Sunday and now I am going for a walk with my dog and get some sleep.

He is not so smart to tell me what he wants by

waiting in the bathroom, I put him there for the picture...

Programming the FPGA of the bladeRFis not so simple

Programming the FPGA of bladeRF

The last few weeks I am trying to program the FPGA.
With the marvelous software Quartus of Altera I managed to make some *.rbf files and transfer them to the FPGA. Well, first I had .sof files but they can be transformed into .rbf files.

I read many manuals, performed some courses from Altera (very instructive by the way!) and experimented a lot. I used some files from the VHDL-bladeRF-subdirectory.

However, I did not manage to blink the LEDs.
Tonight I finally saw why.

C4_CLK running at 38.4 MHz

There is a clock c4_clk running at 38.4 MHz. I tried to reduce that clock with a counter to about 1 Hz. Yes: the blinking LED, the equivalent of Hello World!

I did not have a clock! There was no 38.4 MH signal. Why not?

The clock-signal is derived from a chip, the Si5338. That chip has to be initialised, programmed, setup.

So, there is not a simple way to use the FPGA in an easy way. With easy I mean, a few simple files that grow in time, become more complicated as I understand more and more of the system.

Nuand: use revisions

OK, there is no simple way to experiment with the FPGA. What I want to accomplish is:

• get the rf-output from the LMS-chip
• digitize that signal
• make some hardware with VHDL (correlator, Costas-loop)
• downconvert the signal
• send to the laptop with USB
• do the final processing on the laptop

The FPGA firmware is compiled without the Quartus GUI. I think the nuand-people do that with linux. I want to use Quartus with the GUI in Windows 7 on my laptop.

I have to find out how to tell Quartus to use the files in the master-bladerf-directory.
Nuand advises to use the standard bladerf software and make revisions to that software. Yes, but how to compile the FPGA firmware?

There is no easy way

So, I will do a lot more reading, some more Quartus tutorials and then figure out how to set up a bladerf-Windows-Quartus-project resulting in a nice hostedx115.rbf file!

Please help me

I am not desperate yet, I have to study a lot. But I would love to get some hints, some help, some advice! In the meantime, I will read, think, experiment, and blog about my findings.

Finally I can see GPS-satellites with bladeRF!

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!

bladeRF-car designed with DesignSpark Mechanical

BladeRF-car

A bladerf-car? Yes, somehow. In between my bladerf-project I played a lot with DesignSpark Mechanical, free 3D software.  And now in Holland there are fablabs. What is fablab? Fabrication-laboratory. A workshop with machines like 3D-printers and a laser cutter. Somewhere in the corner is a reprap busy reproducing itself… I am excited.

In the beginning of my career, before being a retired engineer, I was at the start of the computer age. My first job was maintenance of a PDP8 minicomputer. Later, many years later, there was my first PC with a 10 MB hard disk! 10 MB. Incredibly that an average hard disk now is at least 1TB, that is 1000 GB and 1 GB is 100 * 10 MB! So an average hard disk in a laptop (!) has a capacity of 100.000 times my first 10 MB hard disk. The same with memory, the same with transfer speed, the same with…

These times are incredible times. A few years ago I saw the first regular use of LEDs. Now I have a little pocket flashlight with an enormous amount of light. Ok, my old big maglite with 4 D-cells is still in my car if I need to impress someone. The 4 D-cells in an aluminium case impress more than the light of the small lamp inside that maglite however…

There is a new phenomenon: 3D printing. It is not new at all, but there is a breakthrough now. I feel that again I am at the start of something new, something that will change the world.

So, the last few weeks I was playing around with DesignSpark Mechanical. It is free 3D software. With it you can draw 3D objects and produce .STL files. Those files contain the info that a 3D printer needs to 3D print what you did draw with your 3D software. You can upload the file to a printer shop and some time later you get the stuff in real hardware in your postbox. Or you buy or build a 3D printer and print your creative work at home!

Have a look at http://www.thingiverse.com and see what others designed!

About 35 years ago I designed a car made from plywood. It was a kind of puzzle. But once you slide all the pieces in place you have a car to play with. I redesigned that same car with DesignSpark Mechanical. Just for fun and just because you can. I think I am gonna take a piece of Trespa and take it to the fablab. Then I might use a laser cutter to make this car into a real thing.

The car:

And a more or less exploded view:

How to put this all together?

First you take the left side (linkerzijde)

Now first take the bakdeksel (cover of the the rear area) and stick it in the left side:

Not take the right side (rechterzijde) and put it in place:

 Now is the time for the front plate (frontschot), you slide it upwards in place:

At last: the bottom plate (bodemplaat). This is a little tricky: you have to place it at the underside:

 of the car, push it upwards:

 and then slide it forwards. The result:

Or in the original perspective:

Now we can slide the window (voorruit) on its place:

And the roof (cabinedak) with result:

Finally the last two planes called tussenschot and achterschot:

As you can see these last two planes lock the bottom plate, and finally, there it is, the bladerf-car:

Some remarks:
• I designed the car from a plate with a thickness of 6 mm
• In DesignSpark Mechanical I named the various parts, so they can be separately designed
• I used mostly Select and Pull during the design
• It is very handy to select a surface and then choose Plan View to get a plane to draw upon
• With the pull tool you can remove parts by dragging backwards
• If you want to design something like this: study, do and perform all the tutorials and videos you can find on the Internet
it is really difficult at the start, but at the end you draw with ease and confidence!
• Use points, lines, rectangles and construction lines: you can delete them selectively when not needed anymore
• This car is a very simple example with flat planes, DesignSpark Mecanical can do a lot more!
• I did not laser cut or 3D-print this design. I expect some problems: perhaps the gaps should be wider, a fraction of a millimeter. And may be there are errors in this design…
• So, not any warranty, but don’t be afraid: do try this at home yourself!
• I am in no way involved with DesignSpark Mechanical, the software is free and I am just a happy old retired married satisfied customer, really!

The example above was an attempt to design something (the blue ring) to connect a circular LED light (the yellow thing) to the rectangular frame of a microscope (the orangish thing)
I wanted to 3D print the blue adapter but it proved to be a lot easier just to take a piece of plywood and a jigsaw…

Epilogue

At the end of the year, more than ever, I like to look back in (my) history. Because of my age there is more past than future! But the present is interesting enough, so I tend to stay in the present every day as long as I live. The last few weeks I was distracted from my dear bladerf-project because of DesignSpark Mechanical. I took that side path and enjoyed it a lot. I hope you enjoyed it too!
If you want the STL-files of the design of the bladerf-car: just send me an email and I will be more than happy to send it to you.

Happy Christmas, Happy New Year, Happy Whatever!

Que le vaya bien!

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!

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

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.