Digital Beamforming for a Secondary Surveillance Radar Numerical and statistical analysis of the ADS-B reception performance for an SSR system with conventional and adaptive beamforming Master’s thesis in Complex Adaptive Systems MATTIAS WIKLUND DEPARTMENT OF MATHEMATICAL SCIENCES CHALMERS UNIVERSITY OF TECHNOLOGY Gothenburg, Sweden 2024 www.chalmers.se www.chalmers.se Master’s thesis 2024 Digital Beamforming for a Secondary Surveillance Radar Numerical and statistical analysis of the ADS-B reception performance for an SSR system with conventional and adaptive beamforming MATTIAS WIKLUND Department of Mathematical Sciences Division of Applied Mathematics and Statistics Chalmers University of Technology Gothenburg, Sweden 2024 Digital Beamforming for a Secondary Surveillance Radar Numerical and statistical analysis of the ADS-B reception performance for an SSR system with conventional and adaptive beamforming MATTIAS WIKLUND © MATTIAS WIKLUND, 2024. Supervisors: Ida-Maria Svensson, Lars Sundell & Juan Cabello Sánchez, Saab AB Examiner: Adam Andersson, Department of Mathematical Sciences Master’s Thesis 2024 Department of Mathematical Sciences Division of Applied Mathematics and Statistics Chalmers University of Technology SE-412 96 Gothenburg Telephone +46 31 772 1000 Cover: Beam pattern in polar coordinates for adaptive linear constraint minimum power (LCMP) beamforming. A signal of interest is marked with a solid line and an interfering signal with a dashed line. Typeset in LATEX, template by Kyriaki Antoniadou-Plytaria Printed by Chalmers Reproservice Gothenburg, Sweden 2024 iv Digital Beamforming for a Secondary Surveillance Radar Numerical and statistical analysis of the ADS-B reception performance for an SSR system with conventional and adaptive beamforming MATTIAS WIKLUND Department of Mathematical Sciences Chalmers University of Technology Abstract A secondary surveillance radar (SSR) is a surveillance radar system which, rather than detecting reflected signals, relies on interrogations and replies to identify tar- gets. In this master’s thesis, the performance of an SSR system is analyzed. The performance is measured by estimating the probability to detect Automatic De- pendent Surveillance-Broadcast (ADS-B) signals interfered by other signals on the 1090 MHz frequency band. The study is performed using different beamforming techniques, with the aim of enhancing the reception in some directions while re- ducing it in others. The probability of detection is estimated by performing Monte Carlo simulations for different target positions and interference scenarios. Simula- tions are performed for a system using digital beamforming, individually digitizing the received signal in the different antenna elements in the SSR. Both conventional beamforming, using a fixed set of beam positions, and adaptive beamforming, us- ing properties of the received data to optimize the reception, is analyzed. This is benchmarked against a system using analog conventional beamforming, where the signals are decoded in either a sum or a difference channel. The results of this study highlight the potential of digital beamforming, significantly increasing the probabil- ity of detection compared to the analog system, especially in an environment with large amounts of interfering signals. Among the adaptive algorithms tested, linear constraint minimum power (LCMP) beamforming, imposing null constraints on in- terfering signals, exhibited the most promising results. There is potential to build on this work, further analyzing the applicability of the different beamforming tech- niques. Future research areas could include analyzing more adaptive algorithms, as well as optimizing the ones studied in this thesis. Keywords: Secondary surveillance radar, ADS-B, Conventional beamforming, Adap- tive beamforming, Statistical signal processing, Monte Carlo simulations. v Acknowledgements First and foremost, I would like to express my deepest gratitude to my supervisors Ida-Maria Svensson, Lars Sundell and Juan Cabello Sánchez. Ida, for your immense commitment and continuous support throughout the whole project. Your mentor- ship and feedback helped me navigate this thesis. Lars, for all rewarding and ever so interesting discussions about theoretical concepts and ideas. Juan, for always being ready to answer questions and offer advice. Moreover, I would like to thank my examiner Adam Andersson, for contributing with valuable input, and my Saab manager Johan Lindström, for giving me the opportunity to work on this project. To my fellow master’s thesis students at Saab, thank you for enjoyable and much- needed breaks during many long and at times challenging days. Lastly, I would like to express my gratitude to friends and family for supporting me not only during this project, but throughout my Engineering Physics studies at Chalmers. And not least − thank you, Mom and Dad, for your ceaseless love and support. Mattias Wiklund, Gothenburg, June 2024 vii List of Acronyms Below is the list of acronyms that have been used throughout this thesis listed in alphabetical order: ACAS Airborne collision avoidance system ADS-B Automatic Dependent Surveillance-Broadcast AESA Active electronically scanned array CMC Multi-channel system using conventional beamforming DOA Direction of arrival ESPRIT Estimation of signal parameters via rotational invariance techniques FRUIT False replies unsynchronised in time GPS Global Positioning System ICAO International Civil Aviation Organization IFF Identification friend or foe LCMP Linear constraint minimum power LMS Least mean squares LNA Low-noise amplifier LOS Line-of-sight MAE Mean absolute error MDL Minimum detection level MPDR Minimum power distortionless response MUSIC Multiple signal classification PC Principal component PESA Passive electronically scanned array PSR Primary surveillance radar RADAR Radio detection and ranging SIR Signal-to-interference ratio SMI Sample matrix inversion SNR Signal-to-noise-ratio SOI Signal of interest SSR Secondary surveillance radar ULA Uniform linear array ix Contents List of Acronyms ix List of Figures xiii List of Tables xix 1 Introduction 1 1.1 Aims and outline of the thesis . . . . . . . . . . . . . . . . . . . . . . 2 2 Theory 3 2.1 Secondary Surveillance Radar and ADS-B . . . . . . . . . . . . . . . 3 2.2 Interfering signals . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 2.3 Array antennas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.3.1 Uniform linear arrays . . . . . . . . . . . . . . . . . . . . . . . 8 2.4 Properties of a receiving antenna . . . . . . . . . . . . . . . . . . . . 9 2.4.1 Friis transmission equation . . . . . . . . . . . . . . . . . . . . 10 2.4.2 Noise and minimum detectable signal . . . . . . . . . . . . . . 11 2.5 Configuration of SSR systems . . . . . . . . . . . . . . . . . . . . . . 12 2.6 Beamforming . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 2.6.1 Conventional beamforming . . . . . . . . . . . . . . . . . . . . 16 2.6.2 Adaptive beamforming . . . . . . . . . . . . . . . . . . . . . . 18 2.6.2.1 MPDR beamformer . . . . . . . . . . . . . . . . . . 18 2.6.2.2 LCMP beamformer . . . . . . . . . . . . . . . . . . . 20 2.6.2.3 PC beamformer . . . . . . . . . . . . . . . . . . . . . 20 2.7 Estimation of parameters . . . . . . . . . . . . . . . . . . . . . . . . . 21 2.7.1 Covariance matrix . . . . . . . . . . . . . . . . . . . . . . . . 22 2.7.2 Directions of arrival . . . . . . . . . . . . . . . . . . . . . . . . 22 3 Implementation 25 3.1 Environment model . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 3.2 Antenna model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 3.3 Signal model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 3.4 Analysis of adaptive beamforming . . . . . . . . . . . . . . . . . . . . 30 3.5 Probability of detection . . . . . . . . . . . . . . . . . . . . . . . . . . 30 3.6 Monte Carlo simulations . . . . . . . . . . . . . . . . . . . . . . . . . 31 xi Contents 4 Results and analysis 37 4.1 Probability of detection . . . . . . . . . . . . . . . . . . . . . . . . . . 37 4.2 Beamforming examples . . . . . . . . . . . . . . . . . . . . . . . . . . 46 4.2.1 Comparison to adaptive beamforming . . . . . . . . . . . . . . 47 4.2.2 Beamforming in the time domain . . . . . . . . . . . . . . . . 51 4.3 Further analysis of adaptive algorithms . . . . . . . . . . . . . . . . . 53 5 Discussion 61 5.1 Implications of the results . . . . . . . . . . . . . . . . . . . . . . . . 61 5.2 Limitations of the method . . . . . . . . . . . . . . . . . . . . . . . . 62 5.3 Future studies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 6 Conclusion 67 A Appendix I A.1 Complex optimization . . . . . . . . . . . . . . . . . . . . . . . . . . I A.2 Additional results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . II xii List of Figures 2.1 A simple illustration of the principle behind an airborne SSR sys- tem. An aircraft transmits interrogations on the 1030 MHz frequency band. Under certain circumstances, governed by regulations, another transponder-equipped aircraft should then transmit a reply on the 1090 MHz frequency band. . . . . . . . . . . . . . . . . . . . . . . . . 3 2.2 The waveform of a Mode S reply. It begins with a 8 µs long preamble, containing four pulses of length 0.5 µs. The length of the data block is either 56 or 112 µs, depending on whether it is a long or short reply. The information in the data block is encoded using a digital pulse-position modulation, with ”0” being represented by a pulse in the second half of the bit period, whereas ”1” has a pulse in the first half. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 2.3 A simple illustration of the principle behind ADS-B. An aircraft de- termines its position using GPS and transmits the information on the 1090 MHz frequency band. . . . . . . . . . . . . . . . . . . . . . . . . 5 2.4 Illustration of an ADS-B signal of interest (blue) partly overlapped with two interfering signals, one of type Mode S Long (red) and one of type Mode S Short (yellow). The amplitudes are chosen arbitrarily and no noise is added. . . . . . . . . . . . . . . . . . . . . . . . . . . 6 2.5 Three-dimensional illustration of a uniform linear antenna with six elements located on the x axis. The azimuth angle θ is measured from the y axis to the orthogonal projection of the radial line onto the x-y plane. The elevation angle ϕ is measured from the orthogonal projection to the radial line. The shaded region in the x-z plane represents the physical area of the antenna. . . . . . . . . . . . . . . . 8 2.6 Illustration of a uniform linear antenna with six elements with dis- tance d between them in the x-y plane. The wavefront has an angle θ relative to the antenna normal. . . . . . . . . . . . . . . . . . . . . 9 2.7 A simple illustration of an isotropic transmitting antenna and a re- ceiving antenna, where the incoming wavefront has an angle θ relative to the normal of the receiving antenna. The effective radiated power from the transmitter is PtGt and the receiving antenna has a physical area Aphysical. Note that the plane-wave approximation is assumed to hold at the receiver. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 xiii List of Figures 2.8 A two-channel Σ/∆ system with six antenna elements. After shifting the phase of the signals in each element, they are then combined into a left and right part. These signals are then added to and subtracted from each other, respectively, creating the Σ and ∆ channels. These two channels are then digitized. . . . . . . . . . . . . . . . . . . . . . 13 2.9 A multi-channel system where six antenna elements are individually digitized. Moving the phase shifters from hardware to software en- ables the possibility to control multiple independent beams simulta- neously. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 2.10 Example beam patterns with conventional beamforming. Note that this corresponds to the Σ channel in a Σ/∆ system. (a) The case θ0 = 0◦. (b) Contour plot for θ0 ranging from −90◦ to 90◦. . . . . . . 17 2.11 Example beam patterns with conventional beamforming for the ∆ channel in a Σ/∆ system. (a) The case θ0 = 0. (b) Contour plot for θ0 ranging from −90◦ to 90◦. . . . . . . . . . . . . . . . . . . . . . . . 17 3.1 An example square grid. Note that this grid only has 100 pixels, compared to 10,000 pixels used in the simulations. The SSR, marked in blue, is placed at the origin. The green pixel with a ”1” indicates that a target, in this case, has successfully been detected here. . . . . 25 3.2 Beam patterns in polar coordinates for the seven different beam po- sitions, 0◦, ±15◦, ±30◦ and ±45◦, used for conventional beamforming. 28 3.3 (a) The standard deviation of p̂det as a function of the number of Monte Carlo iterations for five randomly selected points. Note that p̂d corresponds to the sample mean Xd of the number of Monte Carlo iterations. (b) The simulation time per pixel as a function of the number of iterations. . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 4.1 Estimated probability of detection for the Σ/∆ system with analog conventional beamforming. (a) The case γ = 5 × 10−3 s−1 km−2. (b) The case γ = 1 × 10−2 s−1 km−2. (c) The case γ = 2 × 10−2 s−1 km−2. Red tones indicate higher probability, while blue tones indicate a lower probability. Note that the maximum line-of-sight distance at 900 km is marked by a white line. . . . . . . . . . . . . . . . . . . . . 38 4.2 Estimated probability of detection for the multi-channel system with digital conventional beamforming (CMC). (a) The case γ = 5 × 10−3 s−1 km−2. (b) The case γ = 1 × 10−2 s−1 km−2. (c) The case γ = 2 × 10−2 s−1 km−2. Red tones indicate higher probability, while blue tones indicate a lower probability. Note that the maximum line- of-sight distance at 900 km is marked by a white line. . . . . . . . . . 39 4.3 Estimated probability of detection for the multi-channel system with adaptive MPDR beamforming. (a) The case γ = 5 × 10−3 s−1 km−2. (b) The case γ = 1 × 10−2 s−1 km−2. (c) The case γ = 2 × 10−2 s−1 km−2. Red tones indicate higher probability, while blue tones indicate a lower probability. Note that the maximum line-of-sight distance at 900 km is marked by a white line. . . . . . . . . . . . . . . . . . . . . 40 xiv List of Figures 4.4 Estimated probability of detection for the multi-channel system with adaptive LCMP beamforming. (a) The case γ = 5 × 10−3 s−1 km−2. (b) The case γ = 1 × 10−2 s−1 km−2. (c) The case γ = 2 × 10−2 s−1 km−2. Red tones indicate higher probability, while blue tones indicate a lower probability. Note that the maximum line-of-sight distance at 900 km is marked by a white line. . . . . . . . . . . . . . . . . . . . . 41 4.5 Estimated probability of detection for the multi-channel system with adaptive PC beamforming. (a) The case γ = 5 × 10−3 s−1 km−2. (b) The case γ = 1 × 10−2 s−1 km−2. (c) The case γ = 2 × 10−2 s−1 km−2. Red tones indicate higher probability, while blue tones indicate a lower probability. Note that the maximum line-of-sight distance at 900 km is marked by a white line. . . . . . . . . . . . . . . . . . . . . 42 4.6 The three different subareas defined by the distance between the tar- get and the SSR. Region 1: 0−300 km. Region 2: 300−600 km. Region 3: 600−900 km. . . . . . . . . . . . . . . . . . . . . . . . . . . 43 4.7 Average probability of detection as a function of the message rate per area γ [s−1 km−2] for the different systems in different subareas. (a) Targets at a distance of 0−300 km. (b) 300−600 km. (c) 600−900 km. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 4.8 Antenna pattern for the Σ/∆ system, where the beam position θ0 steps from −45◦ to 45◦. (a) θ0 = 0◦, with no signal of interest present. (b) θ0 = 15◦, where a signal of interest (solid line) is present at 0◦ together with one interfering signal (dashed line) at 15◦. . . . . . . . 46 4.9 Antenna pattern for a multi-channel system with reception in simul- taneous active beam positions θ0 = 0◦, ±15◦, ±30◦, ±45◦. Each beam is indicated by a unique color. A signal of interest is present at 0◦, together with one interfering signals at 15◦. . . . . . . . . . . . . . . . 46 4.10 An ADS-B signal of interest interfered by a Mode S Long interfering signal, starting 80 µs after the SOI, and a Mode S Short interfering signal, starting 24 µs before the SOI. This corresponds to Cases 2a and 2b in Table 4.3. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 4.11 Antenna patterns for Case 1a in (a), with one interfering signal at +45◦, and for Case 1b in (b), with one interfering signal at +15◦. . . 49 4.12 Antenna patterns for Case 2a in (a), with two interfering signals at ±45◦, and for Case 2b in (b), with two interfering signals at ±15◦. . . 49 4.13 Antenna pattern for Case 3, with four interfering signals at ±45◦ and ±15◦. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 4.14 The incoming signal x(t) before any processing. (a) The actual sim- ulated signal. (b) A simplified image, highlighting the SOI (in blue) and the interfering signal (in red). . . . . . . . . . . . . . . . . . . . . 51 4.15 The signal obtained from applying the MPDR algorithm. (a) The actual simulated signal. (b) A simplified image, highlighting the SOI (in blue) and the interfering signal (in red). The MDL threshold is shown in black and the SIR threshold is shown in cyan. In this case, the SOI is detected. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 xv List of Figures 4.16 The signal obtained from applying the LCMP algorithm. (a) The actual simulated signal. (b) A simplified image, highlighting the SOI (in blue) and the interfering signal (in red). The MDL threshold is shown in black and the SIR threshold is shown in cyan. In this case, the SOI is detected. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 4.17 The signal obtained from applying conventional beamforming (Σ chan- nel with beam position θ0 = 0◦). (a) The actual simulated signal. (b) A simplified image, highlighting the SOI (in blue) and the inter- fering signal (in red). The MDL threshold is shown in black and the SIR threshold is shown in cyan. In this case, the SOI is not detected. 53 4.18 The mean absolute error of the DOA estimation with the ESPRIT algorithm as a function of the direction θSOI for the signal of interest. The results are shown for the three different values of γ. . . . . . . . 54 4.19 The average number of estimated signals in the ESPRIT algorithm as a function of the number of simulated signals. The results are shown for the three different values for γ. The dashed lines in corresponding colors display the distribution of the number of simulated signals. Note that number of simulated signals corresponding to less than 0.1 % of the total number of data points are disregarded. . . . . . . . . . 54 4.20 Average probability of detection as a function of the message rate per area γ [s−1 km−2] for the different systems in different subareas. (a) Targets at a distance of 0−300 km. (b) 300−600 km. (c) 600−900 km. The solid lines are the same as found for MPDR and LCMP in Figure 4.7, where the directions of arrival are estimated. The dashed lines here represent corresponding simulations instead using a priori knowledge of the directions of arrivals. The CMC is included for comparison. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 4.21 Antenna patterns for Case 1a in (a), with one interfering signal at +45◦, and for Case 1b in (b), with one interfering signal at +15◦. . . 58 4.22 Antenna patterns for Case 2a in (a), with two interfering signals at ±45◦, and for Case 2b in (b), with two interfering signals at ±15◦. . . 58 4.23 Antenna pattern for Case 3, with four interfering signals at ±45◦ and ±15◦. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 4.24 The signal-to-interference ratio as a function of the difference in az- imuth ∆θ between the SOI and one interfering signal. Results for 1,000 runs with randomized angles for the interfering signal. The SIR threshold at 6 dB is marked by a black dashed line. . . . . . . . 59 4.25 The signal-to-interference ratio as a function of the proportion of the signal of interest covered by one interfering signal. Results for 1,000 runs with randomized proportions covered. The SIR threshold at 6 dB is marked by a black dashed line. . . . . . . . . . . . . . . . . . . 60 xvi List of Figures 5.1 (a) The discrete Poisson probability density function for the number of interfering signals. The three different cases for the γ parameter are shown. (b) The probability that the number of interfering signals are ≤ 5 as a function of the γ parameter. The three different cases investigated in this work are displayed as black dashed lines. . . . . . 64 A.1 Estimated probability of detection for the multi-channel system with adaptive MPDR beamforming with a priori knowledge of the number of signals and their direction. (a) The case γ = 5 × 10−3 s−1 km−2. (b) The case γ = 1 × 10−2 s−1 km−2. (c) The case γ = 2 × 10−2 s−1 km−2. Red tones indicate higher probability, while blue tones indicate a lower probability. Note that the maximum line-of-sight distance at 900 km is marked by a white line. . . . . . . . . . . . . . . . . . . . . II A.2 Estimated probability of detection for the multi-channel system with adaptive LCMP beamforming with a priori knowledge of the number of signals and their direction. (a) The case γ = 5 × 10−3 s−1 km−2. (b) The case γ = 1 × 10−2 s−1 km−2. (c) The case γ = 2 × 10−2 s−1 km−2. Red tones indicate higher probability, while blue tones indicate a lower probability. Note that the maximum line-of-sight distance at 900 km is marked by a white line. . . . . . . . . . . . . . . . . . . . . III xvii List of Figures xviii List of Tables 3.1 The three different values used for the message rate per area γ, cor- responding to the degree of disturbance. The corresponding Poisson parameters λlong Poi and λshort Poi , describing the average number of received Mode S Long and Mode S Short replies, are also presented together with their sum. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 3.2 Values for the parameters used when modelling the physical envi- ronment experienced by the simulated SSR system. This includes geographical parameters as well as properties of the incoming signals. 27 3.3 Values for the parameters used when modelling the simulated ULA antenna in the SSR system. . . . . . . . . . . . . . . . . . . . . . . . 29 3.4 Coordinates and standard deviations of p̂det for the five randomly selected points. The results are displayed for 10, 100, 1,000 and 10,000 iterations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 4.1 The different algorithms analyzed. The second column describes the number of digitized channels. In the case of Σ/∆, the output Σ and ∆ channel are the two digital channels. In the other cases, all six antenna elements are digitized separately. The third column describes whether analog or digital phase shifters are used, resulting in different types of beamforming. . . . . . . . . . . . . . . . . . . . . . . . . . . 37 4.2 The average probability of detection shown in Figure 4.7 displayed in table form. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 4.3 Directions of arrival for the interfering signals in the five different test cases used for comparing the adaptive algorithms. The subscript ”L” denotes a Mode S Long signal, while ”S” denotes a Mode S Short signal. 48 4.4 The gain evaluated in the directions of the SOI and the interfering signals for the different cases and algorithms. The gain values are displayed in the unit dBi. Note that the numbers presented here correspond to the antenna patterns presented in Figures 4.11, 4.12 and 4.13. The case with conventional beamforming in a Σ channel at beam position θ0 = 0◦ is included as a comparison. . . . . . . . . . . 50 4.5 The average probability of detection shown in Figure 4.20 displayed in table form. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 xix List of Tables xx 1 Introduction Ever since the breakthrough of electromagnetic theory in the 19th century, the ap- plication of radio waves has been a topic of great interest in physics and engineering. Based on the property that surfaces can reflect these waves, it was realized that this phenomenon could be used to detect different objects. The potential for military applications was soon realized, leading to the development of radio detection and ranging (RADAR) systems during World War II [1], [2]. In the decades that fol- lowed, the technology was developed further. The areas of application grew, also including civilian aircraft surveillance, and the acronym ”RADAR” now became a concept of its own. A radar system commonly used in aircraft surveillance is the primary surveillance radar (PSR) [1]. It functions by measuring the reflected waves from aircraft within range. From this, the position of the aircraft can be estimated. One major advantage of the PSR is its independence of the onboard equipment of the targets. However, the system has its limitations. It does not know whether the reflected signals are coming from actual aircraft targets or from other objects. Such unwanted reflections are known as clutter. Moreover, targets close to each other are difficult to identify separately. To address the limitations of the PSR, the secondary surveillance radar (SSR) was developed. Rather than relying on reflected signals, the SSR operates by transmit- ting interrogations to targets. The received reply contains the identity of the target. The SSR requires less power than the PSR, as the signals only has to reach the target rather than also being reflected back. Moreover, the problem with clutter is reduced by transmitting the interrogations and replies on different frequency bands: 1030 MHz for interrogations and 1090 MHz for replies. The SSR has multiple appli- cations related to air traffic control. One such important application is Identification friend or foe (IFF) [3]. This identification system, implemented in both ground- and air-based radars, was originally developed to prevent ”friendly fire” incidents. IFF systems are today used for both military and civilian purposes, and is a requirement for several flying platforms. The SSR is not the only surveillance technology operating on the 1090 MHz band. Automatic Dependent Surveillance-Broadcast (ADS-B) is a technology in which air- craft periodically broadcast their position, identity and other information of interest [1]. Hence, there is no need for interrogations. Since ADS-B signals are transmitted on the same frequency band as SSR replies, 1090 MHz, this enables them to be 1 1. Introduction detected by SSR systems [4]. As air traffic continues to grow, there is an increas- ing risk that ADS-B signals are interfered by other signals on the 1090 MHz band. ADS-B is an important part of air traffic control and flight safety, and missed replies can therefore potentially have dire consequences. Thus, it is of great importance to enhance the ability to decode the information encoded in an ADS-B signal of interest. There are different signal processing techniques which can be applied to detect a signal of interest, despite it being interfered by other signals. One important tech- nique is beamforming, where the receiving antenna in the radar is used in such a way that the reception of signals is enhanced in some directions and reduced in others [5], [6]. There is conventional beamforming, which utilizes the fact that an incom- ing wave has different times of arrival at the different antenna elements building up the antenna. This corresponds to a relative phase shift between the elements. By adding a compensatory phase shift, maximum reception can be obtained in a desired direction. Depending on whether this is performed in hardware or software, the beamforming can be either analog or digital [7]. Moreover, there is adaptive beamforming, in which properties of the received data are used to optimize the reception performance in certain directions [6]. Unlike the conventional case, this enables more possibilities in handling different interference scenarios. 1.1 Aims and outline of the thesis In this thesis, the main aim is to analyze the performance of an airborne SSR system using different types of beamforming techniques. The signal of interest is assumed to be an incoming ADS-B signal. The performance parameter analyzed is the probability of detection, defined here by the amplitude of the signal of interest as well as its relation to the amplitude of possible interfering signals. This parameter is estimated for every point in a simulated spatial grid using Monte Carlo simulations. Simulations are performed for digital systems using both conventional and adaptive beamforming, and their performance is benchmarked against an analog system using conventional beamforming. The outline of the thesis is as follows: In Chapter 2, the underlying theory behind the work is presented. This includes the characteristics of the different signals and systems, antenna theory and mathematical formulations of beamforming. In Chap- ter 3, the implementation is described. This includes all models and assumptions, as well as a description of the Monte Carlo simulations performed. In Chapter 4, the obtained results are presented. This includes main results for the probability of detection, as well as some illustrative beamforming examples and a deeper analysis on the adaptive algorithms. In Chapter 5, the results are interpreted and discussed in relation to the underlying theory and the methods used. Potential limitations of the methods are discussed, as well as an outlook for future studies. In Chapter 6, the findings in the thesis are summarized, thus concluding this work. 2 2 Theory The following chapter contains the theory behind this thesis. This includes the characteristics of SSR replies and ADS-B signals, configuration of SSR systems, fundamental antenna theory, mathematical aspects of conventional and adaptive beamforming, as well as estimation theory. 2.1 Secondary Surveillance Radar and ADS-B A secondary surveillance radar uses an interrogator transmitting signals on the 1030 MHz frequency band [1]. A transmitter-responder, transponder, within range can receive this signal and transmit a reply with a frequency on the 1090 MHz band. According to international rules and regulations, all flying platforms are obliged to have a transponder and antenna for receiving interrogations and replying [4]. A simple illustration of the principle behind an airborne SSR system, which is the focus of this thesis, can be seen in Figure 2.1 below. The reply contains information requested by the interrogation, where different interrogation modes contain different types of information. The different modes are in general characterized by the spacing between the pulses in the transmitted waveform. Figure 2.1: A simple illustration of the principle behind an airborne SSR system. An aircraft transmits interrogations on the 1030 MHz frequency band. Under cer- tain circumstances, governed by regulations, another transponder-equipped aircraft should then transmit a reply on the 1090 MHz frequency band. The most common mode with civilian use is Mode S (Mode Select) [1]. It it ”selec- tive” in the sense that this mode enables selective interrogation of individual aircraft. This is possible thanks to a 24-bit unique code, or address, being assigned to each air- craft. The address is encrypted in the waveform and regulated by the International Civil Aviation Organization (ICAO) [4]. A 24-bit address gives 224 −2 = 16, 777, 214 unique codes (the two codes only containing zeros or ones are not valid). This can 3 2. Theory for example be compared to the older Mode A, where the address instead is given by a four-digit octal number, yielding only 84 = 4, 096 different codes. For increased air traffic, this can result in several aircraft in the same region with identical addresses, with potential safety risks as a consequence. Looking closer on the transmitted waveforms, or codes, a Mode S interrogation starts with two pulses of length 0.8 µs, with a relative time delay of 2 µs [1]. This is known as the preamble, and has the purpose of preventing replies from transponders using other modes. After the preamble follows a longer pulse, the data block, containing the interrogation information. The data block is phase-modulated using differential phase-shift keying modulation, with a phase reversal resulting in bit 1 and no phase reversal resulting in bit 0. There are two types of Mode S interrogations: short with a 16.125 µs data block containing 56 bits and long with a 30.25 µs data block containing 112 bits. The reply from a transponder has a format similar to the interrogation with a preamble and a data block, but with a different composition. The reply preamble has four pulses of length 0.5 µs. The data block that follows contains 56 or 112 bits, with the rate of one bit/µs [1]. The information requested in the interrogation, for example the identity of the aircraft, is encoded in the different parts of the data block. It is encoded using pulse-position modulation, where bit 1 is indicated by a 0.5 µs pulse followed by a flat signal of the same length. Bit 0 is instead indicated by a flat signal followed by a pulse. This is also known as Manchester encoding [8]. Consequently, 0 is coded as [0, 1] and 1 is coded as [1, 0]. An example of this encoding scheme is illustrated in Figure 2.2 below. Figure 2.2: The waveform of a Mode S reply. It begins with a 8 µs long preamble, containing four pulses of length 0.5 µs. The length of the data block is either 56 or 112 µs, depending on whether it is a long or short reply. The information in the data block is encoded using a digital pulse-position modulation, with ”0” being represented by a pulse in the second half of the bit period, whereas ”1” has a pulse in the first half. Different information is encoded in different parts of the waveform. For example, this includes the format number, indicating which type of response it is, the 24-bit 4 2. Theory address and a parity check, used for finding and correcting potential errors in the message [1]. Mathematically, the reply waveform can be described as a(t) = 2treply−1∑ k=0 akrect(t − ktpulse), (2.1) where treply is either 64 or 120 µs, ak is the Manchester encoded data for entry k and rect(t) is a rectangular pulse of width tpulse = 0.5 µs [9]. A subset of Mode S Long replies is ADS-B, Automatic Dependent Surveillance- Broadcast. This is a surveillance system where aircraft broadcasts their position and other relevant flight data periodically, without any need for interrogation [1]. A simple illustration of this principle can be seen in Figure 2.3. ADS-B is ”dependent” in the sense that the transmitted information is dependent on other systems, such as the Global Positioning System (GPS). This is a commonly used surveillance technology and since 2020, civil aircraft in Europe and the United States are required to be compatible with ADS-B. Figure 2.3: A simple illustration of the principle behind ADS-B. An aircraft de- termines its position using GPS and transmits the information on the 1090 MHz frequency band. ADS-B is a Mode S of type ”Extended Squitter”. In this case, a ”squitter” refers to periodic burst transmissions not generated by interrogations [4]. It should be noted that ADS-B messages have different transmissions rates depending on the information transmitted. Messages related to airborne position and velocity are transmitted twice per second. 2.2 Interfering signals In order for the information in the ADS-B waveform to be decoded, the signal must first be detected. There are several factors affecting the possibility to detect a signal of interest (SOI), in this case an ADS-B signal. One such factor is the occurrence of interfering signals. In areas with large amounts of air traffic, there are many signals in air on the 1090 MHz frequency band. Hence, the SOI may be interfered by these 5 2. Theory signals. This could for example be other ADS-B signals or replies from a transponder that has been interrogated by another SSR interrogator. This phenomenon is known as FRUIT (False replies unsynchronized in time) [10]. When such replies overlap with the SOI, this results in ”garbling”. This could potentially result in the target not being detected or detection of non-existing targets. A simple illustration of this phenomenon can be seen in Figure 2.4, where two Mode S replies are interfering with an ADS-B signal of interest. Figure 2.4: Illustration of an ADS-B signal of interest (blue) partly overlapped with two interfering signals, one of type Mode S Long (red) and one of type Mode S Short (yellow). The amplitudes are chosen arbitrarily and no noise is added. The interfering signals overlapping with an SOI can be seen as random points in time following a Poisson point process [11], [12]. This assumption holds as long as the interfering signals occur independently of each other. The number of interfering signals in a time span of limited size can then be seen as a random variable following a Poisson distribution. The discrete probability density function for this distribution, given a discrete random variable, X, is given by P(X = k) = λk Poie−λPoi k! , (2.2) where k is the number of events, with an expected value of λPoi [13]. Translated to the scenario of interest here, the expected value can be written as λPoi = tPoiµ, where tPoi is the time span and µ is the intensity of interfering signals. Analyzing the time span in closer detail, the signal of interest is an ADS-B signal with length tSOI = tADS-B = 120 µs. Now assume that an interfering signal has length tint. Note that the start of a signal can be seen as a Poisson event. If an SOI begins at time t, an interfering signal overlaps with the SOI as long as it begins anywhere in the interval [t−tint, t+tSOI]. Therefore, in order for λPoi to represent the 6 2. Theory expected number of interfering signals, the time span is given by tPoi = tSOI + tint. For example, tint = 120 µs if the interfering signal is of type Mode S Long and tint = 64 µs if it is of type Mode S Short. It is, however, possible to detect an SOI, even though interfering signals are over- lapping with it. By applying different signal processing techniques, the amplitudes of the interfering signals can be suppressed, while the amplitude of the SOI could be amplified. This affects the signal-to-interference-ratio, defined as the quotient be- tween the received signal power and the sum of the power of the interfering signals. For a given time, t, and spatial location, r, this gives SIR(r, t) = PSOI(r, t)∑D i=1 Pi(r, t) , where PSOI(r, t) is received power from the signal of interest and P1(r, t), P2(r, t), . . . , PD(r, t) are the signal power for D interfering signals [14]. Often, this ratio is written in decibel units, according to SIR[dB](r, t) = 10 log10 ( PSOI(r, t)∑D i=1 Pi(r, t) ) . To detect an SOI and decode the information encoded in its waveform, the SIR should typically be above a given threshold. 2.3 Array antennas In the last sections, the focus has been on the characteristics of signals and how interference can impair the possibility to detect them. The focus is now shifted towards the configuration of SSR systems. For an SSR to be able to transfer waves propagating in free space to electrical signals in a circuit, an antenna is required. In this thesis, the focus is on an SSR system based on an array antenna, where the signals in several smaller antennas (antenna elements) interfere with each other constructively and destructively, respectively [5]. In the case of a receiving antenna, this enhances the reception in some directions and reduces it in others, depending on the relative phase between the elements. This process is known as beamforming and is described in more detail in Section 2.6. An incoming signal has different times of arrival at the different elements. For narrowband signals, this time difference can be translated into a phase difference between the elements [6]. In a phased array, the beamforming process is done by utilizing this phase difference and adding a compensatory phase shift at each element. In practice, this is done by having each element connected to a computer- controlled phase shifter. This way, the antenna beam can be electronically steered to a desired direction. There are different types of phased arrays. In a passive electronically scanned array (PESA), all the individual elements are connected to one and the same transmitter/receiver module, while in an active electronically scanned array (AESA), the different elements all have their own transmitter and receiver module [15]. 7 2. Theory 2.3.1 Uniform linear arrays A special case of an array antenna is the uniform linear array (ULA) [6]. This type of array consists of N antenna elements placed along a line with a uniform spacing d between them. For this purpose, it is convenient to introduce a spherical coordinate system defined according to Figure 2.5. A unit vector in this system can be written as â = sin θ cos ϕx̂ + cos θ cos ϕŷ + sin ϕẑ, where θ is the azimuth angle and ϕ is the elevation angle. The figure displays a ULA with six symmetrically placed elements on the x axis, with a physical area in the x-z plane. Figure 2.5: Three-dimensional illustration of a uniform linear antenna with six elements located on the x axis. The azimuth angle θ is measured from the y axis to the orthogonal projection of the radial line onto the x-y plane. The elevation angle ϕ is measured from the orthogonal projection to the radial line. The shaded region in the x-z plane represents the physical area of the antenna. If the elements are symmetrically placed around the origin and on the x axis, the position of the nth element can be written as pn = ( n − N + 1 2 ) d, for n = 1, 2, . . . , N . Note that this work focuses on the horizontal cut ϕ = 0, assuming an airborne radar looking at long distances. An illustration of a uniform linear antenna with N = 6 in the x-y plane (ϕ = 0) is shown in Figure 2.6. This displays an incoming plane wave in the x-y plane, with azimuth angle θ relative to antenna normal on the y axis. 8 2. Theory Figure 2.6: Illustration of a uniform linear antenna with six elements with distance d between them in the x-y plane. The wavefront has an angle θ relative to the antenna normal. Now assume that a plane wave in the x-y plane with azimuth angle θ is received by the ULA. The wave arriving at element n travels an extra distance of ∆x = d sin θ compared to an adjacent element. This corresponds to a time delay of ∆t = ∆x/c, where c is the speed of light. Assuming a narrowband signal, the time delay can in turn be written as a phase shift according to ∆φ = 2π ∆t Tperiod = 2π ∆x λ = 2πd sin θ λ , (2.3) where Tperiod and λ is the period and wavelength of the incoming wave, respectively [6]. How this phase difference relates to the beamforming process is explored more in Section 2.6. 2.4 Properties of a receiving antenna There are other aspects than the occurrence of interfering signals affecting the ability to detect a signal of interest, not least the performance of the receiving antenna in the SSR. One important performance parameter is the gain of the antenna. The gain of a receiving antenna describes how well the antenna can absorb incident power in a given direction and convert it to electrical power [16], [17]. In general, this is described in terms of a ratio of the radiation intensities for a directional antenna and a theoretical isotropic antenna, i.e., an antenna absorbing the same intensity in all directions. Gain can also be defined as G = ηD, where η is the radiation efficiency and D is the directivity. The radiation efficiency is a measure of how efficient an antenna 9 2. Theory can convert between the power of the incoming wave and electrical power. For a lossless antenna, η = 1. The directivity describes the degree of concentration for the absorbed radiation in a given direction. 2.4.1 Friis transmission equation For a given frequency, the gain can also be expressed in terms of effective aperture areas for the directional antenna and a theoretical isotropic antenna. The effective aperture area, Aeff , describes how much power, Pr, from an incoming plane wave which is captured by an antenna and transformed to electrical power [18]. Again focusing on the horizontal cut ϕ = 0, this can be described by the intensity, S(θ), according to Pr = Aeff |S(θ)|. The power absorbed by the antenna can also be related to the incoming plane wave power, Pin, according to Pr = ηPin = ηAphysical|S(θ) · ŷ| = ηAphysical|S(θ)| cos θ, where Aphysical is the physical area of the antenna and ŷ is a unit vector normal to the antenna on the x axis. Note that Pin is free-space power, while Pr is conductive power. This gives Aeff = ηAphysical cos θ. (2.4) Figure 2.7: A simple illustration of an isotropic transmitting antenna and a re- ceiving antenna, where the incoming wavefront has an angle θ relative to the normal of the receiving antenna. The effective radiated power from the transmitter is PtGt and the receiving antenna has a physical area Aphysical. Note that the plane-wave approximation is assumed to hold at the receiver. Relating the effective aperture area to the gain, Gr, of the receiving antenna, this gives Gr = Pr Piso = Aeff Aiso , 10 2. Theory where Aiso is the effective aperture area of an isotropic antenna with output power Piso [16]. For an incident wave with wavelength λ, the effective aperture area of an isotropic antenna is Aiso = λ2 4π . Therefore, the relation between the effective aperture area of the receiving antenna and its gain is given by Aeff = Grλ 2 4π . The effective radiated power from a transmitting source is given by scaling the output power, Pt, with its gain, Gt. Assuming an isotropic transmitting antenna at a distance r from the receiving antenna, its intensity, St = |St(θ)|, is then given by St = PtGt 4πr2 . The available power at the receiver is therefore given by Pr = StAeff = PtGt 4πr2 · Grλ 2 4π = PtGtGr ( λ 4πr )2 = PtGtGr Lp , (2.5) where Lp = (4πr/λ)2 is the path loss. This result is known as Friis transmission equation, and is used to calculate the amount of power received at a certain distance given the characteristics of the transmitting and receiving antennas [16]. An illus- tration of a transmitting and a receiving antenna can be seen in Figure 2.7. Using decibel units instead, Friis transmission equation can be written as P [dBW] r = P [dBW] t + G [dBi] t + G[dBi] r − L[dB] p . Note that power is expressed in dBW (decibel-Watt) and gain in dBi (decibel- isotropic), since this quantity is defined in relation to an isotropic antenna. 2.4.2 Noise and minimum detectable signal To detect a signal, the received power, Pr, should be sufficiently high so that it can be distinguishable from background noise in the system. The noise can for example be thermal Johnson−Nyquist noise, generated by random motion of electrons in the components of the receiving antenna [19]. This sensitivity of the system can be described as the system minimum detection level, MDL, defined according to MDL[dBW] = Noise floor[dBW] + SNR[dB] min . SNRmin is the minimum signal-to-noise-ratio required for the receiver to detect the signal [20]. The noise floor level is given by properties of the system such as band- width and temperature. It is convenient to describe the incoming signals as complex numbers, compactly capturing both magnitude and phase. Assuming the noise is a complex random 11 2. Theory variable, n, its variance is given by the sum of the variance of its real and imaginary parts according to Var[n] = Var[Re(n)] + Var[Im(n)]. Given that both the real and imaginary parts of the noise have variance σ2 n, this gives Var[n] = 2σ2 n [21]. Note that the standard deviation of the noise corre- sponds to its amplitude, meaning that the variance, the square of the standard deviation, corresponds to its power. Assuming a complex normal distribution with zero mean, the noise in each of the N antenna elements, n ∈ CN , is then given by n ∼ CN (0, 2σ2 nIN), where IN is the N -dimensional identity matrix. The proba- bility density function for an N -dimensional complex normal distribution is given by f(z) = 1 πNdet(Σ) exp ( −(z − µ)HΣ−1(z − µ) ) , (2.6) where ”H” is the Hermitian transpose, µ ∈ CN is the mean and Σ ∈ CN×N is the non-singular covariance matrix [21]. 2.5 Configuration of SSR systems There are different ways to utilize the antenna elements in the phased array in the SSR system. One alternative is to use a two-channel system, or a Σ/∆ system [17]. An illustration of such as system can be seen in Figure 2.8. Focusing on a receiving antenna, the power of the incoming signal is amplified using low-noise amplifiers (LNA). The signal in each element is then electronically shifted to the desired phase. Following this, the signal can be divided into two halves, left and right. The division into a left and right half can be described in terms of beam patterns, B(θ). B(θ) describes how well the antenna can receive a signal from certain direc- tion. As before, the horizontal cut ϕ = 0 is assumed. Note that the beam pattern is related to the gain parameter introduced in Section 2.4. The gain is directly pro- portional to the power pattern, P (θ), which is the squared magnitude of the beam pattern [6]. Consequently, G(θ) ∝ P (θ) = |B(θ)|2. 12 2. Theory Figure 2.8: A two-channel Σ/∆ system with six antenna elements. After shifting the phase of the signals in each element, they are then combined into a left and right part. These signals are then added to and subtracted from each other, respectively, creating the Σ and ∆ channels. These two channels are then digitized. The normalized beam pattern for element n can be denoted Bn(θ). Assuming N is an even number of antenna elements, the beam patterns for the left and right parts are given by  BLeft(θ) = ∑N 2 n=1 Bn(θ), BRight(θ) = ∑N n= N 2 +1 Bn(θ). From these two halves, the beam pattern for the so-called sum and difference chan- nels Σ and ∆ can then be obtained according to{ BΣ(θ) = BRight(θ) + BLeft(θ), B∆(θ) = BRight(θ) − BLeft(θ). (2.7) These two channels are then digitized, thus moving the processing from hardware to software. For a receiving antenna, it is then possible to detect a signal in both the Σ and ∆ channels independently. 13 2. Theory Figure 2.9: A multi-channel system where six antenna elements are individually digitized. Moving the phase shifters from hardware to software enables the possibil- ity to control multiple independent beams simultaneously. There are, however, other alternatives. Instead of only having two digital chan- nels, there is a possibility to digitize all antenna elements individually [7], [22]. An illustration of such a multi-channel system can be seen in Figure 2.9. Notably, this moves the phase shifters to a software environment. The enables digital beam- forming, where the antenna can control multiple independent beams simultaneously. Since the beams can have separate directions, this allows for simultaneous coverage of larger areas. This is in contrast to the analog beamforming used in the Σ/∆ system, where the different beam positions must be stepped through one by one. Note that the multiple simultaneous Σ channels makes reception in a ∆ channel redundant in this case. Instead of using the digital phase shifters to obtain maximum constructive interfer- ence in a predetermined direction, the signals in the antenna elements can be com- bined in different linear combinations. The optimal linear combination, or weighting, of the elements can be obtained be using statistics from the received data [6]. This is known as adaptive beamforming and is explored in more detail in the following section. 2.6 Beamforming After describing the technical aspects of an SSR system, the beamforming concept introduced earlier is now analyzed in greater detail. Assume a plane wave arriving at a ULA with N elements. The wavenumber vector, 14 2. Theory k, corresponds to the spatial frequency of an incoming wave. It can be defined as k = 2π λ û, where λ is the wavelength of the wave and û = −(sin θ cos ϕx̂+cos θ cos ϕŷ+sin ϕẑ) is a unit vector in the direction of arrival [5], [6]. For a specific wavenumber, the steering vector, v(k) ∈ CN , describes how the plane wave propagates across the ULA, resulting in a phase shift between the signals in the different elements. The nth element of the steering vector is given by [v(k)]n = exp(−ikTpn), where pn is the position for element n. Moreover, the weight vector, w ∈ CN , describes how the signal is processed after reaching the array. It adjusts the phase and amplitude of the signals received by each antenna element. The choice of w depends on the type of beamforming. From the definition of the steering vector and the weight vector, the beam pattern first mentioned in Section 2.5 can be defined as B(k) = wHv(k). (2.8) Note that the beam pattern can be seen as the Fourier transform of the array response, b(t) [6]. This can be seen by analyzing the time-domain signal b(t) = N∑ n=1 wH n δ(t − tn), where tn is the relative time delay for the incoming wave at element n and δ(t) is the Dirac delta function or unit impulse. Shifting to the frequency domain with a Fourier transform gives B(ω) = F [b(t)] = N∑ n=1 wH n exp(−iωtn). Note that the time delay can be written as tn = ûTpn c , where c is the velocity of propagation for the wave. Moreover, the angular frequency can be written as ω = 2πc/λ. Combining these expressions yields ωtn = 2πc λ ûTpn c = kTpn. The beam pattern can now be written as B(ω) = B(k) = N∑ n=1 wH n exp(−ikTpn) = wHv(k), 15 2. Theory as stated in Equation (2.8). Recall that for a ULA located on the x axis, the position vector can be written as pn = ( n − N + 1 2 ) dx̂. Assuming the horizontal cut ϕ = 0, this results in the steering vector [v(θ)]n = exp [ i 2π λ sin θ ( n − N + 1 2 ) d ] . Observe that the phase shift between the adjacent elements is in accordance with Equation (2.3). The beam pattern can then be written as a function of azimuth according to B(θ) = wHv(θ). As mentioned in Section 2.5, the beam pattern is a useful characteristic for an array antenna as it directly relates to the gain by G(θ) ∝ |B(θ)|2. 2.6.1 Conventional beamforming There are different ways to choose the weight vector w. In conventional beam- forming, the weights are fixed and do not depend on any information from the received signal. The weights are chosen to generate maximum constructive interfer- ence between the elements in a given beam position θ0. In a phased array, this is achieved by shifting the phase for the signals in each element by a proper amount [6]. Consequently, the signals received in each element can be added coherently. This corresponds to defining the weight vector as the steering vector for the antenna, evaluated at the beam position, according to w = v(θ0). (2.9) For a ULA, the beam pattern is then given by B(θ) = wHv(θ) = vH(θ0)v(θ) = = N∑ n=1 exp [ −i 2π λ sin θ0 ( n − N + 1 2 ) d ] exp [ i 2π λ sin θ ( n − N + 1 2 ) d ] = = N∑ n=1 exp [ i 2π λ (sin θ − sin θ0) ( n − N + 1 2 ) d ] = = exp[i(φ1(θ) − φ1(θ0))] + exp[i(φ2(θ) − φ2(θ0))] + · · · + exp[i(φN(θ) − φN(θ0))], where φ1, φ2, . . . , φN are the phases in the different elements. As expected, the magnitude of the sum above has its maximum value for θ = θ0. Figure 2.10 displays normalized example beam patterns illustrating this, with the case θ0 = 0◦ shown in Figure 2.10(a) and a contour plot for θ0 ranging from −90◦ to 90◦ shown in Figure 2.10(b). It can be noted that the beam pattern not only has a peak for θ = θ0, the main lobe, but also smaller peaks symmetrically placed around it, the sidelobes. Sidelobes 16 2. Theory are a characteristic feature of directional antennas and a direct consequence of the antenna geometry [6]. (a) (b) Figure 2.10: Example beam patterns with conventional beamforming. Note that this corresponds to the Σ channel in a Σ/∆ system. (a) The case θ0 = 0◦. (b) Contour plot for θ0 ranging from −90◦ to 90◦. Note that in the case of a Σ/∆ system described in Section 2.5, the definition of w in Equation (2.9) only holds for the Σ channel. For the ∆ channel, the beam pattern of first half of the elements, BLeft, should be subtracted from the second half, BRight, according to Equation (2.7). Therefore, the nth component of the weight vector can in this case be defined according to [w]n = −[v(θ0)]n if n = 1, . . . , N/2 +[v(θ0)]n if n = N/2 + 1, . . . , N = (−1)n≤N/2[v(θ0)]n. (a) (b) Figure 2.11: Example beam patterns with conventional beamforming for the ∆ channel in a Σ/∆ system. (a) The case θ0 = 0. (b) Contour plot for θ0 ranging from −90◦ to 90◦. 17 2. Theory Figure 2.11 displays normalized example beam patterns for the ∆ channel, similar to ones shown in Figure 2.10. Note that the beam pattern has its minimum for the beam position θ = θ0. Note that conventional beamforming can be either analog or digital. As described in Section 2.5, the Σ/∆ system applies analog beamforming. The multi-channel system, on the other hand, digitizes the phase shifters and consequently applies digital beamforming. 2.6.2 Adaptive beamforming In the case of conventional beamforming described above, the beam positions, and therefore weights w, are fixed. Despite its simplicity, conventional beamforming has its limitations. One significant disadvantage is that it does not account for interfering signals. If an interfering signal is present at a main lobe or sidelobe in the beam pattern, there is a risk of this signal being amplified at the expense of the SOI. Moreover, the coverage is limited given a finite number of beam positions. Even in the case of digital beamforming, where it is possible to control multiple beams with different θ0 simultaneously, there is a decrease in amplitude between the different beams. In a scenario where an SOI is located in such a gap, this could impair the probability to detect it. For data-dependent beamforming, or adaptive beamforming, the complex weights are instead determined by optimizing a certain quantity based on the input data [6]. As indicated by the name, these beamformers can adapt to different signal scenarios, including accounting for interfering signals. As previously stated, the optimal weights are calculated using input data to the sys- tem. However, before analyzing the incoming data, optimal beamformers must be derived. In this thesis, three different optimal beamformers are analyzed: The min- imum power distortionless response (MPDR) beamformer, followed by two different modifications of it: the linear constraint minimum power (LCMP) beamformer and the principal component (PC) beamformer. 2.6.2.1 MPDR beamformer In the MPDR beamformer, the goal is to minimize the total received power, while maintaining a distortionless response in the direction of the SOI [5], [6]. In the time domain, the received signal in each element can be written as x(t) ∈ CN . After processing this signal with weights, w, the output signal, y(t), is given by y(t) = wHx(t). The total output power can then, per definition, be obtained as the expected value of the squared magnitude of y(t) according to Pout = E[|y|2] = E[|wHx|2] = E[wHxxHw] = wHE[xxH]w. 18 2. Theory It is now convenient to define the signal covariance matrix between the different antenna elements as Rx = E[xxH]. This matrix captures the statistical properties of the received signal, and plays an important role in adaptive beamforming. The output power can therefore be written as Pout = wHRxw. To maintain a distortionless response in the desired direction of arrival, θSOI, a constraint of unit gain in this direction can be imposed, according to wHv(θSOI) = 1. This leads to the optimization problemmin w wHRxw s.t. wHv(θSOI) = 1. This problem can be solved using complex Lagrangian relaxation. Introducing La- grangian parameter λL, this gives L(w, λL) = wHRxw + λL(wHv(θSOI) − 1) + λH L(vH(θSOI)w − 1). Note that the conjugate of the constraint is included, since this complex optimiza- tion problem effectively has two constraints: one for the real part and one for the imaginary. Also note that w and wH can be treated as two independent variables. More details about the complex optimization are explained in Appendix A.1. Taking the complex gradient with respect to w gives ∇wL(w, λL) = wHRx + λH LvH(θSOI) = 0. Solving for wH gives wH = −λH LvH(θSOI)R−1 x . Note that the constraint now can be written as wHv(θSOI) = −λH LvH(θSOI)R−1 x v(θSOI) = 1, resulting in λH L = − 1 vH(θSOI)R−1 x v(θSOI) . Note that vH(θSOI)R−1 x v(θSOI) is a quadratic form and therefore a scalar. The optimal weight vector can now be written as wH = vH(θSOI)R−1 x vH(θSOI)R−1 x v(θSOI) . (2.10) By applying this weight, interfering signals can be suppressed using the information encoded in the covariance matrix Rx, while the imposed constraint enables unit gain in the direction of the signal of interest [6]. Note that the numerator in Equation (2.10) describes the phase relation between the signals in the different elements when applying the weights. The denominator is effectively a constant scaling factor. 19 2. Theory 2.6.2.2 LCMP beamformer There are several ways to generalize and develop the MPDR beamformer. One alternative is to include more linear constraints in the optimization problem. This leads to the linear constraint minimum power (LCMP) beamformer [6]. For Nc number of constraints, the problem can be written asmin w wHRxw s.t. wHC = gH, where C ∈ CN×Nc is a constraint matrix and g ∈ CNc is a vector representing the signal gains set by the constraints. Following the same steps as in previous section, the Lagrangian can now be written as L(w, λL) = wHRxw + λL(wHC − gH) + λH L(CHw − g). The complex gradient with respect to w gives ∇wL(w, λL) = wHRx + λH LCH = 0. Solving for wH gives wH = −λH LCHR−1 x , resulting in wHC = −λH LCHR−1 x C = gH, and λH L = −gH(CHR−1 x C)−1. The optimal LCMP weights are therefore given by wH = gH(CHR−1 x C)−1CHR−1 x . In the case of imposing direct constraints on the gain for an SOI and D interfering signals, the constraint matrix can be written as C = [ v(θSOI) v(θ1) v(θ2) . . . v(θD) ] , where θSOI is the direction of arrival for the SOI and θ1, θ2, . . . , θD are the directions of arrival for the D interfering signals. For example, imposing a unit gain constraint for the SOI and null constraints for the interfering signals gives g = [ 1 0 0 . . . 0 ]T . 2.6.2.3 PC beamformer Another variation of the MPDR is the principal component (PC) beamformer [6]. The idea here is to perform an eigendecomposition of Rx according to Rx = UΛUH = N∑ i=1 λiΦiΦH i , 20 2. Theory where U is an N×N matrix with Rx eigenvectors Φi as columns, and Λ is a diagonal matrix with corresponding eigenvalues λi. Since Rx is symmetric, its inverse is given by R−1 x = UΛ−1UH = N∑ i=1 1 λi ΦiΦH i . Using principal component analysis, it is then possible to reduce the dimensionality of the data by identifying the components with the largest variance contribution. This can be done by analyzing the eigenvalues λi. The eigenvectors with the largest corresponding eigenvalues represent the directions in which the data has its largest variance. Note that the variance reflects the signal power of the different compo- nents. For the case with an SOI and D interfering signals, the principal components can be chosen as the eigenvectors with the D + 1 largest corresponding eigenvalues. From this, the signal-plus-interference subspace can be defined as the vector space spanned by the eigenvectors of the principal components, according to US+I = [ Φ1 Φ2 . . . ΦD+1 ] . Similarly, the noise subspace can be defined as the span of the remaining eigen- vectors, corresponding to the N − (D + 1) smallest eigenvalues. Note that for the noise subspace to be non-trivial, this requires D + 1 < N . Assuming the noise is uncorrelated to the signals, the eigendecomposition of R−1 x can now be written as R−1 x = US+IΛ−1 S+IU H S+I + UnoiseΛ−1 noiseU H noise. Because of the symmetry of Rx, its eigenvectors are orthogonal. This implies that the signal-plus-interference subspace and the noise subspace are orthogonal comple- ments. Since the signal steering vector v(θSOI) belong to the signal-plus-interference subspace, vH(θSOI)Φi = 0 for i = D + 2, . . . , N . Consequently, the MPDR weights in Equation (2.10) can be reduced to wH = vH(θSOI)US+IΛ−1 S+IU H S+I vH(θSOI)US+IΛ−1 S+IU H S+Iv(θSOI) = ( vH S+I(θSOI)Λ−1 S+I vH S+I(θSOI)Λ−1 S+IvS+I(θSOI) ) UH S+I, where vS+I(θSOI) = UH S+Iv(θSOI) is steering vector projected to the lower-dimensional signal-plus-interference subspace. In conclusion, the PC beamformer can be seen as a projection of the incoming signals into a (D + 1)-dimensional subspace containing sufficient information to obtain the optimal weights w. This reduction in dimen- sionality can be computationally beneficial. 2.7 Estimation of parameters In the previous section, optimal beamformers are derived. In order to implement these in practice, several parameters need to be estimated from the input data. This includes the signal covariance matrix, Rx, and the directions of arrival, θ. 21 2. Theory 2.7.1 Covariance matrix In practical applications, the signal covariance matrix Rx has to be replaced by an estimate R̂x. To this end, the sample matrix inversion (SMI) method can be applied [5], [6]. Assume K number of N -dimensional discrete time samples x1, x2, . . . , xK , independently sampled from a complex normal distribution with zero mean and covariance matrix Rx, x ∼ CN (0, Rx). The joint probability density, corresponding to the likelihood function, is then given by L(Rx) = K∏ k=1 1 πNdet(Rx) exp(−xH k R−1 x xk), in accordance with Equation (2.6). This results in the log-likelihood log L = −K log det(Rx) − K∑ k=1 xH k R−1 x xk + const = = −K log det(Rx) − tr ( R−1 x K∑ k=1 xkxH k ) + const. Note that the cyclic property of the trace is applied here. Using derivative operations for matrices, the matrix gradient of the log-likelihood is then given by ∇Rx log L = −K(R−1 x )T + ( R−1 x K∑ k=1 xkxH k R−1 x )T = 0. Solving for Rx then gives the maximum likelihood estimator R̂x = 1 K K∑ k=1 xkxH k . Assuming an ergodic stochastic process, the time sample average above converges to the ensemble average E[xxH] as the number of samples increases. Therefore, R̂x converges to Rx as K → ∞. 2.7.2 Directions of arrival The derived formulas for the weights not only depend on the covariance matrix, but also on the steering vectors for the SOI and the interfering signals. Given that the antenna geometry and the wavelength of the incident waves are known, the only unknown parameter in the steering vector is the direction of arrival (DOA). One algorithm commonly used to estimate DOAs is ESPRIT (Estimation of signal parameter via rotational invariance techniques). Like the PC beamformer derived in Section 2.6.2.3, this algorithm exploits the fact that the received signals can be decomposed into a signal-plus-interference subspace and a noise subspace [6]. Following this, the signal-plus-interference subspace can be divided into two subarray subspaces, U 1, containing the first N − 1 antenna elements, and U 2, containing the last N − 1 elements. A similar division can be made for the steering vectors. If 22 2. Theory V denotes the matrix containing the steering vectors, V 1 and V 2 represent the steering vector matrices for the two subarrays. Since the steering vectors span the signal-plus-interference subspace, this gives{ U 1 = V 1T , U 2 = V 2T , where T ∈ C(D+1)×(D+1) is a non-singular matrix. In other words, the steering vectors can be expressed as a linear combination of the eigenvectors of the covariance matrix, which span the signal-plus-interference subspace. Moreover, the fact that there is a constant phase shift between the elements can be utilized by writing V 2 = V 1∆Φ, where ∆Φ is a diagonal matrix containing the phase shift between the elements for the different signals according to ∆Φ = diag [ exp(i∆φ1) exp(i∆φ2) . . . exp(i∆φD+1) ] . Recall that the phase shift between adjacent elements is given by ∆φi = 2πd sin θi λ . Combining these expressions, this results in U 2 = V 2T = V 1∆ΦT = U 1T −1∆ΦT . Defining P = T −1∆ΦT , eigencomposition gives that ∆Φ contains the eigenvalues of P . By solving U 2 = U 1P and computing the eigenvalues of P , the directions of arrival, θi, can therefore be estimated. The equation can for example be solved using the method of least squares, yielding P = (UH 1 U 1)−1UH 1 U 2. The full ESPRIT algorithm can be seen in Algorithm 1 below. Algorithm 1 ESPRIT algorithm for DOA estimation 1: Collect measurements x1, x2, . . . , xK . 2: Estimate the covariance matrix R̂x. 3: Perform eigendecomposition on R̂x according to R̂x = UΛUH. 4: Find the signal-plus-interference subspace US+I. 5: Define U 1 and U 2 as the first and last N − 1 rows of US+I, respectively. 6: Solve U 2 = U 1P for P . 7: Calculate eigenvalues λ1, λ2, . . . , λD+1 of P . 8: Obtain the phase shift from λi = exp(i∆φi) ⇒ ∆φi = arg λi. 9: Calculate the DOA from ∆φi = 2πd sin θi/λ ⇒ θi = arcsin[λ∆φi/(2πd)]. Note that the implementation of this algorithm requires knowledge of the number of signals. This could for example be estimated by analyzing the eigenvalues of the covariance matrix, corresponding to the signal power, and require that their magnitude should be above some threshold. 23 2. Theory 24 3 Implementation This chapter describes the different methods used in this thesis, based on the theory described in the previous chapter. It should be mentioned that no real data has been used in this project, but it has rather been generated using the theoretical models and approximations described below. Recall that the aim here is to estimate the probability of detection for ADS-B signals for different interference scenarios. This is done for digital multi-channel systems applying conventional as well as adaptive beamforming. The results are benchmarked against a two-channel Σ/∆ system using analog conventional beamforming. 3.1 Environment model The simulations are performed in a square grid of area Agrid = 106 km2. The spatial resolution is set to Apixel = 10×10 km2, hence a grid consisting of 10,000 pixels. An example, with fewer pixels, is seen Figure 3.1. In the simulations, all pixels in the grid are systematically iterated over and a target of interest, i.e., an ADS-B signal, is assumed to be located in each pixel. The airborne SSR radar is assumed to be located at the origin. Figure 3.1: An example square grid. Note that this grid only has 100 pixels, compared to 10,000 pixels used in the simulations. The SSR, marked in blue, is placed at the origin. The green pixel with a ”1” indicates that a target, in this case, has successfully been detected here. 25 3. Implementation The maximum distance limited by the horizon is set to 900 km. This is based on an approximation that the radar is placed on a platform flying at an altitude of h = 12 km. If two aircraft are at altitude h, the maximum possible line-of-sight distance, RLOS, is given by ( RLOS 2 )2 + (karREarth)2 = (karREarth + h)2, where REarth ≈ 6371 km is the radius of the Earth [1]. Note that due to atmospheric refraction, radio waves are bent towards to surface of the Earth, thus increasing the effective radius by a factor kar. Assuming normal weather conditions, this factor is set to kar = 4/3 [23]. This gives RLOS = 2 √ 2karREarthh + h2 ≈ 900 km. ADS-B signals are assumed to be transmitted at a rate of fsquitter = 2 Hz. The transponder output power is set to P ′ t,SOI = 21 dBW ≈ 126 W. Note that this is the effective power transmitted from the transponders, including possible gains and internal losses. Similarly, the power from the interfering signals is set to P ′ t,int = 24 dBW ≈ 251 W. This assumption is chosen to account for ”pessimistic scenarios”, such that the performance of the system is rather underestimated than overesti- mated. Regardless, the scenarios investigated remain realistic. The numbers are partly based on ICAO standards, regulating that the lowest allowed transponder output power is 24 dBW [4]. Hence, the effective power for the targets is chosen to include losses of 3 dB relative to the interfering signals. In this study, the interfering signals are assumed to be Mode S Long and Mode S Short replies. As explained in Section 2.2, the number of interfering signals can be seen as a Poisson random variable. In this case, the Poisson parameters are λlong Poi = tlong Poi µlong and λshort Poi = tshort Poi µshort. The time parameters are given by tlong Poi = tADS-B +tint,long = 120+120 µs = 240 µs and tshort Poi = tADS-B +tint,short = 120+64 µs = 184 µs. Analyzing the intensity of interfering signals, µ, the total number of interfering signals per second, this can be written as µ = µtρAgrid, where µt is the number of messages per individual transponder per second and ρ is the transponder density (number of transponders per km2). Introducing γ = µtρ, this parameter describes the number of messages per second per km2. Effectively, since γ is proportional to both the density of transponders and how often they transmit a message, this parameter can be thought of as a degree of disturbance. In the simulations, the γ parameter ranges from γ = 5 × 10−3 s−1 km−2 to γ = 2 × 10−2 s−1 km−2. This interval is chosen based on estimations on ρ and µt. Using Flightradar24, ρ is estimated to be in the order of 1 × 10−3 transponders per km2 [24]. Note that this depends on the geographical location as well the air traffic at the given time. It is a bit challenging to estimate µt accurately, as it depends on several factors such as regulations for message rates and the number of interrogators covering the area. It is also affected by the fact that aircraft can interrogate each other through an airborne collision avoidance system (ACAS). Based on earlier 26 3. Implementation work, µt is estimated to be in the order of 10 messages per transponder per second [11], [12]. Consequently, γ is in the order of 1 × 10−2 s−1 km−2. For simplicity, an equal distribution for Mode S Long and Mode S Short is assumed, meaning γlong = γshort = γ/2. From this, three different values for γ are chosen: 5 × 10−3, 1 × 10−2 and 2 × 10−2 s−1 km−2. The corresponding Poisson parameter for Mode S Long replies are given by λlong Poi = tlong Poi µlong = tlong Poi γlongAgrid. The calculation for Mode S Short is performed analogous to the one for Mode S Long. The different values for γ with corresponding Poisson parameters are shown in Table 3.1. Table 3.1: The three different values used for the message rate per area γ, corre- sponding to the degree of disturbance. The corresponding Poisson parameters λlong Poi and λshort Poi , describing the average number of received Mode S Long and Mode S Short replies, are also presented together with their sum. γ [s−1 km−2] λlong Poi λshort Poi λtot Poi 5 × 10−3 0.60 0.46 1.06 1 × 10−2 1.20 0.92 2.12 2 × 10−2 2.40 1.84 4.24 The other parameters related to the environment model described above are sum- marized in Table 3.2 below. Table 3.2: Values for the parameters used when modelling the physical environment experienced by the simulated SSR system. This includes geographical parameters as well as properties of the incoming signals. Parameter Notation Value Grid size Agrid 1000 × 1000 km2 Spatial resolution Apixel 10 × 10 km2 Altitude h 12 km Maximum line-of-sight distance RLOS 900 km Target effective output power P ′ t,SOI 21 dBW Interfering target effective output power P ′ t,int 24 dBW ADS-B squitter frequency fsquitter 2 Hz Duration of ADS-B signal tADS-B 120 µs Duration of Mode S Long reply tint,long 120 µs Duration of Mode S Short reply tint,short 64 µs 3.2 Antenna model A uniform linear array with N = 6 antenna elements is used in this model. The distance between the elements is set to d = 80 mm. For the cases where conven- tional beamforming is applied, a fixed set of seven beam positions are used, namely 27 3. Implementation [ 0◦ ±15◦ ±30◦ ±45◦ ] . Figure 3.2 displays beam patterns in polar coordinates for the different beam positions. Figure 3.2: Beam patterns in polar coordinates for the seven different beam posi- tions, 0◦, ±15◦, ±30◦ and ±45◦, used for conventional beamforming. The bandwidth of the system is set to W = 10 MHz. The simulated time is set to tsim = 1 s. Note that in the case of analog beamforming, this corresponds to the sweep time, i.e., the time it takes for the antenna to step through the beam positions. The physical area of the antenna is approximated as Aphysical = ((N − 1)d + d)2 = (Nd)2 ≈ 0.23 m2, assuming a spacing of d between N antenna elements and a spacing of d/2 at both ends. The effective aperture area is approximated according to Aeff = Aphysical cos θ in accordance with Equation (2.4). This results in a maximum antenna gain G(θ) ≤ Gmax(θ) = 4π λ2 Aeff(θ) = 4π λ2 (Nd)2 cos θ. Note that the wavelength λ is given by λ = c/f ≈ 27.5 cm, where c is the speed of light and f = 1090 MHz is the frequency of the wave. In order to relate the assumed maximum gain to the beam pattern, B(θ) = wHv(θ), the gain as a function of azimuth is approximated as G(θ) = |B(θ)|2 max θ |B(θ)|2 Gmax(θ). (3.1) From this, the maximum possible detection range can be calculated using Friis trans- mission equation (2.5), and setting the received power at the minimum detection level (MDL). This gives Rmax(θ) = λ 4π √ P ′ tGmax(θ) MDL = λ 4π 10(P ′ t [dBW]+G [dBi] max (θ)−MDL[dBW])/20. The minimum signal-to-interference threshold is set to SIRmin = 6 dB. The mini- mum detection level set to MDL = −115 dBW, based on the ICAO requirement of minimum detection level for interrogator systems [4]. The minimum signal-to-noise is set to SNRmin = 10 dB, which gives the total noise variance Pn = 2σ2 n = MDL[dBW] − SNR[dB] min = −125 dBW. 28 3. Implementation The parameters related to the antenna model described are summarized in Table 3.3. Table 3.3: Values for the parameters used when modelling the simulated ULA antenna in the SSR system. Parameter Notation Value Number of antenna elements N 6 Distance between elements d 80 mm Physical area Aphysical 0.23 m2 Beam positions θ0 [ 0◦ ±15◦ ±30◦ ±45◦ ] Bandwidth W 10 MHz Simulated time tsim 1 s Wavelength λ 27.5 cm Maximum gain in antenna normal Gmax(0◦) 15.8 dBi Maximum range in antenna normal Rmax(0◦) 850 km Signal-to-interference threshold SIRmin 6 dB Minimum detection level MDL −115 dBW Signal-to-noise threshold SNRmin 10 dB Noise power Pn −125 dBW 3.3 Signal model The ADS-B signal of interest is generated together with potential interfering signals. The preamble is simulated given the ICAO standard for Mode S. The data block that follows is then randomized with ones and zeros. Recall that the signal waveform a(t) is given by Equation (2.1). A signal reaching the antenna is then modelled as s(t) = √ Pra(t) exp(i2πft), where Pr is the received power calculated using Friis transmission equation (2.5) and f is the frequency of the signal. As per ICAO standards, f = 1090 ± 1 MHz [4]. For simplicity, a downconversion to f = 0 ± 1 MHz is assumed here. Assuming D interfering signals, there are D + 1 signals in total. The total incoming signal can then be written as the vector s(t) ∈ CD+1. The received signal in each antenna element is then given by x(t) = As(t) + n(t), where x(t) ∈ CN is the received signal [6]. n(t) ∈ CN is noise, which is assumed to be sampled from a complex normal distribution as n ∼ CN (0, 2σ2 nIN). A ∈ CN×(D+1) is a matrix containing the steering vectors for the different signals according to A = [ v(θSOI) v(θ1) v(θ2) . . . v(θD) ] . Note that the incoming signals are assumed to be in the far field relative to the receiving antenna. Hence, the plane wave approximation is assumed to hold. The 29 3. Implementation validity of this can be analyzed by looking at dF > 2D2 F/λ, where dF is the Fraunhofer distance, describing the limit between near and far field [16]. DF is the largest dimension of the receiving antenna, which in this case is Nd. Using the parameters defined above, this gives dF ≈ 1.7 m. Since the signals are coming from targets located tens or hundreds of kilometers away, they are consequently in the far field relative to the receiving antenna. Moreover, the signals are assumed to be narrowband signals. This assumption holds if W∆Tmax ≪ 1, where W is the bandwidth of the system and ∆Tmax is the maximum travel time between two elements in the ULA [6]. W = 10 MHz and ∆Tmax = (N − 1)d/c gives W∆Tmax ≈ 0.01. Therefore, the narrowband assumption can be regarded as valid. 3.4 Analysis of adaptive beamforming Implementing adaptive beamforming, the three different optimal beamformers de- scribed in Section 2.6.2 are analyzed: the MPDR, the LCMP and the PC. For the LCMP, a unit gain constraint for the SOI is imposed, with null constraints for the interfering signals. To obtain the statistically optimal weights, the covariance matrix R̂x is estimated. The sample matrix inversion method described in Section 2.7.1 is applied here. Note that a sample frequency of 10 MHz is assumed, corresponding to the bandwidth of the system. This is in accordance with the Nyquist−Shannon sampling theorem, stating that fsampling > 2W , where W is the bandwidth [5]. Since both a real and an imaginary data point are sampled in this case, the effective sampling rate is twice as high, thus satisfying the theorem. Estimates for the directions of arrival, and consequently the steering vectors, for the different signals are obtained by implementing the ESPRIT algorithm described in Section 2.7.2. To obtain the dimension of the signal-plus-interference subspace, used for both the ESPRIT algorithm and in the implementation of the PC beamformer, the total number of signals, D + 1, has to be estimated. This is done by requiring that the eigenvalues of R̂x should be 1 dB above the noise floor. Note that the MDL (10 dB over the noise floor) would be too strict a requirement here, as the SOI could still be amplified after obtaining the optimal weights. This method of determining the subspace dimension effectively leads to the assumption that the total number of signals is at most the number of elements, D + 1 ≤ N . The plausibility of this assumption is discussed further in Section 5.2. 3.5 Probability of detection The probability of detecting an ADS-B signal of interest is modelled based on two main criteria. In order to successfully detect the signal, and therefore be able to decode the information it contains, the following should hold: 30 3. Implementation 1. The signal power of the target of interest at the receiving antenna should be strong enough to reach SNRmin above the noise floor, i.e., reach the minimum detection level (MDL). 2. The signal power of the target at the receiving antenna should be SIRmin above the maximum of the sum of the signal power of the interfering signals. By taking these two criteria into account, both the signal power of the SOI and its relation to the strength of the interfering signals are considered. From these two criteria, a detection binary random variable Xd can be defined as Xd = 1 if target deteced, 0 if target not detected. Mathematically, this assumption can, in linear units, be expressed as the joint prob- ability P(Xd = 1) = P PSOI(rSOI, θSOI) > MDL, PSOI(rSOI, θSOI) max t [∑D k=1 Pk(rk, θk)](t) > SIRmin  . Expressed more explicitly, in terms of the effective radiated power P ′ t from the transmitting antennas in the transponders, the path loss Lp(r) and the gain G(θ) of the receiving antenna, this can be written as P(Xd = 1) = P P ′ t,SOIG(θSOI) Lp(rSOI) > MDL, P ′ t,SOIG(θSOI)/Lp(rSOI) max t [∑D k=1 P ′ t,kG(θk)/Lp(rk)](t) > SIRmin  . From this, the probability of detection can be defined as pd = P(Xd = 1). This probability is estimated for every point (x, y) in the spatial grid, using a sample mean p̂d from Monte Carlo simulations. This is descried in more detail in Section 3.6 below. 3.6 Monte Carlo simulations To estimate the probability of detection for every point in the grid, Monte Carlo simulations in MATLAB are performed. This method is preferred given the consid- erable complexity of the problem, with several random variables involved, making analytical solutions impractical to achieve. The random variables involved in the simulations are: • The number of interfering signals: A Poisson distribution is assumed according to Equation (2.2). This gives Dlong ∼ Poi(λlong Poi ) interfering Mode S Long signals and Dshort ∼ Poi(λshort Poi ) interfering Mode S Short signals. • Spatial position of interfering signals: Uniformly distributed points in a half- circle with radius RLOS, centered at the SSR in the origin. • Temporal position of interfering signals: Uniformly distributed in the interval [120 − tint, 240] µs, where tint is the duration of the interfering signal. Note that this assumes that the SOI begins at t = 120 µs. 31 3. Implementation • Downconverted signal frequency: Normal distribution assumed according to f ∼ N (0, σ2 f ). Setting σf = 1/3 MHz assures the requirement of 0 ± 1 MHz with 99.7 % certainty. • Data block in the waveforms: Uniformly randomized order of ones and zeros. • Noise: Complex normally distributed according to n ∼ CN (0, 2σ2 nIN). For the two-channel Σ/∆ system, there is an additional random variable, namely the beam position which is active for the first ADS-B squitter. If there are Nbeams beam positions and Nsquitters squitters during the sweep time, the first squitter occurs somewhere during the first ∆tsquitter = 1/fsquitter = 0.5 seconds, corresponding to one of the first ⌈Nbeams∆tsquitter/tsim⌉ = ⌈Nbeams/Nsquitters⌉ beam positions. The active beam position for the following squitter is then determined from this using the fixed squitter frequency of 2 Hz. Pseudocode for the Monte Carlo simulations for the Σ/∆ system can be found in Algorithm 2 below. Algorithm 2 Monte Carlo simulation, Σ/∆ system for all iterations do for all positions in the grid do Assume a target is located at the given position. Calculate rSOI and θSOI. Set Xd = 0 if rSOI > Rmax(θSOI). Draw beam position for first ADS-B squitter. Choose the beam positions for the next squitters accordingly. for all beam positions θ0 with a squitter do Draw Dlong ∼ Poi(λlong Poi ) interfering Mode S Long signals. Draw Dshort ∼ Poi(λshort Poi ) interfering Mode S Short signals. Draw positions in space and time for the D = Dlong + Dshort interfering signals. Calculate GΣ(θ) ∝ |vH(θ0)v(θ)|2 and G∆(θ) ∝ |(−1)n≤N/2vH(θ0)v(θ)|2 for SOI and interfering signals. Calculate received power PSOI from SOI and Pk from interfering signals. Calculate SIR. if PSOI > MDL and SIR > SIRmin for Σ or ∆ channel then Target detected, Xd = 1. break end if end for end for end for Pseudocode for the Monte Carlo simulations for the multi-channel system with con- ventional beamforming can be found in Algorithm 3. Compared to Algorithm 2, there is no ∆ channel here, as explained in Section 2.5. Since there are multiple in- dependent beams active simultaneously, there is no random variable for which beam 32 3. Implementation position is active when the first squitter from the ADS-B target arrives. Thus, for every squitter, there is a chance to detect the signal in all beam positions. Algorithm 3 Monte Carlo simulation, multi-channel system with conventional beamforming for all iterations do for all positions in the grid do Assume a target is located at the given position. Calculate rSOI and θSOI. Set Xd = 0 if rSOI > Rmax(θSOI). for all squitters do Draw Dlong ∼ Poi(λlong Poi ) interfering Mode S Long signals. Draw Dshort ∼ Poi(λshort Poi ) interfering Mode S Short signals. Draw positions in space and time for the D = Dlong + Dshort interfering signals. for all beam positions θ0 do Calculate G(θ) ∝ |vH(θ0)v(θ)|2 for SOI and interfering signals. Calculate received power PSOI from SOI and Pk from interfering signals. Calculate SIR. if PSOI > MDL and SIR > SIRmin then Target detected, Xd = 1. break end if end for end for end for end for Pseudocode for the Monte Carlo simulations for the multi-channel system with adap- tive beamforming can be seen in Algorithm 4. The noticeable difference compared to Algorithms 2 and 3 is that there are no predetermined beam positions θ0. In- stead, the optimal weight vector, w, is determined using the input signal x(t). Even though the digital multi-channel system allows for several simultaneous beam posi- tions, only one is used here, as the obtained w should be the statistically optimal one. Note that a DOA estimation θ̂ is obtained for all distinguishable signals, each yield- ing a separate w. Instead of looping through all w, the one resulting in the smallest absolute error relative to the true value is chosen. That is, θ̂SOI = min |θ̂i − θSOI|, where i counts the number of estimated signals. 33 3. Implementation Algorithm 4 Monte Carlo simulation, multi-channel system with adaptive beam- forming for all iterations do for all positions in the grid do Assume a target is located at the given position. Calculate rSOI and θSOI. Set Xd = 0 if rSOI > Rmax(θSOI). for all squitters do Draw Dlong ∼ Poi(λlong Poi ) interfering Mode S Long signals. Draw Dshort ∼ Poi(λshort Poi ) interfering Mode S Short signals. Draw positions in space and time for the D = Dlong + Dshort interfering signals. Simulate the input signal x(t). Estimate the sample covariance matrix R̂x and directions of arrival. Determine the optimal weight vector w. Calculate G(θ) ∝ |wHv(θ)|2 for SOI and interfering signals. Calculate received power PSOI from SOI and Pk from interfering signals. Calculate SIR. if PSOI > MDL and SIR > SIRmin then Target detected, Xd = 1. break end if end for end for end for When performing Monte Carlo simulations, it is of great importance to determine the number of iterations sufficient to obtain good statistics. The simulation has to be run long enough such that a large part of the probability space can be explored. On the other hand, the computational time must also be considered. To estimate the proper number of iterations, tests are performed for five randomly selected points in the grid, which can be found in Table 3.4. If Xd,1, Xd,2, . . . , Xd,n, obtained from each run in the Monte Carlo simulation, are assumed to be independent and identically distributed random variables, the law of large numbers gives the sample average p̂d = Xd = 1 n (Xd,1 + Xd,2 + · · · + Xd,n). As n → ∞, the sample average converges to the actual expected value of Xd, according to Xd → E(Xd,1) = E(Xd,2) = · · · = pd [13]. Therefore, the sample average can be used as an estimator for pd. Analyzing the variance, it is given by Var(Xd) = Var(p̂d) = σ2 d n , 34 3. Implementation where Var(Xd,i) = σ2 d. It can be noted that Var(p̂d) → 0 as n → ∞. As the number of Monte Carlo samples increases, the estimation for pd becomes less uncertain. (a) (b) Figure 3.3: (a) The standard deviation of p̂det as a function of the number of Monte Carlo iterations for five randomly selected points. Note that p̂d corresponds to the sample mean Xd of the number of Monte Carlo iterations. (b) The simulation time per pixel as a function of the number of iterations. Figure 3.3(a) displays the standard deviation of p̂d, √ Var(p̂d), as a function of the number of Monte Carlo iterations for the five randomly selected points. Note that the standard deviation is approximated based on ten separate runs for each point. Figure 3.3(b) displays the average simulation time per pixel as a function of the number of iterations. Table 3.4: Coordinates and standard deviations of p̂det for the five randomly se- lected points. The results are displayed for 10, 100, 1,000 and 10,000 iterations. Standard deviation of p̂d Point Coordinate [km] 101 102 103 104 1 (170, 40) 8.4 × 10−2 7.4 × 10−2 1.2 × 10−2 3.8 × 10−3 2 (330, 660) 1.0 × 10−1 3.1 × 10−2 1.1 × 10−2 2.8 × 10−3 3 (10, 180) 4.2 × 10−2 2.8 × 10−2 4.1 × 10−3 2.0 × 10−3 4 (−10, 240) 1.3 × 10−1 2.4 × 10−2 6.0 × 10−3 3.2 × 10−3 5 (−40, 500) 1.8 × 10−1 4.8 × 10−2 1.1 × 10−2 5.5 × 10−3 Table 3.4 displays the coordinates for the different points, together with the approx- imated standard deviations for the different number of Monte Carlo iterations. From Figure 3.3 and Table 3.4, it can be noted that the standard deviation is in the order of around 10−2 for n = 103 iterations. For n = 104 iterations, the standard deviation is reduced with around 60 %, while the simulation time increases with nearly 900 %. 35 3. Implementation The standard deviation can be used to estimate a confidence interval [13]. Assuming a standard deviation of the sample mean p̂d of 0.01, a 95 % confidence interval is given by [ p̂d − 1.96 σd√ n , p̂d + 1.96 σd√ n ] ≈ p̂d ± 0.02. For this study, this is regarded as a sufficient level of accuracy for estimating the probability in each point of the grid. With regards to this, as well as the simulation time, 1,000 Monte Carlo iterations are deemed sufficient. 36 4 Results and analysis In this chapter, the results obtained from implementing the methods described in the previous chapter are presented. 4.1 Probability of detection Figures 4.1, 4.2, 4.3, 4.4 and 4.5 display the estimated probability of detection p̂d for each pixel in the grid for the different algorithms implemented: analog conventional Σ/∆, digital conventional multi-channel (CMC), adaptive minimum power distor- tionless response (MPDR) multi-channel, adaptive linear constraint minimum power (LCMP) multi-channel and adaptive principal component (PC) multi-channel. A summary of these algorithms can be seen in Table 4.1 below. Table 4.1: The different algorithms analyzed. The second column describes the number of digitized channels. In the case of Σ/∆, the output Σ and ∆ channel are the two digital channels. In the other cases, all six antenna elements are digitized separately. The third column describes whether analog or digital phase shifters are used, resulting in different types of beamforming. Algorithm Digitized channels Phase shifters Beamforming type Σ/∆ 2 Analog Conventional CMC 6 Digital Conventional MPDR 6 Digital Adaptive LCMP 6 Digital Adaptive PC 6 Digital Adaptive The configuration of analog and digital systems is described in Section 2.5. Con- ventional and adaptive beamforming are described in Section 2.6. The results are presented for three different values for the disturbance parameter γ: 5 × 10−3, 1 × 10−2 and 2 × 10−2 s−1 km−2. As seen in Table 3.1, this corresponds to expected values of approximately 1.1, 2.1 and 4.2 number of interfering signals, respectively. 37 4. Results and analysis (a) (b) (c) Figure 4.1: Estimated probability of detection for the Σ/∆ system with analog conventional beamforming. (a) The case γ = 5 × 10−3 s−1 km−2. (b) The case γ = 1 × 10−2 s−1 km−2. (c) The case γ = 2 × 10−2 s−1 km−2. Red tones indicate higher probability, while blue tones indicate a lower probability. Note that the maximum line-of-sight distance at 900 km is marked by a white line. 38 4. Results and analysis (a) (b) (c) Figure 4.2: Estimated probability of detection for the multi-channel system with digital conventional beamforming (CMC). (a) The case γ = 5×10−3 s−1 km−2. (b) The case γ = 1 × 10−2 s−1 km−2. (c) The case γ = 2 × 10−2 s−1 km−2. Red tones indicate higher probability, while blue tones indicate a lower probability. Note that the maximum line-of-sight distance at 900 km is marked by a white line. 39 4. Results and analysis (a) (b) (c) Figure 4.3: Estimated probability of detection for the multi-channel system with adaptive MPDR beamforming. (a) The case γ = 5 × 10−3 s−1 km−2. (b) The case γ = 1 × 10−2 s−1 km−2. (c) The case γ = 2 × 10−2 s−1 km−2. Red tones indicate higher probability, while blue tones indicate a lower probability. Note that the maximum line-of-sight distance at 900 km is marked by a white line. 40 4. Results and analysis (a) (b) (c) Figure 4.4: Estimated probability of detection for the multi-channel system with adaptive LCMP beamforming. (a) The case γ = 5 × 10−3 s−1 km−2. (b) The case γ = 1 × 10−2 s−1 km−2. (c) The case γ = 2 × 10−2 s−1 km−2. Red tones indicate higher probability, while blue tones indicate a lower probability. Note that the maximum line-of-sight distance at 900 km is marked by a white line. 41 4. Results and analysis (a) (b) (c) Figure 4.5: Estimated probability of detection for the multi-channel system with adaptive PC beamforming. (a) The case γ = 5 × 10−3 s−1 km−2. (b) The case γ = 1 × 10−2 s−1 km−2. (c) The case γ = 2 × 10−2 s−1 km−2. Red tones indicate higher probability, while blue tones indicate a lower probability. Note that the maximum line-of-sight distance at 900 km is marked by a white line. 42 4. Results and analysis To numerically analyze how the probability of detection depends on the γ parameter, as well as the distance to the SSR in the origin, the probability can be calculated in different subareas. Three subareas are defined by the distance from the SSR to the target: 0−300 km, 300−600 km and 600−900 km. The three regions are displayed in Figure 4.6. Figure 4.6: The three different subareas defined by the distance between the target and the SSR. Region 1: 0−300 km. Region 2: 300−600 km. Region 3: 600−900 km. As mentioned in Section 3.5, the probability of detection in a point (x, y) is given by pd = P(Xd = 1|(x, y)), where Xd = 1 if target deteced, 0 if target not detected. Denote the probability of detection, given a target located in Region k, as pd(k) = P(Xd = 1|k). From the law of total probability, this can be calculated according to pd(k) = P(Xd = 1|k) = ∑ (xi,yi) P(Xd = 1|k, (xi, yi))P((xi, yi)|k) = = ∑ (xi,yi) ∈ Region k P(Xd = 1|(xi, yi))P((xi, yi)|k). Note that for a uniform distribution of targets, P((xi, yi)|k) = 1/(Number of pixels in Region k). Hence, the total probability of detection given a target in a specific region can be seen as an average probability [13]. Figure 4.7 displays the average probability of detection as a function of the γ pa- rameter for the different systems in the three different subareas. The probability 43 4. Results and analysis values are also displayed in Table 4.2. As