Implementation and testing of a GNSS system consisting of a RF front-end and a software GNSS receiver

Table of contents

Abstract

Thanks

1. Introduction and targets

2 Theory
2.1 Software-Defined Radio - SDR
2.2 Software GNSS Receiver - SGR
2.2.1 Spectral analysis
2.2.2 Acquisition of GPS C/A signals
2.2.3 Tracking of GPS C/A signals
2.2.4 Navigation message and PVT calculations

3 Methods and results
3.1 Analysis and documentation of the code
3.1.1 M2HTML and GraphViz for Borre’s Matlab code
3.1.2 Adaptions of “initSettings.m” before the start of programs
3.1.3 Call sequence in the Borre code
3.1.4 NSL Stereo hardware
3.1.5 NSL Stereo C-programs
3.1.6 NSL Stereo Matlab programs
3.2 Solutions for missing functionality, quality, and performance
3.2.1 Determination of receiver velocities
3.2.2 Signal quality C/No and stopping tracking of bad channels
3.2.3 Handling of records of more than 37 seconds
3.2.4 Statistics and quality of navigation results
3.2.5 Performance improvements
3.3 New signals (non GPS L1)

4 Evaluation & Outlook
4.1 Functionality and speed of new Matlab Code
4.2 Interpretation of NSL Stereo Measurement of non GPS L1 signals
4.3 Proposals for future work

Appendix
A1 Literature
A2 Figures and tables
A3 Abbreviations
A4 DVD contents
A5 List of functions, structures, and files of current Borre code
A6 NSL Stereo handbook

Abstract

An introduction into the theory of software defined receivers and especially in such for detecting GNSS signals, acquiring and tracking GNSS satellites, calculating pseudo ranges, positions, velocity and time (PVT) is presented.

Basis of the practical work was the open source project SoftGPS, programmed in Matlab and published by (Borre 2007). The Radio Frequency front end (RF-FE) used in this project was no longer available and was replaced by one with different behavior: NSL Stereo (amplifier, mixer, sampler, and A/D converter in two chains). Adaptations, corrections and extensions to the Matlab code were ne­ces­sary to work with the new front end and to get new functions.

With Stereo came also new Matlab- and C/C++ code that did not work properly. Parallel to the projected working environment – Ubuntu 16.04 Linux with Matlab 2016a – also Windows 10-64bit and a Windows XP-64bit beta-soft­ware from NSL from January 2013 had to be used due to long delays at NSL to provide updated / wor­king Linux ver­sions: the original software from 2012 for Ubuntu 10 was not working in any newer Linux distribution. Finally a version for Ubuntu 14.04-64bit from Jan 2016 was provided after most of the grabbing of different GNSS-signals was already done.

Code of (Borre 2007) and of NSL for Stereo RF-FE were thoroughly analyzed and documented. Besides own descriptions also the M2HTML documentation generator and GraphViz (for generating dependency graphs) were used. The software was also changed and expanded to archive demands for more modularity, performance, quality and functionality (C/No calculation, output of correct velocities in UTM coordinates, statistics about positions and velocities, continuous processing, ...). As code release tool, Git was used for a complete change history and to be able to recover old versions of the code. With the Git-Bash, identical (UNIX-like) behavior was achieved on both Linux and Windows platforms. Git is more modern than the system used in (Borre 2007) and integrated in Matlab. Even with only 4 parallel processes (in a notebook) and a processing con­ditioned by signal to noise ratios C/No the most time consuming tracking was reduced to about a quarter of the initial processing time.

Optional tasks were to change the software in the direction of real time processing and utilization of a 2nd channel of Stereo RF-FE to receive a GNSS signal other than GPS L1. The continuous processing simu­lating real time operation was achieved by new additional routines for tracking and navigation after evaluating a complete frame during the first 37 seconds, making the tracking functions auto­no­mous (e. g. with complete own file handling) and storing results of every step and feeding them as start va­lues into the next step. The last optional task was canceled as many tests with NSL Stereo RF-FE and attached Software with non GPS L1 signals revealed many deficiencies.

Thanks

I thank Prof. Mathias Hinderer for the hint that my thesis for MSc TropHEE would also be possible in other fields than geology. He advertised also topics of other earth sciences like geodesy where I finished modules in physical and satellite geodesy during my BSc studies in Applied Geo Science. He en­couraged me to look for the topic of software (Matlab)-defined receivers for GNSS in the geodesy department.

I thank Prof. Mathias Becker, the head of the Department of Physical and Satellite Geodesy of the Institute of Geodesy for offering this thesis to me and for the thorough support by him and his staff.

Dr. Stefan Leinen was my every day partner during the work on the thesis. He was always a very well informed and a very helpful adviser. He provided hardware, software, books, scientific papers, advice, moti­vation and weekly feedback.

Dr. David Becker gave some valuable support for using the Matlab “Parallel Computing Toolbox” and did also test runs of the thesis’ software on the institutes server with lots of parallel Matlab workers.

Thanks to the “Matlab Student program” of the company MathWorks I got Matlab, Simulink and some toolboxes (DSP, Parallel Processing, …) needed for the thesis for a low price.

My family deserves my thanks for coping with my reduced availability during the work on this thesis and temporarily conversion of a room into a GNSS lab with open window and GNSS antennas at the end of a beam protruding out of the window.

1. Introduction and targets

There is a strong trend everywhere to replace hardware by software to save material and to get more flexibility: Since many years cheaper printers for example need processing in the CPU and operating system of an attached PC to operate. TV receiver sticks and boxes contain only the essential hardware and most of the processing is done by software (in a computer). In more and more applications this trend is followed to reduce cost for hardware and hardware adaptions and to implement changes to the processing faster by changing high level software. One prominent example are DVB-T USB sticks with the RTL2832U chip (and preferably R820T tuner) and open source projects like RTL-SDR (Laufer 2015), GNU-radio, #SDR, ... providing software to turn the sticks into universal receivers of all frequencies between 30 and 1800 MHz for about only 20 Euros.

Also GNSS receivers in some special fields started to replace specialized hardware (ASICs, FPGAs) and processing in hardware (Verilog, VHDL) by higher level software programs (written in C/C++ and even in Matlab). Using software written in high level languages makes adaptions for changing pro­ce­dures (new frequencies, codes, modulations, …) of GNSS and analysis according to individual require­ments easier and faster. One drawback of this trend is the need for much higher processing power of freely programmable general purpose CPUs compared to specialized hardware (ASICs, FPGAs) and much longer time needed to perform the processing.

Borre (2007) describes how to do the processing of a Software Defined GNSS-Receiver – SDR with Matlab programs for GPS L1 signals. The programs of this “SoftGPS” project of the danish GPS Center at the University of Aalborg and the Colorado Center for Astrodynamics Research – CCAR at the University of Colorado are open source. After 2007 there were some software updates. To make them working with old data (2006/07) on the books DVD Dr. Leinen had to correct some errors. Meanwhile the hardware used for this programs is no longer available. The project is not continued, related web sites are closed.

The Institute of Physical and Satellite Geodesy of TU-DA acquired “Stereo” as Two-channel-Radio Frequency Frontend RF-FE (with related grabbing software from 2012/13) from Nottingham Scientific Limited (NSL) for this MSc thesis that started with the version “v4_vll” (GSM_M_SG/v4_vll) of the Matlab code with corrections by Dr. Leinen. Stereo is no GNSS-receiver but contains just pre-amplifiers, mixer, sampler, 2- and 3-bit analog to digital converters (ADC) and some FPGA-logic to control these chips and to merge the two digital signals into one single data stream for a USB 2 con­nection. For the configuration of the Stereo device and grabbings (= recordings) of GNSS signals soft­ware from NSL without all sources available is necessary. This software was buggy, not sufficiently maintained and initially not working in current Linux distributions at all.

The first goals of this thesis were to make the Stereo RF-FE and the Matlab software of (Borre 2007) working together, to do own grabbings of (different) GNSS signals and process in Matlab the files produced with Stereo. Further targets were a quality analysis of signals (signal to noise ratio) and of navigation results (means, medians, standard deviations, and root mean squares of positions and velo­cities and the deviations from a reference) and the improvement of the code in the fields of modulari­zation and processing speed especially by parallelization. The complete work flow from hardware setup over call sequence of functions and algorithms to the presentation and storage of navigation results was to be thoroughly documented.

Optional/additional tasks were to prepare the code for real time processing by continuous processing of data with different treatment of samples after evaluating at least one complete frame with the navi­ga­tion message to get the ephemerides and the processing of other GNSS signals as GPS L1. The very last task had to be canceled after NSL Stereo and the related software showed strong deficiencies in pro­viding such signals.

The work was done on a Notebook Acer Aspire Timeline X 4830 TG with CPU i5-2430 M (3273 units of http://www.cpubenchmark.net) with 16 GB memory, 1 TBßHD, Matlab 2016a-64bit in Windows 10-64 bit and Ubuntu 16.04 LTS 64 bit in a multiple boot configuration.

2 Theory

2.1 Software-Defined Radio - SDR

Still most of radio receivers and GNSS receivers are implemented completely in hardware. For small and cheaper devices most of the components are packed into a single hardware chip, an Application Specific Integrated Circuit – ASIC. In Contrast to a general purpose and freely programmable CPU of computer like a PC this ASIC was tailor-made for a single task and cannot be used for a completely different task just by calling a different software program. In a Software-Defined Radio – SDR – only the absolutely necessary parts – antenna, filters, preamplifier, (down-)mixer and Analog/Digital Con­ver­ter (ADC) - are implemented in hardware. The faster the general purpose freely programmable hard­ware becomes, the less receiver hardware is necessary. Generally hardware filters and (down-) mixers are not really necessary – the original High Frequency signal could be converted directly by an ADC. Sampling the original signal means to get more samples than necessary and high computational de­mands to handle them.

So today SDR-projects use more hardware than the absolutely necessary parts and transfer lower Inter­mediate Frequency (IF) or base band (IF=0) signals quantized with a few bits per sample to the software. Sampling and quantization of the received signals delivers a stream of numbers that is further treated by software programs. By changing the software, the signals can be treated in strongly varying ways. Examples are GNU-Radio and SDR-project for very cheap receiver hardware (DVB T receiver in USB-sticks for converting PCs into TVs but run in debug mode) receive e. g. TV, Radio, Air Traffic Communi­cation (ACARS, ADS-B), NOAA & Meteor-M weather reports from satellites, Marine Auto­matic Identifi­cation & weather fax (AIS, WEFAX), GSM signals… depending on the software on a PC for deciphering the stream of numbers.

The software of the SDR provides Digital Signal Processing (DSP). DSP offers for example numerical filtering (Finite and Infinite Impulse Response filters – FIR & IIR filters) and control loops like Phase Lock Loops (PLL) or Delay Lock Loops (DLL) to produce certain frequencies and phases.

2.2 Software GNSS Receiver - SGR

The related hardware delivers again a stream of encoded samples, i. e. a stream of numbers. In GNSS re­ceivers the software has to analyze the stream of numbers from the ADC and to recognize for a se­lected GNSS the transmitting satellites (their coding and/or frequencies), their Doppler shift in fre­quen­cies caused by the relative motion between satellite and receiver, the navigation message with orbit parameters, time, variation coefficients, time needed to travel to the receiver, pseudo ranges, position, velocity, time (PVT).

The discovery of satellites and the coarse determination of frequencies and phases during a short period is called acquisition. The more exact determination of frequencies and phases during a long period, that also enables the software defined receiver to demodulate the sample stream exactly by removing carrier and identifying code, is called tracking.

From the tracking results the navigation message is assembled. Its start is found by identifying sub frames. Orbit parameters (ephemerides) and times are decoded from the message and allow further navigation calculations like satellite and receiver positions – see figure 1:

Abbildung in dieser Leseprobe nicht enthalten

Figure 1: Software Defined GNSS Receiver (Tsui 2005)

2.2.1 Spectral analysis

Every GNSS signal has its specific frequency spectrum. The signal intensity (the power or squared signal amplitudes) around the intermediate frequency (IF) or more general for a frequency range between 0 and sampling frequency (that must be at least twice IF) is calculated by a procedure pres­ented by P. D. Welch in 1967 with overlapped portions of samples for averaged peridiograms. Non-over­lapped sample portions were already used by M. S. Bartletts in 1948/50 for periodiograms. See (Proakis & Manolakis 2007) and (Beucher 2015) for details. Matlab provides a pwelch() function. The function “probeData.m” from (Borre 2007) for example uses 100 ms of samples – that is 100 ∙ samplesPerCode = 2600000 samples for the default sampling frequency of 26 MHz of NSL Stereo. Function pwelch() forms 8 chunks of these samples overlapping by default by 50%. The maximum power is expected near IF – that is 6.5 MHz for GPS L1 in NSL Stereo. Figure 2 shows the power spec­trum of a GPS L1 signal:

Abbildung in dieser Leseprobe nicht enthalten

Figure 2: Power Spectral Density - PSD - from pwelch()

Acquisition and tracking of the GNSS signal can only be successful if the Power Spectrum has the ex­pec­ted appearance.

2.2.2 Acquisition of GPS C/A signals

The target of satellite acquisition is to identify visible satellites - in case of GPS L1 signals - by their pseudo random noise (PRN) code because they are all sending with the same carrier frequency only changed by the Doppler shift caused by relative movements between satellites and receiver. The signal looks like noise. Acquisition tries to find within the noise some of the up to 32 satellite identifying Codes by correlation of receiver generated codes with the received signal. Unknown Doppler shift and unknown code phase require a multiple of 32 correlation tests.

The code sent by the satellites modulated onto the L1 carrier is called Coarse/Acquisition- or PRN-code. Such Codes are members of the larger family of Gold-codes: Robert Gold introduced the production of pseudo codes by linear feedback shift registers. The carrier with a frequency of 1575 MHz is modulated with the satellite identifying PRN-code at 1.023 MHz repeating 1023 bit long sequences of the same code each 1 milliseconds long (=1023 / 1.023 MHz) and the even much lower frequent navigation message at 50 bits/second. PRN-code and navigation message are not directly detectable from the digitized IF that comes from the RF-FE (6.5 MHz default of NSL Stereo for GPS L1). The navigation message can be deciphered not before tracking was success­ful for at least 36 seconds.

The PRN codes of up to 32 GPS satellites must be generated by the program equivalent of some shift registers. As the codes are 1023 bits long at least all possible 1023 integer code phases (in units of entire chips = “code bits”) combined with many possible Doppler shifts – in (Borre 2007) ±7 kHz (or up to ±10 kHz for very fast moving receivers) in 0.5 kHz steps are used (= 2 ∙ 2 ∙ 7 +1 = 29 Doppler variants) in acquisition. Of course the code phase can be determined more precisely (e. g. in steps of one data sample = 1/26 code bit for default NSL Stereo sampling rate of 26 MHz). For details see (Borre 2007) and (Tsui 2005).

The carrier is removed from the signal by multiplying the signal with the generated carrier (as) which produces the two sums “signal + carrier” and “signal – carrier” (= PRN- code + navi­gation message). The latter is processed further. It has a real part called ‘I’ (for In-phase) and an ima­gi­nary part called ‘Q’ (for Quadrature).

One very desired property of the PRN codes is that the cross correlation between different codes is very low compared to the auto correlation (cross correlation with itself). So the high number of gene­rated PRN code combinations can be correlated with the received signal for searching high correlation values which identify the satellite (its PRN) and delivers coarse values of its code phase and Doppler shift. The Matlab code of (Borre 2007) does this for two ranges of data samples each one 1 millisecond long. Two are used because the change of a navigation bit (20 ms long) can fall into one of these two sample ranges and could disturb the correlation with generated codes. The correlations of the “better” millisecond are used to determine PRN and first coarse values of PRN code phase and Doppler shift.

The cross correlation or con­vo­lu­tion in the time domain is replaced by the equivalent sequence of (Fast) Fourier Trans­for­ma­tion (FFT), complex multiplication in the frequency domain, and Inverse (Fast) Fourier Transfor­ma­tion (IFFT) back to the time domain. The results are squared because only signal power instead of amplitude is relevant for the search. The highest and the second-highest peak are compared. If their ratio is larger than a certain threshold of e. g. 1.75 then the related satellite counts as acquired. Figure 3 shows the correlation above the frequency and code phase area for a single PRN (here 21) with the IF = 9.55 MHz of a RF-FE used in (Borre 2007), that is no longer available. Now the NSL Stereo RF-FE has IF = 6.5 MHz instead. Figure 4 shows the corre­la­tion / acquisition results plotted by function “plotAcquisition.m” as a bar chart for all 32 possible PRN codes:

Abbildung in dieser Leseprobe nicht enthalten

In a second step another FFT for 10 ms of samples (260 000 samples for the default sampling fre­quen­cy of 26 MHz of NSL Stereo for GPS L1) is used to improve the precision of the carrier fre­quen­cy (more exactly of the Doppler shifted IF). PRN code phase (also called Code Offset) can be reported as number of samples (0..26000 for NSL Stereo for GPS L1).

2.2.3 Tracking of GPS C/A signals

Tracking must demodulate the signal from acquired satellites over longer periods and must do this always exactly. Demodulation i. e. removing carrier and PRN-code from the signal to get the navi­ga­tion message bits, is based on mixing the signal with receiver generated versions of carrier and PRN-code that are in agreement with the received ones.

Acquisition only delivers first coarse values of carrier frequency (more precisely of IF) and PRN code phase. These are the start values for tracking. Frequencies and phases change con­stantly due to satellite (and receiver) motion and the software generated versions of carrier and PRN-code need continuous adjustment. Only if received and software generated carrier and PRN-codes are in agreement for frequency and phase the demodulation works correctly. Tracking uses there­fore two interwoven control circuits or Lock Loops (implemented in DSP software) with a nume­rical controlled oscillator (NCOs) for the carrier and a numerical controlled PRN code generator to get more exact and current frequencies and phases of both carrier (IF) and PRN code. These control circuits / Look Loops are only able to adjust small devia­tions. They must be fed with nearly correct values at the start to be able to work properly.

During tracking all samples are read and processed in chunks of the PRN code period of 1 ms (roughly block size = 26000 samples in NSL Stereo) and the tracking results are also saved in millisecond distances. To decipher later the navigation message, successful tracking of at least one GPS C/A frame (= 5 sub frames à 6 seconds and à 6 ∙50 = 300 navigation bits) is necessary. The start sample can belong to every sub frame, also to incomplete ones. Therefore at least 6 sub frames = 36 seconds must be tracked. To be able to skip some samples at the start of sampling files another second is added.

Processing at least in Matlab is not fast enough to handle single samples arriving in real time and a RF-FE like NSL-Stereo normally produces first a big file, then export it to different data formats and finally the software can read samples from that file. Normally more than 36 seconds (e. g. 37 seconds) are used as the first tracking period. (Borre 2007) uses by default a structure with 37000 elements of its components for frequen­cies, phases, real (I) and imaginary (Q) amplitudes of samples and deciphering of navigation bits does not start before these 37000 elements are successfully filled with data. As in acquisition also in tracking but not only for a few milliseconds but for all 37000 milliseconds the carrier (more exactly its IF) is removed from the signal by multiplication with the synthetic carrier (as ) generated by program code. The result is a base band signal (around IF=0 Hz) containing PRN code and navigation message. To remove also the PRN code from the signal this signal must be multiplied with a software generated PRN code (correctly in frequency and phase).

Both is achieved by two interwoven control circuits: a PLL (as bit switch insensitive Costas loop) for the carrier (its IF) and a DLL for the PRN code phase. The PLL has a NCO (as adjustable sinus ge­ne­ra­tor) and a filter. The DLL has a numerically controlled PRN code generator (rectangle signal gene­ra­tor). Related coefficients are calculated for every millisecond. For the DLL the generated code is shifted e. g. half a chip in both directions to get an early, prompt and late version of the 1023 bit PRN code. The code phase is refined during tracking by looking on early and late squared amplitudes (= powers) of real (=I) and imaginary (=Q) parts of the signal in the DLL – see (Tsui 2005) for details. The PRN code generator for the DLL is controlled by the variable “codeNco” influenced by the “codeError” cal­culated from early and late IE, IL and QE, QL powers. The correct code is generated when early and late powers are of same size. Because this thesis has an emphasis on programming (instead of mathe­ma­tics) mathematical symbols are often replaced by variable names also used in the Matlab code for better recognition:

Abbildung in dieser Leseprobe nicht enthalten

The PLL is controlled by the (arcus tangens of the) ratio of prompt values of Q and I-amplitudes – see (Stewart 2015) for details. The correct carrier frequency is achieved when the prompt Q-power is 0:

The following figure 5 shows the combined Lock Loops to generate exact carrier (IF) frequencies and PRN code phases to remove carrier and PRN code from the signal (= demodulation) to get the basis for the navigation message:

Abbildung in dieser Leseprobe nicht enthalten

Figure 5: Phase (Costas) and Delay Lock Loops (PLL & DLL) from (Borre 2007)

With properly working control circuits carrier (IF) and PRN code are removed from the sampled signal (the signal is demodulated) and the prompt values of I- amplitudes contain the navi­ga­tion message.

2.2.4 Navigation message and PVT calculations

The start of the navigation message must be found before the message bits can be recovered and de­ciphered. Preambles help to find the start point.

Navigation data is received with a boud rate of 50 (50 bps, every 20 ms a new bit). It is organized in frames (1500 bits/30sec) of 5 sub frames of 300 bits/6sec each. Sub frame 1-3 is constantly repeated: sub frame 1 contains satellite clock and health data, sub frame 2 and 3 contain satellite ephemeris data. Sub frames 4 and 5 are different in every frame and contain almanac data for other satellites. After 25 frames (25 ∙30 sec = 750 sec = 12.5 min) 25 ∙2 = 50 sub frames of number 4 or 5 are sent with data from foreign satellites and the sequence is repeated. A sub frame is organized in 10 words à 30 bits (24 data + 6 parity bits). The first word of every sub frame is the telemetry word (TLM). The second word of every sub frame is the Handover Word (HOW). The TLM contains 8 bits of the preamble (fixed bit pattern) for frame synchronization. The preamble is followed by 16 reserved bits and parity infor­ma­tion. All words end with 6 bit parity for data validation. Its values are calculated by multiple XORing of the 24 data bits of the word and the last 2 parity bits from the last word. For frame/ sub frame structures for GPS see (Borre 2007) or (Tsui 2005).

Preambles help to find the start of navigation message bits and (sub) frames. The preamble bit sequence [1 -1 -1 -1 1 -1 1 1] is extended (inflated) bit-wise to the 20-fold because during a single bit of navigation data (and therefore also during a single bit of the preamble) 20 complete PRN codes are send. This is correlated with the I_P values (e. g. with 2-bit encoding they have the values -3, -1, 1 or 3 multiplied by the number of samples per millisecond) changed to 1 for positive and to -1 for negative values. The “inflated” preamble sequence is moved along the 37000 binary values derived from IP and the correlation between both (maximal 8 ∙20 ∙12=160) is calculated. High correlation values of more than 153 (= 95% of the correlation maximum) are received at positions were the sample bits are equal to the preamble bits. Often more than the expected 6 positions (from 6 sub frames) are found. Such positions of sub frame start points must be verified. 6 seconds = 6000 ms (6000 elements in a structure trackResults) later the preamble must be repeated i. e. produce the next high correlation value.

If this preamble repetition test is passed for a position, then parity checks are done for the first two words (TLM + HOW) of a sub frame. The parity values are calculated and compared with transmitted values. Because also 2 parity bits from the predecessor word are needed for parity calcu­lation (see (Borre 2007), 2+ 2 ∙30 = 62 bits must be read. Every bit is 20 ms = 20 PRN cycles long. Related to the start of a preamble we must read I_P (In- phase prompt) values from position -2 ∙ 20 = -40 to position +60 ∙20 -1 = +1199. Theo­reti­cally 20 consecutive I_P values should have the same value because they belong to the same navi­gational bit, but that is seldom the case. Chunks of 20 values (equivalent to 20 ms) of I_P are therefore summed. If these sums are positive the represented bit is a ‘+1’. Zero or negative sums represent a ‘-1’. 62 of these bits are the input of a function to calculate the parity bits by XORing the bits.

If all these correlation-, repetition-, and parity tests for preambles are passed, then the navigation bit sequence is generated. A whole frame of 1500 navigation bits (5 sub frames ∙ 300 navigation bit/sub frame, equates to 1500 ∙ 20 complete PRN code periods à 1ms = 30000 ms) plus 2 navigation bits in front of the preamble (= last two parity bits from last word of predecessor frame = 2 ∙ 20 ms = 40 ms) of I_P (In-phase prompt) data is read into “navBitsSamples”. Chunks of 20 are summed up because 20 PRN code periods are contained in one navigation bit. If the sum is positive the 20 are represented by a single 1 otherwise by a 0.

The ephemeris-function is called to decipher data like time and orbit characteristics from the navi­ga­tion bits (1 complete frame + 2 last parity bits from pre­decessor). One gets the Time Of Week – TOW – also called z-count, the time since GPS time rollover (mid­night Saturday→Sunday) in units of 1.5 se­conds with number truncated by 2 bits (22 units = 4 units = 4 ∙1.5 sec = 6 sec). Z-count is Time Of Week with a reduced precision of 6 seconds from the first part of the second word (HOW) of a sub frame. Parity checks are repeated for all sub frames. The sub frame ID = 3 bits = position 50:52 in sub frame = position 20:22 in second word (HOW) of a sub frame. See (Borre 2007) and (Tsui 2005).

Figure 6: First two words of every sub frame of GPS L1 C/A from (Borre 2007)

For every sub frame different data is extracted and written to an ephemerides structure (“eph”). The sub frames 1-3 are described in detail in (Tsui 2005). Time Of (GPS-)Week is read from a sub frame from bit position 31-47 (= position 1-17 in the second word = HOW). That is the by 2 bit truncated code for the number of units of 1.5 seconds each with a precision of only 4 units = 4 ∙1.5 sec = 6 seconds. This number is multiplied by 6 (seconds) and 30 (seconds) are subtracted because it is the time for the end of a 30-second-frame but the time at the start of the frame is needed.

The ephemerides describe for the reference time the form of the elliptic satellite orbit (semi major axis A and eccentricity e), its orientation in space (inclination ‘i’ between orbit and equatorial plane and longi­tude Ω of ascending node), satellites position on the orbit (argument of perigee ω and mean ano­maly M), and parameters for change rates (doted parameters) and corrections. The elements of satellite orbi­tal motion are described in (Seeber 2003) – see figure 7:

Figure 7: Kepler's orbit parameters from (Seeber 2003)

Every millisecond of tracking data has also a related last sample and gets an absolute time (since start of the GPS week) assigned. The transmission time is interpolated between these records (one per ms). Two interpolations are calculated in function “findTransTime.m” - see chapter 3.1.3. The receiving time is at least 68.802 ms later because this time is needed for the wave to travel the shortest distance of about 20000 km from a GPS satellite to the Earth’ surface.

PVT stands for position, velocity and time (of the receiver). The transmission time of sub frame starts in GPS seconds since GPS week start and the times of tracking records (spaced by milliseconds) were already calculated.

The satellite positions are calculated for the transmission times based on the ephemerides. Details are found in (Borre 2007), (Tsui 2005), and (Seeber 2003). Receiver positions X, Y, and Z in earth centered earth fixed ECEF-coordinates are calculated. The positions are results of the least square method fed with the pseudo ranges called distances between satellites and receiver as product of speed of light ‘c’ and difference of transmission and receiving time of the signals. The velocities are derived from the Doppler shifts of the signals. The three positions and c ∙ dt are components of a vector x in the equation of motion with m≥4 satellites. Matlab solves this easily for the x-vector with ECEF positions and time error using the method of least squares.

Similarly works the calculation of the 3 velocities that is based on Doppler shifts.

The “Inaccuracy” or Dilution Of Precision (DOP) of this calculation depends on the A matrix of the least squares method. Details of the Geometrical-, Horizontal-, Vertical-, Positional- and Time- compo­nent of DOP are found in (Bauer 2002) and (Kahmen 2006).

Multiplying ECEF coordinates (for positions and velocities) with rotation matrices changes them to UTM coordinates.

3 Methods and results

Normally methods and results are covered in separate chapters, but here the results are quite clear because they are explicitly requested and the methods are strongly related to the requested results. The methods have to provide a way to get these results. Therefore both are put close together and the main divi­sion of the material is driven by topics of documentation, new functionality and tuning for speed and modularity.

3.1 Analysis and documentation of the code

The first version of the code on the book DVD of (Borre 2007) is roughly documented in the books appen­dix A and in a few comment lines each starting with ‘%’ at least at the top of the Matlab files. The code needed for the NSL Stereo RF-FE was from 2012/13 and was to be be downloaded from the NSL website. The NSL Stereo booklet describes mainly the hardware and C-programs to configure it and to record grabbing files. Of the accompanying Matlab code only the start-program “visibility_mgui.m” is docu­mented. More and better documentation was requested. A better docu­men­ta­tion has not only a high value for personal starting to work with NSL Stereo RF-FE and related pro­grams but was also necessary to find the positions in the code were in later steps corrections and extensions where to made.

A good means to produce complete (in the sense that all files are treated and nothing is ignored) and con­sistent documentation is a documentation generator like “M2HTML” (Matlab to HTML) that generates text and graphical documentation to read with web browsers. This static view does not inform about the time sequence of function calls – the dynamic aspect of the Matlab code to document. The call sequence was investigated manually by reading the code and following the function calls. Also break points for the debugger were set to stop processing of Matlab code, check the contents of variables and to step line by line through the code processing only small parts of the code between checks of the program state.

M2HTML is free Matlab code that investigates other Matlab code and generates two types of docu­men­tation in HTML, i. e. viewable with a web browser: text documentation (list of files and some few comment lines from the top of each file) and call or dependency graphs visualizing which file calls which others. The later is mainly done by GraphViz a free tool to draw different styles of graphs for a text file containing just the source and target nodes of a directed graph (a digraph) in so called dot-language. GraphViz positions nodes and draws edges between them. GraphViz programs are normally called automatically by M2HTML generating inver­ted tree graphs (root points up). But the dot-com­mand of GraphViz not only offers different styles of graphs (inverted tree becoming wider and wider further down, graphs developing from a center into all directions, …) but also exporting them to different file formats like GIF, JPEG, portable network graphic (*.png), scalable vector graphic (*.svg) and much more. Looking for commands within the GraphViz package and calling them directly to draw varying graphs according to node relationship files (*.dot) generated by M2HTML can produce much nicer and easier to handle diagrams than the ones generated by defaults. In Linux the GraphViz path was added automatically, in Windows the current version 2.38 of GraphViz no longer sets the path (install into “C:\Program Files (x86)\Graphviz238\”) and the website gene­rated by M2HTML as docu­men­tation of some Matlab code does no longer show the graph.

Free Matlab code of M2HTML (m2html.zip) and some documentation can be downloaded from the site http://www.artefact.tk/software/matlab/m2html/tutorial.php. It investigates Matlab files in indi­ca­ted folders (and sub folders) and reads some lines at the top of Matlab files. Matlab comment lines start with a ‘%’, such for sections with ‘%%’ and complete blocks of comments are between ‘%{‘ and ‘%}’ (both alone on their lines). GraphViz as graphical extension to visualize edges between nodes comes from http://www.graphviz.org/, where also the dot-language is explained. Both support multiple operating systems including Linux and Windows. In Windows the search path for GraphViz must be set manually e. g. to: C:\Program Files (x86)\Graphvix238\bin.

In Windows M2HTML doesn’t show the graph in image file “graph.png”. It generates a file “graph.html” for every investigated folder that contains amongst other HTML code an empty map-section i. e. no area shape definitions:

Abbildung in dieser Leseprobe nicht enthalten

This means M2HTML expects that GraphViz had generated from node-edge-definition in the file “graph.dot” an image file “graph.png” and a map file “graph.map” with area shape definitions also to be copied to the map section in “graph.html”. If the png-file does not exist, no graph is displayed. M2HTML generates the file “graph.dot” (the digraph definition) that is transferred to GraphViz’ “dot.exe” to produce the related png-file for an inverted tree graph and a related image map so that a click into the graph calls the hyper referenced HTML documentation file. The command “dot” is called with the following options and parameters:

dot -Tcmap -Tpng -o graph.map -o graph.png graph.dot

If this does not work properly, or one wants the graph in other file formats or if one wants other graphs than the default inverted tree, the *.dot-file (see figure 9 and 11 with all edges of a directed graph – digraph) can be loaded into GraphViz’s “GVEdit” (see figure 10) with path “C:\Program Files (x86)\Graphviz328\bin\gvedit.exe”. There one can select under ‘Graph→Settings’ the layout engine (‘dot’ for inverted tree, ‘neato’ for …, ‘fdp’ for concentric graphs), the output file type (png, gif, jpg, svc, ….), and the output file name (see figure 8):

Abbildung in dieser Leseprobe nicht enthalten

Figure 8: GraphViz' GVEdit graphic export

Abbildung in dieser Leseprobe nicht enthalten

Figure 9: First part of a dot-file

Abbildung in dieser Leseprobe nicht enthalten

Figure 10: GraphViz’ GVEdit with directed graph

Abbildung in dieser Leseprobe nicht enthalten

Figure 11: Last part of a dot-file

As soon as the file “graph.png” is in the folder with the other documentation generated by M2HTML the graph is displayed when the related link within the generated documentation is selected. Also the HTML-map-file that defines selectable areas in the graph to branch to the text documentation of a selected object of the graph is not generated automatically in the Windows version of GraphViz. One could try the dialog of figure 8 with the values: dot, ismap, “graph.map”, …

GraphViz’ version for Linux works much better. But the Windows version of M2HTML can be “re­paired” by editing “m2html.m” line 874 specifying the absolute path to “dot” or to rela­tives like “fdp”. This can also change permanently the graphs style (tree or circular or …):

Abbildung in dieser Leseprobe nicht enthalten

Instead of calling “m2html.m” in Matlab, the graphical user interfaces to it can be called with “wizard.m” or “wizard2.m”.

As the Borre Matlab Code was changed and extended many times during the work on this thesis the use of a code version control tool was recommendable. Also signal recordings (grabbings) and their analysis results from Borre or NSL Stereo code benefit from code version control systems. The old Borre code used a code version system (CVS), that placed special comments like the following with file name, version numbers, date and time in the header of the Matlab files (here in “init.m”):

- CVS record:
- $Id: init.m,v 1.14.2.21 2006/08/22 13:46:00 dpl Exp $

For this thesis the system was changed to Git, because Git is more modern, more widespread, inte­gra­ted in Matlab, … Git is a distributed code version control tool. It keeps complete files in different ver­sions instead only the latest version of files and collections of changes (deltas) to preceding versions. It is distributed, i. e. all users have a full set of files (also serving as a backup) without constant con­nection to a central server. The server database with meta information (e. g. SHA-hashes of all files, pro­gram­mer name, date, time, description) = Git directory = Git repository is used seldom and writing access can be postponed without hindering the work on the local copy = working directory. Between these two locations and states there is a third one called “staging area” or “index” where bookkeeping is done about files that are changed in the working directory and are to be committed next time to the git directory – see figure 12.

Abbildung in dieser Leseprobe nicht enthalten

Figure 12: State transitions in Code Version Control System Git from

All former versions, responsible person, the date and time they were committed, descriptions e. g. as commit message, … can be restored every time from Git. One consequence of hashing (check-sum­ming with a one-way-function like SHA1) everything is the loss of special comment lines with ver­sion information in the code files. Every file and folder is identified by a hash. If one inserts com­ment lines (with version information) the hash is changed completely. There­fore there should be no such spe­cial version comments for Git. Only one kind of comment with the hash value itself as com­ment is accep­ted and could be activated with some effort but this was not used. Git was used without using a server. Instead a sub folder “.git” was the repository and placed in the working directory.

Git has command line (text) and graphical user interfaces. Using Git-bash in the work for this thesis pro­vided a common UNIX-like environment with a Borne-Again-Shell in Linux and Windows. The main storage for Git was an “invisible” sub folder “.git” in the working directory (GSR_M_SG/v4_ll for Borre code, GSR_M_SG/20130108_stereo_v25-win64-beta/bin and 20160104_stereo_v25-linux64/bin for NSL Stereo code):

Abbildung in dieser Leseprobe nicht enthalten

Git is even used for the complicated development of the Linux kernel (thousands of files, hundreds of programmers, fast changing and appearing new hardware to be supported and dozens of years to work). Nevertheless there was a severe problem installing Git in Ubuntu 16.04 also discussed in bulletin boards for many month. If someone suffers such problems (or problems with Matlab crashes, installing and repairing the complete software environment of the thesis in Linux or Windows), the diaries (two volumes with more than 300 pages together) written during the work on this thesis should be consulted to get the solution.

Also for C-Programs there are documentation generators but there are only two programs for NSL Stereo: one (stereo_app) uses proprietary libraries not available in source code so that a proper docu­mentation was not possible and the other (stereo_exp) was analyzed line by line only with an editor and no other tools.

3.1.1 M2HTML and GraphViz for Borre’s Matlab code

The first version of Matlab code comes with the book DVD of (Borre 2007). The version GSR_M_SG\v4_vll contains some improvements and corrections published on web sites after 2007. The global settings are in file “initSettings.m”. There is a graphical program “setSettings.m” to change the contents of “initSettings.m” that should no longer be used: the code file “initSettings.m” has com­ments and is there­fore better understandable, it is easily changed by every editor (including that of Matlab) and it got new variables for the extended code written during this thesis that are not accounted for in “setSettings.m”. So if some settings must be changed, then “initSettings.m” must be edited di­rect­ly with an editor!

The program itself is started via “init.m”, that sets path to the subfolders “include” and “geoFunctions” of GSR_MSG\v4_vll and reads the variables in “initSettings.m”.

For the documentation generator the m2html-folder is put at the side of “v4_ll” i. e. to the same folder level. Also a target directory for the generated documentation (here it will be “doc”) can be put here. Before starting documentation generation the parent folder of “v4_ll” (and “m2html” and “doc”) is made current directory in Matlab. With Matlab command “addpath” the m2html-folder is added and then “m2html” with the following parameters for source (mfiles = v4_vll) and target (htmldir = doc) directories … is started:

addpath m2html

m2html ('mfiles','v4_vll', 'htmldir','doc', 'recursive','on', 'graph','on', 'global','on')

The parameters of the function call to “m2html” are explained in the following table 1:

Table 1: Parameters of the m2html function

Abbildung in dieser Leseprobe nicht enthalten

The directory with the generated documentation (here “doc”) contains the file “index.html” that starts viewing the documentation. See figure 13 with “index.html”:

Abbildung in dieser Leseprobe nicht enthalten

Figure 13: “index.html” in M2HTML generated doc-folder for v4_vll

Every entry is a link that can be followed into sub directories or to documentation of a single file with some lines from comments at the files beginning and (cross reference) lists of functions called from this one and list of files from where this one is called. Following the link of “v4_vll” shows figure 14 with list of files contained in this folder:

Abbildung in dieser Leseprobe nicht enthalten

Figure 14: Contents of single folder of M2HTML generated documentation

Here the files “init.m” and “postProcessing.m” have no adequate description because the comments at the head of the files were not placed correctly. A dependency graph can be generated for the whole pro­ject and for every (sub) folder of it. The processing starts with “init” that in turn calls “initSettings”, “probeData”, and “postProcessing”. Then “postProcessing” calls most of the other files/functions. See figure 15 and figure 16 with two halves of the dependency graph:

Abbildung in dieser Leseprobe nicht enthalten

Figure 15: M2HTML dependency graph for v4_vll, left part

Abbildung in dieser Leseprobe nicht enthalten

Figure 16: M2HTML dependency graph for v4_vll, right part

Selecting a single file name delivers description (figure 17) and cross reference (figure 18), here for “init” with the still incorrect comment line (% –--------------------------):

Abbildung in dieser Leseprobe nicht enthalten

Figure 17: M2HTML description of 'init'

Abbildung in dieser Leseprobe nicht enthalten

Figure 18: M2HTML cross references for 'init'

The complete M2HTML generated documentation of the final state (code changed during the work on the thesis) is found on the DVD accompanying this thesis.

The generated documentation eases understanding the processing but the time sequence of calling functions is not completely clear only visiting cross reference lists and dependency graphs. Therefore the code was read, every function call followed and the sequence documented on the following pages.

3.1.2 Adaptions of “initSettings.m” before the start of programs

Before the processing of an “exported” grabbing file can be started with “init.m” the file “initSettings.m” is edited at least to name the file to process. The most important variables (= compo­nents or fields of struct “settings”) are explained in the following table 2. Column ‘New?’ indicates variables created during the work on this thesis with “yes”:

Table 2: Variables (components/fields of struct ‘settings’) in file “initSettings.m”

3.1.3 Call sequence in the Borre code

When the variables in “initSettings.m” are all correctly set, then the program “init.m” can be started. The sequence of functions (and scripts) called is the following for VLLen=0 and for only one cycle of pro­cessing that delivers all ephemerides once. The new code added during this thesis can also do several cycles of acquisition, tracking, and navigation calculations and also shorter cycles of short tracking and reduced navigation calculations after a full first cycle with new functions “tracking2.m” and “postNavigation2.m”, explained in chapter 3.2.

init.m-program

It adds paths to sub folders “include” and “geoFunctions” of “v4_vll”, and calls the functions “initSettings.m”, “probeData.m”, and “postProcessing.m”. The last one is only called after the user enters a ‘1’. If the user enters a ‘0’ the program stops here.

initSettings.m-function

It is called by “init.m” and sets variables “settings.XXX=…” like “fileName” of the file to read, ...

probeData.m-function

It is called by “init.m”, gets the settings defined in “initSettings.m” as parameter, extracts the path and name of the file and reads from that file a short piece of data. Start position (skipNumberOfBytes) and length depend on values in “initSettings.m”. A time domain plot for the data read from file is plotted. A frequency domain graph (Power Spectral Density – PSD) is calculated with the “pwelch()” function. Finally a histogram of sample values is plotted.

The one-sided (= real input) function pwelch(data, 32758, 2048, 16368, samplingFreq/1e6) uses chunks of 32758 samples (= segment size = Hamming window size; default is so that the complete data is split into 8 chunks with 50% overlap between each other), the next chunk overlaps the former one by 2048 samples (default is 50% of window), 16368 DFT points are used to calculate each estimate. It returns the values of PSD for the frequency range from 0 to “samplingFreq”/2=26e6 /2=13 MHz. Maximum power is expected at IF=6.5 MHz. The two-sided (= complex I/Q-input) pwelch() returns the Power Spectral Density - PSD for a frequency range from 0 to “samplingFreq”. After plotting the frequency domain plot, a histogram (number of samples per value or range of values) for the data samples is calculated with A/D-converted signal amplitudes and plotted.

postProcessing.m-program

It is called by “init.m” if the user enters a ‘1’ to initiate GNSS processing (0 = end of program). It calls several functions – in particular “acquisition.m”, “tracking.m”, and “postNavigation.m” - to process 37000 ms of the data file. It identifies the satellites (aquisition of space vehicles =sv or rather pseudo ran­dom numbers = PRN) with the data for a few milliseconds. PRN, fre­quency (norm value of GPS L1 + Doppler shift) and phase is determined. Then a frame (= 5 sub frames of 6 seconds each with navigation data) is tracked. Due to a start not at the beginning of a sub frame one more (6 sub frames ∙ 6 sec/sub frame + 1sec = 37000 ms) is tracked. The program “postProcessing.m” opens the data file (rb) for read of big-endian (most significant bit first) data. The endianess can change with the receiver (mixer to IF + A/D-converter) to little endian (most significant bit last) which is standard e. g. in Intel/ PC equipment. Function fread(fileID, sizeA, precision) is called with some milliseconds of samples and precision = “settings.dataType = schar” (signed char from -27 to +27). In case of I/Q samples I0,Q0, I1,Q1, I2,Q2, … real_n,imag_n,... complex numbers are build as I0+i Q0, I1+i Q1, I2+i Q2, … With this data function “acquisition.m” is called.

acquisition.m-function

Is called by “postProcessing.m” with some ms of data (= “longSignal”) and the settings from “initSettings.m” as parameters. It calculates samplesPerCode=26000 (= 1ms for GPS C/A, see “probeData.m”) and creates 2 signal vectors “signal1” and “signal2” of this length filled with the 1st respectively 2nd millisecond of data from “longSignal”. “Signal0DC” is a data vector without DC i. e. with average 0. “phasePoints” is a vector [0:26000] ∙ 2 pi 1/samplingFreq. NumberOfFrqBins = round (acqSearchBand ∙ 2) +1 = 2 ∙ 14 (kHz) + 1 = 29 frequency bins for searching in 500 Hz steps within the range of ±14 kHz around standard fre­quen­cy of L1 carrier down mixed to IF to find the Doppler-shifted frequency. Then the function “makeCaTable.m” from sub folder “include” is called.

makeCaTable.m-function

Is found in sub folder “include” and called ones by “acquisition.m” with settings from “initSettings.m” as parameter to fill “caCodesTable” a matrix of dimension 32 ∙ samplesPerCode = 32 ∙ 26000. It calls 32 times (for PRN = 1 : 32, i. e. for all satellites) the function “generateCAcode.m” from the sub folder “include” with a PRN as parameter. In each loop one (PRN-)line of the caCodesTable-matrix is filled. “codeValueIndex” is a vector with samplesPerCode = 26000 elements (see probeData.m) but with only 1023 different values (1, 2, 3, …. 1023). Several adjacent elements have the same value. The 1023 different values of caCode are taken many times and copied to caCodesTable in the manner prescribed by variable “codeValueIndex”.

generateCAcode.m-function

Is located in sub folder “include” and called by “makeCAtable.m” for a certain PRN to generate the GPS C/A code for that PRN. It uses two shift register each 10 bit long. The 1st one G1 has polynomial 1+x3 +x10, the 2nd G2 has polynomial 1+x2+x3+x⁶+x⁸+x⁹+x10 . G2 is shifted depending on PRN, then the 1023 bit long Code of G1 and shifted G2 are multiplied to produce the PRN dependent 1023 bits of C/A code for one satellite. Cross correlation of different C/A codes is very low compared to an auto corre­lation (if code phase is also the same). A vector of length of 1023 (= 210 bits) is returned to variable “caCode” in function “makeCaTable.m”.

aquisition.m-function, continued 1

Having received the 32 ∙26000-matrix caCodesTable each row filled according to the codeValueIndex vector from “makeCaTable.m” with C/A-(PRN-) codes generated for a satellite, some variables are ini­tialized (with zeroes) to receive results later:

- results is a numberOfFrqBins ∙ samplesPerCode = 29 ∙ 26000 matrix,
- frqBins is a vector of numberOfFrqBins = 29 elements,
- acqResults.carrFreq, acqResults.codePhase and acqResults.peakMetric are vectors with 32 elements each for up to 32 satellites/PRNs.

In a large loop ’for PRN = acqSatelliteList‘ the complex conjugate of the FFT of a row of the matrix “caCodesTable” is saved in vector “caCodeFreqDom”.

In a sub-loop ”for frqBinIndex = 1: numberOfFrqBins” (stepping through possible frequency range around intermediate frequency IF) the data from the 1st ms (signal1-vector) and from the 2nd ms (signal2-vector) is multiplied element-wise by a vector sigCarr with values. I1 and I2 are the real parts whereas Q1 and Q2 are the imaginary parts of these products. Complex number vectors I1+j∙Q1 and I2+j∙Q2 are FFT‘ed to vectors IQfreqDom1/2. The inverse FFTs of element-wise product of IQfreqDom1/2 (fourier-transformed samples from 1st and 2nd millisecond) and caCodeFreqDom vec­tor (fourier-transformed PRN-codes from generator) is equivalent to the direct correlation/ convo­lu­tion of sampled and generated Codes in the time domain. The correlation results are squared and the milli­second (=26000 values) with the larger maximum (1st or 2nd ms with 26000 values) is saved in matrix results(frqBinIndex,:). Investigating 2 ms is necessary because 1 bit of navigational data (boud rate 50/s = every 20 ms a bit might change) can begin in one of the 2 ms and disturb the correla­tion/ con­volution. A millisecond without any change in navigation data is needed. Of the selected milli­se­cond there are 29 versions for different frequencies and the correct code phase (0 - 1022 possible inte­ger values for a code length of 1023 or finer fractional values down to one single sample out of 26000) is unknown.

After finishing the sub-loop “for frqBinIndex = 1 : numberOfFrqBins” value and position/index of maxima in the results-matrix are calculated: frequencyBinIndex= row number where the maximum of the values in that row is higher than the maxima of all other rows. CodePhase = column number (out of 1:26000) where the maximum of the values in that column is higher than the maxima of all other co­lumns. “CodePhase” is later saved to “acqResults.codePhase(PRN)” if the satellite shows a large first to second peak ratio. Beside the maximum peak also a second to maximum peak at least one chip away is searched within the same row (frequency is fix) of the results-matrix. The quotient “peakSize”/ “secondPeakSize” is saved in “acqResults.peakMetric(PRN)”. If it is larger than acqThreshold = 1.75 the sa­tellite/PRN is further processed and its PRN printed. Not acquired satellites are represented with a dot:

Acquiring satellites...

(. . 03 . 05 06 . . 09 . 11 . . . 15 . . 18 . . 21 22 . . . 26 . . 29 . . . )

Before calling function “generateCAcode.m” again, a new “codeValueIndex”-vector is calculated. There are 10230 different values for 260000 elements. Many elements have the same value. The remainder (= 0, 1, 2, … 1022) of the division by 1023 is calculated to determines how the 1023 different values in caCode=generateCAcode(PRN) are to be distributed to the 260000 elements of vector “longCaCode”. This vector is element-wise multiplied with DC-free signal samples from “signal0DC”-variable „blown up“ from 26000 to 260000 (representing 10 ms) elements by index codePhase:[codePhase + 10 ∙samplesPerCode -1]. Each of the 26000 data samples is used 10-times. Reason for so many data points: more points in the time domain result in finer frequency resolution in the frequency domain. The product is saved in vector “xCarrier”. The FFT of “xCarrier” is calculated for data points (i. e. more than elements in “xCarrier” so that zero padding appears) and stored in “fftxc”. More points in the time domain are again rela­ted to a finer resolution in the frequency domain. The maximum in “fftxc” and where it is located is stored in “fftMax” and “fftMaxIndex”.

postProcessing.m-program, continued 1

Back from “acquision.m”, the PRNs of satellites with strong signals, and their frequencies (of carrier mixed down to IF), and codePhase are stored in “acqResults.*”. Function “plotAcquisition (acqResults)” is called.

plotAcquisition.m-function

It is called by “postProcessing.m” with acqResults-structure (carrFreq, codePhase and peakMetric) as para­meter and plots a bar-chart of “acqResults.peakMetric” over PRN (1-32). Stronger Signals (acquired ones) are marked in green. The related bar-chart was already shown in the theory chapter.

postProcessing.m-program, continued 2

Back from “plotAcquisition.m” the functions “preRun.m” and “showChannelStatus.m” (both from sub fol­der “include”) are called if at least one satellite/PRN has a “peakMetric” greater than “acqThreshold” (first maximum divided by second maximum peak > 1.75).

preRun.m-function

It is located in sub folder “include” and called by the program “postProcessing.m” with acqResults- and settings- structure as parameter. It creates a column vector of numberOfChannels=12 (for sa­tellites/ PRN with strongest signals – one per line) elements each a structure with the components “PRN”, “acquiredFreq”, “codePhase”, and “status”. The elements are sorted for descending “peakMetric” (strongest on top, weakest at bottom, e. g. PRN = 21, 22, 15, 18 26, 6, ...). It returns to “postProcessing.m” where imme­diately the function “showChannelStutus.m” is called.

showChannelStatus.m-function

It is located in sub folder “include” and called by program “postProcessing.m” with channel- and settings- structure as parameter and prints the contents of “channel” as table on the screen. This is data from the old Borre code with a grabbing file from the Book DVD (not from NSL Stereo). Therefore the (inter­mediate) frequency is around IF=9.548 MHz (Stereo has 6.5 MHz) and the Code Offset in samples can be up to 38192 for the sampling frequency of 38.192 MHz (Stereo has sampling frequency 26 MHz and therefore a Code Offset of maximal 26000):

Abbildung in dieser Leseprobe nicht enthalten

postProcessing.m-program, continued 3 (start tracking)

Back from “showChannelStatus.m” tracking starts. Acquisition only delivered PRN and rough values of frequency and phase. These values must be determined more exactly and recalculated all the time be­cause they are changing with moving satellites (and receivers). If VLLen=0, then classical SLL (serial lock loop) instead of VLL (vector lock loop) is done. A file “trackingResults.mat” (the saved trackResults-struct) that is needed for later VLL-tracking is produced. VLL is not treated in this thesis. Because the SLL-tracking needs a lot of time, “startTime” is taken for later time measure­ment and printed on the screen. Function “tracking.m” is called with channel- and settings-structure as parameter. “tracking.m” returns the trackResults-structure:

tracking.m-function

It is called from “postProcessing.m” with the channel-struct of selected satellites and settings. It creates the structure “trackResults” with the components/fields:

- status (‘T’ for probably trackable or ‘–‘ for not acquired and not trackable satellites),
- absoluteSample (row vector with msToProcess = 37000 elements with sample numbers initialized with 0; later filled with growing numbers from 0 to 37 ∙ 26000000),
- codeFreq (row vector with msToProcess=37000 elements initialized with Infinity) and
- remCodePhase (row vector, 37000 elements each with a remainder between 0 and 1) for C/A-code,
- carrFreq and remCarrPhase for carrier (down mixed to IF),
- I_P, I_E, I_L, Q_E, Q_P and Q_L for prompt, early and late correlators for real In-phase and ima­ginary Quadrature values for,

- delay (DLL) and phase (PLL) lock loop discriminators (all vectors with msToProcess = 37000 ele­ments),…

The structure “trackResults” is duplicated a “numberOfChannels” (=12) times to keep the results for all tracked satellites.

Carrier (frequency and phase) and Code-phase tracking must update results from acquisition because the values do not stay constant and exact. Replicas of the carrier and the C/A code are necessary to re­move both from the signal by multiplication (= demodulation) in order to get the navigation data. Mul­ti­plication of 2 signals creates new signals: one with the sum and another with difference of the fre­quencies. The signal with the sum of frequencies (double frequency) is filtered away by a low pass filter. The signal with the difference of the 2 frequencies is the de­modu­lated one.

The function “calcLoopCoef.m” from the sub folder “include” is called twice with parameter xxxNoiseBandwidth = LBW = Bl , xxxDampingRatio = zeta = ζ, and xxxLoopGain = k. Where xxx is dll for code tracking with a delay lock loop with values 2, 0.7, 1.0 and xxx is pll for carrier tracking with phase lock loop in a Costas loop with values 25, 0.7, 0.25.

calcLoopCoeff.m-function

It is located in sub folder “include” and called by “tracking.m” to calculate PLL- and DLL- loop coeffi­cients, natural frequency, tau1 and tau2. The tau’s are returned to the caller. See (Borre 2007) and (Tsui 2005) for details. Both tau’s are returned to “tracking.m” to calculate (later) the commands for the numerically controlled oscillator (NCO) in the carrier lock loop and the numerical controlled PRN code generator in the delay lock loop to minimize their errors. It let the generated signals lock into measured signals, i. e. get the same frequency and phase.

tracking.m-function, continued 1

Back from the two calls to “calcLoopCoeff.m” and having received two filter coefficients tau1/2 for the filter of code DLL and two filter coefficients tau1/2 for the filter of carrier PLL (many lines “later” they are used to control NCO and PRN code generator) the function “waitbar(fractional_length_0_1, Text_displayed, para­meter_name_value_string_pairs)” is called to create a progress bar showing later how many milli­se­conds (of the 37000) are already processed for the current satellite/PRN/Channel-No.

In a large loop “for channelNr = 1:numberOfchannels” (1:12) the channels are processed if their PRN is not zero (if the satellite showed a decent signal). Function “generateCAcode.m” returns the C/A-code for the current PRN to “caCode” (vector with 1023 elements – each +1 or -1 instead of +1 and 0). Then caCode(1023) is pre-padded and caCode(1) is post padded to get a new caCode-vector with 1025 ele­ments. “oldCodeNco” and “oldCodeError” are initialized to 0 for control of the PRN code generator in the code tracking loop (DLL with Early, Prompt and Late generated Codes). “oldCarrNco” and “oldCarrError” are initialized to 0 for control of the NCO (Numerical Controlled Oscillator) in the carrier tracking loop (PLL, Costas Loop with I and Q check). For PLLs see also (Stewart 2015).

An inner for-loop (for loopCnt = 1: codePeriods) all 37000 milliseconds (1 ms = 1 complete C/A code) are read from the file ms for ms and processed. loopCnt/codePeriods = 0.x= fractional length of the waitbar (a progress bar), updated each time another chunk of 400 ms (original value was 50 ms that wasted much more processing time) is processed if setting showWaitbar=1.

Abbildung in dieser Leseprobe nicht enthalten

Figure 19: Waitbar displayed during tracking

After implementing a new function “CNoVSM.m” the waitbar also displays the current signal to noise ratio - SNR (or carrier to noice ratio C/No) in db. But with parallel processing – enabled by new code written during this thesis – the waitbar must be switched off by setting showWaitbar=0.

One complete C/A-Code has 1023 bits (with 1sec/1.023 e6 per bit) and needs 1 ms. One C/A-Code bit is represented by circa 26 data samples. In every call to fread() blksize bytes (schar) of real data (S0, S1, S2, …, S26000) or up to 2 ∙blksize bytes (schar) of complex data (I0,Q0,I1,Q1,… ,I26000, Q26000) is read. In case of complex data I’s and Q’s are put together as I+i ∙Q. Early and late code vectors of “blksize” elements are constructed from generated “caCode” by starting half a chip (earlyLateSpc = 0.5) earlier and later than the prompt code. The related phase series are also vectors of “blksize” values with step size of “codePhaseStep”.

Multiplying with the generated carrier removes the carrier and multiplying further with the right ge­ne­rated Code removes also the PRN code from the signal. The baseband data is multiplied by the early, late and prompt versions of generated C/A code for the delay lock loop – DLL and summed up for one Code length (1ms = PDIcarr = PDIcode = 0.001 sec).

The discriminator for the PLL/Costas loop is “arctan()” of prompt values of Q and I = carrError. This changes “carrNco” controlling the numerically controlled oscillator NCO. The NCO for the DLL for the code phase is controlled by variable “codeNco” influenced by the “codeError” calculated from early and late I and Q energies. If setting Cno.enableVSM=1, then the newly developed function “CNoVSM.m” is called

CNoVSM.m-function

It was developed during work on this thesis. It is called from “tracking.m”, gets 400 ms of the structure “trackResults” and calculates the carrier to noise ratio. See chapter 3.2 for details.

postProcessing.m-program, continued 4

Back from tracking.m-function the time needed for processing the data of 37000 ms is printed together with the message “Tracking is over ( elapsed time ...)”. Structures “trackResults”, “settings”, “acqResults” and “channel” are saved in file “trackingResults.mat”. The new code developed during this thesis has an addi­tio­nal “tracking2.m” that can be called several times for smaller amounts of data after a 37000 ms chunk is processed (with tracking.m and postNavigation.m) and had already de­li­vered the ephe­merides.

postNavigation.m-function

It is called from “postProcessing.m” with structures “trackResults” and “settings” as input. It computes the structures “navSolutions” (pseudoranges, receiver clock error, receiver coordinates) and “eph” (ephe­me­rides). A stop watch is started with “tic” at the beginning of the processing and stopped with “toc” at the end of the file. Processing is only done if “msToProcess” is at least 36000 (5 complete sub frames / frame each 6 seconds long + one potentially incomplete/unusable sub frame = 6 sub frames ∙ 6 sec / sub frame = 36 seconds and a grace second for optional skipping of samples at the start of the data file = 37 seconds) and at least 4 satellites were trackable (‘T’ in trackResults).

Navigation data is received with a bout rate of 50 (50 bits / sec., every 20 ms a new navigation bit). It is organized in frames (1500 bits/30sec) of 5 sub frames of 300 bits/6sec each. Sub frame 1-3 are con­stant­ly repeated: sub frame 1 contains satellite clock and health data, sub frame 2 and 3 contain satellite ephemeris data. Sub frames 4 and 5 are different in every frame and contain almanac data for other satellites. After 25 frames (25 ∙30sec = 750 sec = 12.5 min) 25 ∙2=50 sub frames of number 4 or 5 are sent with data from foreign satellites and the sequence is repeated. A sub frame is organized in 10 words a 30 bits (24 data + 6 parity bits). The first word of every sub frame is the telemetry word (TLM). The second word of every sub frame is the Handover Word (HOW). The TLM contains 8 bits of the preamble (fixed bit pattern) for frame synchronization. The preamble is followed by 16 reserved bits and parity information. All words end with 6 bit parity for data validation. Its values are calculated by multiple XORing of the 24 data bits of the word and the last 2 parity bits from the last word. For the structure of frames and sub frames of GPS L1 see (Borre 2007) and (Tsui 2005).

The first sub frame (of a frame) and its preamble for every satellite/PRN in “activeChnList” is found by the function “findPreamles.m”.

findPreambles.m-function

It is called from “postNavigation.m” with structures “trackResults” and “settings” as input parameters. Returned is an array “firstSubframe” with the positions of the first sub frame for all the satellites found in the “activeChnList” with all satellites with valid preambles that is also returned. The preamble is the fixed bit sequence [1 -1 -1 -1 1 -1 1 1]. With Matlab’ “kron()” the preamble is extended bit-wise to the 20-fold into vector “preamble_ms” because during a single bit of navigation data (and therefore also during a single bit of the preamble) 20 complete C/A codes are send:

preamble_ms = [ 1 1 1 … 1 -1 -1 -1 …-1 1 1 1 … 1]

The vector contains element 1-20 = 1, element 21 – 80 = -1, element 81-100 = 1, element 101-120 = -1, and finally ele­ment 121-160 = 1.

Variable “activeChnList” is filled with the sequential numbers of trackable satellites from “trackResults.status”. In a for-loop over the channel numbers contained in “activeChnList” a vector called ‘bits’ with the 37000 elements (ms or C/A periods) of I_P (in phase prompt) data is read from “trackResults(channelNr)”. Where this vector contains a negative value (or zero) a “-1” is set and where it contains a positive value a “+1” is set by:

bits(bits > 0) = 1; bits(bits <=0) = -1;

The bits-vector and the preamble_ms-vector are cross correlated with Matlab function xcorr() to “tlmXcorrResult”, a vector with 2 ∙37000 elements. One half of the elements are values of cross corre­lation and the other half are the related lags (displacements). Within the 2nd half elements of “tlmXcorrResult” the positions of values larger than 153 are searched and saved into variable “index”. The maximal cross correlation is 8 ∙ 20 = 160 (-160 for inverted preambles) theoretically. E. g. 24 posi­tions with such values are saved in “index” and not only 6 that are expected from 6 frames = 36000 ms. Because only data of a selected channel = a single satellite is analyzed – preambles should be one sub frame = 6 sec = 6000 ms be apart from each other. That is normally not the case.

In the loop “for i=1:size(index)” the distances (= 6000?) between a preamble and its forerunner are checked. Valid candidates for preambles must pass further checks. The parity bits for the first two words (TLM + HOW) of a sub frame are calculated and compared with transmitted values. Because also two parity bits from the predecessor word are needed for parity calculation 2+ 2 ∙30 = 62 bits must be read. Every bit is 20ms = 20 C/A-cycles long. Related to the start of a preamble we must read I_P (in phase prompt) values from “trackResults(channelNr)” from position -2 ∙20 = -40 to (60 ∙ 20) -1 = 1199. Theoretically 20 consecutive samples should have the same value because they belong to the same navigation bit, but that is not the case. Chunks of 20 samples are summed. If these sums are posi­tive the represented bit is a “+1”. Zero or negative sums represent a “-1”. The variable “bits” has 62 elements each a 1 or a -1. This is the input of function “include/navPartyChk.m” (two calls with bit 1-32 and with bit 31 - 62).

navPartyChk.m-function

It is found in folder “include”, called by “findPreamples.m” and returns ≠0 (corresponding to TRUE) if the parity check for a word (+ 2 last bits from forerunner parity) is successful, otherwise it returns a 0 (corresponding to FALSE). More precisely successful tests return a negated ndat(2) (D30* = last parity bit from predecessor word that causes a negation of the data bits before the XORing). Table 3 lists the bits used for parity calculations with names in (Borre 2007) and in the Matlab code:

Table 3: Navigation bits used for parity calculations

Abbildung in dieser Leseprobe nicht enthalten

findPreambles.m-function – continuation 1

Back from “navPartyChk.m” (two calls for the two words TLM + HOW) for all channels in variable “activeChnList” the start points (at which ms) the first sub frame (preamble) are received in variable “subFrameStart”.

postNavigation.m-function – continuation 1

Back from “findPreambles.m” the variable subFrameStart = ms informs where the preambles start for all active channels. In a loop „for channelNr=activeChnList“ a whole frame of 1500 bits (5 sub frames ∙ 300 bit/sub frame, equates to 1500 ∙ 20 code bits a 1ms = 30000 ms) plus 2 bits in front of the preamble (= last two parity bits from last word of predecessor frame = 2 ∙20 ms = 40 ms) of I_P (in phase prompt) data is read from “trackResults(channelNr)” into “navBitsSamples”. Chunks of 20 bits are summed up because 20 C/A-code bits are contained in 1 navigation bit. If the sum is positive the 20 are represented by a single 1 otherwise by a 0. Matlab internal function “dec2bin.m” returns a string of the binary representation to variable “navBitsBin”.

Function “ephemeris.m” is called with “navBitsBin” as input (1 complete frame + 2 last parity bits from predecessor) to get the Time Of Week – TOW – also called z-count, the time since GPS time rollover (mid­night Saturday→Sunday) in units ot 1.5 seconds with the number truncated by 2 bits (22 units = 4 units = 4 ∙1.5 sec = 6 sec). Z-count is Time Of Week with a reduced precision of 6 seconds from the first part of the second word (HOW) of a sub frame.

ephemeris.m-function

It is found in sub folder “include”, called by “postNavigation.m” and returns the Time Of Week – TOW from HOW. It checks length and data type of inputs and investigates in loop “for i=1:5” five sub frames of 300 navigational bits each. In the loop “for j = 1:10” a parity check for all ten words à 30 bits of the current sub frame is done by calling the function “parityCheck.m” from the include-folder.

parityCheck.m-function

[...]