**Simulation Framework of Wireless Sensor Network (WSN) Using MATLAB/SIMULINK Software**

Qutaiba I. Ali

262 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 2

2007, vol., no., pp.342-347, 25-28 Sept. 2007. ISBN: 978-0-7695-2974-5

Applications, Second Edition, Wiley, ISBN 0-471-33341-7.

and Applications. Springer. USA. ISBN 978-0-387-21948-6.

no.8, pp.1209-1218, ISSN : 0733-8716.

Inc. ISBN 0-07-046213-5.

1-4244-0504-1.

0-12-598062-0.

Nelson, Barry. (2002). Stochastic Modeling: Analysis & Simulation. USA. Dover Publications

Parsons, J.D. (1996). The Mobile Radio Propagation Channel, USA, Wiley, ISBN 0471964158 Rezaeian, M. (July 2006), Symmetric Characterization of Finite State Markov Channels, Information Theory, 2006 IEEE International Symposium on, pp.2734-2738, 9-14, E-ISBN

Ross, S. (2007) Introduction to Probability Models, 9th edition, USA, Academic Press, ISBN

Sanchez-Salas, D.A.; Cuevas-Ruiz, J.L.; (2007), N-states Channel Model using Markov Chains, Electronics, Robotics and Automotive Mechanics Conference, 2007. CERMA

Trivedi, K. S. (2002). *Probability, Statistics with Reliability, Queueing and Computer Science*

Vucetic, B.; Du, J.; (Oct 1992), Channel modeling and simulation in satellite mobile communication systems, Selected Areas in Communications, IEEE Journal on, vol.10,

Yin, George. Zhang, Qing. (2005). Discrete-Time Markov Chains. Two Time-Scale Methods

Additional information is available at the end of the chapter

http://dx.doi.org/10.5772/46467

## **1. Introduction**

A wireless sensor network consists of spatially distributed autonomous sensors to cooperatively monitor physical or environmental conditions, such as temperature, sound, vibration, pressure, motion or pollutants. The development of wireless sensor networks was motivated by military applications such as battlefield surveillance. They are now used in many industrial and civilian application areas, including industrial process monitoring and control, machine health monitoring, environment and habitat monitoring, healthcare applications, home automation, and traffic control [1-2].

A smart sensor node is a combination of sensing, processing and communication technologies. Figure 1 shows the basic architectural components of a sensor node. The sensing unit senses the change of parameters, signal conditioning circuitry prepares the electrical signals to convert to the digital domain, the sensed analog signal is converted and is used as the input to the application algorithms or processing unit, the memory helps processing of tasks and the transceiver is used for communicating with other sensors or the base stations or sinks in WSN[3], see figure 1.

Sensors can monitor temperature, pressure, humidity, soil makeup, vehicular movement, noise levels, lighting conditions, the presence or absence of certain kinds of objects or substances, mechanical stress levels on attached objects, and other properties. Their mechanism may be seismic, magnetic, thermal, visual, infrared, acoustic, or radar. A smart sensor is also capable of self-identification and self-diagnosis. The mechanisms of smart sensors work in one of three ways: by a line of sight to the target (such as visual sensors), by proximity to target (such as seismic sensors), and by propagation like a wave with possible bending (such as acoustic sensors)[4,5].

© 2012 Ali, licensee InTech. This is an open access chapter distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. © 2012 Ali, licensee InTech. This is a paper distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Simulation Framework of Wireless Sensor Network (WSN) Using MATLAB/SIMULINK Software 265

and energy), updates to current models as well as performance improvements. Furthermore, a new QT based GUI was added providing a scenario designer, a visualizer to view network scenarios (2D and 3D), a packet tracer for debugging, an analyzer for statistics and a file

OPNET Modeler Wireless Suite [8]–[10] is a commercial modeling and simulation tool for various types of wireless networks. It is developed by developed by OPNET Technologies, Inc. and based on the well-known product OPNET Modeler. The simulation environment uses a fast discrete event simulation engine operating with a 32-bit/ 64-bit fully parallel simulation kernel, which is available for Windows and Linux. The OPNET Modeler provides an object-oriented modeling approach and a hierarchical modeling environment. Although there are no special routing protocols for wireless sensor network available, at least different propagation and modulation techniques as well as a ZigBee (802.15.4) MAC layer are provided. Additional modules have to be customized or developed from the scratch. The simulations of wireless networks can be run as discrete event, hybrid or analytical, encompassing terrain, mobility and path-loss models. Due to the open interface external object files, libraries as well as other simulators can be integrated to the OPNET Modeler. Optional a System-in-the-Loop is available to interface simulations with live systems. Furthermore, the OPNET Modeler Wireless Suite provides grid computing support so that simulations can be executed in a distributed

TOSSIM (TinyOS mote simulator) [12]–[15] is a discrete event simulator for TinyOS sensor networks that is part of the official TinyOS package. TOSSIM takes advantage of the component based architecture of TinyOS by integrating it transparently by providing a new hardware resource abstraction layer that simulates the TinyOS network stack at the bit level for normal PCs. Due to this approach low-level protocols up to top-level applications can be simulated with TOSSIM. TOSSIM has an external communication system so that transmitted packets can be monitored and even new packets can be injected to network. Furthermore, the configuration of the debug options is fine grained providing the desired debug output at runtime. TOSSIM offers three network connectivity models: simple connectivity, static connectivity and space connectivity. The running simulations can be visualized and

OMNeT++ [17]–[20] is an object-oriented discrete network simulation framework. The architecture is rather generic so that various problem domains can be simulated such as protocol modeling, validation of hardware architectures and modeling of wired and wireless communication networks. OMNeT++ is not a simulator, but it rather provides a framework and tools to write simulations. It is highly portable so that it can be run on the most common operating systems such as Windows, Linux and Mac OSX. There are a couple

editor to edit the scenarios directly.

manner[11]. **c. TOSSIM** 

**d. OMNeT++** 

**b. OPNET Modeler Wireless Suite** 

controlled by the Java-based GUI TinyViz[16].

**Figure 1.** Basic architectural components of a smart sensor

## **2. Review of existing simulation environments for WSNs**

In this section, a selection of existing simulation environments for WSNs is discussed. Basically, the investigated simulation environments can be divided into two major types: adaptive development and new development. The adaptive development covers simulation environments that already existed before the idea of WSNs emerged. These simulation environments were then extended to support wireless functionality and were then adapted for the use with WSNs. In contrast, new developments cover new simulators, which were created solely for simulating WSNs, considering sensor specific characteristics from the beginning. Both types have advantages and disadvantages, but basically it can be stated that while the evolutionary adaptation has some advantages in reusing well-tested ideas and source code as well as the bigger user and developer basis, the new developments have their advantages in focusing on the special characteristics and the functioning of sensor nodes.

### **a. GloMoSim/QualNet**

GloMoSim [6] is a scalable simulation environment for wireless and wired network systems, which uses the parallel discrete-event simulation capability provided by Parsec [7], a cbased simulation language for sequential and parallel execution of discrete-event simulation models. Both, GloMoSim as well as Parsec, were developed by the Parallel Computing Lab. at UCLA. GloMoSim offers basic functionality to simulate wireless networks, even for ad hoc networks (e.g. AODV, DSR). However, the current version of GloMoSim does not offer any sensor network specific features in the default package so that without any further efforts no WSNs can be simulated meaningfully. In 2000 QualNet [6], [7], a commercial derivate of GloMoSim, was created and with GloMoSim 2.0, the last version of GloMoSim, was released under an academic license. From this point in time, no further improvements to GloMoSim were made, whereas the development of QualNet expedited. In October 2009, version 5.0 of QualNet was released including enhancements such as a new sensor network library for ZigBee, new network security library, parallel updates, new models (e.g. battery and energy), updates to current models as well as performance improvements. Furthermore, a new QT based GUI was added providing a scenario designer, a visualizer to view network scenarios (2D and 3D), a packet tracer for debugging, an analyzer for statistics and a file editor to edit the scenarios directly.

## **b. OPNET Modeler Wireless Suite**

OPNET Modeler Wireless Suite [8]–[10] is a commercial modeling and simulation tool for various types of wireless networks. It is developed by developed by OPNET Technologies, Inc. and based on the well-known product OPNET Modeler. The simulation environment uses a fast discrete event simulation engine operating with a 32-bit/ 64-bit fully parallel simulation kernel, which is available for Windows and Linux. The OPNET Modeler provides an object-oriented modeling approach and a hierarchical modeling environment. Although there are no special routing protocols for wireless sensor network available, at least different propagation and modulation techniques as well as a ZigBee (802.15.4) MAC layer are provided. Additional modules have to be customized or developed from the scratch. The simulations of wireless networks can be run as discrete event, hybrid or analytical, encompassing terrain, mobility and path-loss models. Due to the open interface external object files, libraries as well as other simulators can be integrated to the OPNET Modeler. Optional a System-in-the-Loop is available to interface simulations with live systems. Furthermore, the OPNET Modeler Wireless Suite provides grid computing support so that simulations can be executed in a distributed manner[11].

### **c. TOSSIM**

264 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 2

**Figure 1.** Basic architectural components of a smart sensor

nodes.

**a. GloMoSim/QualNet** 

**2. Review of existing simulation environments for WSNs** 

In this section, a selection of existing simulation environments for WSNs is discussed. Basically, the investigated simulation environments can be divided into two major types: adaptive development and new development. The adaptive development covers simulation environments that already existed before the idea of WSNs emerged. These simulation environments were then extended to support wireless functionality and were then adapted for the use with WSNs. In contrast, new developments cover new simulators, which were created solely for simulating WSNs, considering sensor specific characteristics from the beginning. Both types have advantages and disadvantages, but basically it can be stated that while the evolutionary adaptation has some advantages in reusing well-tested ideas and source code as well as the bigger user and developer basis, the new developments have their advantages in focusing on the special characteristics and the functioning of sensor

GloMoSim [6] is a scalable simulation environment for wireless and wired network systems, which uses the parallel discrete-event simulation capability provided by Parsec [7], a cbased simulation language for sequential and parallel execution of discrete-event simulation models. Both, GloMoSim as well as Parsec, were developed by the Parallel Computing Lab. at UCLA. GloMoSim offers basic functionality to simulate wireless networks, even for ad hoc networks (e.g. AODV, DSR). However, the current version of GloMoSim does not offer any sensor network specific features in the default package so that without any further efforts no WSNs can be simulated meaningfully. In 2000 QualNet [6], [7], a commercial derivate of GloMoSim, was created and with GloMoSim 2.0, the last version of GloMoSim, was released under an academic license. From this point in time, no further improvements to GloMoSim were made, whereas the development of QualNet expedited. In October 2009, version 5.0 of QualNet was released including enhancements such as a new sensor network library for ZigBee, new network security library, parallel updates, new models (e.g. battery TOSSIM (TinyOS mote simulator) [12]–[15] is a discrete event simulator for TinyOS sensor networks that is part of the official TinyOS package. TOSSIM takes advantage of the component based architecture of TinyOS by integrating it transparently by providing a new hardware resource abstraction layer that simulates the TinyOS network stack at the bit level for normal PCs. Due to this approach low-level protocols up to top-level applications can be simulated with TOSSIM. TOSSIM has an external communication system so that transmitted packets can be monitored and even new packets can be injected to network. Furthermore, the configuration of the debug options is fine grained providing the desired debug output at runtime. TOSSIM offers three network connectivity models: simple connectivity, static connectivity and space connectivity. The running simulations can be visualized and controlled by the Java-based GUI TinyViz[16].

### **d. OMNeT++**

OMNeT++ [17]–[20] is an object-oriented discrete network simulation framework. The architecture is rather generic so that various problem domains can be simulated such as protocol modeling, validation of hardware architectures and modeling of wired and wireless communication networks. OMNeT++ is not a simulator, but it rather provides a framework and tools to write simulations. It is highly portable so that it can be run on the most common operating systems such as Windows, Linux and Mac OSX. There are a couple

of simulation frameworks that enable OMNeT++ to be used for wireless sensor networks[21]. The most common of these frameworks are discussed in the following subsections.

Simulation Framework of Wireless Sensor Network (WSN) Using MATLAB/SIMULINK Software 267

Avrora [41]–[43] is a set of simulation and analysis tools for programs written for AVR micro-controllers. It has support for different sensor platforms, such as Mica2 and MicaZ, allowing wireless network simulation, dynamic instrumentation and static analysis. Since 2004, Avrora is developed in a research project of the UCLA compiler group. The special characteristic of Avrora is that it operates on the instruction-level, i.e. actual microcontroller programs can be run in the simulator, instead of just simulating software

J-Sim [44]–[46] is a component-based compositional simulation environment based on the autonomous component architecture (ACA). The basic entities of ACA are components, which communicated with each other by sending and receiving data using their ports. Application specific models can be defined by sub-classing the specified classes of the WSN simulation framework and adapting them to the desired behavior. At the moment, 802.11 is

ATEMU [47], [48] is one of the first instruction-level software emulators for AVR based systems. Additionally peripheral devices of the MICA2 sensor node platform such as radio is supported. Although at the moment only the MICA2 hardware is supported, ATEMU can be easily extended to support other sensor node platforms. Although ATEMU is the most accurate instruction-level emulator for wireless sensor network research, it lacks from

EmStar [49]–[51] is an environment for WSNs built from Linux-class devices, so called micro servers. In comparison to motes, micro servers are much less constrained in computational power and data storage size so that they can handle more complex tasks such as image and audio processing. EmStar consists of simulation and emulation tools, which utilize a modular, but not strictly layered, architecture. EmStar provides different simulation modes: a pure simulation mode, an emulation mode, a real mode and a hybrid mode. EmStar provides various services that are used and combined to provide network functionality for wireless embedded systems. This includes link drivers for the lowest-layer interfaces to network resources, pass-through modules that implement various types of filter and passive processing, and routing modules such as Flooding, DSR, Sink, StateSync

SENS [52], [53] is an application-oriented simulator for WSNs. It has a modular, layered architecture so that components for applications, network communication and the physical environment can be easily interchanged and extended. Due to different component

used as MAC Layer and AODV is provided as routing protocol.

simulation speed, being 30 times slower than TOSSIM, for example.

**f. Avrora** 

models.

**g. J-Sim** 

**h. ATEMU** 

**i. EmStar** 

and Centroute.

**j. SENS** 

Mobility Framework: The Mobility Framework [22]–[24], developed in the Telecommunication Networks Group (TKN) at the Technical University of Berlin, provides only basic support for mobile and wireless networks. It includes some basic layers such as MAC layers (Aloha, CSMA) and network layers (flooding) as well as some basic mobility functionality and some basic application layer.

MiXiM: MiXiM [25], [26] is a merger of several OMNeT++ frameworks to support mobile and wireless simulations. It uses the mobility support, the connection management, and the general structure from the Mobility Framework (MF); the radio propagation models from the CHannel SIMulator (ChSim); and the protocol library from the MAC simulator, the Positif frame- work [27], and the Mobility Framework.

Castalia: Castalia [28], [29] is a simulator for WSNs (WSN), body area networks (BAN) and generally networks of low-power embedded devices that is based on OMNeT++. It is developed at the National ICT Australia since 2006 and made public as open source under the Academic Public License in 2007.

INET Framework: The INET Framework [30], [31] is a framework for OMNeT++ that contains various implementations of common protocols, such as IPv4, IPv6, TCP, UDP etc., as well as several application models. The INET Framework is not specialized on mobile and wireless networks, but has some support for it.

NesCT: NesCT [32] is not a real framework, but rather a translator from the programming language NesC to C++ classes for OMNet++.

### **e. NS-2**

NS (the Network Simulator) [33], [34] is an object-oriented discrete event simulator targeting at networking research. NS-2 is written in C++ and OTcl, an object-oriented version of Tcl. A huge amount of contributed protocol source codes can be found on the website http://nsnam.isi.edu/nsnam/index.php/Contributed Code. among them there are also some for WSNs interesting wireless protocols such as different variations of 802.11, 802.16, IR-UWB, BlueTooth and 802.15.4. Despite the great number of contributing researchers the support for wireless sensor network specific protocols is rather low. As special wireless sensor network framework the Mannasim Framework [36] should be highlighted that provides sensor network specific protocols such as LEACH and Directed Diffusion. Also the extension NS2-MIUN [37] provides some wireless sensor network specific contributions with the focus on intrusion detection. SensorSim: SensorSim [38]–[40] is a simulation framework for modeling sensor networks that built up on NS-2. It provides additional features for modeling sensor networks such as sensor channel models, power models (battery and radio), lightweight protocol stacks for wireless micro-sensors, scenario generation and hybrid simulation.

## **f. Avrora**

266 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 2

subsections.

**e. NS-2** 

functionality and some basic application layer.

the Academic Public License in 2007.

Positif frame- work [27], and the Mobility Framework.

and wireless networks, but has some support for it.

language NesC to C++ classes for OMNet++.

generation and hybrid simulation.

of simulation frameworks that enable OMNeT++ to be used for wireless sensor networks[21]. The most common of these frameworks are discussed in the following

Mobility Framework: The Mobility Framework [22]–[24], developed in the Telecommunication Networks Group (TKN) at the Technical University of Berlin, provides only basic support for mobile and wireless networks. It includes some basic layers such as MAC layers (Aloha, CSMA) and network layers (flooding) as well as some basic mobility

MiXiM: MiXiM [25], [26] is a merger of several OMNeT++ frameworks to support mobile and wireless simulations. It uses the mobility support, the connection management, and the general structure from the Mobility Framework (MF); the radio propagation models from the CHannel SIMulator (ChSim); and the protocol library from the MAC simulator, the

Castalia: Castalia [28], [29] is a simulator for WSNs (WSN), body area networks (BAN) and generally networks of low-power embedded devices that is based on OMNeT++. It is developed at the National ICT Australia since 2006 and made public as open source under

INET Framework: The INET Framework [30], [31] is a framework for OMNeT++ that contains various implementations of common protocols, such as IPv4, IPv6, TCP, UDP etc., as well as several application models. The INET Framework is not specialized on mobile

NesCT: NesCT [32] is not a real framework, but rather a translator from the programming

NS (the Network Simulator) [33], [34] is an object-oriented discrete event simulator targeting at networking research. NS-2 is written in C++ and OTcl, an object-oriented version of Tcl. A huge amount of contributed protocol source codes can be found on the website http://nsnam.isi.edu/nsnam/index.php/Contributed Code. among them there are also some for WSNs interesting wireless protocols such as different variations of 802.11, 802.16, IR-UWB, BlueTooth and 802.15.4. Despite the great number of contributing researchers the support for wireless sensor network specific protocols is rather low. As special wireless sensor network framework the Mannasim Framework [36] should be highlighted that provides sensor network specific protocols such as LEACH and Directed Diffusion. Also the extension NS2-MIUN [37] provides some wireless sensor network specific contributions with the focus on intrusion detection. SensorSim: SensorSim [38]–[40] is a simulation framework for modeling sensor networks that built up on NS-2. It provides additional features for modeling sensor networks such as sensor channel models, power models (battery and radio), lightweight protocol stacks for wireless micro-sensors, scenario Avrora [41]–[43] is a set of simulation and analysis tools for programs written for AVR micro-controllers. It has support for different sensor platforms, such as Mica2 and MicaZ, allowing wireless network simulation, dynamic instrumentation and static analysis. Since 2004, Avrora is developed in a research project of the UCLA compiler group. The special characteristic of Avrora is that it operates on the instruction-level, i.e. actual microcontroller programs can be run in the simulator, instead of just simulating software models.

## **g. J-Sim**

J-Sim [44]–[46] is a component-based compositional simulation environment based on the autonomous component architecture (ACA). The basic entities of ACA are components, which communicated with each other by sending and receiving data using their ports. Application specific models can be defined by sub-classing the specified classes of the WSN simulation framework and adapting them to the desired behavior. At the moment, 802.11 is used as MAC Layer and AODV is provided as routing protocol.

## **h. ATEMU**

ATEMU [47], [48] is one of the first instruction-level software emulators for AVR based systems. Additionally peripheral devices of the MICA2 sensor node platform such as radio is supported. Although at the moment only the MICA2 hardware is supported, ATEMU can be easily extended to support other sensor node platforms. Although ATEMU is the most accurate instruction-level emulator for wireless sensor network research, it lacks from simulation speed, being 30 times slower than TOSSIM, for example.

### **i. EmStar**

EmStar [49]–[51] is an environment for WSNs built from Linux-class devices, so called micro servers. In comparison to motes, micro servers are much less constrained in computational power and data storage size so that they can handle more complex tasks such as image and audio processing. EmStar consists of simulation and emulation tools, which utilize a modular, but not strictly layered, architecture. EmStar provides different simulation modes: a pure simulation mode, an emulation mode, a real mode and a hybrid mode. EmStar provides various services that are used and combined to provide network functionality for wireless embedded systems. This includes link drivers for the lowest-layer interfaces to network resources, pass-through modules that implement various types of filter and passive processing, and routing modules such as Flooding, DSR, Sink, StateSync and Centroute.

### **j. SENS**

SENS [52], [53] is an application-oriented simulator for WSNs. It has a modular, layered architecture so that components for applications, network communication and the physical environment can be easily interchanged and extended. Due to different component

implementations, which varies in the degree of realism, application-specific environments can be created and simulated. Due to the chosen approach, SENS enables application portability because the same source code can be run with in a simulation or deployed on actual sensor nodes.

Simulation Framework of Wireless Sensor Network (WSN) Using MATLAB/SIMULINK Software 269

including matrix algebra, complex arithmetic, linear systems, differential equations, signal processing, optimization, nonlinear systems, and many other types of scientific computations. The most important feature of *MATLAB* is its programming capability, which is very easy to learn and to use, and which allows user-developed functions. It also allows access to Fortran algorithms and C codes by means of external interfaces. There are several optional toolboxes written for special applications such as signal processing, control systems design, system identification, statistics, neural networks, fuzzy logic, symbolic computations, and others. *MATLAB* has been enhanced by the very powerful Simulink

Simulink is a software package for modeling, simulating, and analyzing dynamical systems. It supports linear and nonlinear systems, modeled in continuous time, sampled time, or a hybrid of the two. Systems can also be multi-rate, i.e., have different parts that are sampled or updated at different rates. For modeling, Simulink provides a graphical user interface (GUI) for building models as block diagrams, using click-and-drag mouse operations. With this interface, you can draw the models just as you would with pencil and paper (or as most textbooks depict them). Simulink includes a comprehensive block library of sinks, sources, linear and nonlinear components, and connectors. You can also customize and create your own blocks Models are hierarchical. This approach provides insight into how a model is organized and how its parts interact. After you define a model, you can simulate it, using a choice of integration methods, either from the Simulink menus or by entering commands in MATLAB's command window. The menus are particularly convenient for interactive work, while the command-line approach is very useful for running a batch of simulations (for example, if you are doing Monte Carlo simulations or want to sweep a parameter across a range of values). Using scopes and other display blocks, you can see the simulation results while the simulation is running. In addition, you can change parameters and immediately see what happens, for "what if" exploration. The simulation results can be put in the MATLAB workspace for post processing and visualization. And because MATLAB and Simulink are integrated, you can simulate, analyze, and revise your models in either environment at any

In order to demonstrate the concepts of the suggested simulation methodology, a simple WSN model was built as shown in figure 2[60]. This network consisted of three sensors (slaves) sending their measured data samples to a master node. In this chapter, MATLAB Simulink communication block set was used to build a complete WSN system. Simulation procedure includes building the hardware architecture of the transmitting nodes, modeling both the communication channel and the receiving master node architecture. Bluetooth was chosen to undertake the physical layer communication with respect to different channel parameters (i.e., Signal to Noise ratio, Attenuation and Interference). The simulation model was examined using different topologies under various conditions and numerous results

program[59].

point. [59]**.** 

were collected.

**3.1. Simulating a simple WSN in Simulink MATLAB** 

### **k. SENSE**

SENSE (Sensor Network Simulator and Emulator) [54], [55] is a simulator for WSNs that is based on a novel component-oriented simulation methodology, which promotes extensibility and reusability. At the same time, the simulation efficiency and the scalability was considered. In the component repository of SENSE there are already different components available from the application to the physical layer including IEEE 802.11, AODV, DSR, SSR, SHR as well as Battery Models and a Power Model. At the moment, there does not seem to be any further tools included in SENSE so that, for example, a visualization tool to analyze the network behavior graphically is missing.

#### **l. Shawn**

Shawn [56]–[58] is customizable sensor network simulator based on an algorithmic approach. The primary design goals of Shawn are: simulate the effect caused by a phenomenon, scalability and support for extremely large networks and free choice of the implementation model. Models are the interfaces that are used by Shawn to control the simulation without knowing the exact implementation. Each implementation of a model specifies the actual behavior. Currently there are several implementations for the transmission model provided such as Pure CSMA & CSMA/CA, (Slotted) Aloha, Random Drop etc. Additionally to the three core models, there are several other models provided by Shawn for random variables, distance estimation and mobility. The simulation environment provides a sort of virtual world in which the different parts of the simulation are located. The simulated nodes are located in a single world instance and the nodes themselves are containers for processors. The application logic is implemented as instances of processors.

In this chapter, Simulink MATLAB was adopted to be the simulation tool of wireless sensor network (WSN). The main advantage of the suggested method is to determine the effect of the different channel parameters (i.e., Signal to Noise ratio, Attenuation and Interference) on the system behavior.

## **3. The proposed WSN simulation methodology**

The environment in which we build our simulation model was MATLAB. The name MATLAB stands for matrix laboratory. *MATLAB*, developed by MathWorks Inc., is a software package for high performance numerical computation and visualization. The combination of analysis capabilities, flexibility, reliability, and powerful graphics makes *MATLAB* the premier software package for scientific researchers. *MATLAB* provides an interactive environment with hundreds of reliable and accurate built-in mathematical functions. These functions provide solutions to a broad range of mathematical problems including matrix algebra, complex arithmetic, linear systems, differential equations, signal processing, optimization, nonlinear systems, and many other types of scientific computations. The most important feature of *MATLAB* is its programming capability, which is very easy to learn and to use, and which allows user-developed functions. It also allows access to Fortran algorithms and C codes by means of external interfaces. There are several optional toolboxes written for special applications such as signal processing, control systems design, system identification, statistics, neural networks, fuzzy logic, symbolic computations, and others. *MATLAB* has been enhanced by the very powerful Simulink program[59].

268 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 2

tool to analyze the network behavior graphically is missing.

**3. The proposed WSN simulation methodology** 

actual sensor nodes.

**k. SENSE** 

**l. Shawn** 

of processors.

the system behavior.

implementations, which varies in the degree of realism, application-specific environments can be created and simulated. Due to the chosen approach, SENS enables application portability because the same source code can be run with in a simulation or deployed on

SENSE (Sensor Network Simulator and Emulator) [54], [55] is a simulator for WSNs that is based on a novel component-oriented simulation methodology, which promotes extensibility and reusability. At the same time, the simulation efficiency and the scalability was considered. In the component repository of SENSE there are already different components available from the application to the physical layer including IEEE 802.11, AODV, DSR, SSR, SHR as well as Battery Models and a Power Model. At the moment, there does not seem to be any further tools included in SENSE so that, for example, a visualization

Shawn [56]–[58] is customizable sensor network simulator based on an algorithmic approach. The primary design goals of Shawn are: simulate the effect caused by a phenomenon, scalability and support for extremely large networks and free choice of the implementation model. Models are the interfaces that are used by Shawn to control the simulation without knowing the exact implementation. Each implementation of a model specifies the actual behavior. Currently there are several implementations for the transmission model provided such as Pure CSMA & CSMA/CA, (Slotted) Aloha, Random Drop etc. Additionally to the three core models, there are several other models provided by Shawn for random variables, distance estimation and mobility. The simulation environment provides a sort of virtual world in which the different parts of the simulation are located. The simulated nodes are located in a single world instance and the nodes themselves are containers for processors. The application logic is implemented as instances

In this chapter, Simulink MATLAB was adopted to be the simulation tool of wireless sensor network (WSN). The main advantage of the suggested method is to determine the effect of the different channel parameters (i.e., Signal to Noise ratio, Attenuation and Interference) on

The environment in which we build our simulation model was MATLAB. The name MATLAB stands for matrix laboratory. *MATLAB*, developed by MathWorks Inc., is a software package for high performance numerical computation and visualization. The combination of analysis capabilities, flexibility, reliability, and powerful graphics makes *MATLAB* the premier software package for scientific researchers. *MATLAB* provides an interactive environment with hundreds of reliable and accurate built-in mathematical functions. These functions provide solutions to a broad range of mathematical problems Simulink is a software package for modeling, simulating, and analyzing dynamical systems. It supports linear and nonlinear systems, modeled in continuous time, sampled time, or a hybrid of the two. Systems can also be multi-rate, i.e., have different parts that are sampled or updated at different rates. For modeling, Simulink provides a graphical user interface (GUI) for building models as block diagrams, using click-and-drag mouse operations. With this interface, you can draw the models just as you would with pencil and paper (or as most textbooks depict them). Simulink includes a comprehensive block library of sinks, sources, linear and nonlinear components, and connectors. You can also customize and create your own blocks Models are hierarchical. This approach provides insight into how a model is organized and how its parts interact. After you define a model, you can simulate it, using a choice of integration methods, either from the Simulink menus or by entering commands in MATLAB's command window. The menus are particularly convenient for interactive work, while the command-line approach is very useful for running a batch of simulations (for example, if you are doing Monte Carlo simulations or want to sweep a parameter across a range of values). Using scopes and other display blocks, you can see the simulation results while the simulation is running. In addition, you can change parameters and immediately see what happens, for "what if" exploration. The simulation results can be put in the MATLAB workspace for post processing and visualization. And because MATLAB and Simulink are integrated, you can simulate, analyze, and revise your models in either environment at any point. [59]**.** 

### **3.1. Simulating a simple WSN in Simulink MATLAB**

In order to demonstrate the concepts of the suggested simulation methodology, a simple WSN model was built as shown in figure 2[60]. This network consisted of three sensors (slaves) sending their measured data samples to a master node. In this chapter, MATLAB Simulink communication block set was used to build a complete WSN system. Simulation procedure includes building the hardware architecture of the transmitting nodes, modeling both the communication channel and the receiving master node architecture. Bluetooth was chosen to undertake the physical layer communication with respect to different channel parameters (i.e., Signal to Noise ratio, Attenuation and Interference). The simulation model was examined using different topologies under various conditions and numerous results were collected.

Simulation Framework of Wireless Sensor Network (WSN) Using MATLAB/SIMULINK Software 271

generated in the Simulink model by a baseband MFSK block set to 79 symbols and a separation of 1MHz. If a hop frequency value 0 is input, a -39MHz complex sinusoid is generated. If a 1 is entered, a -38 MHz complex sinusoid is generated and so on. In the model, the hop sequences are generated by a simple random number generator, not using the actual method specified in the standard. The transmitter is turned off after 366 bits using a Gain block to multiply the frame with a mask of 36600 ones and 26500 zeros.

 AWGN Channel: The AWGN Channel block adds white Gaussian noise to a real or complex input signal. When the input signal is real, this block adds real Gaussian noise and produces a real output signal. When the input signal is complex, this block adds

 Path Loss: This block reduces the amplitude of the input signal by an amount specified. The loss can be specified directly using the "Decibels" mode, or indirectly using the "Distance and Frequency" mode. The reciprocal of the loss is applied as a gain, e.g., a loss of +20 dB, which reduces the signal by a factor of 10 corresponds to a gain value of 0.1. 802.11b interferer: This block adds signals that have the same frequency of the data signal to make interference between the data signal and other signals( i.e. a Wireless

 Multiport Switch: In order to simulate the multiple access and multiplexing functions of the channel, this block was used. It chooses between a number of inputs. The first input is called the control input, while the rest of the inputs are called data inputs. The value of the control input determines which data input is passed through to the output port.

 Demodulation and decoding: This block is used to extract the original information-bearing signal from a modulated carrier wave, and to recover the information contends in it. Zero-Order Hold: This block samples and holds its input for the specified sample period. The block accepts one input and generates one output, both of which can be scalar or vector. If the input is a vector, all elements of the vector are held for the same

 Un-buffer: This block un-buffers an Mi-by-N frame-based input into a 1-by-N samplebased output. That is, inputs are un-buffered row-wise so that each matrix row becomes an independent time-sample in the output. The rate at which the block receives inputs

Down-sampling to 8ksamples/s: This block down-samples the input to a lower rate by

Scope RX**:** It was used to display the received signal and compare it with the original

As known, a piconet can includes up to seven slaves and one master. In this example three signals were sent from three sensors (slaves) to the receiving component (master) representing one piconet , the information obtained by the sensors are used to estimate the

is generally less than the rate at which the block produces outputs.

**2. The medium which consists of the following blocks** 

Local Area Network (WLAN) transmission).

**3. The receiver consists of the following blocks:** 

sample period.

deleting the repeating samples.

signal to discover the system behavior.

Hop Sequence Generator: same as mentioned earlier.

complex Gaussian noise and produces a complex output signal.

**Figure 2.** Simple WSN model

The architecture of the system could be explained as follows:

#### **1. The transmitter**

This system was based on Bluetooth technology that is considered as the backbone of transmission operation. Bluetooth is a short-range radio link technology that operates in the 2.4 GHz Industrial, Scientific, and Medical (ISM) band[60]. In this system we modulated the signal using Gaussian frequency shift keying (GFSK) over a radio channel with maximum capacity of 1 Mbps. The transmitter consists of the following blocks:


generated in the Simulink model by a baseband MFSK block set to 79 symbols and a separation of 1MHz. If a hop frequency value 0 is input, a -39MHz complex sinusoid is generated. If a 1 is entered, a -38 MHz complex sinusoid is generated and so on. In the model, the hop sequences are generated by a simple random number generator, not using the actual method specified in the standard. The transmitter is turned off after 366 bits using a Gain block to multiply the frame with a mask of 36600 ones and 26500 zeros.

## **2. The medium which consists of the following blocks**

270 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 2

**Figure 2.** Simple WSN model

using 256 quantization level.

between samples.

**1. The transmitter** 

The architecture of the system could be explained as follows:

capacity of 1 Mbps. The transmitter consists of the following blocks:

This system was based on Bluetooth technology that is considered as the backbone of transmission operation. Bluetooth is a short-range radio link technology that operates in the 2.4 GHz Industrial, Scientific, and Medical (ISM) band[60]. In this system we modulated the signal using Gaussian frequency shift keying (GFSK) over a radio channel with maximum

 Sensor signal stage: It is represented by a sensor to sense the physical signals such as temperature, pressure…etc, then transducing them into an electrical signal. In addition, this stage includes the A/D convertor which converts the signal from Analog to Digital

Up-sampling to 64ksamples/s: Up-samples the input to a higher rate by inserting zeros

Payload FEC encode: Encodes the data to enable error correction(an FEC encoder may

 Bluetooth Clock: Each Bluetooth device has a free-running 28-bit Bluetooth clock. The clock ticks 3,200 times per second or once every 312.5 µsec, representing a clock rate of 3.2 KHz. Hop Sequence Generator: For devices to communicate with each other, they must transmit and receive on the same frequency at the same time. The hop sequence generator generates a sequence of hop frequencies in the range 0 to 78. It can generate

either the connection state hop sequence, a random white sequence, or be fixed. Encoder and modulator: The 366 data bits are transmitted at 1 Mbps and modulated using Gaussian frequency shift keying (GFSK). GFSK effectively transmits +150 kHz signal relative to the carrier for a 1bit, and a 150 kHz signal for a 0 bit. The carrier signal is

include a binary convolutional encoder followed by a puncturing device).


### **3. The receiver consists of the following blocks:**


As known, a piconet can includes up to seven slaves and one master. In this example three signals were sent from three sensors (slaves) to the receiving component (master) representing one piconet , the information obtained by the sensors are used to estimate the

Bluetooth performance as well as to study the media effect. Noise and interference are added to the signals in order to simulate the channel effect and measure Bit Error Rate (BER) and Frame Error Rate (FER). The following figures shows the system performance under different working conditions.

**Figure 3.** Signals sent from the three sensors

Simulation Framework of Wireless Sensor Network (WSN) Using MATLAB/SIMULINK Software 273

(a) SNR & BER (b) SNR & FER

(a) Average Rate=6 (b) Average Rate=12

(c) Average Rate=25 (d) Average Rate=50 (e)Average Rate=100

**Figure 5.** Relationship between SNR & (BER, FER)

**Figure 6.** Received signals with different rate of interference

(a) SNR=20dB (b) SNR=15dB

(c) SNR=12dB (d) SNR=10dB

**Figure 4.** Received signals with different SNR

**Figure 5.** Relationship between SNR & (BER, FER)

different working conditions.

**Figure 3.** Signals sent from the three sensors

**Figure 4.** Received signals with different SNR

Bluetooth performance as well as to study the media effect. Noise and interference are added to the signals in order to simulate the channel effect and measure Bit Error Rate (BER) and Frame Error Rate (FER). The following figures shows the system performance under

(a) SNR=20dB (b) SNR=15dB

(c) SNR=12dB (d) SNR=10dB

**Figure 6.** Received signals with different rate of interference

Simulation Framework of Wireless Sensor Network (WSN) Using MATLAB/SIMULINK Software 275

**Figure 9.** WSN model with masters & slaves

The following figures shows the system performance under different working conditions.

**Figure 7.** Received signals with different rate of interference & different SNR

**Figure 8.** Relationship between interference & (BER, FER)

#### **3.2. More complex example**

As known from Blutooth operation, each piconet consists of one master and seven slaves and each master of a specific piconet may acts as a slave for another piconet which means the ability to expand the network to respond to more than seven sensors.

In this example two piconets are connected, so that the first piconet consists of three sensors connected to the master, and the later is connected as a slave to the second piconet. The second piconet consists of two slaves and one master as shown in Figure 9 below:

**Figure 9.** WSN model with masters & slaves

**Figure 7.** Received signals with different rate of interference & different SNR

(a) SNR=15dB, Average Rate=6 (b) SNR=12dB, Average Rate=25

**Figure 8.** Relationship between interference & (BER, FER)

As known from Blutooth operation, each piconet consists of one master and seven slaves and each master of a specific piconet may acts as a slave for another piconet which means

 (a) Interference & BER (b) Interference & FER

In this example two piconets are connected, so that the first piconet consists of three sensors connected to the master, and the later is connected as a slave to the second piconet. The second piconet consists of two slaves and one master as shown in Figure 9

the ability to expand the network to respond to more than seven sensors.

**3.2. More complex example** 

below:

The following figures shows the system performance under different working conditions.

Simulation Framework of Wireless Sensor Network (WSN) Using MATLAB/SIMULINK Software 277

(a) Average Rate=6 (b) Average Rate=12

(c) Average Rate=25 (d) Average Rate=50

**Figure 12.** Received signals with different rate of interference

**Figure 13.** Relationship between interference & (BER, FER)

From the above results, it is obvious that the behavior of the system was successfully described using the suggested simulation methodology. It is also important to mention that this simulation method provides the ability to change the different system parameters to create new environment and hence, new simulation scenarios. This new simulation

(a) Interference & BER (b) Interference & FER

**Figure 10.** Received signals with different Signal to Noise Ratio (SNR)

**Figure 11.** Relationship between SNR & (BER, FER)

**Figure 12.** Received signals with different rate of interference

(a) SNR=20dB (a) SNR=20dB

(a) SNR=20dB (a) SNR=20dB

(a) SNR & BER (b) SNR & FER

**Figure 10.** Received signals with different Signal to Noise Ratio (SNR)

**Figure 11.** Relationship between SNR & (BER, FER)

**Figure 13.** Relationship between interference & (BER, FER)

From the above results, it is obvious that the behavior of the system was successfully described using the suggested simulation methodology. It is also important to mention that this simulation method provides the ability to change the different system parameters to create new environment and hence, new simulation scenarios. This new simulation

methodology proves the ability of the Simulink MATLAB to be a useful and flexible approach to study the effect of different physical layer parameters on the performance of wireless sensor networks.

Simulation Framework of Wireless Sensor Network (WSN) Using MATLAB/SIMULINK Software 279

Includes sensor network package containing models such as propagation, battery, radio model and sen-sor protocol stack including AODV and 802.11 → project seems to be abandoned

Complete emulation of the AVR instruction set with partial Mica2 support; TinyOS based code can be run → very special application area; slow simulation speed; project seems to be abandoned

Provides network functionality for wireless embedded systems; EmTOS can be used to run TinyOS applications as EmStar module → project seems be abandoned (download links broken)

Provides very basic network and physical layer sup-port. Source can be compiled for TinyOS. → project does not seem to be developed any further

Includes battery and power models, MAC layers (802.11) as well as network protocols (AODV, DSR, SSR, SHR) → does not seem to be developed any further

Algorithmic approach that concentrates on lower layers, no special WSN protocols → very active project -lot of recent changes in SVN

Detailed simulation of the end nodes and their architecture, Physical layer parameters, different modulation & encoding techniques, communication channel modeling(SNR, effect of different Noise schemes, Interference, distance, etc..), various methods to monitor and record results, making use of the rich library of Matlab/Simulink.

J-Sim Java; configuration

ATEMU AVR micro-

EmStar C

SENS C++

SENSE C++

Shawn C++

SIMULINK C , Java

**Table 1.** The features of the different simulation methods of WSN

MATLAB

Tcl/Java

controller binaries

Table (1) below, summarize the features of the different simulation methods of WSN including the suggested one in this chapter, MATLAB SIMULINK.


including the suggested one in this chapter, MATLAB SIMULINK.

Programming

configuration by GUI; internals C++

basic modules C++; larger structures NED

OTcl

controller binaries

wireless sensor networks.

Simulation Environment

GloMoSim/

OPNET Mod-eler Wireless Suite

TOSSIM (part of

OMNeT++

TinyOS) nesC

NS-2 C++; configu-ration

Avrora AVR micro-

QualNet C and Parsec

methodology proves the ability of the Simulink MATLAB to be a useful and flexible approach to study the effect of different physical layer parameters on the performance of

Table (1) below, summarize the features of the different simulation methods of WSN

Language **WSN support** 

GloMoSim: basic mobility and radio propagation models; 802.11; QualNet: additionally battery and energy model; ZigBee → GloMoSim seems to be outdated; QualNet seems to be more up-to-date, but commercial

Different propagation models; 802.11, ZigBee; some MANET protocols, but no special WSN support → powerful tool with a nice GUI, but expensive

All TinyOS-based WSN protocols can be simulated with TOSSIM without modifications → good ap-proach especially if implementation should also be used with TinyOS-based nodes

Several frameworks that add WSN functionality to OMNet++ such as MiXiM, Castalia, etc. → active project with a huge user base; Eclipse-based IDE for development

Huge amount of protocols available contributed by NS-2 users → complex configuration; unclear situa-tion due to large number of different user contributed implementations

Particularly for programs written for AVR micro-controller with support for support for Mica2 and Mi-caZ → very special application area; project seems to be still active -still changes in CVS


**Table 1.** The features of the different simulation methods of WSN

## **4. Conclusions**

In this chapter, a new simulation methodology of wireless sensor networks (WSN) was presented. MATLAB/Simulink was used as the tool to build the simulation environment. The strength of this simulation method falls in the ability to study the effect of different physical layer parameters (channel noise and interference, Signal to noise ratio…etc.) on the system behavior. The other advantage of this method is its flexibility in building the end nodes and sensors. This simulation methodology could be used to build different WSN types and opens the doors to use the MATLAB in this new field.

Simulation Framework of Wireless Sensor Network (WSN) Using MATLAB/SIMULINK Software 281

http://www.opnet.com/support/des model library/images/MANET scrnsht.jpg, 2012.

[13] Levis, P., Lee, N.,Welsh, M., Culler, M., "TOSSIM: Accurate and scalable simulation of entire TinyOS applications," in Proceedings of the 1st international conference on Embedded networked sensor systems. ACM New York, NY, USA, 2003, pp. 126–

[14] Levis, P., Lee, N., "Tossim: A simulator for tinyos networks," UC Berkeley, September,

[15] Notani, S., "Performance Simulation of Multihop Routing Algorithms for Ad-Hoc Wireless Sensor Networks Using TOSSIM," in Advanced Communication Technology,

[16] Computer Science Division at UC Berkeley. Visualisation of a TOSSIM simulation with TinyViz. [Online]. Available: http://www.tinyos.net/tinyos-1.x/doc/tutorial/imgs/

[18] Varga, A., et al., "The OMNeT++ discrete event simulation system," in Proceedings of

[19] Varga, A., "OMNeT++ Discrete event simulation system. User Manual," Technical

[20] Varga, A., Hornig, R., "An overview of the OMNeT++ simulation environment," in Proceedings of the 1st international conference on Simulation tools and techniques for communications, networks and systems & workshops table of contents. ICST (Institute for Computer Sciences, Social-Informatics and Telecommunications Engineering) ICST,

[22] Mobility Framework for OMNeT++ Community. [Online]. Available: http://mobility-

[23] Drytkiewicz, W., Sroka, S., Handziski, V., Koepke, A., Karl, H., "A mobility framework

[24] L¨obbers, M., Willkomm, D., K¨opke, A., Karl, H., "Framework for Simulation of

[25] MiXiM developers. MiXiM project. [Online]. Available: http://mixim.sourceforge.net/,

[26] K¨opke, A., Swigulski, M., Wessel, K., Willkomm, D., Haneveld, P., Parker, T., Visser, O., Lichte, H. Valentin, S., "Simulating wireless and mobile networks in OMNeT++ the MiXiM vision," in Proceedings of the 1st international conference on Simulation tools and techniques for communications, networks and systems & workshops table of contents. ICST (Institute for Computer Sciences, Social-

the European Simulation Multiconference (ESM'2001), 2001, pp. 319–324.

http://omnetpp.org/doc/omnetpp40/ide-overview/pictures/img1.png, 2012.

[11] OPNET Technologies, Inc. [Online]. Available:

137.

2003.

tinyviz-screenshot1.gif, 2012.

Brussels, Belgium, Belgium, 2008.

fw.sourceforge.net, 2012.

2012.

http://www.omnetpp.org/

[12] Computer Science Division at UC Berkeley. [Online]. Available: http://www.cs.berkeley.edu/\_pal/research/tossim.html, 2012

2008. ICACT 2008. 10th International Conference on, vol. 1, 2008.

[17] OMNeT++ Community. (2010, May) OMNeT++. [Online]. Available:

University of Budapest, Dept. of Telecommunications, 2006.

[21] OMNeT++ Community. OMNeT++ 4.0 IDE.[Online]. Available:

for omnet++," in 3rd International OMNeT++ Workshop, 2003.

Mobility in OMNeT++(Mobility Framework)," 2004.

## **Author details**

Qutaiba I. Ali *Mosul University, Computer Engineering Department, Iraq*

## **Acknowledgement**

Many thanks to my students Akram & Hussien for their assistance.

## **5. References**


[11] OPNET Technologies, Inc. [Online]. Available:

280 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 2

*Mosul University, Computer Engineering Department, Iraq*

Many thanks to my students Akram & Hussien for their assistance.

2007 IEEE IMTC conference. Poland, Warsaw 2007; 1(3): 1-6.

http://www.scalable-networks.com/products/qualnet/,2012.

[8] OPNET Technologies, Inc. [Online]. Available: http://www.opnet.com/,2012.

networks," Jisuanji Gongcheng/ Computer Engineering, vol. 33, no. 4, 2007.

[7] S. Technologies, "Qualnet v. 3.9. 5 user's guide," 2006.

In this chapter, a new simulation methodology of wireless sensor networks (WSN) was presented. MATLAB/Simulink was used as the tool to build the simulation environment. The strength of this simulation method falls in the ability to study the effect of different physical layer parameters (channel noise and interference, Signal to noise ratio…etc.) on the system behavior. The other advantage of this method is its flexibility in building the end nodes and sensors. This simulation methodology could be used to build different WSN types and opens the doors to use the MATLAB in this

[1] Lewis, F "Wireless Sensors Networks, Smart Environments: Technologies, Protocols, and Applications', ed. Cook DJ, Das SK, John Wiley, New York, 2004; 1-18. [2] Akyildiz, IF, Sankarasubramaniam, F, Cayirci, E, "A Survey on Sensor Networks", IEEE

[3] Rabaey, J, Ammer, M, da Silva, J.L., D. Patel, and S. Roundy, "Picoradio supports ad hoc ultra-low power wireless networking," *Computer*, vol. 33, no. 7, pp. 42–48, July

[4] Gupta, G, Mukhopadhyay, SC, Sutherland, M., Demidenko, S., "Wireless Sensor Network for Selective Activity Monitoring in a home for the Elderly", Proceedings of

[5] Callaway, E., Gorday, P. , Hester, L., "Home Networking with IEEE 802.15.4: A Developing Standard for Low-Rate Wireless Personal Area Networks", IEEE Commun

[9] Jiang, H., Wang, P. ,Liu, H., "Research on OPNET simulation model in wireless sensor

[10] Jurc´k, P., Koubˆaa,A., "The IEEE 802.15. 4 OPNET Simulation Model: Reference

**4. Conclusions** 

new field.

**Author details** 

**Acknowledgement** 

Commun Mag 2002; 102-114.

Qutaiba I. Ali

**5. References** 

2000.

Mag 2002; 69-77.

Guide v2. 0," 2007.

[6] QualNet. [Online]. Available:

http://www.opnet.com/support/des model library/images/MANET scrnsht.jpg, 2012.


Informatics and Telecommunications Engineering) ICST, Brussels, Belgium, Belgium, 2008.

Simulation Framework of Wireless Sensor Network (WSN) Using MATLAB/SIMULINK Software 283

[42] Titzer, B., "Avrora: The AVR simulation and analysis framework," Master's thesis,

[43] Titzer, B., Lee, D. Palsberg, J., "Avrora: Scalable sensor network simulation with precise timing," in Proceedings of the 4th international symposium on Information processing

[44] Department of Computer Science at University of Illinois at Urbana-Champaign). J-Sim.

[45] Sobeih, A., Chen, W., Hou, J., Kung, L., "J-sim: A simulation environment for wireless sensor networks," in Proceedings of the 38th annual Symposium on Simulation. IEEE

[46] Hou, J., Kung, L., "J-Sim: A Simulation and emulation environment for wireless sensor networks," IEEE Wireless Communications Magazine, vol. 13, no. 4, pp. 104–119,

[47] Center for Satellite and Hybrid Communication Networks (CSHCN) at University of Maryland. Atemu. [Online]. Available: http://www.cshcn.umd.edu/research/atemu/,

[48] Polley, J., Blazakis, D., "Atemu: A fine-grained sensor network simulator," in Sensor and Ad Hoc Communications and Networks, 2004. IEEE SECON 2004. 2004 First

[49] Laboratory for Embedded Collaborative Systems (LECS) at UCLA. [Online]. Available:

[50] Elson, J., Bien, S., "Emstar: An environment for developing wireless embedded systems software," Center for Embedded Networked Sensing (CENS) Technical Report, vol. 9,

[51] Girod, L., Ramanathan, N.," A Software Environment for Developing and Deploying Heterogeneous Sensor Actuator Networks," Center for Embedded Network Sensing, p.

[52] Open Systems Laboratory at University of Illinois at Urbana-Champaign. [Online].

[53] Sundresh, S., Kim, W., Agha, G., "SENS: A sensor, environment and network simulator," in Proceedings of the 37th annual symposium on Simulation. IEEE

[54] Computer Science Department at Rensselaer Polytechnic Institute (RPI). [Online].

[55] Chen, G., Branch, J., "Sense: A sensor network simulator," Advances in Pervasive

[57] Kroeller, A., Pfisterer, D., "Shawn: A new approach to simulating wireless sensor

[58] Fekete, S., Kroller, A., Fischer, S., Pfisterer, D., "Shawn: The fast, highly customizable sensor network simulator," in Networked Sensing Systems, 2007. INSS'07. Fourth

[56] SwarmNet project. [Online]. Available: http: //shawn.sourceforge.net/, 2012.

Annual IEEE Communications Society Conference on, 2004, pp. 145–152.

University of California, Los Angeles, 2004.

http://www.lecs.cs.ucla.edu/emstar/, 2012.

Available: http://osl.cs.uiuc.edu/sens/, 2012.

Computer Society Washington, DC, USA, 2004.

Available: http://www.ita.cs.rpi.edu/sense/, 2012.

Computing and Networking, pp. 249–267, 2004.

networks," Arxiv preprint cs/0502003, 2005.

International Conference on, 2007, pp.299–299. [59] MATLAB Web Site: http://www.mathworks.com/

2006.

2012.

2003.

101, 2007.

in sensor networks. IEEE Press Piscataway, NJ, USA, 2005.

Computer Society Washington, DC, USA, 2005, pp. 175–187.

[Online]. Available: http://sites.google.com/site/jsimofficial, 2012.


[42] Titzer, B., "Avrora: The AVR simulation and analysis framework," Master's thesis, University of California, Los Angeles, 2004.

282 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 2

Available: http://www.consensus.tudelft.nl/software.html, 2012.

systems. ACM New York, NY, USA, 2007, pp. 407–408.

[34] Downard I., DC, N., "Simulating sensor networks in ns-2,"2004.

http://apachepersonal.miun.se/\_qinwan/resources.htm, 2012.

http://www.isi.edu/nsnam/nam/nambig.gif, 2012.

http://nesl.ee.ucla.edu/projects/sensorsim/, 2012.

[41] UCLA Compilers Group). Avrora. [Online]. Available:

http://compilers.cs.ucla.edu/avrora/, 2012.

http://inet.omnetpp.org/, 2012.

http://nesct.sourceforge.net/, 2012.

edu/nsnam/ns/, 2012.

dcc.ufmg.br/, 2012.

2000, pp. 104–111.

2006.

2008.

2012.

2008.

Informatics and Telecommunications Engineering) ICST, Brussels, Belgium, Belgium,

[27] University of Twente and TU Delft. Positif, MAC Simulator and T-MAC. [Online].

[28] National ICT Australia. Castalia. [Online]. Available: http://castalia.npc.nicta.com.au/,

[29] Boulis, A., "Castalia: revealing pitfalls in designing distributed algorithms in WSN," in Proceedings of the 5th international conference on Embedded networked sensor

[30] OMNeT++ Community. INET framework for the OMNeT++. [Online]. Available:

[31] Ariza-Quintana, A., Casilari, E., Cabrera, A., "Implementation of MANET routing protocols on OMNeT++," in Proceedings of the 1st international conference on Simulation tools and techniques for communications, networks and systems & workshops table of contents. ICST (Institute for Computer Sciences, Social-Informatics and Telecommunications Engineering) ICST, Brussels, Belgium, Belgium,

[32] OMNeT++ Community. NesCT for the OMNeT++. [Online]. Available:

[33] NS-2 developers. The Network Simulator – ns-2. [Online]. Available: http://www.isi.

[35] NS-2 developers. Visualisation of a ns-2 simulation with NAM. [Online]. Available:

[36] Departamento de Ciˆencia da Computac¸ ˜ao, Universidade Federal deMinas Gerais. Mannasim Framework for ns-2. [Online]. Available: http://www.mannasim.

[37] Computer Science, Mid Sweden University , Sweden. NS2-MIUN. [Online]. Available:

[38] Networked and Embedded Systems Laboratory (NESL) at the University of California at Los Angeles (UCLA). SensorSim framework. [Online]. Available:

[39] Park, S., Savvides, A., Srivastava, M., "SensorSim: a simulation framework for sensor networks," in Proceedings of the 3rd ACM international workshop on Modeling, analysis and simulation of wireless and mobile systems. ACM New York, NY, USA,

[40] "Sensor Sim: A Simulation Framework for Networks Sensors", Electrical Engineering Department, University of California, Los Angeles, Retrieved October, vol. 16,

	- [60] Ali, Q., Abdulmaojod, A., Ahmed, H.," Simulation & Performance Study of Wireless Sensor Network (WSN) Using MATLAB, IJEEE Journal,2010.

**Chapter 0**

**Chapter 13**

**Semi-Analytic Techniques for**

Additional information is available at the end of the chapter

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 Bit Error Rate (BER) are difficult to determine analytically for transmission schemes involving correction codes and communication channels with Inter-Symbol Interference (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, BERs of the order of 10−<sup>8</sup> <sup>−</sup> <sup>10</sup>−<sup>9</sup> require at least 109 <sup>−</sup> 1010 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.

powerful tool for the analysis of complex systems.

cited.

In order to overcome the limitation of Monte Carlo and analytical techniques, semi-analytic approaches have been previously implemented [1, 2]. 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

> ©2012 Borio and Cano, licensee InTech. This is an open access chapter distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly

©2012 Borio and Cano, licensee InTech. This is a paper distributed under the terms of the Creative CommonsAttribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use,

distribution, and reproduction in any medium, provided the original work is properly cited.

**Fast MATLAB Simulations**

Daniele Borio and Eduardo Cano

http://dx.doi.org/10.5772/46470

**1. Introduction**

**Chapter 0 Chapter 13**

## **Semi-Analytic Techniques for Fast MATLAB Simulations**

Daniele Borio and Eduardo Cano

Additional information is available at the end of the chapter

http://dx.doi.org/10.5772/46470

## **1. Introduction**

284 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 2

Sensor Network (WSN) Using MATLAB, IJEEE Journal,2010.

[60] Ali, Q., Abdulmaojod, A., Ahmed, H.," Simulation & Performance Study of Wireless

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 Bit Error Rate (BER) are difficult to determine analytically for transmission schemes involving correction codes and communication channels with Inter-Symbol Interference (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, BERs of the order of 10−<sup>8</sup> <sup>−</sup> <sup>10</sup>−<sup>9</sup> require at least 109 <sup>−</sup> 1010 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 [1, 2]. 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.

©2012 Borio and Cano, licensee InTech. This is an open access chapter distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. ©2012 Borio and Cano, licensee InTech. This is a paper distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

A classical example is the model of a transmission chain where a Traveling Wave Tube Amplifier (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

The TWTA is a memory-less device and can be characterized using *AM-AM conversion* and *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

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 Energy per bit (*Eb*) is computed. Since the noise properties are

erfc

where *Eb*,*<sup>i</sup>* has been obtained by simulating the transmission chain (including the non-linear device) in the absence of noise and transmitting the symbol *si*. The parameter *N*<sup>0</sup> is the noise power spectral density and it is a known value of the system. Finally, the BER of the system

> *Ns*−1 ∑ *i*=0

2

where *pi* is the probability that the symbol *si* will be transmitted and *Ns* is the number of

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

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.

When considering the previous example it is possible to identify three building blocks that

*Eb*,*<sup>i</sup> N*<sup>0</sup>

*pi*erfc

*Eb*,*<sup>i</sup> N*<sup>0</sup>

, (1)

Semi-Analytic Techniques for Fast MATLAB Simulations 287

, (2)

*BERi* <sup>=</sup> <sup>1</sup>

*BER* <sup>=</sup> <sup>E</sup> [*BERi*] <sup>=</sup> <sup>1</sup>

2

non-linear amplifier.

non-linear devices.

is obtained from

**1.1. Building blocks**

symbols of the signal constellation.

transmission of the signal component is simulated.

play different roles in the semi-analytic framework:

known, the BER for the *i*th symbol, *si*, is given by

**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.

**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.

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.

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 MATLAB 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.

A classical example is the model of a transmission chain where a Traveling Wave Tube Amplifier (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 conversion* and *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 Energy per bit (*Eb*) is computed. Since the noise properties are known, the BER for the *i*th symbol, *si*, is given by

$$BER\_{l} = \frac{1}{2} \text{erfc}\left(\sqrt{\frac{E\_{b,i}}{N\_0}}\right),\tag{1}$$

where *Eb*,*<sup>i</sup>* has been obtained by simulating the transmission chain (including the non-linear device) in the absence of noise and transmitting the symbol *si*. The parameter *N*<sup>0</sup> 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 = \operatorname{E}\left[BER\_{i}\right] = \frac{1}{2} \sum\_{i=0}^{N\_{\text{s}}-1} p\_{i} \text{erfc}\left(\sqrt{\frac{E\_{b,i}}{N\_{0}}}\right),\tag{2}$$

where *pi* is the probability that the symbol *si* will be transmitted and *Ns* 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**

2 MATLAB / Book 2

**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)

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

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 MATLAB 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.

communications channel is modeled as a non-linear device, which affects only the signal component,

**Figure 2.** Model adopted for the evaluation of the BER using a semi-analytic approach. The

complexity, combining the strengths of Monte Carlo and analytical approaches.

 

> 

 

> 

 

between analytical models and Monte Carlo simulations.

and the addition of a noise term supposed to be Gaussian.

 

 

> 

> >

> > > -


 

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


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.

 

 

 

analytic components of the system.

advanced tracking loops for new GNSS signals.

closed-loop approach will be detailed in Section 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

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

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

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

metrics such as tracking jitter or Mean Time to Lose Lock (MTLL) [12].

 

 

 

 

Semi-Analytic Techniques for Fast MATLAB Simulations 289

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 Energy per bit to Noise power spectral density ratio (*Eb*/*N*0).

This type of configuration is defined here as *sequential* since 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.

A second type of configuration has been recently considered for the analysis of tracking loops in Direct Sequence Spread Spectrum (DSSS) and Global Navigation Satellite System (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 Integrate and Dump (I&D) blocks that rely on simple operations that can be analytically modeled. For this reason, semi-analytic models exploiting the knowledge of the I&D blocks and simulating only the non-linear parts of the system have been developed [7–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 I&D have been suggested. [7] modeled the correlator outputs evaluated by the I&D blocks as linear combinations of independent Gaussian random variables. Correlation among the different correlators was obtained by using, for the generation of different I&D 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

4 MATLAB / Book 2

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

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

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 Energy per bit

This type of configuration is defined here as *sequential* since there is no feedback between the different blocks and each element of the chain is run sequentially. The principle of this

A second type of configuration has been recently considered for the analysis of tracking loops in Direct Sequence Spread Spectrum (DSSS) and Global Navigation Satellite System (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 Integrate and Dump (I&D) blocks that rely on simple operations that can be analytically modeled. For this reason, semi-analytic models exploiting the knowledge of the I&D blocks and simulating only the non-linear parts of the system have been developed [7–10]. This resulted in efficient analysis tools, which require

Different techniques for modeling the output of the I&D have been suggested. [7] modeled the correlator outputs evaluated by the I&D blocks as linear combinations of independent Gaussian random variables. Correlation among the different correlators was obtained by using, for the generation of different I&D 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

reduced processing time with the applicability of Monte Carlo simulations.

is used to determine *Eb* from the simulated decision variables.

• **simulation block** • **estimation block** • **analytical model**.

and Section 3.

**1.2. Main configurations**

to Noise power spectral density ratio (*Eb*/*N*0).

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.

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 Mean Time to Lose Lock (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 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 Additive White Gaussian Noise (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.

The signal at the input of a GNSS receiver, in a one-path AWGN channel and in the presence

where *yl*(*t*) is the signal transmitted by the *l*th GNSS L1 satellite, *L* is the total number of

• *cl*(·) is the *l*th pseudo-random sequence extracted from a family of quasi-orthogonal codes

• *τ*0,*l*, *fd*0,*<sup>l</sup>* and *ϕ*0,*<sup>l</sup>* are the delay, Doppler frequency and phase introduced by the

It is noted that GNSS signals adopt a DSSS modulation. The pseudo-random sequences, *cl*(*t*), allow the simultaneous transmission of several signals at the same time and in the same band. Moreover, *cl*(*t*) sequences are characterized by sharp correlation functions that allow the precise measurement of the signal travel time. The travel time is then converted into

The pseudo-random sequence, *cl*(*t*), is composed of several terms including a primary

*yl*(*t*) + *i*(*t*) + *η*(*t*), (3)

Semi-Analytic Techniques for Fast MATLAB Simulations 291

2*π*(*fRF* + *fd*0,*l*)*t* + *ϕ*0,*<sup>l</sup>*

*cl*,(*<sup>i</sup>* mod *Nc* )*sb*(*t* − *iTh*). (5)

, (4)

*r*(*t*) =

*L*−1 ∑ *l*=0

satellites in view, *i*(*t*) is the received interference signal and *η*(*t*) is the noise term.

**2.1. The GNSS signal**

where

of RF interference, can be modeled as

Each useful signal, *yl*(*t*), can be expressed as

*yl*(*t*) = <sup>2</sup>*Cldl*

• *Cl* is the power of the *l*th useful signal;

and used for spreading the signal spectrum;

• *fRF* is the centre frequency of the GNSS signal.

distances that allows a GNSS receiver to determine its position.

*cl*(*t*) =

+∞ ∑ *i*=−∞

• *dl*(·) is the navigation message;

communication channel, and

spreading sequence and a subcarrier:

 *t* − *τ*0,*<sup>l</sup> cl t* − *τ*0,*<sup>l</sup>* cos

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 Global Positioning System (GPS) L1 receiver affected by an interference signal caused by third order harmonics of a Digital Video Broadcasting - Terrestrial (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 Radio-Frequency (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 Ultra Wideband (UWB) transmissions [15], Distance Measuring Equipment (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 Additive White Gaussian Noise (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(t) = \sum\_{l=0}^{L-1} y\_l(t) + i(t) + \eta(t),\tag{3}$$

where *yl*(*t*) is the signal transmitted by the *l*th GNSS L1 satellite, *L* is the total number of satellites in view, *i*(*t*) is the received interference signal and *η*(*t*) is the noise term.

Each useful signal, *yl*(*t*), can be expressed as

$$y\_I(t) = \sqrt{2\mathbf{C}\_I} d\_I \left( t - \pi\_{0,l} \right) c\_l \left( t - \pi\_{0,l} \right) \cos \left( 2\pi (f\_{RF} + f\_{d0,l})t + \varphi\_{0,l} \right), \tag{4}$$

where

6 MATLAB / Book 2

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

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

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 Global Positioning System (GPS) L1 receiver affected by an interference signal caused by third order harmonics of a Digital Video Broadcasting - Terrestrial (DVB-T) signal. A comprehensive description of the sequential approach applied to an inter-system interference case is provided

The reception of GNSS signals is challenging due to low signal power, possible severe channel conditions and the presence of Radio-Frequency (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 Ultra Wideband (UWB) transmissions [15], Distance Measuring Equipment (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

**2. Sequential approach: The inter-system interference case**

demanding error counting processes.

behind the sequential semi-analytic approach.

in the following.

loops is discussed in Section 3.


It is noted that GNSS signals adopt a DSSS modulation. The pseudo-random sequences, *cl*(*t*), allow the simultaneous transmission of several signals at the same time and in the same band. Moreover, *cl*(*t*) 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, *cl*(*t*), is composed of several terms including a primary spreading sequence and a subcarrier:

$$c\_l(t) = \sum\_{i = -\infty}^{+\infty} c\_{l,(i \text{ mod } N\_\ell)} s\_b(t - iT\_h). \tag{5}$$

The signal *sb*(*t* − *iTh*) in (5) represents the subcarrier of duration *Th*, which determines the spectral characteristics of the transmitted GNSS signal. The GPS L1 Coarse/Acquisition (C/A) component is Bi-Phase Shift Keying (BPSK) modulated, whereas the Galileo E1 signal adopts a Composite Binary-Offset Carrier (CBOC) scheme. The sequence *cl*,*i*, of length *Nc*, 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(t) = y(t) + i(t) + \eta(t),\tag{6}$$

where *M* is the modulation order, *h* is the subcarrier index, *p* is the symbol index, *Nd* represents the total number of transmitted symbols and *Ip*,*<sup>h</sup>* 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

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

where *an* 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

The terms of order *n* > 1 in (10) model the amplifier non-linearities and the ratios *an*/*a*<sup>1</sup> are expected to be small for *n* > 1. Since only the third harmonics will fall into the GPS L1 band,

3

The signal *i*(*t*) is filtered and down-converted by the receiver front-end and signal *iF*(*t*), the

After signal conditioning, the sequence *rBB*[*n*] is correlated with local replicas of the useful signal code and carrier as shown in Figure 4. Since the code delay, *τ*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 *rBB*[*n*], can be computed in order to produce a decision variable less affected by noise and interference. In this way, the output

where *τ*, *fd* and *ϕ* 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 *Tc* = *NTs* 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,

*i*(*t*) = *a*3*i*

*DVB*−*T*(*t*), (10)

Semi-Analytic Techniques for Fast MATLAB Simulations 293

*DVB*−*T*(*t*). (11)

*rBB*[*n*]*c* (*nTs* − *τ*) exp {−*j*2*π fdnTs* − *jϕ*} , (12)

∞ ∑ *n*=1 *ani n*

*p*(*t*) =

conversion functions discussed in Section 1 for the TWTA case.

filtered version of *i*(*t*), will affect receiver operations.

of the *k*th complex correlator can be expressed as

(*k*+1)*N*−1 ∑ *n*=*kN*

*Sk* <sup>=</sup> <sup>1</sup> *N*

the interference signal at the antenna of a GNSS receiver is given by

Finally, the interference term in (7) can be modeled as *iBB*[*n*] = *iF*(*nTs*).

signal up-conversion.

amplifier input/output function:

**2.3. The acquisition process**

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*(*t*), 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}[n] = y\_{BB}(nT\_s) + i\_{BB}(nT\_s) + \eta\_{BB}(nT\_s) = y\_{BB}[n] + i\_{BB}[n] + \eta\_{BB}[n],\tag{7}$$

where the notation *x*[*n*] is used to denote discrete time sequences sampled at the frequency *fs* = <sup>1</sup> *Ts* . In addition, the index "BB" is used to denote a filtered signal down-converted to baseband. Furthermore, the signal *yBB*[*n*] in (7) can be written as

$$y\_{BB}[n] = \sqrt{\mathbb{C}}d\left(nT\_{\rm s} - \tau\_{0}\right) \varepsilon\left(nT\_{\rm s} - \tau\_{0}\right) \exp\left\{j2\pi f\_{0}nT\_{\rm s} + \varphi\_{0}\right\}.\tag{8}$$

The noise term, *ηBB*[*n*], is AWGN with variance *σ*<sup>2</sup> *BB*. This variance depends on the filtering, down-conversion and sampling strategy applied by the receiver front-end and can be expressed as *σ*<sup>2</sup> *BB* = *N*0*BRX*, where *BRX* is the front-end bandwidth and *N*<sup>0</sup> is the Power Spectral Density (PSD) of the input noise *η*(*t*). The ratio between the carrier power, *C*, and the noise PSD, *N*0, defines the Carrier-to-Noise density power ratio (*C*/*N*0), one of the main signal quality indicators used in GNSS.

#### **2.2. The DVB-T interfering signal**

The interference term in (6), *i*(*t*), 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 Orthogonal Frequency Division Multiplexing (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}(t) = \frac{1}{\sqrt{M}} \sum\_{p=0}^{N\_d - 1} \sum\_{h=0}^{M-1} I\_{p,h} \exp\left(\frac{j2\pi ht}{M}\right),\tag{9}$$

where *M* is the modulation order, *h* is the subcarrier index, *p* is the symbol index, *Nd* represents the total number of transmitted symbols and *Ip*,*<sup>h</sup>* 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(t) = \sum\_{n=1}^{\infty} a\_n i\_{DVB-T}^n(t)\_\prime \tag{10}$$

where *an* 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* > 1 in (10) model the amplifier non-linearities and the ratios *an*/*a*<sup>1</sup> are expected to be small for *n* > 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

$$\dot{a}(t) = a\_3 i\_{DVB-T}^3(t). \tag{11}$$

The signal *i*(*t*) is filtered and down-converted by the receiver front-end and signal *iF*(*t*), the filtered version of *i*(*t*), will affect receiver operations.

Finally, the interference term in (7) can be modeled as *iBB*[*n*] = *iF*(*nTs*).

#### **2.3. The acquisition process**

8 MATLAB / Book 2

The signal *sb*(*t* − *iTh*) in (5) represents the subcarrier of duration *Th*, which determines the spectral characteristics of the transmitted GNSS signal. The GPS L1 Coarse/Acquisition (C/A) component is Bi-Phase Shift Keying (BPSK) modulated, whereas the Galileo E1 signal adopts a Composite Binary-Offset Carrier (CBOC) scheme. The sequence *cl*,*i*, of length *Nc*, defines the primary spreading code of the *l*th GNSS signal. In the following, only the BPSK case is

A GNSS receiver is able to process the *L* useful signals independently since the spreading

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*(*t*), 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

where the notation *x*[*n*] is used to denote discrete time sequences sampled at the frequency

filtering, down-conversion and sampling strategy applied by the receiver front-end and can

Spectral Density (PSD) of the input noise *η*(*t*). The ratio between the carrier power, *C*, and the noise PSD, *N*0, defines the Carrier-to-Noise density power ratio (*C*/*N*0), one of the main

The interference term in (6), *i*(*t*), 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 Orthogonal Frequency Division Multiplexing (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

*rBB*[*n*] = *yBB*(*nTs*) + *iBB*(*nTs*) + *ηBB*(*nTs*) = *yBB*[*n*] + *iBB*[*n*] + *ηBB*[*n*], (7)

. In addition, the index "BB" is used to denote a filtered signal down-converted to

*r*(*t*) = *y*(*t*) + *i*(*t*) + *η*(*t*), (6)

*Cd* (*nTs* − *τ*0) *c* (*nTs* − *τ*0) exp {*j*2*π f*0*nTs* + *ϕ*0} . (8)

*BB* = *N*0*BRX*, where *BRX* is the front-end bandwidth and *N*<sup>0</sup> is the Power

*BB*. This variance depends on the

considered. The results can be easily extended to different subcarriers.

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

baseband. Furthermore, the signal *yBB*[*n*] in (7) can be written as

*iDVB*−*T*(*t*) = <sup>1</sup>

<sup>√</sup>*<sup>M</sup>*

*Nd*−1 ∑ *p*=0

*M*−1 ∑ *h*=0

*Ip*,*<sup>h</sup>* exp

 *j*2*πht M*

, (9)

The noise term, *ηBB*[*n*], is AWGN with variance *σ*<sup>2</sup>

*yBB*[*n*] = <sup>√</sup>

signal quality indicators used in GNSS.

**2.2. The DVB-T interfering signal**

these operations, (6) becomes:

*fs* = <sup>1</sup> *Ts*

be expressed as *σ*<sup>2</sup>

can be represented as

codes are quasi-orthogonal. Therefore, expression (3) can be simplified to

After signal conditioning, the sequence *rBB*[*n*] is correlated with local replicas of the useful signal code and carrier as shown in Figure 4. Since the code delay, *τ*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 *rBB*[*n*], 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

$$S\_k = \frac{1}{N} \sum\_{n=kN}^{(k+1)N-1} r\_{BB}[n] \mathbf{c}\left(nT\_s - \tau\right) \exp\left\{-j2\pi f\_d n T\_s - j\varphi\right\},\tag{12}$$

where *τ*, *fd* and *ϕ* 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 *Tc* = *NTs* 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,

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

non-coherent integration can be implemented as illustrated in Figure 4. More specifically, the impact of the navigation message, *d*(·), is removed through squaring, |*Sk*| 2, and the final decision variable is computed as

$$D = \frac{1}{K} \sum\_{k=0}^{K-1} |\mathbb{S}\_k|^2 \tag{13}$$

*D*<sup>0</sup> *D*<sup>1</sup>

correctly declared False alarm

*H*<sup>1</sup> Missed detection Signal detection

**Table 1.** Confusion matrix describing the four events that can happen in the binary test performed by

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:

The probabilities of missed detection and correct signal absence decision are obtained as 1 − *Pd*(*β*) and 1 − *Pf a*(*β*), respectively. The performance of the acquisition process is characterized in terms of Receiver Operation Curves (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

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*<sup>0</sup> 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 *Sk* are independent and identically distributed (i.i.d.) complex Gaussian random variables with independent real and

> <sup>−</sup> *<sup>β</sup>* 2*σ*<sup>2</sup> *n*

> > √ *Kλ*; *βσ*<sup>2</sup> *n*

<sup>−</sup> *<sup>x</sup>*<sup>2</sup>+*a*<sup>2</sup> 2 

Q-function of order *K*. The function *IK*(·) is the modified Bessel function of first kind and

*<sup>λ</sup>* <sup>=</sup> <sup>|</sup><sup>E</sup> [*Sk*]<sup>|</sup>

*σ*2 *n*

 *K*−1 ∑ *i*=0

2

1 *i*!  *β* 2*σ*<sup>2</sup> *n*

*<sup>n</sup>* is the variance of the real and imaginary part of *Sk*. The parameter

*i*

, (17)

*IK*−1(*ax*)*dx* is the generalized Marcum

(16)

(18)

*Pd*(*β*) = Prob(*D* > *β*|*H*1) Probability of detection (14)

Semi-Analytic Techniques for Fast MATLAB Simulations 295

*Pf a*(*β*) = Prob(*D* > *β*|*H*0) Probability of false alarm. (15)

Signal absence

*H*<sup>0</sup>

and the four events described in Table 1 can occur.

of the different decision thresholds.

**2.4. The semi-analytic approach**

imaginary parts, it is possible to show [23] that

*<sup>b</sup> x x a K*−<sup>1</sup>

*Pf a*(*β*) = exp

*Pdet*(*β*) = Q*<sup>K</sup>*

exp

the acquisition process.

and

and

where *QK*(*a*; *<sup>b</sup>*) = <sup>+</sup><sup>∞</sup>

order *K*. In (16) and (17), *σ*<sup>2</sup>

*λ* is given by

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

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


The null hypothesis, *H*<sup>0</sup> assumes that the correlator outputs, *Sk*, are made of noise alone. Since the pseudo-random sequences, *c*(·), 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, *yBB*[*n*], then the useful signal component is almost completely filtered out at the correlator output. Thus, also in this case, the *H*<sup>0</sup> hypothesis is verified.

Furthermore, *H*<sup>1</sup> is the alternative hypothesis and assumes that the signal is present and the local code and carrier replica are perfectly aligned. If *H*<sup>1</sup> is declared, then rough estimates of *τ*<sup>0</sup> and *f*<sup>0</sup> are also obtained.

Depending on the result of the test, *D* > *β*, two decisions can be taken by the receiver



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

and the four events described in Table 1 can occur.

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(\beta) = \text{Prob}(D > \beta | H\_1) \quad \text{Probability of detection} \tag{14}$$

and

 

2, and the final

 

 

10 MATLAB / Book 2

2

 

· ( ) <sup>1</sup>

*Sk D*

*<sup>K</sup>* ∑ <sup>⋅</sup>

2, (13)

 

( ) <sup>1</sup> *<sup>N</sup>* ∑ <sup>⋅</sup>

**Figure 4.** Schematic representation of the operations performed by the acquisition block of a GNSS

the impact of the navigation message, *d*(·), is removed through squaring, |*Sk*|

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

• *H*1: the signal is present and the local code and carrier replica are aligned.

correlator output. Thus, also in this case, the *H*<sup>0</sup> hypothesis is verified.

*<sup>D</sup>* <sup>=</sup> <sup>1</sup> *K*

non-coherent integration can be implemented as illustrated in Figure 4. More specifically,

*K*−1 ∑ *k*=0 |*Sk*|

where *K* is the total number of correlation samples that are non-coherently integrated. It should be noted that for *K* = 1 only coherent integration is used. In order to determine the signal presence, the receiver compares *D* with a decision threshold, *β*. If *D* is greater than *β*

• *H*0: the signal is not present or it is not correctly aligned with the local code and carrier

The null hypothesis, *H*<sup>0</sup> assumes that the correlator outputs, *Sk*, are made of noise alone. Since the pseudo-random sequences, *c*(·), 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, *yBB*[*n*], then the useful signal component is almost completely filtered out at the

Furthermore, *H*<sup>1</sup> is the alternative hypothesis and assumes that the signal is present and the local code and carrier replica are perfectly aligned. If *H*<sup>1</sup> is declared, then rough estimates of

Depending on the result of the test, *D* > *β*, two decisions can be taken by the receiver

 

Local Code

Local Carrier

decision variable is computed as

then the useful signal is declared present.

[ ] *BB r n*

 

replica .

*τ*<sup>0</sup> and *f*<sup>0</sup> are also obtained.

• *D*0: the signal is declared not present • *D*1: the signal is declared present

receiver.

$$P\_{fa}(\beta) = \text{Prob}(D > \beta | H\_0) \quad \text{Probability of false alarm.} \tag{15}$$

The probabilities of missed detection and correct signal absence decision are obtained as 1 − *Pd*(*β*) and 1 − *Pf a*(*β*), respectively. The performance of the acquisition process is characterized in terms of Receiver Operation Curves (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*<sup>0</sup> 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 *Sk* are independent and identically distributed (i.i.d.) complex Gaussian random variables with independent real and imaginary parts, it is possible to show [23] that

$$P\_{fa}(\beta) = \exp\left\{-\frac{\beta}{2\sigma\_n^2}\right\} \sum\_{i=0}^{K-1} \frac{1}{i!} \left(\frac{\beta}{2\sigma\_n^2}\right)^i \tag{16}$$

and

$$P\_{\det}(\beta) = \mathbf{Q}\_{\mathcal{K}}\left(\sqrt{K\lambda}; \sqrt{\beta \sigma\_n^2}\right),\tag{17}$$

where *QK*(*a*; *<sup>b</sup>*) = <sup>+</sup><sup>∞</sup> *<sup>b</sup> x x a K*−<sup>1</sup> exp <sup>−</sup> *<sup>x</sup>*<sup>2</sup>+*a*<sup>2</sup> 2 *IK*−1(*ax*)*dx* is the generalized Marcum Q-function of order *K*. The function *IK*(·) is the modified Bessel function of first kind and order *K*. In (16) and (17), *σ*<sup>2</sup> *<sup>n</sup>* is the variance of the real and imaginary part of *Sk*. The parameter *λ* is given by

$$\lambda = \frac{|\mathbb{E}\left[\mathbb{S}\_k\right]|^2}{\sigma\_n^2} \tag{18}$$

12 MATLAB / Book 2 296 MATLAB – A Fundamental Tool for Scienti c Computing and Engineering Applications – Volume 2 Semi-Analytic Techniques for Fast MATLAB Simulations <sup>13</sup>

**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.

and defines the Signal-to-Noise Ratio (SNR) at the correlator outputs.

The correlator output can be considered i.i.d. 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 *Sk* 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 *σ*<sup>2</sup> *<sup>n</sup>* and *λ*.

The analytical knowledge of the system can be further exploited to simplify the evaluation of *σ*2 *<sup>n</sup>* and *λ*. In particular, since *iBB*[*n*] and *ηBB*[*n*] 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|\mathbb{E}\left[\mathbb{S}\_{k}\right]\right|^{2} = \mathbb{C}\_{\prime} \tag{19}$$

 

ROC

highlighted in different colors.

determine the system ROC.

and

as [13, 24]

 

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

> <sup>=</sup> *Ci* 2*Tc*

where *Ci* is the interference power and *ka* is the Spectral Separation Coefficient (SSC) defined

The function *Gi*(*f*) in (24) is the normalized PSD of the DVB-T interference signal after

The function *Gc*(*f*) models the effect of the correlation on the interfering signal. Correlation can be modeled as an additional filtering stage and *Gc*(*f*) can be shown to be well approximated by the PSD of the subcarrier used in the despreading process. Also, *Gc*(*f*) is normalized to have a unit integral. It is noted that different subcarriers lead to different *Gc*(*f*) and thus, *iBB*[*n*] 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 *iBB*[*n*] is obtained. The samples of *iBB*[*n*] are used to estimate the normalized PSD, *Gi*(*f*). This can be easily obtained using the MATLAB functions developed for spectral analysis. In this case, the pwelch function is used. The function *Gi*(*f*) is used to compute the SSC, which is then used to

 *BRX* /2 −*BRX* /2

 *BRX* /2 −*BRX* /2

Simulation block

Analytical model -

1 2 Var *Sη*,*<sup>k</sup>* <sup>=</sup> *<sup>N</sup>*<sup>0</sup> 2*Tc*

1 2 Var *Si*,*<sup>k</sup>*

*ka* =

front-end filtering. In addition, *Gi*(*f*) is normalized such that

 

Semi-Analytic Techniques for Fast MATLAB Simulations 297

(22)


*ka*, (23)

*Gi*(*f*)*Gc*(*f*)*d f* . (24)

*Gi*(*f*)*d f* = 1. (25)

 

where *C* is the useful signal power and is one of the known parameters of the system. Thus, *λ* can be derived from *C* and *σ*<sup>2</sup> *n*.

Finally, exploiting the linearity of the correlation process, it is possible to express *Sk* as

$$S\_k = S\_{r,k} + S\_{\eta,k} + S\_{i,k} \tag{20}$$

which is a linear combination of a useful signal term, derived from *yBB*[*n*], a noise term, derived from *ηBB*[*n*], and an interference term derived from *iBB*[*n*]. The variance *σ*<sup>2</sup> *<sup>n</sup>* can be obtained as

$$
\sigma\_n^2 = \frac{1}{2} \text{Var}\left\{ \mathbf{S}\_k \right\} = \frac{1}{2} \text{Var}\left\{ \mathbf{S}\_{\eta,k} \right\} + \frac{1}{2} \text{Var}\left\{ \mathbf{S}\_{i,k} \right\}.\tag{21}
$$

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

**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.

$$\frac{1}{2}\text{Var}\left\{\mathbf{S}\_{\eta,k}\right\} = \frac{N\_0}{2T\_c} \tag{22}$$

and

12 MATLAB / Book 2

 

**Figure 5.** Schematic representation of the full Monte Carlo simulation system for the ROC evaluation of

The correlator output can be considered i.i.d. 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 *Sk* 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

The analytical knowledge of the system can be further exploited to simplify the evaluation of

*<sup>n</sup>* and *λ*. In particular, since *iBB*[*n*] and *ηBB*[*n*] are modeled as zero mean random processes, they only contribute to the variance of the correlator outputs. Thus, neglecting residual errors

where *C* is the useful signal power and is one of the known parameters of the system. Thus,

which is a linear combination of a useful signal term, derived from *yBB*[*n*], a noise term, derived from *ηBB*[*n*], and an interference term derived from *iBB*[*n*]. The variance *σ*<sup>2</sup>


Finally, exploiting the linearity of the correlation process, it is possible to express *Sk* as

2 Var *Sη*,*<sup>k</sup>* + 1 2 Var *Si*,*<sup>k</sup>* 

Var {*Sk*} <sup>=</sup> <sup>1</sup>

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

the acquisition of a GPS L1 signal in the presence of DVB-T third order harmonics.

and defines the Signal-to-Noise Ratio (SNR) at the correlator outputs.

*<sup>n</sup>* and *λ*.

due to delay and frequency partial misalignments, yields

*σ*2 *<sup>n</sup>* <sup>=</sup> <sup>1</sup> 2 *n*.

[ ] *BB r n Dk r n*[ ]

 

<sup>2</sup> = *C*, (19)

*Sk* = *Sr*,*<sup>k</sup>* + *Sη*,*<sup>k</sup>* + *Si*,*k*, (20)

ROC

*<sup>n</sup>* can be

. (21)

 

*i n*[ ]

 

 

  *y n*[ ]

η[ ] *n*

that need to be estimated are *σ*<sup>2</sup>

*λ* can be derived from *C* and *σ*<sup>2</sup>

*σ*2

obtained as


$$\frac{1}{2}\text{Var}\left\{\mathbf{S}\_{i,k}\right\} = \frac{\mathbf{C}\_i}{2T\_c}\mathbf{k}\_{a\nu} \tag{23}$$

where *Ci* is the interference power and *ka* is the Spectral Separation Coefficient (SSC) defined as [13, 24]

$$k\_d = \int\_{-B\_{\rm RX}/2}^{B\_{\rm RX}/2} \mathcal{G}\_i(f)\mathcal{G}\_{\mathcal{C}}(f) df. \tag{24}$$

The function *Gi*(*f*) in (24) is the normalized PSD of the DVB-T interference signal after front-end filtering. In addition, *Gi*(*f*) is normalized such that

$$\int\_{-\mathcal{B}\ge 2}^{\mathcal{B}\_{\text{RX}}/2} G\_l(f) df = 1. \tag{25}$$

The function *Gc*(*f*) models the effect of the correlation on the interfering signal. Correlation can be modeled as an additional filtering stage and *Gc*(*f*) can be shown to be well approximated by the PSD of the subcarrier used in the despreading process. Also, *Gc*(*f*) is normalized to have a unit integral. It is noted that different subcarriers lead to different *Gc*(*f*) and thus, *iBB*[*n*] 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 *iBB*[*n*] is obtained. The samples of *iBB*[*n*] are used to estimate the normalized PSD, *Gi*(*f*). This can be easily obtained using the MATLAB functions developed for spectral analysis. In this case, the pwelch function is used. The function *Gi*(*f*) is used to compute the SSC, which is then used to determine the system ROC.

**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.

<sup>10</sup>-6 <sup>10</sup>-5 <sup>10</sup>-4 <sup>10</sup>-3 <sup>10</sup>-2 <sup>10</sup>-1 <sup>10</sup><sup>0</sup> <sup>10</sup>-3

Coherent integration time, *Tc* 1 ms

Centre frequency difference, Δ*f* 0 Hz Receiver bandwidth, *BRX* 8 MHz

Number of Monte Carlo Simulation runs 10<sup>6</sup>

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

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

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

Interference to signal power ratio, *Ci*

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

the semi-analytic approach described above.

K = 1

Pfa

Parameter Value *C*/*N*<sup>0</sup> 35 dB-Hz

*<sup>C</sup>* 30 dB

Semi-Analytic Monte Carlo

Semi-Analytic Techniques for Fast MATLAB Simulations 299

10-2

computational complexity.

Pd

10-1

10<sup>0</sup>

K = 10

K = 5

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

#### **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 MATLAB 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 MATLAB. 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 *fI* = *fRF* + Δ*f* , where Δ*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 Δ*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

14 MATLAB / Book 2

GPS L1 Centre Frequency

<sup>1560</sup> <sup>1565</sup> <sup>1570</sup> <sup>1575</sup> <sup>1580</sup> -120

**Figure 7.** Representation of a normalized PSD realization of the third harmonic of the DVB-T signal and

The developed semi-analytic approach is shown in Figure 6 where the simulation, estimation

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 MATLAB 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 MATLAB. 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 *fI* = *fRF* + Δ*f* , where Δ*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 Δ*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

Frequency [MHz]


the frequency response of the GPS L1 front-end filter.

and analytic components are clearly highlighted.

**2.5. Performance comparison**



Normalized PSD



0

DVB-T Harmonics Centre Frequency

Front-end Filter

DVB-T 3rd Harmonics PSD

**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.


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

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.

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 Semi-Analytic Tracking Loop Simulations (SATLSim) toolbox is a set of MATLAB 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:

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 Numerically Controlled Oscillator

Efficient tracking loop simulations can be obtained by substituting the I&D blocks with their

• Δ*τ<sup>d</sup>* and Δ*τ<sup>s</sup>* are the code and subcarrier delay errors. The delay Δ*τ<sup>s</sup>* is present only when

• *Tc* = *NTs* is the coherent integration time where *N* is the number of samples used to

• *Rl* (Δ*τd*, Δ*τs*) is the correlation function between incoming and locally generated code and is a function of both code and subcarrier delay errors. When the SLL is not used,

• *η<sup>c</sup>* is a zero-mean noise term whose variance depends on the input noise power, front-end filtering and the correlation process operated by the I&D blocks. More details on the

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

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, Δ*fd*, Δ*ϕ*, Δ*τ<sup>d</sup>* and Δ*τ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 *η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 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

*Rl* (Δ*τd*, Δ*τs*) exp {*j*Δ*ϕ*} + *ηc*, (26)

Semi-Analytic Techniques for Fast MATLAB Simulations 301

(NCO) that is used for generating new local signal replicas.

• Δ*fd* and Δ*ϕ* are the residual frequency and phase errors;

a SLL is used to correctly align the signal subcarrier [27];

*Rl* (Δ*τd*, Δ*τs*) is replaced by the standard code correlation function;

by the simulation scenario) and estimates produced by the NCO.

√

compute a single correlator;

properties of *η<sup>c</sup>* can be found in [9].

digital tracking loops.

where

analytical model. More specifically, a correlator output can be modeled as:

*<sup>C</sup>*sin (*π*Δ*fdTc*) *π*Δ*fdTc*


In the following, the closed-loop approach for the simulation of digital tracking loops is considered and the MATLAB 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 *τ*<sup>0</sup> 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 Delay Lock Loop (DLL) is usually used to provide accurate estimates of the code delay, *τ*0, and track delay variations due to the relative motion between receiver and satellite. The Doppler frequency, *f*0, is recovered using either a Frequency Lock Loop (FLL) or a Phase Lock Loop (PLL). If a PLL is used then the carrier phase, *ϕ*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 Subcarrier Lock Loop (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 I&D 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 Numerically Controlled Oscillator (NCO) that is used for generating new local signal replicas.

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

$$\sqrt{\mathcal{C}} \frac{\sin \left( \pi \Delta f\_d T\_c \right)}{\pi \Delta f\_d T\_c} \mathcal{R}\_l \left( \Delta \tau\_d, \Delta \tau\_s \right) \exp \left\{ j \Delta \varphi \right\} + \eta\_{c\prime} \tag{26}$$

where

16 MATLAB / Book 2

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 Semi-Analytic Tracking Loop Simulations (SATLSim) toolbox is a set of MATLAB 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

In the following, the closed-loop approach for the simulation of digital tracking loops is considered and the MATLAB 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

If the signal is successfully acquired then different tracking loops are used to refine the estimate of the signal parameters. A Delay Lock Loop (DLL) is usually used to provide accurate estimates of the code delay, *τ*0, and track delay variations due to the relative motion between receiver and satellite. The Doppler frequency, *f*0, is recovered using either a Frequency Lock Loop (FLL) or a Phase Lock Loop (PLL). If a PLL is used then the carrier phase, *ϕ*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 Subcarrier Lock Loop (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 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 I&D blocks.

**3. Closed-loop approach: Digital tracking loops**

• http://www.ngs.noaa.gov/gps-toolbox/SATLSim.htm • http://plan.geomatics.ucalgary.ca/publications.php.

In the following, the joint simulation of DLL and SLL is considered.

the following websites:

code delay *τ*<sup>0</sup> and Doppler frequency *f*0.


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.

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, Δ*fd*, Δ*ϕ*, Δ*τ<sup>d</sup>* and Δ*τ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 *η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

18 MATLAB / Book 2 302 MATLAB – A Fundamental Tool for Scienti c Computing and Engineering Applications – Volume 2 Semi-Analytic Techniques for Fast MATLAB Simulations <sup>19</sup>

**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.

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:

$$
\phi\_k = \phi\_{k-1} + \frac{T\_c}{2} \left( \delta \phi\_{k-1} + \delta \phi\_{k-2} \right),
\tag{27}
$$

**Initialization:**

**Main Simulation Loop**

• Initial parameters: sampling frequency, integration time…

• Noise generation (GenerateNoiseVector.m)

6. Loop Filter update (UpdateFilter.m)

**Figure 10.** Structure of the SATLSim toolbox and list of the different MATLAB functions.

3. Error-to-Signal conversion (GenerateSignalCorrelation.m)

5. Discriminator update (UpdateDiscriminator.m)

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 *C*/*N*0. 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

> *<sup>τ</sup>*<sup>ˆ</sup> <sup>±</sup> <sup>1</sup> 2

where *τ*ˆ is the best code delay estimate and *ds* 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 *τ*ˆ 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 *C*/*N*<sup>0</sup> is used to determine the

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

*ds* (28)

Semi-Analytic Techniques for Fast MATLAB Simulations 303

1. NCO update (UpdateNCO.m) 2. Evaluation of the estimation error

4. Signal and Noise combining

**End Loop on the simulation runs**

correlator amplitude and the variance of the noise component, *ηc*.

• Evaluate tracking jitter **End Loop on the C/N0 values**

**End Loop on the Early-minus-Late spacing**

• Plot tracking results

Late correlators, computed for the delays

**3.1. Code structure**

• Loop filter design (FilterDesign.m)

• Input (true) parameters generation

**Loop on the simulation runs:**

**Loop on the Early-minus-Late spacing: Loop on the C/N0 values:**

where *ϕ*ˆ *<sup>k</sup>* denotes the *k*-th estimate of the tracking parameter under consideration and *δϕ*ˆ *<sup>k</sup>* is its estimated rate of change. The rate *δϕ*ˆ *<sup>k</sup>* 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 I&D 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, *δϕ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.

#### **Initialization:**

18 MATLAB / Book 2

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

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

where *ϕ*ˆ *<sup>k</sup>* denotes the *k*-th estimate of the tracking parameter under consideration and *δϕ*ˆ *<sup>k</sup>* is its estimated rate of change. The rate *δϕ*ˆ *<sup>k</sup>* 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 I&D 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

The operations required to convert the correlator outputs into a new estimate of the parameter rate, *δϕ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

*Tc*

 


<sup>2</sup> (*δϕ*<sup>ˆ</sup> *<sup>k</sup>*−<sup>1</sup> <sup>+</sup> *δϕ*<sup>ˆ</sup> *<sup>k</sup>*−2), (27)

 

Simulation block

 

 

feedback NCO [28], characterized by the following update equation:

*<sup>ϕ</sup>*<sup>ˆ</sup> *<sup>k</sup>* = *<sup>ϕ</sup>*<sup>ˆ</sup> *<sup>k</sup>*−<sup>1</sup> +

filter that is used to produce new estimates of the tracking parameters.

Analytical model

 


 

SATLSim toolbox.

using the approach described in [9].

 

 

 


#### **Main Simulation Loop**

#### **Loop on the Early-minus-Late spacing:**


**End Loop on the Early-minus-Late spacing**

**Figure 10.** Structure of the SATLSim toolbox and list of the different MATLAB functions.

### **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 *C*/*N*0. 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

$$
\hat{\mathbf{r}} \pm \frac{1}{2} d\_s \tag{28}
$$

where *τ*ˆ is the best code delay estimate and *ds* 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 *τ*ˆ 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 *C*/*N*<sup>0</sup> is used to determine the correlator amplitude and the variance of the noise component, *ηc*.

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 *C*/*N*<sup>0</sup> 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 *C*/*N*0, 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 I&D 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 = \frac{1}{\mathcal{C}/N\_0 T\_c}.\tag{29}
$$

*Cn* =

⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣

*Rl* � *ds* <sup>2</sup> , *dsc* 2 �

*Rl* � *ds* <sup>2</sup> , 0�

*Rl* � *ds* <sup>2</sup> , *dsc* 2 �

of different tracking algorithms [9].

a) Tracking jitter of the DLL alone b) Tracking jitter of the SLL alone

speed of light.

c) Jitter of the combined delay estimate.

**3.4. Sample results**

1 *Rl*

*Rl* (*ds*, 0) *Rl*

where *dsc* is the subcarrier Early-minus-Late spacing.

the rate of change of both code and subcarrier delay.

are shown in Figure 11. The figure is divided into three parts:

� *ds* <sup>2</sup> , *dsc* 2 � *Rl* � *ds* <sup>2</sup> , 0� *Rl* � *ds* <sup>2</sup> , *dsc* 2 �

� *ds* <sup>2</sup> , *dsc* 2 � *Rl* � *ds* <sup>2</sup> , 0� *Rl* � *ds* <sup>2</sup> , *dsc* 2 �

*Rl* � 0, *dsc* 2 �

1 *Rl*

*Rl* (0, *dsc*) *Rl*

� 0, *dsc* 2 �

� 0, *dsc* 2 �

The NCO update (UpdateNCO) is performed on both code and subcarrier loops and the estimated errors, Δ*τ<sup>d</sup>* and Δ*τ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 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

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

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

In addition to this, the curves are shown in Figure 11a) and Figure 11b). More specifically, three different methodologies have been employed for determining the tracking jitter.

1 *Rl*

*Rl* (*ds*, 0)

Semi-Analytic Techniques for Fast MATLAB Simulations 305

⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦

, (31)

� *ds* <sup>2</sup> , *dsc* 2 �

� *ds* <sup>2</sup> , *dsc* 2 �

1

*Rl* (0, *dsc*) *Rl*

1 *Rl*

� 0, *dsc* 2 � *Rl* � *ds* <sup>2</sup> , 0�

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

$$\sqrt{\mathbb{C}} \frac{\sin \left( \pi \Delta f\_d T\_\varepsilon \right)}{\pi \Delta f\_d T\_\varepsilon} \exp \left\{ j \Delta \varphi \right\} + \eta\_\varepsilon \tag{30}$$

where Δ*fd* is obtained by comparing the true Doppler frequency against the loop filter output. Δ*ϕ* 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 Double Estimator (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 × *Nsim* matrix, where *Nsim* 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

$$\mathbf{C}\_{n} = \begin{bmatrix} 1 & R\_{l}\left(\frac{d\_{s}}{2}, \frac{d\_{s}}{2}\right) & R\_{l}\left(\frac{d\_{s}}{2}, 0\right) & R\_{l}\left(\frac{d\_{s}}{2}, \frac{d\_{s}}{2}\right) & R\_{l}\left(d\_{s}, 0\right) \\ R\_{l}\left(\frac{d\_{s}}{2}, \frac{d\_{s}}{2}\right) & 1 & R\_{l}\left(0, \frac{d\_{s}}{2}\right) & R\_{l}\left(0, d\_{s}\right) & R\_{l}\left(\frac{d\_{s}}{2}, \frac{d\_{s}}{2}\right) \\ R\_{l}\left(\frac{d\_{s}}{2}, 0\right) & R\_{l}\left(0, \frac{d\_{s}}{2}\right) & 1 & R\_{l}\left(0, \frac{d\_{s}}{2}\right) & R\_{l}\left(\frac{d\_{s}}{2}, 0\right) \\ R\_{l}\left(\frac{d\_{s}}{2}, \frac{d\_{s}}{2}\right) & R\_{l}\left(0, d\_{s}\right) & R\_{l}\left(0, \frac{d\_{s}}{2}\right) & 1 & R\_{l}\left(\frac{d\_{s}}{2}, \frac{d\_{s}}{2}\right) \\ R\_{l}\left(d\_{s}, 0\right) & R\_{l}\left(\frac{d\_{s}}{2}, \frac{d\_{s}}{2}\right) & R\_{l}\left(\frac{d\_{s}}{2}, 0\right) & R\_{l}\left(\frac{d\_{s}}{2}, \frac{d\_{s}}{2}\right) & 1 \end{bmatrix},\tag{31}$$

where *dsc* is the subcarrier Early-minus-Late spacing.

The NCO update (UpdateNCO) is performed on both code and subcarrier loops and the estimated errors, Δ*τ<sup>d</sup>* and Δ*τ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**

20 MATLAB / Book 2

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 *C*/*N*<sup>0</sup> 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 *C*/*N*0, 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

In the code provided, theoretical formulae for the computation of the tracking jitter are also

The simulation of a standard PLL requires the generation of the Prompt correlator alone (GenerateSignalCorrelation). The Prompt correlator is the output of the I&D 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

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

where Δ*fd* is obtained by comparing the true Doppler frequency against the loop filter output. Δ*ϕ* is the phase error obtained as the difference between the true phase and the phase estimate

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

In the Double Estimator (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 × *Nsim* matrix, where *Nsim* 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

. (29)

exp {*j*Δ*ϕ*} + *ηc*, (30)

*σ*2 *<sup>n</sup>* <sup>=</sup> <sup>1</sup> *C*/*N*0*Tc*

*<sup>C</sup>*sin (*π*Δ*fdTc*) *π*Δ*fdTc*

√

**3.3. Double estimator (DoubleEstimator.m)**

implemented and used as a comparison term for the simulation results.

quantities of interest such as the tracking jitter.

real and imaginary parts with variance [9]

**3.2. Standard PLL (PLL.m)**

simplifies to

produced by the NCO.

by changing this function.

correlation matrix

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:


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.

In addition to this, the curves are shown in Figure 11a) and Figure 11b). More specifically, three different methodologies have been employed for determining the tracking jitter.

Parameter Value Sampling Frequency *fs* = 8 MHz Integration Time *Tc* = 4 ms DLL Early-minus-Late spacing 0.1955 *μ*s (0.2 chips) SLL Early-minus-Late spacing 0.1955 *μ*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)

Semi-Analytic Techniques for Fast MATLAB Simulations 307

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

latest curve is not available for the combined delay estimate.

loop requiring only limited computation resources.

in Table 3.

**4. Conclusions**

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 *C*/*N*0. 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

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

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 *C*/*N*<sup>0</sup> values, a good agreement between theoretical and simulation results is found. However, for *C*/*N*<sup>0</sup> 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 *C*/*N*<sup>0</sup> 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 MATLAB toolbox is able to effectively describe the non-linear behavior of the

In this chapter, the development of fast semi-analytic techniques using MATLAB 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

**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.


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

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 *C*/*N*0. 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.

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 *C*/*N*<sup>0</sup> values, a good agreement between theoretical and simulation results is found. However, for *C*/*N*<sup>0</sup> 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 *C*/*N*<sup>0</sup> 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 MATLAB toolbox is able to effectively describe the non-linear behavior of the loop requiring only limited computation resources.

## **4. Conclusions**

22 MATLAB / Book 2

a) Non-coherent discriminator - DLL

Theory

5

5

5

Tracking jitter [m]

10

Tracking jitter [m]

10

Tracking jitter [m]

10

<sup>20</sup> <sup>25</sup> <sup>30</sup> <sup>35</sup> <sup>40</sup> <sup>0</sup>

b) Non-coherent discriminator - SLL

Theory

<sup>20</sup> <sup>25</sup> <sup>30</sup> <sup>35</sup> <sup>40</sup> <sup>0</sup>

c) Non-coherent discriminator - Combined delay

<sup>20</sup> <sup>25</sup> <sup>30</sup> <sup>35</sup> <sup>40</sup> <sup>0</sup>

[dB-Hz]

C/N0

**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.

[dB-Hz]

Theory

C/N0

[dB-Hz]

Estimated from the actual error Estimated from the loop filter output

Estimated from the actual error Estimated from the loop filter output

Estimated from the actual error

C/N0

In this chapter, the development of fast semi-analytic techniques using MATLAB 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

#### 24 MATLAB / Book 2 308 MATLAB – A Fundamental Tool for Scienti c Computing and Engineering Applications – Volume 2 Semi-Analytic Techniques for Fast MATLAB Simulations <sup>25</sup>

would require. In this way, the strengths of both analytical and Monte Carlo methods are effectively combined.

[8] J. S. Silva, P. F. Silva, A. Fernandez, J. Diez, and J. F. M. Lorga. Factored correlator model: A solution for fast, flexible, and realistic GNSS receiver simulations. In *Proc.*

Semi-Analytic Techniques for Fast MATLAB Simulations 309

[9] D. Borio, P. B. Anantharamu, and G. Lachapelle. Semi-analytic simulations: An extension to unambiguous BOC tracking. In *Proc. of the ION/ITM (International Technical Meeting)*,

[10] Daniele Borio, Pratibha Anantharamu, and Gérard Lachapelle. SATLSim: a semi-analytic framework for fast GNSS tracking loop simulations. *GPS Solutions*, 15:427–431, 2011. [11] J.M. Geist. Computer generation of correlated gaussian random variables. *Proceedings of*

[12] A. R. Golshan. Loss of lock analysis of a firstorder digital code tracking loop and comparison of results to analog loop theory for BOC and NRZ signals. In *Proc. of the ION/NTM (National Technical Meeting)*, pages 299 – 305, San Diego, CA, January 2005. [13] J. W. Betz. Effect of narrowband interference on GPS code tracking accuracy. In *Proc. of*

[14] D. Borio. GNSS acquisition in the presence of continuous wave interference. *IEEE*

[15] T. Van Slyke, W. Kuhn, and B. Natarajan. Measuring interference from a UWB transmitter in the GPS l1 band. In *Proc. of the IEEE Radio and Wireless Symposium*, pages 887 – 890,

[16] A. Simsky, T. De Wilde, D. Mertens, E. Koitsaly, and J.-M. Sleewaegen. First field experience with L5 signals: DME interference reality check. In *Proc. of ION/GNSS*, pages

[17] D. Borio, S. Savasta, and L. Lo Presti. On the DVB-t coexistence with galileo and GPS system. In *Proc. of the* 3*rd ESA Workshop on Satellite Navigation User Equipment Technologies*

[18] M. Wildemeersch, A. Rabbachin, E. Cano, and J. Fortuny. Interference assessment of DVB-t within the GPS l1 and galileo e1 band. In *Proc. of the* 5*th ESA European Workshop on GNSS Signals and Signal Processing (NAVITEC)*, pages 1 – 8, Noordwijk, The Netherlands,

[19] Steven M. Kay. *Fundamentals of Statistical Signal Processing, Volume 2: Detection Theory*,

[20] ETSI. Digital video broadcasting (DVB); framing structure, channel coding and

[21] E. D. Kaplan and C. Hegarty, editors. *Understanding GPS: Principles and Applications*.

[22] H. L. Van Trees. *Detection, Estimation, and Modulation Theory - Part 1*. Wiley-Interscience,

[23] Daniele Borio. *A statistical theory for GNSS signal acquisition*. Phd thesis, Politecnico di

[24] John W. Betz. Effect of partial-band interference on receiver estimation of *C*/*N*0. In *Proc.*

[25] J. Fortuny-Guasch, M. Wildemeersch, and D. Borio. Assessment of DVB-T impact on GNSS acquisition and tracking performance. In *Proc. of the ION/GNSS*, pages 347–356,

[26] P. Anantharamu, D. Borio, and G. Lachapelle. Sub-carrier shaping for BOC modulated GNSS signals. *EURASIP Journal on Advances in Signal Processing*, 2011(1):133, 2011.

*Transactions on Aerospace and Electronic Systems*, 46(1):47 – 60, January 2010.

*of ION/GNSS*, pages 2676 – 2686, Forth Worth, TX, September 2007.

pages 1023–1036, San Diego, CA, January 2010.

*ION/NTM*, pages 16–27, Anaheim, CA, January 2000.

*(NAVITEC)*, Noordwijk, The Netherlands, December 2006.

modulation for digital terrestrial television, 2006. EN 300 744.

*of the ION/NTM*, pages 817 – 828, Long Beach, CA, January 2001.

volume 2. Prentice Hall, 1rt edition, February 1998.

Artech House, 2nd edition, November 2005.

1st edition, September 2001.

San Diego, CA, January 2011.

Torino, April 2008.

*the IEEE*, 67(5):862 – 863, may 1979.

Orlando, FL, March 2008.

December 2010.

29 – 37, Savannah, GA, September 2009.

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 MATLAB toolbox. This toolbox has been developed for the analysis of digital tracking loops and fully exploits the flexibility of the MATLAB programming language. The code has been organized in functions that can be easily replaced by different MATLAB 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 MATLAB language are an effective tool for the analysis of complex communications systems.

## **Author details**

Daniele Borio and Eduardo Cano *EC Joint Research Centre, Institute for the Protection and Security of the Citizen, Italy*

## **5. References**


[8] J. S. Silva, P. F. Silva, A. Fernandez, J. Diez, and J. F. M. Lorga. Factored correlator model: A solution for fast, flexible, and realistic GNSS receiver simulations. In *Proc. of ION/GNSS*, pages 2676 – 2686, Forth Worth, TX, September 2007.

24 MATLAB / Book 2

would require. In this way, the strengths of both analytical and Monte Carlo methods are

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

The second example considered the closed-loop approach and specific focus was devoted to the SATLSim MATLAB toolbox. This toolbox has been developed for the analysis of digital tracking loops and fully exploits the flexibility of the MATLAB programming language. The code has been organized in functions that can be easily replaced by different MATLAB 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 MATLAB language are an effective tool for

[1] W. H. Tranter, K. S. Shanmugan, T. S. Rappaport, and K. L. Kosbar. *Principles of Communication Systems Simulation with Wireless Applications*. Prentice Hall, January 2004. [2] M. C. Jeruchim, P. Balaban, and K. S. Shanmugan. *Simulation of communication systems*.

[3] F.M. Gardner and J. D. Baker. *Simulation Techniques: Models of Communication Signals and*

[4] M. Jeruchim. Techniques for estimating the bit error rate in the simulation of digital communication systems. *IEEE Journal on Selected Areas in Communications*, 2(1):153 – 170,

[5] M. Pent, L. Lo Presti, G. D'Aria, and G. De Luca. Semianalytic BER evaluation by simulation for noisy nonlinear bandpass channels. *IEEE Journal on Selected Areas in*

[6] M.T. Core, R. Campbell, P. Quan, and J. Wada. Semianalytic BER for PSK. *IEEE*

[7] A. R. Golshan. Post-correlator modeling for fast simulation and joint performance analysis of GNSS code and carrier tracking loops. In *Proc. of the ION/NTM (National*

*EC Joint Research Centre, Institute for the Protection and Security of the Citizen, Italy*

Kluwer Academic/Plenum Publishers, new york edition, 2000.

*Transactions on Wireless Communications*, 8(4):1644 –1648, April 2009.

*Technical Meeting)*, pages 312 – 318, Monterey, CA, January 2006.

effectively combined.

alleviating those requirements.

Daniele Borio and Eduardo Cano

**Author details**

**5. References**

January 1984.

the analysis of complex communications systems.

*Processes*. John Wiley & Sons, 2003.

*Communications*, 6(1):34 –41, January 1988.

	- [27] M. S. Hodgart and P. D. Blunt. A dual estimate receiver of binary offset carrier (BOC) modulated signals global navigation satellite systems. *Electronics Letters*, 43(16):877–878, August 2007.
	- [28] S. A. Stephens and J. B. Thomas. Controlled-root formulation for digital phase-locked loops. *IEEE Transactions on Aerospace and Electronic Systems*, 31(1):78 –95, january 1995.
	- [29] A. J. V. Dierendonck, P. Fenton, and T. Ford. Theory and performance of narrow correlator spacing in a GPS receiver. *NAVIGATION: the Journal of The Institut of Navigation*, 39(3):265 – 283, Fall 1992.

26 MATLAB / Book 2

310 MATLAB – A Fundamental Tool for Scienti c Computing and Engineering Applications – Volume 2

[27] M. S. Hodgart and P. D. Blunt. A dual estimate receiver of binary offset carrier (BOC) modulated signals global navigation satellite systems. *Electronics Letters*, 43(16):877–878,

[28] S. A. Stephens and J. B. Thomas. Controlled-root formulation for digital phase-locked loops. *IEEE Transactions on Aerospace and Electronic Systems*, 31(1):78 –95, january 1995. [29] A. J. V. Dierendonck, P. Fenton, and T. Ford. Theory and performance of narrow correlator spacing in a GPS receiver. *NAVIGATION: the Journal of The Institut of Navigation*,

August 2007.

39(3):265 – 283, Fall 1992.

## *Edited by Vasilios N. Katsikis*

This excellent book represents the second part of three-volumes regarding MATLABbased applications in almost every branch of science. The present textbook contains a collection of 13 exceptional articles. In particular, the book consists of three sections, the first one is devoted to electronic engineering and computer science, the second is devoted to MATLAB/SIMULINK as a tool for engineering applications, the third one is about Telecommunication and communication systems and the last one discusses MATLAB toolboxes.

MATLAB - A Fundamental Tool for Scientific Computing and Engineering Applications -

Volume 2

MATLAB

A Fundamental Tool for Scientific Computing

and Engineering Applications - Volume 2

*Edited by Vasilios N. Katsikis*

Photo by monsitj / iStock