Semi-Analytic Techniques for Fast MATLAB Simulations

1. Introduction

Advances in electronics and telecommunications are leading to complex systems able to efficiently use the available resources. Fast electronics, complex modulation schemes and correction codes enable transmissions on channels with unfavorable characteristics, coexistence between different services in the same frequency bands and high transmission rates. However, the complexity of such communications systems often prevents analytical characterizations. For example, figures of merit such as the BER are difficult to determine analytically for transmission schemes involving correction codes and communication channels with ISI and fading. In such cases, the system is characterized through Monte Carlo simulations [1, 2]. The Monte Carlo framework involves simulations of the whole system under analysis. For example, when considering a communications system, the whole transmission-reception chain is simulated. A large number of sequences are sent through the simulated system and the message recovered by the simulated receiver is compared to the original transmitted sequence. This comparison allows one to determine the average number of transmission errors and compute the BER.

Monte Carlo simulations can be applied to almost any system although their implementation and computation requirements can be significantly high. In addition to this, precision problems can arise when the quantity to be estimated is significantly low. For example, BER of the order of $${10}^{-8}-{10}^{-9}$$require at least $${10}^9-{10}^{10}$$simulation runs. Conversely, analytical models have a limited applicability and usually adopt approximations (i.e., model linearization) that can yield a poor description of the system under analysis.

In order to overcome the limitation of Monte Carlo and analytical techniques, semi-analytic approaches have been previously implemented [2, 1]. In a semi-analytic framework, the knowledge of the system under analysis is exploited to reduce the computational load and complexity that full Monte Carlo simulations would require. In this way, the strengths of both analytical and Monte Carlo methods are effectively combined. Semi-analytic techniques are a powerful tool for the analysis of complex systems.

The characteristics and relationships among the three aforementioned methods are shown in Figure 1 : semi-analytic approaches represent a good compromise in terms of applicability and complexity, combining the strengths of Monte Carlo and analytical approaches.

Figure 1.

Different approaches available for the analysis of complex systems. Semi-analytic approaches represent a compromise in terms of applicability and complexity (computational and implementation) between analytical models and Monte Carlo simulations.

The main goal of this chapter is to provide a general overview of semi-analytic techniques for the simulation of communications systems. Specific emphasis is given to their implementation in MT and two examples from the communications context are analyzed in detail.

Despite their potential, semi-analytic techniques have received limited attention from the communications community. Reference books on simulation of communications systems such as [1], [2] and [3] dedicate only a few pages to this kind of techniques. The focus is usually on the computation of the BER, which represents one of the first applications of semi-analytic techniques in communication system analysis [4, 5]. In this case, the model depicted in Figure 2 is adopted. The communications channel is modeled as a non-linear device, which distorts the signal component alone, with the addition of a noise term that is assumed to be Gaussian. This model is quite general and can be used to represent several communications channels.

Figure 2.

Model adopted for the evaluation of the BER using a semi-analytic approach. The communications channel is modeled as a non-linear device, which affects only the signal component, and the addition of a noise term supposed to be Gaussian.

A classical example is the model of a transmission chain where a TWTA is used to amplify the useful signal before transmission. The TWTA is highly non-linear and can lead to signal distortions. Since the signal is injected into the TWTA before transmission, the noise component entering the amplifier is negligible. In this case, the model depicted in Figure 2 is appropriate for describing the transmission chain including a non-linear amplifier.

The TWTA is a memory-less device and can be characterized using AM-AM conversionand AM-PM conversion(AM = Amplitude Modulation, PM = Phase Modulation) curves [3]. When a base-band signal model is used, the amplifier input and output are complex quantities; moreover, the response of the device usually depends only on the amplitude (instantaneous power) of the input signal. AM-AM and AM-PM conversion curves define the relationship between the input/output signal amplitudes and phases as a function of the input amplitude. Using these conversion curves, it is possible to simulate the behavior of the TWTA and other non-linear devices.

In the semi-analytic framework, the additivity of the noise component is exploited to compute the BER. More specifically, only the signal transmission chain is simulated and for each possible signal symbol, the Eb is computed. Since the noise properties are known, the BER for the $$i$$th symbol, $$s_i$$, is given by

$$BE{R}_i\mathrm{<mo>=</mo>}\frac{1}{2}erfc\left(\sqrt{\frac{E_{b\mathrm{,}i}}{N_0}}\right)\mathrm{,}$$E1

where $$E_{b\mathrm{,}i}$$has been obtained by simulating the transmission chain (including the non-linear device) in the absence of noise and transmitting the symbol$$s_i$$. The parameter $$N_0$$is the noise power spectral density and it is a known value of the system. Finally, the BER of the system is obtained from

$$BER\mathrm{<mo>=</mo>}E\left[BE{R}_i\right]\mathrm{<mo>=</mo>}\frac{1}{2}{\displaystyle \sum _{i\mathrm{<mo>=</mo>0}}^{N_s-1}}{p}_{i}erfc\left(\sqrt{\frac{E_{b\mathrm{,}i}}{N_0}}\right)\mathrm{,}$$E2

where $$p_i$$is the probability that the symbol $$s_i$$will be transmitted and $$N_s$$is the number of symbols of the signal constellation.

This simple example clearly illustrates the principles of semi-analytic techniques: the analytical knowledge of the system is exploited to reduce the computation load and complexity that full Monte Carlo simulations would require. In this case, only the transmission of the signal component is simulated.

In the literature, several generalizations of the aforementioned BER computation technique have been proposed. For example, [5] considered the case where the noise term at the input of the non-linear amplifier is not negligible. An equivalent model is proposed where the noise at the input is propagated after the non-linearity. [5] also considered the presence of a bandpass channel. More recently, [6] proposed a methodology for estimating the BER in the presence of ISI. All these examples show the potential and flexibility of the semi-analytic approach.

1.1. Building blocks

When considering the previous example it is possible to identify three building blocks that play different roles in the semi-analytic framework:

• simulation block

• estimation block

• analytical model.

The simulation block corresponds to that part of the system that is actually simulated. In the previous example, this block corresponds to the signal generation, the non-linear amplifier and the correlation receiver simulation. These blocks were used to determine the decision variable employed for recovering the transmitted symbol. The analytical model exploits the properties of the system to determine the quantities of interest. In the previous example, the fact that the noise introduced by the communication channel is Gaussian was exploited to determine the BER as a function of the Eb of each transmitted symbol. The estimation block is used as the interface between the simulated and analytical parts of the system. In BER computation case, the simulation part allows one to generate the different decision variables, whereas the analytical model is expressed as a function of energy per bit. The estimation block is used to determine Eb from the simulated decision variables.

The three functional blocks can be connected according to different configurations leading to different types of semi-analytic approaches. In the next section, two of these configurations are briefly discussed. Examples for each type of semi-analytic system are given in Section 2 and Section 3.

1.2. Main configurations

When considering the BER example, it is possible to note that the simulation, estimation and analytical blocks are connected in series. The simulation block is used at first to compute the different decision variables. The estimation block determines the Eb associated with each variable and finally the analytical model is used to compute the BER from the EN0.

This type of configuration is defined here as sequentialsince there is no feedback between the different blocks and each element of the chain is run sequentially. The principle of this configuration is shown in the upper part of Figure 3.

Figure 3.

Basic blocks and main configurations adopted in semi-analytic approaches. An estimation block is used for determining key signal and system parameters and interfacing the simulation and analytic components of the system.

A second type of configuration has been recently considered for the analysis of tracking loops in DSSS and GNSS receivers. The most computationally demanding operation in a DSSS/GNSS receiver is despreading, i.e, the correlation of the incoming samples with local replicas of the code and carrier. This operation is performed by the ID blocks that rely on simple operations that can be analytically modeled. For this reason, semi-analytic models exploiting the knowledge of the ID blocks and simulating only the non-linear parts of the system have been developed [7, 8, 9, 10]. This resulted in efficient analysis tools, which require reduced processing time with the applicability of Monte Carlo simulations.

Different techniques for modeling the output of the ID have been suggested. [7] modeled the correlator outputs evaluated by the ID blocks as linear combinations of independent Gaussian random variables. Correlation among the different correlators was obtained by using, for the generation of different ID outputs, a subset of the same random variables. This technique becomes complex as the number of correlators increases. Another attempt made in [8] assumed that the correlator outputs were independent which, in general, is not a realistic condition. Finally, [9, 10] suggested the use of a technique based on the Cholesky decomposition detailed in [11]. This approach allows one to easily generate an arbitrary number of correlated Gaussian random variables. In this way, [9, 10] were able to simulate advanced tracking loops for new GNSS signals.

Regardless of the type of approach used for modeling the correlator outputs, the aforementioned semi-analytic configuration can be represented as in the bottom part of Figure 3. In this case, an analytical model is used to generate quantities that will be propagated by simulation. In the tracking loop case, an analytical model is used to generate the correlator outputs that are then processed through simulations. The non-linear parts of the system are fully simulated and quantities such as the loop discriminator and filter outputs are computed. Finally, an estimation block is used to interface the simulation and analytical components of the semi-analytic scheme. A new estimate of the signal parameters is obtained and used to generate new correlator outputs. The estimation block is also used to compute performance metrics such as tracking jitter or MTLL [12].

The three functional blocks described in Section 1.1 are connected in a loop and thus this type of configuration is named closed-loop approach. The technique developed by [9, 10] and the closed-loop approach will be detailed in Section 3.

2. Sequential approach: The inter-system interference case

As shown in Figure 3, the sequential approach, at first, requires an initial simulation block to generate the random processes and system functions that cannot be analytically described. Subsequently, the simulated processes are employed by the estimation unit to obtain key parameters required by the analytical model to compute metrics of interest. In the analytical model, the estimated parameters are plugged into mathematical expressions to obtain the desired final metrics. In the sequential architecture, the gain in computational load mainly depends on the simplifications allowed by the analytical model. The use of such a model allows one to simulate only a part of the system and eventually avoid computationally demanding error counting processes.

The computational complexity of Monte Carlo simulations increases significantly when two or more communication systems coexist within the same environment. In this case, it is necessary to account for the interaction of the different systems and determine potential inter-system interference. In addition, computational requirements of Monte Carlo simulations increase dramatically with the number of random elements included in each block of the communication chain. These requirements can be considerably reduced by adopting semi-analytic techniques in which the evaluated metrics are analytically expressed as a function of parameters estimated through simulations. This principle is the core idea behind the sequential semi-analytic approach.

In order to better illustrate the principles of the semi-analytic sequential approach, an inter-system interference scenario is considered in this section. The case of a satellite navigation receiver affected by interference generated by a communications system is considered. The primary system (i.e., the victim system) considered here is a GPS L1 receiver affected by an interference signal caused by third order harmonics of a DVB-T signal. A comprehensive description of the sequential approach applied to an inter-system interference case is provided in the following.

The reception of GNSS signals is challenging due to low signal power, possible severe channel conditions and the presence of RF interference. The presence of RF interference can be particularly troublesome and the performance of a GNSS receiver can vary significantly depending on the type of interference. For this reason, significant research efforts have been devoted to the characterization of the receiver performance in the presence of different types of interference [13, 14]. Furthermore, the impact of interference originated by specific communication technologies, such as UWB transmissions [15], DME signals [16] and DVB-T harmonics [17, 18], has been thoroughly investigated. It is noted that the interference impact strongly depends on the strategy adopted by a GNSS receiver for processing the useful signals. Moreover, different impacts are expected depending on the receiver operating mode. The first task of a GNSS receiver is to determine the presence of a specific GNSS signal. This task is accomplished by the acquisition block that implements a statistical test for the detection of useful signals. After acquisition, the useful GNSS signals are passed to the tracking stage that refines the estimates of different signal parameters. Since acquisition and tracking implement different processing strategies, RF interference will affect these two receiver blocks differently. In the following, the acquisition stage is considered. A semi-analytic approach for the analysis of GNSS tracking loops is discussed in Section 3.

The acquisition of a GNSS signal can be formulated as a classical detection problem [19], where the signal of interest is buried in noise. The outcome of the acquisition process is twofold. First, a decision relative to the signal presence is provided. If the signal is present, a rough estimate of signal parameters (defined in the following) is also obtained. The received useful GNSS signal, which is impaired by AWGN and interference, is processed by the acquisition block yielding a decision variable. If the decision variable is greater than a decision threshold the signal presence is declared. This decision variable is calculated by using the digital samples provided by the receiver front-end. The signal model and the acquisition process are briefly summarized in the following sections.

2.1. The GNSS signal

The signal at the input of a GNSS receiver, in a one-path AWGN channel and in the presence of RF interference, can be modeled as

$$r\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{\displaystyle \sum _{l\mathrm{<mo>=</mo>0}}^{L-1}}{y}_l\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}+i\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}+\eta \mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>,}$$E3

where $$y_l\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}$$is the signal transmitted by the $$l$$th GNSS L1 satellite, $$L$$is the total number of satellites in view, $$i\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}$$is the received interference signal and $$\eta \mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}$$is the noise term.

Each useful signal, $$y_l\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}$$, can be expressed as

$$y_l\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}\sqrt{2C_l}{d}_l\left(t-{\tau }_{\mathrm{0,}l}\right)c_l\left(t-{\tau }_{\mathrm{0,}l}\right)\cos \left(2\pi \mathrm{<mo\; stretchy="false">(</mo>}{f}_{RF}+f_{d\mathrm{0,}l}\mathrm{<mo\; stretchy="false">)</mo>}t+{\phi }_{\mathrm{0,}l}\right)\mathrm{,}$$E4

where

• $$C_l$$$$l$$
• $$d_l\mathrm{<mo\; stretchy="false">(</mo>}\cdot \mathrm{<mo\; stretchy="false">)</mo>}$$
• $$c_l\mathrm{<mo\; stretchy="false">(</mo>}\cdot \mathrm{<mo\; stretchy="false">)</mo>}$$$$l$$
• $${\tau }_{\mathrm{0,}l}$$$$f_{d\mathrm{0,}l}$$$${\phi }_{\mathrm{0,}l}$$
• $$f_{RF}$$

It is noted that GNSS signals adopt a DSSS modulation. The pseudo-random sequences, $$c_l\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}$$, allow the simultaneous transmission of several signals at the same time and in the same band. Moreover, $$c_l\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}$$sequences are characterized by sharp correlation functions that allow the precise measurement of the signal travel time. The travel time is then converted into distances that allows a GNSS receiver to determine its position.

The pseudo-random sequence, $$c_l\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}$$, is composed of several terms including a primary spreading sequence and a subcarrier:

$$c_l\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{\displaystyle \sum _{i\mathrm{<mo>=</mo>}-\infty }^{+\infty }}{c}_{l\mathrm{,<mo\; stretchy="false">(</mo>}i\mathrm{mod}{N}_c\mathrm{<mo\; stretchy="false">)</mo>}}{s}_b\mathrm{<mo\; stretchy="false">(</mo>}t-i{T}_h\mathrm{<mo\; stretchy="false">)</mo>.}$$E5

The signal $$s_b\mathrm{<mo\; stretchy="false">(</mo>}t-i{T}_h\mathrm{<mo\; stretchy="false">)</mo>}$$in (5) represents the subcarrier of duration$$T_h$$, which determines the spectral characteristics of the transmitted GNSS signal. The GPS L1 CA component is BPSK modulated, whereas the Galileo E1 signal adopts a CBOC scheme. The sequence$$c_{l\mathrm{,}i}$$, of length$$N_c$$, defines the primary spreading code of the $$l$$th GNSS signal. In the following, only the BPSK case is considered. The results can be easily extended to different subcarriers.

A GNSS receiver is able to process the $$L$$useful signals independently since the spreading codes are quasi-orthogonal. Therefore, expression (3) can be simplified to

$$r\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}y\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}+i\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}+\eta \mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>,}$$E6

where the index $$l$$has been dropped for ease of notation.

The received signal in (6) is filtered and down-converted by the receiver front-end. Filtering is of particular importance since it determines which portion of the interfering signal, $$i\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}$$, will effectively enter the receiver. After down-conversion and filtering, the input signal is sampled and quantized. In this analysis, the impact of quantization and sampling is neglected. After these operations, (6) becomes:

$$r_{BB}\mathrm{<mo\; stretchy="false">[</mo>}n\mathrm{<mo\; stretchy="false">]</mo><mo>=</mo>}{y}_{BB}\mathrm{<mo\; stretchy="false">(</mo>}n{T}_s\mathrm{<mo\; stretchy="false">)</mo>}+i_{BB}\mathrm{<mo\; stretchy="false">(</mo>}n{T}_s\mathrm{<mo\; stretchy="false">)</mo>}+{\eta }_{BB}\mathrm{<mo\; stretchy="false">(</mo>}n{T}_s\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{y}_{BB}\mathrm{<mo\; stretchy="false">[</mo>}n\mathrm{<mo\; stretchy="false">]</mo>}+i_{BB}\mathrm{<mo\; stretchy="false">[</mo>}n\mathrm{<mo\; stretchy="false">]</mo>}+{\eta }_{BB}\mathrm{<mo\; stretchy="false">[</mo>}n\mathrm{<mo\; stretchy="false">]</mo>,}$$E7

where the notation $$x\mathrm{<mo\; stretchy="false">[</mo>}n\mathrm{<mo\; stretchy="false">]</mo>}$$is used to denote discrete time sequences sampled at the frequency$$f_s\mathrm{<mo>=</mo>}\frac{1}{T_s}$$. In addition, the index “BB” is used to denote a filtered signal down-converted to baseband. Furthermore, the signal $$y_{BB}\mathrm{<mo\; stretchy="false">[</mo>}n\mathrm{<mo\; stretchy="false">]</mo>}$$in (7) can be written as

$$y_{BB}\mathrm{<mo\; stretchy="false">[</mo>}n\mathrm{<mo\; stretchy="false">]</mo><mo>=</mo>}\sqrt{C}d\left(n{T}_s-{\tau }_0\right)c\left(n{T}_s-{\tau }_0\right)\exp \left\{j2\pi f_{0}n{T}_s+{\phi }_0\right\}\mathrm{.}$$E8

The noise term, $${\eta }_{BB}\mathrm{<mo\; stretchy="false">[</mo>}n\mathrm{<mo\; stretchy="false">]</mo>}$$, is AWGN with variance$${\sigma }_{BB}^2$$. This variance depends on the filtering, down-conversion and sampling strategy applied by the receiver front-end and can be expressed as$${\sigma }_{BB}^2\mathrm{<mo>=</mo>}{N}_0{B}_{RX}$$, where $$B_{RX}$$is the front-end bandwidth and $$N_0$$is the PSD of the input noise$$\eta \mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}$$. The ratio between the carrier power, $$C$$, and the noise PSD, $$N_0$$, defines the CN, one of the main signal quality indicators used in GNSS.

2.2. The DVB-T interfering signal

The interference term in (6), $$i\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}$$, originates from DVB-T emissions. The DVB-T system is the European standard for the broadcasting of digital terrestrial television signals and has been adopted in many countries, mainly in Europe, Asia and Australia. The standard employs an OFDM-based modulation scheme operating in the VHF III ($$174-230$$MHz), UHF IV ($$470-582$$MHz) and UHF V ($$582-862$$MHz) bands [20]. It is noticeable that none of these bands fall within the bands allocated for GNSS signals. However, the second harmonics of UHF IV and third order harmonics of UHF V could coincide with the GPS L1 band, and, thus cause harmful interference. The case of the third order harmonics of a DVB-T signal is considered here. The DVB-T transmitted signal can be represented as

$$i_{DVB-T}\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}\frac{1}{\sqrt{M}}{\displaystyle \sum _{p\mathrm{<mo>=</mo>0}}^{N_d-1}}{\displaystyle \sum _{h\mathrm{<mo>=</mo>0}}^{M-1}}{I}_{p\mathrm{,}h}\exp \left(\frac{j2\pi ht}{M}\right)\mathrm{,}$$E9

where $$M$$is the modulation order, $$h$$is the subcarrier index, $$p$$is the symbol index, $$N_d$$represents the total number of transmitted symbols and $$I_{p\mathrm{,}h}$$models the $$h$$th constellation point of the $$p$$th symbol. Here, the term “subcarrier” should not be confused with the subcarrier used to modulate the GNSS signals in (5). In the OFDM context, several components are transmitted in parallel on different overlapping frequency bands. The term subcarrier denotes each individual transmitted component. In GNSS, the subcarrier is an additional component that modulates the transmitted signal and plays a role analogous to the carrier used for the signal up-conversion.

Third order harmonics are the consequence of the malfunctioning of the transmitter electronics. In particular, the presence of these harmonics are due to the non-linearities of an amplifier. The output of an amplifier can be modeled using a polynomial expansion of the amplifier input/output function:

$$p\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{\displaystyle \sum _{n\mathrm{<mo>=</mo>1}}^{\infty }}{a}_n{i}_{DVB-T}^n\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>,}$$E10

where $$a_n$$are the polynomial coefficients of the Taylor series expansion of the amplifier input/output function. This type of model is an alternative to the AM-AM and AM-PM conversion functions discussed in Section 1 for the TWTA case.

The terms of order $$n\mathrm{<mo>></mo>1}$$in (10) model the amplifier non-linearities and the ratios $$a_n\mathrm{<mo>/</mo>}{a}_1$$are expected to be small for$$n\mathrm{<mo>></mo>1}$$. Since only the third harmonics will fall into the GPS L1 band, the interference signal at the antenna of a GNSS receiver is given by

$$i\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{a}_3{i}_{DVB-T}^3\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>.}$$E11

The signal $$i\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}$$is filtered and down-converted by the receiver front-end and signal$$i_F\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}$$, the filtered version of$$i\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}$$, will affect receiver operations.

Finally, the interference term in (7) can be modeled as$$i_{BB}\mathrm{<mo\; stretchy="false">[</mo>}n\mathrm{<mo\; stretchy="false">]</mo><mo>=</mo>}{i}_F\mathrm{<mo\; stretchy="false">(</mo>}n{T}_s\mathrm{<mo\; stretchy="false">)</mo>}$$.

2.3. The acquisition process

After signal conditioning, the sequence $$r_{BB}\mathrm{<mo\; stretchy="false">[</mo>}n\mathrm{<mo\; stretchy="false">]</mo>}$$is correlated with local replicas of the useful signal code and carrier as shown in Figure 4. Since the code delay, $${\tau }_0$$, and the Doppler frequency, $$f_0$$, of the useful signal in (8) are unknown to the receiver, several delays and frequencies are tested by the acquisition block. In addition to this, several correlators, computed using subsequent portions of the input signal$$r_{BB}\mathrm{<mo\; stretchy="false">[</mo>}n\mathrm{<mo\; stretchy="false">]</mo>}$$, can be computed in order to produce a decision variable less affected by noise and interference. In this way, the output of the $$k$$th complex correlator can be expressed as

Figure 4.

Schematic representation of the operations performed by the acquisition block of a GNSS receiver.

$$S_k\mathrm{<mo>=</mo>}\frac{1}{N}{\displaystyle \sum _{n\mathrm{<mo>=</mo>}kN}^{\mathrm{<mo\; stretchy="false">(</mo>}k+\mathrm{1<mo\; stretchy="false">)</mo>}N-1}}{r}_{BB}\mathrm{<mo\; stretchy="false">[</mo>}n\mathrm{<mo\; stretchy="false">]</mo>}c\left(n{T}_s-\tau \right)\exp \left\{-j2\pi f_{d}n{T}_s-j\phi \right\}\mathrm{,}$$E12

where$$\tau$$, $$f_d$$and $$\phi$$are the code delay, Doppler frequency and carrier phase tested by the receiver. The parameter $$N$$is the number of samples used for computing a single correlation output and $$T_c\mathrm{<mo>=</mo>}N{T}_s$$defines the coherent integration time. It is noted that the computation of correlation outputs is essential for the proper functioning of a GNSS receiver and they are both used in acquisition and tracking modes [21]. To further improve the acquisition performance, non-coherent integration can be implemented as illustrated in Figure 4. More specifically, the impact of the navigation message, $$d\mathrm{<mo\; stretchy="false">(</mo>}\cdot \mathrm{<mo\; stretchy="false">)</mo>}$$, is removed through squaring, $$\mathrm{<mo>|</mo>}{S}_k{\mathrm{<mo>|</mo>}}^2$$, and the final decision variable is computed as

$$D\mathrm{<mo>=</mo>}\frac{1}{K}{\displaystyle \sum _{k\mathrm{<mo>=</mo>0}}^{K-1}}\mathrm{<mo>|</mo>}{S}_k{\mathrm{<mo>|</mo>}}^2\mathrm{,}$$E13

where $$K$$is the total number of correlation samples that are non-coherently integrated. It should be noted that for $$K\mathrm{<mo>=</mo>1}$$only coherent integration is used. In order to determine the signal presence, the receiver compares $$D$$with a decision threshold,$$\beta$$. If $$D$$is greater than $$\beta$$then the useful signal is declared present.

It is noted that, as in any binary test, two hypotheses are possible:

• $$H_0$$
• $$H_1$$

The null hypothesis, $$H_0$$assumes that the correlator outputs, $$S_k$$, are made of noise alone. Since the pseudo-random sequences, $$c\mathrm{<mo\; stretchy="false">(</mo>}\cdot \mathrm{<mo\; stretchy="false">)</mo>}$$, are selected to have good autocorrelation properties, if the code delay and Doppler frequencies tested by the receiver do not match the parameters of the input signal, $$y_{BB}\mathrm{<mo\; stretchy="false">[</mo>}n\mathrm{<mo\; stretchy="false">]</mo>}$$, then the useful signal component is almost completely filtered out at the correlator output. Thus, also in this case, the $$H_0$$hypothesis is verified.

Furthermore, $$H_1$$is the alternative hypothesis and assumes that the signal is present and the local code and carrier replica are perfectly aligned. If $$H_1$$is declared, then rough estimates of $${\tau }_0$$and $$f_0$$are also obtained.

Depending on the result of the test, $$D\mathrm{<mo>></mo>}\beta$$, two decisions can be taken by the receiver

• $$D_0$$
• $$D_1$$

and the four events described in Table 1 can occur.

	$$D_0$$	$$D_1$$
$$H_0$$	Signal absence correctly declared	False alarm
$$H_1$$	Missed detection	Signal detection

Table 1.

Confusion matrix describing the four events that can happen in the binary test performed by the acquisition process.

The off-diagonal events in Table 1 correspond to the different errors that the acquisition block can commit. The following probabilities are usually associated with the events in Table 1 :

$$P_d\mathrm{<mo\; stretchy="false">(</mo>}\beta \mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}Prob\mathrm{<mo\; stretchy="false">(</mo>}D\mathrm{<mo>></mo>}\beta \mathrm{<mo>|</mo>}{H}_1\mathrm{<mo\; stretchy="false">)</mo>}Probabilityofdetection$$E14

and

$$P_{fa}\mathrm{<mo\; stretchy="false">(</mo>}\beta \mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}Prob\mathrm{<mo\; stretchy="false">(</mo>}D\mathrm{<mo>></mo>}\beta \mathrm{<mo>|</mo>}{H}_0\mathrm{<mo\; stretchy="false">)</mo>}Probabilityoffalsealarm\mathrm{.}$$E15

The probabilities of missed detection and correct signal absence decision are obtained as $$1-P_d\mathrm{<mo\; stretchy="false">(</mo>}\beta \mathrm{<mo\; stretchy="false">)</mo>}$$and$$1-P_{fa}\mathrm{<mo\; stretchy="false">(</mo>}\beta \mathrm{<mo\; stretchy="false">)</mo>}$$, respectively. The performance of the acquisition process is characterized in terms of ROC [22], which plots the detection probability as a function of the false alarm rate. ROC curves capture the behavior of a detector as a function of the different decision thresholds.

2.4. The semi-analytic approach

The goal of the sequential semi-analytic approach considered in this section is the evaluation of the ROC in the presence of DVB-T interference. The full Monte Carlo approach would consist of simulating the full transmission/reception scheme shown in Figure 5 and generating several realizations of $$D$$both under $$H_0$$and$$H_1$$. Probabilities of detection and false alarm would then be determined through error counting techniques. This approach is computationally demanding and does not exploit the analytical knowledge of the system. More specifically, under the hypothesis that the correlator outputs $$S_k$$are iid complex Gaussian random variables with independent real and imaginary parts, it is possible to show [23] that

$$P_{fa}\mathrm{<mo\; stretchy="false">(</mo>}\beta \mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}\exp \left\{-\frac{\beta }{2{\sigma }_n^2}\right\}{\displaystyle \sum _{i\mathrm{<mo>=</mo>0}}^{K-1}}\frac{1}{i\mathrm{<mo>!</mo>}}{\left(\frac{\beta }{2{\sigma }_n^2}\right)}^i$$E16

and

$$P_{det}\mathrm{<mo\; stretchy="false">(</mo>}\beta \mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}Q{}_K\left(\sqrt{K\lambda }\mathrm{<mo>;</mo>}\sqrt{\beta {\sigma }_n^2}\right)\mathrm{,}$$E17

where $$Q_K\mathrm{<mo\; stretchy="false">(</mo>}a\mathrm{<mo>;</mo>}b\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{\displaystyle {\int }_b^{+\infty }}x{\left(\frac{x}{a}\right)}^{K-1}\exp \left\{-\frac{x^2+a^2}{2}\right\}{I}_{K-1}\mathrm{<mo\; stretchy="false">(</mo>}ax\mathrm{<mo\; stretchy="false">)</mo>}dx$$is the generalized Marcum Q-function of order$$K$$. The function $$I_K\mathrm{<mo\; stretchy="false">(</mo>}\cdot \mathrm{<mo\; stretchy="false">)</mo>}$$is the modified Bessel function of first kind and order$$K$$. In (16) and (17), $${\sigma }_n^2$$is the variance of the real and imaginary part of$$S_k$$. The parameter $$\lambda$$is given by

$$\lambda \mathrm{<mo>=</mo>}\frac{{\left|E\left[S_k\right]\right|}^2}{{\sigma }_n^2}$$E18

and defines the SNR at the correlator outputs.

Figure 5.

Schematic representation of the full Monte Carlo simulation system for the ROC evaluation of the acquisition of a GPS L1 signal in the presence of DVB-T third order harmonics.

The correlator output can be considered iid complex Gaussian random variables even in the presence of DVB-T interference. More specifically, the large number of terms in the sum performed in (12) allows one to invoke the central limit theorem and assume $$S_k$$is Gaussian. Independence derives from the down-sampling performed by the correlators. Since only one correlator is produced every $$N$$samples, the statistical correlation between subsequent correlators is significantly reduced. The lack of correlation translates into independence for Gaussian random variables. Thus, models (16) and (17) can be used and the only parameters that need to be estimated are $${\sigma }_n^2$$and$$\lambda$$.

The analytical knowledge of the system can be further exploited to simplify the evaluation of $${\sigma }_n^2$$and$$\lambda$$. In particular, since $$i_{BB}\mathrm{<mo\; stretchy="false">[</mo>}n\mathrm{<mo\; stretchy="false">]</mo>}$$and $${\eta }_{BB}\mathrm{<mo\; stretchy="false">[</mo>}n\mathrm{<mo\; stretchy="false">]</mo>}$$are modeled as zero mean random processes, they only contribute to the variance of the correlator outputs. Thus, neglecting residual errors due to delay and frequency partial misalignments, yields

$${\left|E\left[S_k\right]\right|}^2\mathrm{<mo>=</mo>}C\mathrm{,}$$E19

where $$C$$is the useful signal power and is one of the known parameters of the system. Thus, $$\lambda$$can be derived from $$C$$and$${\sigma }_n^2$$.

Finally, exploiting the linearity of the correlation process, it is possible to express $$S_k$$as

$$S_k\mathrm{<mo>=</mo>}{S}_{r\mathrm{,}k}+S_{\eta \mathrm{,}k}+S_{i\mathrm{,}k}\mathrm{,}$$E20

which is a linear combination of a useful signal term, derived from$$y_{BB}\mathrm{<mo\; stretchy="false">[</mo>}n\mathrm{<mo\; stretchy="false">]</mo>}$$, a noise term, derived from$${\eta }_{BB}\mathrm{<mo\; stretchy="false">[</mo>}n\mathrm{<mo\; stretchy="false">]</mo>}$$, and an interference term derived from$$i_{BB}\mathrm{<mo\; stretchy="false">[</mo>}n\mathrm{<mo\; stretchy="false">]</mo>}$$. The variance $${\sigma }_n^2$$can be obtained as

$${\sigma }_n^2\mathrm{<mo>=</mo>}\frac{1}{2}Var\left\{S_k\right\}\mathrm{<mo>=</mo>}\frac{1}{2}Var\left\{S_{\eta \mathrm{,}k}\right\}+\frac{1}{2}Var\left\{S_{i\mathrm{,}k}\right\}\mathrm{.}$$E21

Figure 6.

Schematic representation of the semi-analytic approach adopted for the evaluation of the ROC in the presence of DVB-T interference. The three functional elements of the semi-analytic approach are highlighted in different colors.

Using the results derived in [24], [13] and [23], it is possible to show

$$\frac{1}{2}Var\left\{S_{\eta \mathrm{,}k}\right\}\mathrm{<mo>=</mo>}\frac{N_0}{2T_c}$$E22

and

$$\frac{1}{2}Var\left\{S_{i\mathrm{,}k}\right\}\mathrm{<mo>=</mo>}\frac{C_i}{2T_c}{k}_a\mathrm{,}$$E23

where $$C_i$$is the interference power and $$k_a$$is the SSC defined as [24, 13]

$$k_a\mathrm{<mo>=</mo>}{\displaystyle {\int }_{-B_{RX}\mathrm{<mo>/</mo>2}}^{B_{RX}\mathrm{<mo>/</mo>2}}}{G}_i\mathrm{<mo\; stretchy="false">(</mo>}f\mathrm{<mo\; stretchy="false">)</mo>}{G}_c\mathrm{<mo\; stretchy="false">(</mo>}f\mathrm{<mo\; stretchy="false">)</mo>}df\mathrm{.}$$E24

The function $$G_i\mathrm{<mo\; stretchy="false">(</mo>}f\mathrm{<mo\; stretchy="false">)</mo>}$$in (24) is the normalized PSD of the DVB-T interference signal after front-end filtering. In addition, $$G_i\mathrm{<mo\; stretchy="false">(</mo>}f\mathrm{<mo\; stretchy="false">)</mo>}$$is normalized such that

$${\displaystyle {\int }_{-B_{RX}\mathrm{<mo>/</mo>2}}^{B_{RX}\mathrm{<mo>/</mo>2}}}{G}_i\mathrm{<mo\; stretchy="false">(</mo>}f\mathrm{<mo\; stretchy="false">)</mo>}df\mathrm{<mo>=</mo>1.}$$E25

The function $$G_c\mathrm{<mo\; stretchy="false">(</mo>}f\mathrm{<mo\; stretchy="false">)</mo>}$$models the effect of the correlation on the interfering signal. Correlation can be modeled as an additional filtering stage and $$G_c\mathrm{<mo\; stretchy="false">(</mo>}f\mathrm{<mo\; stretchy="false">)</mo>}$$can be shown to be well approximated by the PSD of the subcarrier used in the despreading process. Also, $$G_c\mathrm{<mo\; stretchy="false">(</mo>}f\mathrm{<mo\; stretchy="false">)</mo>}$$is normalized to have a unit integral. It is noted that different subcarriers lead to different $$G_c\mathrm{<mo\; stretchy="false">(</mo>}f\mathrm{<mo\; stretchy="false">)</mo>}$$and thus, $$i_{BB}\mathrm{<mo\; stretchy="false">[</mo>}n\mathrm{<mo\; stretchy="false">]</mo>}$$will have different effects depending on the type of modulation considered.

The only unknown parameter in the previous equation is the SSC, which needs to be estimated using Monte Carlo simulations. Also, the interfering DVB-T signal is fully simulated. The resulting signal is filtered by the receiver front-end and the sequence $$i_{BB}\mathrm{<mo\; stretchy="false">[</mo>}n\mathrm{<mo\; stretchy="false">]</mo>}$$is obtained. The samples of $$i_{BB}\mathrm{<mo\; stretchy="false">[</mo>}n\mathrm{<mo\; stretchy="false">]</mo>}$$are used to estimate the normalized PSD,$$G_i\mathrm{<mo\; stretchy="false">(</mo>}f\mathrm{<mo\; stretchy="false">)</mo>}$$. This can be easily obtained using the MT functions developed for spectral analysis. In this case, the pwelch function is used. The function $$G_i\mathrm{<mo\; stretchy="false">(</mo>}f\mathrm{<mo\; stretchy="false">)</mo>}$$is used to compute the SSC, which is then used to determine the system ROC.

The developed semi-analytic approach is shown in Figure 6 where the simulation, estimation and analytic components are clearly highlighted.

Figure 7.

Representation of a normalized PSD realization of the third harmonic of the DVB-T signal and the frequency response of the GPS L1 front-end filter.

2.5. Performance comparison

A comparison between a full Monte Carlo simulation and a semi-analytic technique, implemented for the evaluation of the acquisition performance of a GPS L1 receiver impaired by third order harmonics of a DVB-T signal, is presented in this section. Initially, the DVB-T interfering signal in time domain is programmed in MT by following the DVB-T standard and the non-linear amplifier model, as illustrated in Figure 5. Note that the simulation of the interfering signal is required for both Monte Carlo and semi-analytic techniques. Subsequently, the estimated PSD of the interfering signal, needed for the estimation of the SSC in the semi-analytic method, is obtained by applying the pwelch function of MT. A realization of the normalized PSD of the interfering signal is depicted in Figure 7. The centre frequency of the interfering signal is set to$$f_I\mathrm{<mo>=</mo>}{f}_{RF}+\Delta f$$, where $$\Delta f$$is the frequency shift of the interference signal with respect to the centre frequency of the GPS L1 signal. The impact of selecting different values of $$\Delta f$$on the acquisition performance of a GPS L1 receiver is analyzed in [25]. Furthermore, the frequency response of the GPS L1 front-end filter is also plotted in Figure 7. In this case, the selected filter bandwidth is $$8$$MHz.

Sample results comparing ROC curves obtained using semi-analytic and Monte Carlo simulations are shown in Figure 8. The parameters used for the analysis are reported in Table 2. From Figure 8, it can be observed that the Monte Carlo and semi-analytic approaches provide similar results and the curves obtained using the two methods overlap. However, the semi-analytic approach provides increased precision, particularly when small values need to be estimated, and a significant reduction in terms of computational complexity. Full Monte Carlo simulations require the implementation of the full transmission/reception chain and the evaluation of the ROC with a computational complexity significantly higher than that of the semi-analytic approach described above.

Figure 8.

Comparison between ROC curves obtained using semi-analytic and Monte Carlo simulations. The semi-analytic framework considered provides increased precision and requires a lower computational complexity.

Parameter	Value
CN	$$35$$dB-Hz
Coherent integration time,$$T_c$$	$$1$$ms
Interference to signal power ratio,$$\frac{C_i}{C}$$	$$30$$dB
Centre frequency difference,$$\Delta f$$	$$0$$Hz
Receiver bandwidth,$$B_{RX}$$	$$8$$MHz
Number of Monte Carlo Simulation runs	$${10}^6$$

Table 2.

Parameters used for the evaluation of the ROC curves shown in Figure 8.

Additional results relative to the impact of DVB-T interference on GNSS can be found in [25].

3. Closed-loop approach: Digital tracking loops

As anticipated in Section 1, a second configuration, called closed-loop approach, has been recently proposed for the simulation of digital tracking loops in DSSS/GNSS receivers. The SATLSim toolbox is a set of MT functions implementing the semi-analytic closed-loop approach for the analysis of digital tracking loops. The SATLSim toolbox has been developed by [9, 10] and can be downloaded from the following websites:

• http://www.ngs.noaa.gov/gps-toolbox/SATLSim.htm

• http://plan.geomatics.ucalgary.ca/publications.php.

In the following, the closed-loop approach for the simulation of digital tracking loops is considered and the MT code developed in the SATLSim toolbox is briefly analyzed.

A description of the correlator model used for reducing the computational complexity of the system is at first provided. The samples given by (7) at the input of a GNSS receiver are processed by the different functional blocks with different objectives. The acquisition process described in Section 2 is the first stage of a GNSS receiver and has the goal of determining the signal presence and provide a rough estimate of its parameters. These parameters include the code delay $${\tau }_0$$and Doppler frequency$$f_0$$.

If the signal is successfully acquired then different tracking loops are used to refine the estimate of the signal parameters. A DLL is usually used to provide accurate estimates of the code delay, $${\tau }_0$$, and track delay variations due to the relative motion between receiver and satellite. The Doppler frequency, $$f_0$$, is recovered using either a FLL or a PLL. If a PLL is used then the carrier phase, $${\phi }_0$$, is also estimated. The code delay and carrier phase allow the receiver to determine its position whereas the Doppler frequency can be used for computing the user velocity.

As indicated in Section 2, a subcarrier can be used for shaping the spectrum of the transmitted GNSS signal and improving its robustness against multipath. The presence of a subcarrier makes code tracking more complex since the correlation function of the transmitted signal may have multiple peaks. More specifically, fine delay estimation is obtained by maximizing the correlation between input signal and local code: the correlation function is maximized only when the delay of the locally generated code matches the delay of the input signal. The presence of several peaks in the correlation function may cause the DLL to converge to a local maximum causing biases in the delay estimation. For this reason, several solutions have been proposed to avoid lock on secondary correlation peaks [26, 27]. An effective solution is represented by the SLL proposed by [27]. In this case, the subcarrier is seen as a periodic waveform that further modulates the transmitted signal. The delays of code and subcarrier are decoupled and estimated separately. In this way, the ambiguous one-dimensional signal correlation is projected in an unambiguous bi-dimensional function. In the following, the joint simulation of DLL and SLL is considered.

In a GNSS tracking loop, the incoming signal is correlated with several locally generated code and carrier replicas and different correlator outputs are produced. This process is analogous to the correlation operations described in Section 2 and is performed by the ID blocks.

Each correlator output is a function of the input signal and the parameters previously estimated by the tracking loop. The correlator outputs are passed to the non-linear discriminator that produces a first estimate of the tracking error that the loop is trying to minimize. The tracking error is filtered and passed to the NCO that is used for generating new local signal replicas.

Efficient tracking loop simulations can be obtained by substituting the ID blocks with their analytical model. More specifically, a correlator output can be modeled as:

$$\sqrt{C}\frac{\sin \left(\pi \Delta f_d{T}_c\right)}{\pi \Delta f_d{T}_c}{R}_l\left(\Delta {\tau }_d\mathrm{,}\Delta {\tau }_s\right)\exp \left\{j\Delta \phi \right\}+{\eta }_c\mathrm{,}$$E26

where

• $$\Delta f_d$$$$\Delta \phi$$
• $$\Delta {\tau }_d$$$$\Delta {\tau }_s$$$$\Delta {\tau }_s$$
• $$T_c\mathrm{<mo>=</mo>}N{T}_s$$$$N$$
• $$R_l\left(\Delta {\tau }_d\mathrm{,}\Delta {\tau }_s\right)$$$$R_l\left(\Delta {\tau }_d\mathrm{,}\Delta {\tau }_s\right)$$
• $${\eta }_c$$$${\eta }_c$$

From (26), it is possible to reconstruct the correlator outputs given the estimation errors generated by the tracking loops. Thus, the correlation process does not need to be simulated and only the estimation errors are determined using a Monte Carlo approach. Based on this principle, the simulation scheme shown in Figure 9 can be adopted for the fast simulation of digital tracking loops.

Figure 9.

Semi-analytic scheme adopted for the simulation of GNSS tracking loops. Each element of the scheme proposed for the analysis of tracking loops has been implemented in a different function of the SATLSim toolbox.

The functional elements in Figure 9 have been grouped to form the simulation block, the analytical model and the estimation part. The analytical model is used to convert the signal parameter errors, $$\Delta f_d$$, $$\Delta \phi$$, $$\Delta {\tau }_d$$and$$\Delta {\tau }_s$$, into the signal components of the correlator outputs. At the same time, the analytical model is used to determine the variance and correlation of the different noise terms used to simulate$${\eta }_c$$. Since the noise components are simulated using parameters determined by the analytical model, the “noise generation” block is shared between the analytical and simulation parts. The remaining parts of the loop, including the non-linear discriminator, loop filter and NCO, are fully simulated. Finally, the estimation block determines the residual signal parameter errors by comparing true values (determined by the simulation scenario) and estimates produced by the NCO.

By modifying these functional blocks, it is possible to simulate different tracking loops. In the simulation scheme implemented in the SATLSim toolbox, a new estimate of the tracking parameters (Doppler frequency, carrier phase and code and subcarrier delays) is generated by an NCO model. This model accounts for the integration process performed by a real NCO and different update equations can be used [28]. A commonly used model is the rate-only feedback NCO [28], characterized by the following update equation:

$${\widehat{\phi }}_k\mathrm{<mo>=</mo>}{\widehat{\phi }}_{k-1}+\frac{T_c}{2}\left(\delta {\widehat{\phi }}_{k-1}+\delta {\widehat{\phi }}_{k-2}\right)\mathrm{,}$$E27

where $${\widehat{\phi }}_k$$denotes the $$k$$-th estimate of the tracking parameter under consideration and $$\delta {\widehat{\phi }}_k$$is its estimated rate of change. The rate $$\delta {\widehat{\phi }}_k$$is generally provided by the loop filter. It is noted that when several parameters are considered, equation (27) is used to update each term independently. The new parameter estimate is compared to the true value and a new estimation error is computed. This error is then used for the generation of the signal component at the output of the ID block using equation (26). The noise term, generated separately, is then added to the signal component. When several correlators are required, the correlation among the different noise components has to be accounted for. This is simulated using the approach described in [9].

The operations required to convert the correlator outputs into a new estimate of the parameter rate, $$\delta {\phi }_k$$, are fully simulated and correspond to the functional blocks that can be found in a real tracking loop. For instance, the correlator outputs are used to update the nonlinear discriminator and the loop filter. It is noted that a similar simulation scheme can be used for analyzing Kalman filter based tracking. In this case, the correlator outputs are fed to a Kalman filter that is used to produce new estimates of the tracking parameters.

3.1. Code structure

The structure of the code developed in the SATLSim toolbox is provided in Figure 10. In this case, the code is used to estimate the tracking jitter of the loop as a function of different parameters, such as the Early-minus-Late spacing and the input CN. In particular, the non-linear discriminator may use several correlators to compute the cost function that the loop is trying to minimize. A DLL usually requires at least two correlators, named Early and Late correlators, computed for the delays

$$\widehat{\tau }\pm \frac{1}{2}{d}_s$$E28

where $$\widehat{\tau }$$is the best code delay estimate and $$d_s$$is the Early-minus-Late spacing. Early and Late correlators are computed symmetrically with respect to the best delay estimate and the non-linear discriminator computes a cost function proportional to the misalignment between these two correlators. Since the code correlation function is symmetric, the output of the discriminator is minimized when $$\widehat{\tau }$$corresponds to the delay of the input signal. The SLL works using similar principles. The performance of DLL and SLL depends on the Early-minus-Late spacing that is a simulation parameter. The CN is used to determine the correlator amplitude and the variance of the noise component,$${\eta }_c$$.

Figure 10.

Structure of the SATLSim toolbox and list of the different MT functions.

The parameters required for initializing the simulation procedure are accessible through the function InitSettings. These parameters include the sampling frequency, the loop bandwidth and the coherent integration time that are used to design the loop filters through the function FilterDesign. In the code provided, standard formulae from [21] are used. However, FilterDesign can be modified in order to adopt a different approach, such as the controlled-root formulation proposed by [28]. During the initialization phase, the true input parameters are also generated. The simulation core consists of three nested loops, on the Early-minus-Late spacing, for different CN values and for the number of simulation runs. The loop on Early-minus-Late spacing can be absent if, for example, only a PLL is considered. For each Early-minus-Late spacing and for a fixed CN, a noise vector containing the noise components of the correlator outputs is generated. The vector length is equal to the number of simulation runs and all the noise components are generated at once for efficiency reasons.

All intermediate results, such as the discriminator and loop filter outputs, are stored in auxiliary vectors and are used at the end of the loop on the simulation runs to evaluate quantities of interest such as the tracking jitter.

In the code provided, theoretical formulae for the computation of the tracking jitter are also implemented and used as a comparison term for the simulation results.

3.2. Standard PLL (PLL.m)

The simulation of a standard PLL requires the generation of the Prompt correlator alone (GenerateSignalCorrelation). The Prompt correlator is the output of the ID block computed with respect to the best delay estimate provided by the loop [21]. For this reason, the noise generation (GenerateNoiseVector) simply consists of simulating a one dimensional complex Gaussian white sequence with independent and identically distributed real and imaginary parts with variance [9]

$${\sigma }_n^2\mathrm{<mo>=</mo>}\frac{1}{C\mathrm{<mo>/</mo>}{N}_0{T}_c}\mathrm{.}$$E29

When simulating a standard PLL alone, perfect code synchronization is assumed and (26) simplifies to

$$\sqrt{C}\frac{\sin \left(\pi \Delta f_d{T}_c\right)}{\pi \Delta f_d{T}_c}\exp \left\{j\Delta \phi \right\}+{\eta }_c\mathrm{,}$$E30

where $$\Delta f_d$$is obtained by comparing the true Doppler frequency against the loop filter output. $$\Delta \phi$$is the phase error obtained as the difference between the true phase and the phase estimate produced by the NCO.

In SATLSim, the function UpdateDiscriminator implements a standard Costas discriminator. Different phase discriminators, as indicated in [21], can be easily implemented by changing this function.

3.3. Double estimator (DoubleEstimator.m)

In the DE case, i.e. when DLL and SLL are jointly used, the function GenerateNoiseVector, responsible for the generation of the correlator noise, produces a $$5\times N_{sim}$$matrix, where $$N_{sim}$$is the number of simulation runs. The five rows of this matrix correspond to the five correlators required by the DE that are characterized by the following correlation matrix

$$\mathrm{,}{C}_n=\left[\begin{array}{ccccc}1& R_l\left(\frac{d_s}{2},\frac{d_{sc}}{2}\right)& R_l\left(\frac{d_s}{2},0\right)& R_l\left(\frac{d_s}{2},\frac{d_{sc}}{2}\right)& R_l\left(d_s,0\right)\\ R_l\left(\frac{d_s}{2},\frac{d_{sc}}{2}\right)& 1& R_l\left(0,\frac{d_{sc}}{2}\right)& R_l\left(0,d_{sc}\right)& R_l\left(\frac{d_s}{2},\frac{d_{sc}}{2}\right)\\ R_l\left(\frac{d_s}{2},0\right)& R_l\left(0,\frac{d_{sc}}{2}\right)& 1& R_l\left(0,\frac{d_{sc}}{2}\right)& R_l\left(\frac{d_s}{2},0\right)\\ R_l\left(\frac{d_s}{2},\frac{d_{sc}}{2}\right)& R_l\left(0,d_{sc}\right)& R_l\left(0,\frac{d_{sc}}{2}\right)& 1& R_l\left(\frac{d_s}{2},\frac{d_{sc}}{2}\right)\\ R_l\left(d_s,0\right)& R_l\left(\frac{d_s}{2},\frac{d_{sc}}{2}\right)& R_l\left(\frac{d_s}{2},0\right)& R_l\left(\frac{d_s}{2},\frac{d_{sc}}{2}\right)& 1\end{array}\right],$$E31

where $$d_{sc}$$is the subcarrier Early-minus-Late spacing.

The NCO update (UpdateNCO) is performed on both code and subcarrier loops and the estimated errors, $$\Delta {\tau }_d$$and$$\Delta {\tau }_s$$, are used to compute new correlator signal components (GenerateSignalCorrelation). Two nonlinear discriminators (UpdateDiscriminator) and loop filters (UpdateFilter) are run in parallel to determine the rate of change of both code and subcarrier delay.

The DE provides an example of how several tracking loops, operating in parallel, can be easily coupled in order to provide more realistic simulations accounting for the interaction of different tracking algorithms [9].

3.4. Sample results

In this section, sample results obtained using the SATLSim toolbox are shown for the DE case. Results for the analysis of the PLL can be found in [10]. Specific focus is devoted to the analysis of the tracking jitter, which is one of the main metrics used for the analysis of digital tracking loops. The tracking jitter quantifies the amount of noise transferred by the tracking loop to the final parameter estimate [29]. The tracking jitter is the standard deviation of the final parameter estimate normalized by the discriminator gain. The non-linear discriminator is usually a memoryless device characterized by an input/output function relating the parameter estimation error to the discriminator output. The discriminator gain is the slope of this function in the neighborhood of zero (hypothesis of small estimation error).

Tracking jitter results obtained using non-coherent discriminators [21] for both DLL and SLL are shown in Figure 11. The figure is divided into three parts: [a)]

1. Tracking jitter of the DLL alone

2. Tracking jitter of the SLL alone

3. Jitter of the combined delay estimate.

This is due to the fact that the DE jointly uses a DLL, for estimating the code delay, and a SLL, for determining the subcarrier delay. Subcarrier and code delay are then combined to obtain the final estimate of the travel time of the transmitted signal [27]. Thus, three different jitters are evaluated for the different estimates produced by the system. Tracking jitter has been expressed in meters by multiplying the standard deviation of the delay estimates by the speed of light.

Figure 11.

Tracking jitter obtained using the SATLSim toolbox. a) Tracking jitter of the DLL alone. b) Tracking jitter of the SLL alone. c) Jitter of the combined delay estimate.

In addition to this, the curves are shown in Figure 11 a) and Figure 11 b). More specifically, three different methodologies have been employed for determining the tracking jitter. The theoretical curve corresponds to approximate formulas obtained by linearizing the input/output function of the non-linear discriminator. These formulas are valid only for small tracking errors or equivalently for high CN. The jitter obtained from the actual error has been obtained by evaluating the variance of the code phase error. It is noted that in a real tracking loop the code phase error is not directly accessible since the true code phase is unknown. Thus, the tracking error can be evaluated by measuring the error at the loop filter output, which is an observable point, and propagating its variance through the loop. The tracking jitter obtained by propagating the variance at this measurable point corresponds to the curve denoted by “Estimated from the loop filter output”. The relationship between the variances of the discriminator output and the true tracking error can easily be evaluated when the loop is working in its linear region. The measured curve was introduced to further validate the theoretical model and test the correctness of the simulation methodology. This latest curve is not available for the combined delay estimate.

The parameters used for the evaluation of the tracking jitter, shown in Figure 11, are provided in Table 3.

Parameter	Value
Sampling Frequency	$$f_s\mathrm{<mo>=</mo>8}$$MHz
Integration Time	$$T_c\mathrm{<mo>=</mo>4}$$ms
DLL Early-minus-Late spacing	$$0.1955$$$$\mu$$s (0.2 chips)
SLL Early-minus-Late spacing	$$0.1955$$$$\mu$$s (0.2 chips)
DLL Loop Order	$$1$$
SLL Loop Order	$$1$$
DLL Loop Bandwidth	$$0.5$$Hz
SLL Loop Bandwidth	$$0.5$$Hz
Modulation type	BOC(1, 1)

Table 3.

Parameters used for the evaluation of the tracking jitter shown in Figure 11.

From the results shown in Figure 11, it is observed that the developed semi-analytic technique is able to effectively capture the behavior of the system. For high CN values, a good agreement between theoretical and simulation results is found. However, for CN lower than $$22$$dB-Hz theoretical and simulation results start diverging. This is more clear in parts a) and c) of the figures. For such low CN values, the loop is no longer working in the linear region of the discriminator input/output function. Thus, the theoretical model is unable to capture the behavior of the loop that is losing lock. The semi-analytic technique implemented in the SATLSim MT toolbox is able to effectively describe the non-linear behavior of the loop requiring only limited computation resources.

4. Conclusions

In this chapter, the development of fast semi-analytic techniques using MT has been analyzed. In the semi-analytic framework, the knowledge of the system under analysis is exploited to reduce the computational load and complexity that full Monte Carlo simulations would require. In this way, the strengths of both analytical and Monte Carlo methods are effectively combined.

Two examples of semi-analytic techniques have been thoroughly analyzed and used to illustrate the two main configurations developed within the semi-analytic framework. The first example illustrates the sequential configuration where simulations and the analytical model are used sequentially. This type of configuration provides increased precision with respect to full Monte Carlo simulations, particularly when the quantities to be estimated assume small values. In addition to this, a significant reduction in terms of computational complexity is achieved. In the example considered, full Monte Carlo simulations require the implementation of a full transmission/reception chain including the interaction between two different systems, DVB-T and GNSS. This requirement led to a significant computational and development complexity. The considered semi-analytic approach is an effective solution for alleviating those requirements.

The second example considered the closed-loop approach and specific focus was devoted to the SATLSim MT toolbox. This toolbox has been developed for the analysis of digital tracking loops and fully exploits the flexibility of the MT programming language. The code has been organized in functions that can be easily replaced by different MT modules. In this way, different loop components such as discriminators, loop filters and NCO models can be integrated in the SATLSim toolbox. The efficiency of semi-analytic techniques and the reduced development time enabled by the MT language are an effective tool for the analysis of complex communications systems.