Radio signal DOA estimation : Implementing radar signal direction estimation on an FPGA.

Linköpings universitet SE–581 83 Linköping

Linköping University | Department of Electrical Engineering

Master thesis, 30 ECTS | Elektroteknik

2019 | LiTH-ISY-EX--19/5199--SE

Radio signal DOA estimation

Implementing radar signal direction estimation on an FPGA.

Alfred Patriksson

Supervisor : Francis Görmarker Examiner : Kent Palmkvist

Upphovsrätt

Detta dokument hålls tillgängligt på Internet – eller dess framtida ersättare – under 25 år från publiceringsdatum under förutsättning att inga extraordinära omständigheter uppstår. Tillgång till dokumentet innebär tillstånd för var och en att läsa, ladda ner, skriva ut enstaka kopior för enskilt bruk och att använda det oförändrat för ickekommersiell forskning och för undervisning. Överföring av upphovsrätten vid en senare tidpunkt kan inte upphäva detta tillstånd. All annan användning av dokumentet kräver upphovsmannens medgivande. För att garantera äktheten, säkerheten och tillgängligheten finns lösningar av teknisk och admin-istrativ art. Upphovsmannens ideella rätt innefattar rätt att bli nämnd som upphovsman i den omfattning som god sed kräver vid användning av dokumentet på ovan beskrivna sätt samt skydd mot att dokumentet ändras eller presenteras i sådan form eller i sådant sam-manhang som är kränkande för upphovsmannens litterära eller konstnärliga anseende eller egenart. För ytterligare information om Linköping University Electronic Press se förlagets hemsida http://www.ep.liu.se/.

Copyright

The publishers will keep this document online on the Internet – or its possible replacement – for a period of 25 years starting from the date of publication barring exceptional circum-stances. The online availability of the document implies permanent permission for anyone to read, to download, or to print out single copies for his/hers own use and to use it unchanged for non-commercial research and educational purpose. Subsequent transfers of copyright cannot revoke this permission. All other uses of the document are conditional upon the con-sent of the copyright owner. The publisher has taken technical and administrative measures to assure authenticity, security and accessibility. According to intellectual property law the author has the right to be mentioned when his/her work is accessed as described above and to be protected against infringement. For additional information about the Linköping Uni-versity Electronic Press and its procedures for publication and for assurance of document integrity, please refer to its www home page: http://www.ep.liu.se/.

c

Abstract

This master’s thesis covers the design and implementation of a monopulse direction of arrival (DOA) estimation algorithm on an FPGA. The goal is to implement a complete system that is capable of estimating the bearing of an incident signal. In order to determine the estimate quality both a theoretical and practical noise analysis of the signal chain is performed. Special focus is placed on the statistical properties of the transformation from I/Q demodulated signals with correlated noise to a polar representation. The pros and cons for three different methods of calculating received signal phasors are also covered.

The system is limited to two receiving channels which constrains this report to a 2D analysis. In addition the used hardware is limited to C-band signals.

We show that an FPGA implementation of monopulse techniques is definitely viable and that an SNR higher than ten dB allows for a gaussian approximation of the polar rep-resentation of an I/Q signal.

Glossary

LP/BP/HP Low- Band- and High-pass. CW Continuous Wave

ADC Analog to Digital Converter. LNA Low-Noise Amplifier.

RF Radio Frequency. (Signal before demodulation). BB Baseband. (Signal after demodulation).

SNR Signal to Noise Ratio X Stochastic variable

x Vector variable ˆx Variable estimate xn Quantity x of channel n.

xi In-phase component of signal x. xq Quadrature component of signal x.

j Complex unit: j2=´1 c Speed of light in a vacuum.

ρ Polar coordinate distance from origin

θ Polar coordinate angle - counter-clockwise from positive x-axis φ Signal phase

fs ADC sampling rate fn Nyquist frequency

fc Carrier frequency of modulated signal fb Baseband frequency of demodulated signal x(t) Transmitted signal to estimate direction of. y(t) Complex Signal in(t) +j ¨ qn(t)at receiver n. z(t) Polar conversion of y(t) =ρexp(jφ)

e(t) Additive noise.

G(θ) Antenna gain - θ relative main lobe

HB_LPtu Ideal low-pass filter with cutoff frequency B H_HPB tu Ideal high-pass filter with cutoff frequency B

Contents

Abstract iii

Glossary iv

Contents v

List of Figures vii

List of Tables ix 1 Introduction 1 1.1 Motivation . . . 1 1.2 Aim . . . 1 1.3 Research questions . . . 1 1.4 Delimitations . . . 1 2 Theory 2 2.1 Monopulse Techniques . . . 2 2.2 Phasor Calculations . . . 10 2.3 Noise Analysis . . . 12 3 Hardware 31 3.1 MikTran . . . 31 3.2 Platform . . . 32 3.3 Test Equipment . . . 33 4 Method 34 4.1 ADC Calibration . . . 34 4.2 Performance measurements . . . 35

4.3 Antenna directivity measurement . . . 36

4.4 Calibration Routine . . . 36 5 Implementation 37 5.1 FPGA . . . 39 5.2 MATLAB . . . 40 6 Results 41 6.1 ADC Calibration . . . 41 6.2 Noise Measurements . . . 42

6.3 Antenna Directivity Measurements . . . 43

6.4 Calibration . . . 45

6.5 Direction Estimates . . . 47

7.1 Results . . . 57 7.2 Method . . . 57 7.3 The work in a wider context . . . 57

8 Conclusion 59

List of Figures

2.1 Wave properties in near and far-field. . . 3

2.2 Squinted antennas. . . 4

2.3 Ideal case incident signal amplitude. . . 5

2.4 Ideal case amplitude sum and difference signals. . . 6

2.5 Ideal case amplitude ratio. . . 7

2.6 Difference in incident signal travel distance. . . 7

2.7 Inverse sine ambiguity with Dλ=1.25 and∆φ=´ 3 4π. . . 9

2.8 Complex signal representation of I/Q data. . . 10

2.9 Comparison of Görtzel and DTFT binning errors. . . 11

2.10 Görtzel Algorithm relative error for N=1024 . . . 12

2.11 Mathematical representations of elliptical noise. . . 14

2.12 Monte Carlo simulation of uncorrelated pure WGN in polar coordinates. . . 17

2.13 Monte Carlo simulation of correlated pure WGN in polar coordinates. . . 18

2.14 Example low SNR WGN signal distribution. . . 19

2.15 Correlation induced amplitude bias. . . 20

2.16 Correlation induced maximum amplitude bias. . . 21

2.17 Correlation induced phase bias. . . 22

2.18 Maximum correlation induced phase bias. . . 23

2.19 Comparison of polar and cartesian WGN. SNR = {3,6,10} dB . . . 25

2.20 Overlap between polar and Cartesian normal distributions. . . 26

2.21 PDF of discretized ˆN . . . 29

2.22 Probability of correctly resolving phase ambiguity. . . 30

3.1 Picture of the MikTran circuit board. . . 31

3.2 System with mounted horn antennas. . . 32

3.3 System with mounted spiral antennas. . . 33

3.4 Antenna used as transmitter during measurements. . . 33

5.1 High level system architecture. . . 37

5.2 LVDS packet structure . . . 38

5.3 Network packet structure . . . 38

5.4 FPGA Module Structure . . . 39

6.1 Voltage calibration factor γ. . . . 41

6.2 Voltage calibration factor γ at maximum input power. . . . 42

6.3 Inherent noise spectrum of system. . . 42

6.4 Measured gain and amplitude ratio functions for horn antenna pair. . . 44

6.5 Measured gain and amplitude ratio functions for spiral antenna pair. . . 45

6.6 Horn antenna phase difference measurements and separation distance estimate. . 46

6.7 Spiral antenna phase difference measurements and separation distance estimate. . 46

6.8 Horn antenna direction estimate at -10dB SNR . . . 47

6.10 Horn antenna direction estimate at 5dB SNR . . . 49

6.11 Horn antenna direction estimate at 10dB SNR . . . 50

6.12 Horn antenna direction estimate at 20dB SNR . . . 51

6.13 Spiral antenna direction estimate at -5dB SNR . . . 52

6.14 Spiral antenna direction estimate at 0dB SNR . . . 53

6.15 Spiral antenna direction estimate at 5dB SNR . . . 54

6.16 Spiral antenna direction estimate at 10dB SNR . . . 55

List of Tables

1

Introduction

1.1

Motivation

The radio spectrum is a finite resource shared by all. Engineers are constantly striving for higher data rates which in turn causes a need to utilize this resource more effectively. One way to achieve this is to limit signal spatial overlap through methods such as beamforming. In order to effectively do this one needs to be able to efficiently and accurately determine what direction a radio signal is being transmitted from.

Being able to cheaply estimate DOA of radio signals also makes localization tools more available for things like emergency services. It can also allow easy location of dangerous in-terference sources like GPS jammers which have become increasingly commercially available for low costs.

1.2

Aim

The goal of this thesis project is to implement a 2D direction of arrival (DOA) estimator of RF pulses in the C-band. The project is to be implemented on an existing platform known as "MikTran". The resulting implementation should avoid relying on specific protocol properties of the incoming radio signals and instead aim to be as general as possible.

1.3

Research questions

1. Is it possible to combine amplitude and phase method estimates in a synergistic way? 2. How do the implementation restrictions limit the estimation quality?

3. What geometric position of the antennas is optimal? (Distance and squint angle).

1.4

Delimitations

Because of restrictions in the implementation platform this project is limited to a two-antenna array. This in turn means that only 2D DOA estimation is possible and we therefore assume that any elevation differences between transmitters and receivers are negligible. The analysis is also limited to a single non-moving target.

2

Theory

2.1

Monopulse Techniques

It is possible to unambiguously estimate signal direction from a single received pulse if some criteria are met. This section discusses some of the standard methods, their drawbacks as well as how they can be combined.

The basic structure of the problem is that two antennas measure some transmitted sig-nal at two different points in space. In order to make the problem more tractable this report assumes that the distance between the antennas is far smaller than the distance to the trans-mitter. This allows for the plane wave signal assumption as shown in fig. 2.1.

This report covers pulsed sinusoidal signals. This means that the general structure of the transmitted signals follow the format described by eq. (2.1) where P(t)is a square pulse train with unknown parameters.

x(t) =C ¨ sin(ωct+φx)¨P(t) (2.1) None of the signal parameters C, ωcand φxare known. Neither is the behaviour of the pulsing function P(t).

The signal measured at the receivers is dependent on the incident angle θ and the unkown distance dnbetween the transmitter and receiver n according to eq. (2.2).

τn= dn c yn(t) =Gn(θ)¨xn(t ´ τn) yn(t) =Gn(θ)¨C sin(ωc(t ´ τn) +φx)¨P(t ´ τn) yn(t) =Gn(θ)¨C sin(ωct+ (φx´ ωcτ))¨P(t ´ τn) (2.2)

The signal quantities Anand φnmeasurable at the receivers is provided by the definitions in eq. (2.3).

φn =φx´ ωcτn+2πNn , NnPZ An =C ¨ Gn(θ)

yn(t) =An¨sin(ωct+φn)

2.1. Monopulse Techniques

The properties measurable at the receivers are received sinusoid amplitudes A ¨ Gn(θ)and

phases φn. These properties depend on the unknown quantities θ and dn and are therefore not immediately useful without further processing. Note that the phase measurements are ambiguous! This is because of the 2π-periodicity of the sin function.

Figure 2.1: Wave properties in near and far-field.

Amplitude Method

By using two antennas one can calculate the received signal direction through the received signal amplitudes. The primary problem is eliminating the unknown source signal amplitude A to make the problem solely dependent on the antenna gain functions G(θ). Some criteria

have to be fulfilled in order for this problem to be tractable.

The antenna gain functions Gn(θ)are assumed to be known along with some estimate of

the measured signals amplitudes An from eq. (2.3). This is a system of two equations with two unknowns: θ and Ax.

A1=Ax¨G1(θ)ô θ=G₁´1 A1 Ax  A2=Ax¨G2(θ)ô θ=G₂´1 A2 Ax  (2.4)

This system is solvable when G´1n exist, are unambiguous, and not identical.

In most practical situations two identical antennas are used. At first glance this makes the problem unsolvable since if G1”G2then the equations in eq. (2.4) are identical and the system is degenerate. This can be overcome by fixing the antennas at an angle relative to the boresight axis as in fig. 2.2.

2.1. Monopulse Techniques

θs

θs

D

Figure 2.2: Squinted antennas.

This squint angle θsreshapes eq. (2.4) to eq. (2.5).

A1=Ax¨G(θ ´ θs) A2= Ax¨G(θ+θs)

(2.5)

The idealized case with full sidelobe suppression and a completely noise-free signal re-sults in amplitudes similar to those shown in fig. 2.3. The peaks of the antennas correspond to the squint angles in fig. 2.2.

2.1. Monopulse Techniques -4/₈ -3/₈ -2/₈ -1/₈ 1/₈ 2/₈ 3/₈ 4/₈ [Rad] 0.2 0.4 0.6 0.8 1 Signal amplitude _A 1 A 2

Figure 2.3: Ideal case incident signal amplitude.

A naive approach is to immediately divide the amplitudes A1and A2in order to eliminate the unknown Ax. This leaves inverting the ratio between the gain functions. While this is theoretically doable in the ideal noise free case the nonlinearity of the equation makes it unsuitable for estimation purposes.

Instead one can form a ratio of the sum and difference signals as described by eq. (2.6) [8]. D=A1´A2 S=A1+A2 R(θ) = D S = A1´A2 A1+A2 = Ax Ax¨ D1(θ)´D2(θ) D1(θ) +D2(θ) (2.6)

The quantity R(θ)behaves much better since dividing by the sum normalizes the quantity

to the range[´1, 1]as shown in fig. 2.5. This function is approximately linear around θ =0 with the slope k being solely dependent on the antenna gain functions. If the antennas have been characterized it is possible to calculate the slope and as a consequence one can invert R(θ)like in eq. (2.7).

ˆθR=R´1(A1, A2)«k´1¨

A1´A2

A1+A2 (2.7)

The linear estimate tends to underestimate the incident signal angle. However, is not neces-sary to use a linear approximation as long as the function R(θ)is invertible. This is generally

2.1. Monopulse Techniques -4/₈ -3/₈ -2/₈ -1/₈ 1/₈ 2/₈ 3/₈ 4/₈ [Rad] 0.5 1 1.5 Sum Signal -4/ 8 -3_/ 8 -2_/ 8 -1_/ 8 1_/ 8 2_/ 8 3_/ 8 4_/ 8 [Rad] -0.8 -0.6 -0.4 -0.2 0.2 0.4 0.6 0.8 Difference Signal

2.1. Monopulse Techniques -4/₈ -3/₈ -2/₈ -1/₈ 1/₈ 2/₈ 3/₈ 4/₈ [Rad] -1 -0.8 -0.6 -0.4 -0.2 0.2 0.4 0.6 0.8 1 Ratio R( )

Figure 2.5: Ideal case amplitude ratio.

Phase Method

The phase method exploits the phase delay caused by the travel time τnbeing longer to one of the antennas. With a known distance D between the receivers it is a simple matter of finding the distance r marked in fig. 2.6 for any incident angle θ.

D

θ θ

r

2.1. Monopulse Techniques

Forming a right angle triangle with the distance D between the receivers as the hy-potenuse one gets expression (2.8) for the difference in travel distance. With the distance known it is easy to calculate the travel time difference∆τ since r=c ¨∆τ.

r=D sin(θ)ô∆τ= D

c ¨sin(θ) (2.8)

By introducing the difference between the receiver phases through the definition in eq. (2.3) one gets the expression in eq. (2.9).

∆φ=φ1´ φ2 (2.9)

Carrying this through and introducing∆τ one can make the following simplifications. ∆τ =τ2´ τ1ñ

∆φ= (φx´ ωcτ1+2πN1)´(φx´ ωcτ2+2πN2)ô ∆φ=ωc¨(τ2´ τ1) +2π(N2´N1)ô

This ends up with a different way of expressing∆φ in equation eq. (2.10).

∆φ=ωc¨∆τ+2πN , N PZ (2.10)

This neatly eliminates the unknown absolute phase term φxonly leaving the phase dif-ference∆φ. Combining equations (2.8) and (2.10) provides a complete formula for the phase difference in eq. (2.11), expressed in known quantities.

∆φ=ωc¨D

c ¨sin(θ) +2πN (2.11)

Some notational simplifications can be made at this stage. Since wc=2π fcand c = fc¨ λ one can write:

ωc¨D c =2π D fcλ =2π D λ ”2πDλ (2.12)

Here the quantity D_λis introduced to help simplify notation. This is the distance between the antennas expressed in units of wavelengths. Using this definition and eq. 2.11 one gets an expression eq. (2.13) for θ expressed in known quantities.

θ=sin´1 _{∆φ ´ 2πN} 2πD_λ  =sin´1  1 D_λ¨ _∆φ 2π ´N  , N PZ (2.13)

The inverse sine function is valid (i.e. real-valued) for input values in the domain[´1, 1] while the quantity∆φ lies in the range[´π, π]. In the worst case this means that solutions become ambiguous once | ˘ π

2π´N| ă Dλ ô | ˘0.5 ´ N| ă Dλsince eq. (2.13) becomes solvable for N other than zero.

2.1. Monopulse Techniques -4 -3 -2 - 2 3 4 Signal phase -2/₄ -1/₄ 1_/ 4 2_/ 4 Incident angle N = -1 N = 0 N = 1 Ambiguity point Example Solutions

Figure 2.7: Inverse sine ambiguity with Dλ=1.25 and∆φ=´34π.

Method Combination

The amplitude and phase methods have complimentary properties that combine well. The ambiguity inherent in the phase method can be resolved using the unambiguous amplitude method. With the ambiguity resolved one can use the more precise and unbiased direction estimates from the phase method to cover the entire incident signal angle domain. Using the quantities θAand∆φ one needs to resolve the ambiguity N in eq. (2.13).

It takes some work to robustly merge the measurements while handling the discontinu-ities in the phase information since they are noisy. When∆φ « ˘π there are two values of N that can solve the equation.

The easiest way to combine the angle θAand phase difference∆φ is to take a ’step back-wards’ by calculating sin(θA)which leads to eq. (2.14).

θA=θPô θA=sin´1  1 Dλ ¨ _∆φ 2π ´N  ô sin(θA) = 1 D_λ ¨ _∆φ 2π ´N  ô sin(θA) D_λ = ∆φ 2π ´N N= ∆φ 2π ´ sin(θA) D_λ (2.14)

This also provides an easy heuristic for estimating the trustworthiness of a given mea-surement. Since N has to be integer valued the fractional part of ˆN should be small. If this is

2.2. Phasor Calculations

not true then we risk an error in our ambiguity resolution, with a large error in the direction estimate as a result.

With N determined one can then perform the phase direction calculations from eq. (2.13) without ambiguity.

2.2

Phasor Calculations

The received signal I/Q representation is not immediately useful when it comes to calculating the incident signal angle. Rather; the derived signal characteristics of amplitude and phase are needed to estimate the direction of the transmitter. These therefore need to to estimated from the I/Q sample data. This report compares three different techniques for performing this estimation.

Instantaneous

The instantaneous method can provide amplitude and phase information from I/Q data on a per-sample basis and has a low inherent latency because of it. The first step is forming the complex signal defined in eq. (2.15).

yn(t) =in(t) +j ¨ qn(t) (2.15) By performing this step the signal amplitude and phase quantites can be derived by trans-forming the quantity ynto its polar representation as shown in (2.16).

An(t) = b in(t)2_+qn(t)2 φn(t) =atan2(qn(t), in(t)) zn(t) =An(t)¨exp(j ¨ φn(t)) (2.16)

The major drawback of this method is that it isn’t frequency selective nor does it suppress any noise.

φ

A

I Q

Figure 2.8: Complex signal representation of I/Q data.

Rather than calculate the phase and amplitude separately it is possible to use the CORDIC algorithm to calculate both at the same time which simplifies implementation. The CORDIC algorithm is a method for calculating vector rotations through an iterative process. In princi-ple this is done by finding the rotation angle α that aligns a vector with the positive real axis. One can then determine the phase as φ = ´αand amplitude A = Re(zR)from the rotated sample value zr.[13]

2.2. Phasor Calculations

Fourier Transform

The discrete time fourier transform is a very common tool used for signal analysis. Rather than immediately calculating the phasor characteristics one can calculate Yn[f] =Ftyn[t]u. If the carrier frequency fc is known one calculates the amplitude and phase from Yn[fc]as in eq. (2.16). If the carrier frequency isn’t known the best strategy is probably to iterate through Yn[f]and select the AC component with the highest amplitude.

The primary benefit of using the fourier transform method is that both noise and interfer-ing signals on other frequencies are heavily suppressed. White gaussian noise is by definition evenly spread over the frequency spectrum. This means that if an N-sample DTFT is used the noise energy is reduced by a factor N.

The main drawbacks of using the fourier transform is that it is computationally expensive and delays the phasor estimates by at least the time it takes to collect all the samples. The Vivado VHDL IP implementations of the FFT generally also separate input, computation, and output steps introducing additional large delays.[10] The frame lengths N are also fixed to powers of two.

Görtzel Algorithm

The Görtzel algorithm is a less well-known version of the Fourier Transform. Rather than calculate the complete DTFT it allows the user to only calculate a specific component Yn[f]. This in turn implies that the carrier frequency fcmust be known. In return the computational complexity is reduced. If only log(N)components are needed the Görtzel algorithm is more computationally efficient. Implementing the Görtzel algorithm also reduces the phasor cal-culation delay since the phasor can be output as soon as the last sample has been collected.[3] The Enhanced Görtzel algorithm can also circumvent the binning errors caused by the fourier transform bin frequencies not necessarily exactly matching the carrier frequency. This has the added benefit of allowing the operator to pick window lengths that are not powers of two.[11] 0 200 400 600 800 1000 Number of samples 1.81 1.812 1.814 1.816 1.818 1.82 1.822 1.824

Phase estimate [Rad]

FFT Görtzel

2.3. Noise Analysis

The primary drawback of the Görtzel method is that the equivalent filter implementing the algorithm is metastable. In practice this means that numerical errors tend to accumulate over multiple iterations.[11] As such the number of fractional bits used to store the interme-diate calculation will have dramatic effects on the output quality, as shown in fig. 2.10.

0 5 10 15 20 25 30 35

Number of fractional bits

-80 -70 -60 -50 -40 -30 -20 -10 0 Amplitude error [dB]

Figure 2.10: Görtzel Algorithm relative error for N=1024

These calculations use double precision floating point calculations as a baseline truth when comparing the results of fixed point calculations.

2.3

Noise Analysis

The noise analysis is separated into three main sections corresponding to the major compo-nents of the device. The analog front-end noise, the transformed noise when switching to polar representation, and finally the noise in the direction estimate.

Sampling Noise

The sampling noise is likely composed of two major components. The first is thermal noise, also known as Johnson-Nyquist noise. Any electrical system at temperatures above 0K are affected by noise due to the random movement of electrons. This is essentially gaussian white noise at the frequencies we’re interested in. The thermal noise power PT is proportional to both temperature and bandwidth as described by eq. (2.17).[5]

PT=4 ¨ kb¨T ¨∆ f W (2.17)

Working at room temperature with a bandwidth of around 10 MHz gives an inherent thermal noise level of -104 dBm. This noise exists in the pre-amplifier stage and is therefore amplified along with the signal.

2.3. Noise Analysis

Noise is also introduced through quantization errors in the ADC once the signal has been amplified. Even without sampling error a 12-bit number can not perfectly represent the ana-log values. There’s an average error depending on the sample width and input voltage range. The maximum input voltage amplitude is 0.625 Volts.[12] This yields the voltage resolution ∆Vdefined in eq. (2.18).

∆V= 2 ¨ 0.625

212 =0.3052 mV (2.18)

The quantization can then be viewed as additive white gaussian noise with a RMS voltage VQudescribed by (2.19).[2]

VQu= ∆2

V

12 (2.19)

With the value of∆V from (2.18) and a resistance of 50Ω this gives a quantization power level equivalent to approximately -68 dBm. In the worst case scenario with the maximum LNA amplification (62 dB) applied the quantization and thermal noises should result in a terminated signal noise power of approximately -42 dBm. The AD9361 Data-sheet claims a noise figure of approximately 4 dB at maximum RX gain so an experimentally measured noise power should be around -38 dBm at worst.[12]

The number of bits it takes to represent this noise is calculated in eq. (2.20). PW =0.001 ¨ 10 PdBm 10 W V=?P ¨ R Nsteps= V ∆V Nbits=log2 Nsteps

Nbits=log₂ ?P_∆VW¨R



(2.20) Inserting P=´38dBm and R=50Ω gives Nbits«3.2bits so one should expect that in the worst case the input noise covers up at least the lowest three bits.

Since the noise components are normally distributed and the sum of normally distributed stochastic variables is itself normally distributed one can reasonably model the input noise as normally distributed. The noise distribution of the sampled signal is thus described by eq. (2.21). X „ N(µ,Σ) pX(x) = 1 2π|Σ|¨exp  ´1 2¨(x ´ µ)Σ(x ´ µ) T (2.21)

Because the AD9361 circuit has DC rejection functionality the mean values µ are assumed to be zero. The input noise can therefore be described by three parameters σi, σq and ρ as shown in eq. (2.22). The parameter rho is the noise correlation coefficient between the chan-nels while σiand σqare the noise variances for the respective channels.

e(t) = ei(t) eq(t)  =N  0,  σ_i2 ρσiσq ρσiσq σq2  (2.22) Visualizing this distribution is an important part of understanding the noise analysis. While the parameters in eq. (2.22) map directly to the measurable Cartesian signal noise distribution the matrixΣ is unwieldy when it comes to analysis.

Rather than describe the variances along the coordinate axes one can use the parameters of the resulting ellipse to describe the distribution. The new parameters λ1, λ2, and φe are shown in fig. 2.11.

2.3. Noise Analysis λ1 λ2 λI λQ φe

Figure 2.11: Mathematical representations of elliptical noise.

Any non-degenerate bivariate normal distribution X can be transformed into a uncorre-lated distribution Y. This is done by rotating the coordinate system to align with the eigenvec-tors e. For a normally distributed variable X „ N(µ,Σ)one can therefore diagonalizeΣ. This

yields a uncorrelated normal distribution with variance equal to the eigenvalues λ1and λ2in the new coordinate space. The angle of the major axis eigenvector is equal to tan´1(e1,2/e1,1). A new quantity referred to as the elliptical variance is introduced in eq. (2.23). This is an expression for the radius of the elliptical distribution as a function of the angle relative to the major axis and it turns out to be useful in later computations.

σ_e2(φ) =λe(φ) = cos(φ) 2 σ₁2 +sin(φ) 2 σ₂2 !´1 (2.23)

Phasor transformation

The phasor calculators take sample values as input and produce phasors as output. We’re interested in how they in the ideal case would transform a particular input distribution into output distributions.

How does a particular I/Q signal with possibly correlated AWGN e(t)transform into to phasors? In other words: what is the polar representation of a bivariate normal distribution

2.3. Noise Analysis

with arbitrary mean µ and covariance matrixΣ? The polar transformation of an IQ signal is described by eq. (2.24).

i=ρcos(φ), µi=µρcos(µφ) q=ρsin(φ), µq=µρsin(µφ)

di dq=ρdρ dφ

(2.24)

Of primary interest is describing to what extent the shape and qualities of the output distribution depends on the input SNR. The problem is that the full expression for this distri-bution (shown in eq. (2.25) is essentially intractable.

2π|Σ|pZ(z) =exp ´1 2¨ "  ρcos(φ)´ µρcos(µφ) σ1 2 + ρsin(φ)´ µρsin(µφ) σ2 2#! = exp ´1 2¨ "  ρ σe(φ) 2 +  µρ σe(µφ) 2 ´2ρµρ cos(φ)cos(µφ) σ₁2 +sin(φ)sin(µφ) σ₂2 !#! (2.25) By setting the variances to unity the polar gaussian PDF in eq. (2.25) is greatly simplified into eq. (2.26). pZ(z) = 1 2πexp  ´1 2¨ h ρ2+µ2_ρ´2ρµρcos(φ ´ µφ) i (2.26) Even for the simplified case of eq. (2.26) the expression is too complicated to be of immediate theoretical use. The expression in (2.25) simplifies when µρÑ0 or µρÑ 8corresponding to the very low and very high SNR cases. These are covered in the following chapters.

Pure Noise

The integral of any probably density function over its entire domain is by definition equal to one. This combined with eq. (2.24) provides eq. (2.27).

1= 8 ż ´8 8 ż ´8 pY(y)dy1dy2= 8 ż 0 π ż ´π pZ(z)ρdρ dφ (2.27)

Setting µρin eq. (2.25) to zero one can express eq. (2.27) as eq. (2.28).

1= 1 2πσ1σ2 8 ż 0 π ż ´π exp  ´1 2¨ ρ2 σe(φ)2  ρdρ dφ (2.28)

The variable K defined in eq. (2.29) is introduced for the sake of readability. K(φ) =exp  ´1 2¨ 1 σe(φ)2  0 ă K ă 1 (2.29)

Introducing K makes it obvious that one can evaluate the integral for ρ and φ separately as is done in eq. (2.30). 1= 1 2πσ1σ2 8 ż 0 π ż ´π K(φ)ρ 2 ¨ ρ dr dφ= 1 2πσ1σ2 π ż ´π " 1 2¨ K(φ)ρ2 ln(K(φ)) #8 0 dφ (2.30)

2.3. Noise Analysis

Because K(φ)is less than one the upper limit of this integral converges to zero. This leaves

eq. (2.31). 1 4πσ1σ2 π ż ´π 1 ln(K(φ)) dφ= ´1 4πσ1σ2 π ż ´π 1 cos2_(φ) σ₁2 + sin2(φ) σ₂2 dφ= ´1 4πσ1σ2 π ż ´π σ_e2(φ)dφ (2.31)

From the definition of a CDF we finally get an expression for the angular PDF. Combin-ing this expression with the original coordinate system rotation φe provides a closed form expression eq. (2.32) for pΦin the original coordinate system.

p_Φ(φ) =´ 1 4πλ1λ2

¨ σ_e2(φ ´ ˆφ) (2.32)

The only thing that is left to determine is pρ(ρ, φ). For a given value of θ we know that the bivariate distribution along to this direction is a one-dimensional normal distribution Nφwith standard deviation λφ(φ). The issue is that the distance ρ is the absolute value of this distribution. However, the absolute value of a zero-mean normal distribution is the Chi distribution! One can therefore express pρas eq. (2.33).

pρ(ρ) = #

λφ(φ)¨ χ(ρ) ρ ě0

2.3. Noise Analysis

-4

-2

2

4

I-component

-2

-1

1

2

3

Q-component

2.3. Noise Analysis

-4

-2

2

4

I-component

-3

-2

-1

1

2

3

Q-component

Figure 2.13: Monte Carlo simulation of correlated pure WGN in polar coordinates.

The figures include a 95% confidence interval marked in red. Signal SNR sensitivity

The signal distribution undergoes a complicated transformation as the SNR rises from ´8. The generalized expression for the noise distribution in eq. (2.25) is essentially intractable and so some heuristics need to introduced in order to estimate the phase and amplitudes distributions.

Since the received signal phase from . (2.2) is from a travelling sinusoidal signal the mea-sured phases are continuously "wandering" around the phasor circle. However; the noise distribution ellipse is pointed in a specific direction depending on the noise covariance ma-trixΣ. This means that in the presence of correlated or uneven component noise the phasor noise distribution eZbecomes elliptical. This means that the phasor distribution Z becomes

2.3. Noise Analysis

-2

2

4

6

I-component

-3

-2

-1

1

2

3

Q-component

Figure 2.14: Example low SNR WGN signal distribution.

Four times per period the polar noise ellipse major and minor axes are at a right angle to the phasor vector. At these times the noise is uncorrelated (but not necessarily evenly distributed) and as a result the amplitude distribution follows a Generalized Non-central Chi distribution. At these points it is therefore possible to analytically express the expected dis-tribution.

By comparing Monte Carlo simulations of correlated noise to the Generalized Chi distri-bution one can therefore determine the bias error induced by the signal components being correlated as shown in fig. 2.15.

2.3. Noise Analysis

Figure 2.15: Correlation induced amplitude bias.

Rather than show the phase varying amplitude error one can plot the maximum error as a function of SNR and signal correlation as is done in fig. 2.16.

2.3. Noise Analysis

Figure 2.16: Correlation induced maximum amplitude bias.

The phase error induced by correlated signal components is harder to describe analyt-ically. Monte Carlo simulations and numerical integrations of the phase distributions are shown in fig. 2.17 and 2.18.

2.3. Noise Analysis

2.3. Noise Analysis

Figure 2.18: Maximum correlation induced phase bias.

In the high SNR limit the noise correlation essentially disappears and the distributions are normal. A comparison between Cartesian and polar normal distributions for varying SNR is shown in fig. 2.19. Since the noise can be assumed to be uncorrelated the distribution is rotated to be aligned with the x-axis (φ=0) without loss of generality.

When calculating the parameters for the polar distribution one needs to estimate the polar variances σρand σφfrom σiand σq.

When making the gaussian polar approximation we assume that σidepends solely on σρ and σq depends solely on σφ. This assumption makes some intuitive sense since the radial vector and tangential vector point along i and q, respectively. If one makes this assumption then σρis simply equal to σi.

2.3. Noise Analysis

Calculating σφis harder. We introduce the identities in eq. (2.24) which can be derived from the expected mean of the log-normal distribution.

E[exp(iφ)] =exp ´

σ_φ2 2 ! ñ E[cos(φ)] =exp ´ σ_φ2 2 ! E[sin(φ)] =0 (2.34)

Using (2.34) one can calculate how the variances of i and q depend on σφwhich yields eq. (2.35).

Var[i] =E[i2]´E[i]2=ρ2¨E[cos(φ)2]´exp(´σ_φ2)ô Var[i] =ρ2¨1 2E[1+cos(2φ)]´ ρ 2_¨exp(´σ2 φ) Var[q] =E[q2]´E[q]2=ρ2¨E[sin(φ)2]´0 ô Var[q] =ρ2¨1 2E[1 ´ exp(2φ)] σ_i2= ρ 2 2 ¨  1 ´ exp(´σ_φ2)2 σ_q2= ρ 2 2 ¨  1 ´ exp(´2σ_φ2) (2.35)

Equation (2.35) shows that σiis far less dependent on σφthan σqis because of the higher exponent. This fact lends credence to the validity of the simplifying assumption that σφonly depends on σq.

By solving eq. (2.35) for σφone gets eq. (2.36)

σ_φ2=´1 2 ¨log 1 ´ 2 ¨ σ_q2 ρ2 ! (2.36) Using these approximations and plotting the Cartesian Gaussian distribution ellipse and the polar Gaussian distribution ellipse transformed into the Cartesian coordinate space one gets fig. 2.19.

2.3. Noise Analysis 2 4 6 -3 -2 -1 0 1 2

3 _{Polar Gaussian Approximation}

Cartesian Gaussian Distribution

1 2 3 4 5 6 7 -3 -2 -1 0 1 2

3 _{Polar Gaussian Approximation}

Cartesian Gaussian Distribution

8 10 12 -3 -2 -1 0 1 2

3 _{Polar Gaussian Approximation}

Cartesian Gaussian Distribution

Figure 2.19: Comparison of polar and cartesian WGN. SNR = {3,6,10} dB

The overlap between the polar and Cartesian normal distributions is shown in fig. 2.20. This measure is calculated by taking the ratio of the intersection and union of the sets formed by the 99% confidence intervals.

2.3. Noise Analysis 0 2 4 6 8 10 12 14 16 18 20 SNR [dB] 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Distribution overlap

Figure 2.20: Overlap between polar and Cartesian normal distributions.

These figures show that the polar Gaussian approximation works very well for higher SNR values but breaks down in the low SNR range. This is not unexpected because the fun-damental distribution actually changes in the lower SNR ranges. In the very low SNR range the amplitude is chi-distributed and the phase is drawn from the wrapped normal distribu-tion. As such a Gaussian approximation in the polar domain will necessarily fail at lower SNR.

2.3. Noise Analysis

Direction Estimation

With knowledge of the phasor distributions one can analyze the distribution of the amplitude ratio R from eq. (2.6) and phase difference∆φ from (2.9). These are immediately connected to the incident signal angle and therefore provide a measure of the quality of our direction estimate.

This analysis assumes that the signal amplitudes and phase are uncorrelated and follow a normal distribution, as shown in eq. (2.37). The quality of this assumption can be seen in fig. 2.20. In practice this means that the results in this chapter can only really be considered valid for SNR greater than 10 dB.

An „N(µAn, λAn)

φn „N(µφn, λφn)

(2.37) For a given pulse we sample a set of values A1, A2, φ1, φ2. The parameters Dλand k are assumed to be known without error.

The first step is to calculate the statistical properties of the amplitude direction estimate

θAand phase difference∆φ.

θA=k´1¨ A1 ´A2 A1+A2 ∆φ=φ1´ φ2

Working from the defining equations (2.7) and (2.9) and using the normal distribution addition rules one gets the distribution parameters in eq. (2.39) and (2.38) for θAand∆φ.

∆φ „ N µ∆φ, λ∆φ
µ_∆φ=µφ1´ µφ2 λ_∆φ=λφ1+λφ2 (2.38) θA„N µθA, λθA
µ_θA « µA1´ µA2 µA1+µA2 λ_θA «k´2¨ λA1+λA2 µA₁+µA2
2 (2.39)

With∆φand θA determined the next step is to calculate the distribution for N using eq. (2.14). Making the simplifying assumption that sin(θA)« θA one gets the distribution in eq. (2.40). N „ N(µN, λN) µN « µ_∆φ 2π ´ µθA Dλ λN« λ_∆φ 4π2 + λθA D_λ2 (2.40)

With the distribution of N determined one can finally analyze the expression for θPin eq. (2.13). θP=sin´1  _∆φ 2πDλ ´ N D_λ 

The distribution of θPends up fairly complicated if one attempts to solve it exactly. By linearizing the inverse sine function around the point _DN

λ one gets a much more manageable

expression, shown in eq. (2.41). This linearization better represents the underlying distribu-tion as Dλincreases, since N can attain more discrete values and we linearize around N.

2.3. Noise Analysis

Note that N is assumed to be constant in the following analysis steps. This is a necessary simplifying assumption to avoid calculating the bivariate distribution of θ and N. Since N is rounded towards the closest integer this is a fairly reasonable assumption as long as λNis not too large. f(x+c)« f(c) + d dxf(c)¨(x ´ c) +O(x 2₎ d dφsin ´1 N Dλ  = 1 Dλ¨ b 1 ´_DN λ N ă Dλ θ=sin´1  _∆φ 2πD_λ´ N D_λ  « ´sin´1 N D_λ  +  _∆φ 2πD_λ + N D_λ  ¨ 1 Dλ¨ b 1 ´_DN λ (2.41)

Which finally yields an expression for the distribution of the direction estimate θ shown in eq. (2.42). θ „N(µ_θ, λ_θ) µ_θ =sin´1  _µ ∆φ 2πDλ ´µN Dλ  λ_θ«   1 2πD2_λ¨ b 1 ´ _DN λ   2 ¨ λ_∆φ (2.42) Disambiguation Failure

When attempting to resolve the value of N a rounding strategy is used where the estimate ˆN is rounded to the nearest integer. This transforms the continuous distribution fC(n)into the discrete equivalent fD[˜n], shown in eq. (2.43).

P(˜n= N) = żN+0.5

N´0.5

fC(n)dn (2.43)

2.3. Noise Analysis -2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5 N 0.1 0.2 0.3 0.4 0.5 0.6 Probability _Continuous Discrete Figure 2.21: PDF of discretized ˆN

What is most important is not the specific distribution, but rather the probablility of failure Pfailwhich is equivalent to P(N ‰ Nˆ ) =1 ´ P(N). This is shown in fig. 2.22 over a range of variances and bias errors.

2.3. Noise Analysis

3

Hardware

3.1

MikTran

The SDR platform MikTran is a proprietary platform built from commercially available com-ponents and is shown in fig. 3.1. The platform was designed at FOI and consists of a Zynq ZC7020 SoC and an AD9361 IC.

The Zynq ZC7020 SoC is a chip with an FPGA and a dual-core ARM Cortex-9 CPU run-ning Arch Linux. The AD9361 IC is an RF transceiver capable of demodulating IQ data with a carrier frequency between 80 MHz and 6 GHz with a bandwidth of up to 56 MHz.[12]

3.2. Platform

The AD9361 chip operates independently from the BBP. Settings are written through SPI commands implemented in an operating-system independent library known as the AD9361 No-OS library. This library is written in C and runs on the PS. Samples are read from the AD9361 into the PL through an LVDS interface.

3.2

Platform

In order to perform angular measurements the system was mounted on a servo-powered turntable. The turntable waws powered by a HS-645MG 180 degree servo controlled by an Arduino Uno. The position of the turntable could be controlled through MATLAB. The turntable was in turn mounted on a Manfrotto 161MK2 camera tripod with built in spirit levels.

In order to allow for different antennas and antenna spacings a plank with holes drilled at 5 cm intervals was mounted onto the turntable along with the MikTran module. The entire setup can be seen in figures 3.2 and 3.3.

3.3. Test Equipment

Figure 3.3: System with mounted spiral antennas.

The system requires three separate cables to control. Ethernet and power cables for Mik-Tran along with a USB cable for the Arduino. The servo control performed by the Arduino could in the future be implemented on the MikTran card in order to simplify the design.

3.3

Test Equipment

In order to evaluate the system and measure the antenna gain functions a separate transmitter was necessary. This was simply a horn antenna mounted on a home-made aluminium stand as shown in figure 3.4. The output signal was generated by an Agilent N5181A pulsed sine-wave generator.

4

Method

Servo lag compensation

While servos are generally fairly accurate the servo used had some issues properly catching up to the requested angles, often undershooting by a degree or two. In order to remove this measurement error the servo was swept from its lower limit to the upper limit and then back again when performing measurements.

Working under the assumption that the servo lag was equally severe in both directions it is then possible to calculate the correlation delay between the positive and negative direction sweeps. By shifting the measurements towards the center by half of this delay it is possible to eliminate the impact of the servo on the measurements.

This method was used in all measurements that required sweeping the antennas. Using this method greatly increased the horn antenna measurement quality due to their great size and inertia. The spiral antenna setup was not as affected by servo lag due to their small size and placement closer to the center of rotation resulting in a smaller moment of inertia.

4.1

ADC Calibration

The AD9361 circuit is not by default calibrated to match absolute power units. By connecting a high quality signal generator to the system using cables with known attenuation it is possi-ble to match the sample values received from the system with real units. The signal generator used was an Agilent N5181A and the cables used were Pasternack 36" precision cables with a loss of Gcable =´0.9 dB at 6 GHz.[9] The AD9361 has a preamplifier before the ADC with some gain GLNA.

A pure sinusoid with known power PindBm is applied to the system. The actual power at the system input is then Pin+GcabledBm because of the cable losses. The power at the ADC is then defined by eq. (4.1).

PADC=Pin+GLNA+Gcable dBm PADC=exp₁₀ Pin+GLNA+Gcable

10 ´3



4.2. Performance measurements

The RF equipment has a standard impedance of 50Ω. This means that a sine with power PADCW has a RMS voltage of VRMS=

?

50 ¨ PADCV and a peak voltage of VMAX= ?

2 ¨ VRMS. The ADC saturation point in the AD9361 IC is at 0.625V.[12] The ADC sample width is 12 bits signed which means that for a given peak voltage sample value S the corresponding measured voltage is VS=S ¨0.625₂11 V.

Since the measured sample voltages VSand the actual input voltage VMAXshould be the same one can calculate a correction factor γ for the voltages such that eq. (4.2) holds.

VS=γ ¨ VMAX (4.2)

By expanding until the known terms Pinand S are reached and solving that for γ one gets the expression in eq. (4.3).

VS=γ ¨ VMAX S ¨0.625 211 =γ ¨ ? 2 ¨ VRMS S ¨0.625 211 =γ ¨ ? 2 ¨a50 ¨ PADC S ¨0.625 211 =γ ¨10 ¨ a PADC S ¨0.0625 211 =γ ¨ d

exp₁₀ Pin+GLNA+Gcable

10 ´3  S ¨0.0625 211 =γ ¨exp10  Pin+GLNA+Gcable 20 ´ 3 2  γ=S ¨0.0625 211 ¨exp10 

´Pin+GLNA+Gcable

20 +

3 2



(4.3) From eq. (4.3) one can then determine γ if the signal generator output power and the measured sample values are known!

4.2

Performance measurements

System performance was measured at three distinct junctions in order to evaluate the theo-ries.

• After sampling. This measures thermal and ADC quantization noise. • After phasor calculation. Determine method robustness.

• After direction estimation. Determine information merge efficiency.

Of primary interest is how input signal SNR propagates through the system and affects the direction estimates. It is also necessary to assess the validity of the assumptions made in this thesis. Is the noise sufficiently close to WGN? To what degree are the two channels actually uncorrelated?

By terminating the input ports and measuring signal input in the absence of an applied signal a reasonable measurement of the signal-independent noise present in the signal chain can be obtained. By varying system parameters such as LO frequency and sample rate one can perform a regression analysis in order to build a model of approximate noise levels given a set of parameters.

There might also be noise at specific frequencies due to things such as clock leakage or nearby high-powered transmitters so a full sweep of the available LO frequencies was per-formed to check for inherent sources of non-WGN in the system. This is done by combining

4.3. Antenna directivity measurement

the fourier transform spectra for measurements made using varying LO frequencies. Since complex signals are measured it is possible to get single-sided spectra which makes it simple to stitch together multiple DTFT spectra.

4.3

Antenna directivity measurement

In order to use amplitude monopulse techniques one needs to know the antenna gain char-acteristics. While datasheets are available it is generally a good idea to measure the antenna gain functions before use. Ideally one would measure these in both an anechoic chamber as well as in the field but gaining access to such a testing room was not feasible.

Instead the antenna gain functions were field-tested. This was done by placing a horn antenna (shown in fig. 3.4) at a distance of approx. 20 meters pointed towards the system.

By rotating the system while measuring the amplitude of a signal transmitted from this antenna the amplitude function can be estimated. By applying a pure sine-wave and using the fourier phasor transform method one can reduce the amount of external noise influencing this experiment. Maintaining a high SNR is also important.

With the antenna gain functions measured it is a simple matter of performing a curve fit to the region between the main lobes to get an amplitude ratio model. Two models were tested: linear and third degree polynomial. While more sophisticated models are likely possible to implement the polynomial approximation offers a good tradeoff between accuracy and implementation complexity. Curve fitting was performed in MATLAB and the calculated polynomial coefficients then transferred to the platform.

4.4

Calibration Routine

In the ideal case both of the receiver channels have identical characteristics. In reality this is not quite the case even with closely matched antennas and cables. The channels may have both different gains due to antenna differences as well as different phase delays because of signal path length differences.

In the case of different signal path lenghts eq. 2.10 is founded on the incorrect assumption that both channels are affected by the same absolute phase delay φx. Rather; there will be some constant difference between the phase delay because of the added length to the signal path.

With a transmitter placed along the boresight axis the phase difference between the chan-nels should be zero. By measuring the actual phase difference one gets a measure of the phase bias φB. This bias is recorded and later subtracted from any future measurements in order to effectively ’zero’ the system.

The final primary physical parameter left to estimate is the distance between the antennas, Dλ. While this can in principle be measured physically there is quite a bit of margin for error. This error is mainly because of the fact that the antennas are not point receivers and so it is difficult to determine what points to measure the distance between.

Rather than measuring the distance it is possible to calculate it. By differentiating eq. (2.11) and inserting eq. (2.12) one gets eq. (4.4).

d(∆φ)

dθ =2πDλ¨cos(θ) (4.4)

This means that through measuring the slope of ∆φ in the region around zero (where cos(θ)«1) one gets an estimate of the distance between the antennas by solving eq. (4.4) for Dλ.

5

Implementation

The goal of this project was to provide a usable VHDL module that is synthesizable on the MikTran FPGA platform and can be used to estimate the direction of an incident signal. In principle this final product does not require any surrounding infrastructure and can be run on the FPGA directly. However, in practice this project could not have been executed with-out the surrounding software infrastructure in order to analyze data, configure settings, and visualize results. As such the final product ends up having to be a self-contained fully usable system rather than a simple software module.

The overall high level system architecture is shown in figure 5.1.

Figure 5.1: High level system architecture.

LVDS interface

The LVDS standard transmits data using low-voltage differential signals. This reduces both noise and power draw compared to regular logic gates.[6]

The clock of the LVDS interface is controlled by the AD9361 and is proportional to the sample rate. It is a DDR interface which means that it transmits data on both the rising and falling edges of the bus clock.

With both channels active the channel needs to transfer 12 ¨ 4 = 48 bits of data for each set of samples. Since the AD9361 has twelve output channels and uses differential signalling it can transmit six bits per edge. This means that each sample requires a total of four clock cycles to transfer to the FPGA. The format of this transfer is shown in fig. 5.2.

Figure 5.2: LVDS packet structure

The samples are then reconstructed in the FPGA.

AXI Interface

The PL and PS subsystems on the Zynq ZC702 SoC communicate through an AXI Interface. A codebase for this interface was provided by FOI and only slightly modified to fit the project. The AXI interface enables DMA and interrupt propagation between the FPGA and proces-sors. It is used to implement a number of asynchronous registers as well as two synchronous 64-bit wide data channels that can be used to transmit synchronized sequences of data.

The asynchronous registers are primarily used for writing parameters and control regis-ters while the data channel are used to transmit RX and TX channel data. The TX channel was not used for this project. The RX channel was used to receive a number of different types of data from the direction estimation module such as raw samples, phasors, debug information and direction estimates.

Network Server

The signal processing was implemented in the PL. However, in order to control the system settings and stream data to a workstation there needed to be a network server running on the MikTran board. Since the library to control the AD9361 IC was written for the PS the server could be written in C and did not have to be implemented on the FPGA.

The network server implements a simple but effective protocol on top of TCP/IP. The TCP protocol has strong guarantees regarding packet loss and ordering which simplifies the design immensely. The server is entirely passive and only responds to commands from the client. The packet structure is shown in fig. 5.3.

IDENTIFIER COMMAND PACKET LENGTH

COMMAND SPECIFIC DATA FIELD

8 16 32

5.1. FPGA

The identifier field is a randomly chosen constant which exists to validate packets and make sure the server and client send and receive data in sync. The command field is a simple command mapping where the functions that are to be available over the network are enumer-ated by unique IDs. The command data field is used to provide function arguments to these functions. If more bytes of data are needed the Packet Length field defines how many bytes of raw data will follow this packet header.

The network interface exposes all of the AD9361 control functions defined in the No-OS library to the client. It also allows for user defined functions such as parameter writes and sample reading to be performed over the network. By using the AXI interface the server can also write and read both registers and streamed data from the PL. Once the server had been implemented only minor adjustments and protocol extensions had to be made to it over the course of the project. Implementing the server and a MATLAB client were quite significant up-front time investments but the choice sped up the design cycle immensely and was ultimately well worth it.

While the resulting interface works well it is not robust. It relies on good-faith actors on both sides of the network channel to remain stable and would need to be extended or replaced if employed in an open network.

5.1

FPGA

The actual direction estimation algorithm and signal preprocessing was implemented on the FPGA. The design was subdivided into a number of modules as shown in fig. 5.4.

phasor

(Görtzel)

phasor

(Inst.)

axi_ fo

Settings

Samples

Phasors

Directions

data_packer

rf_rx

axi_reg

rx_in

rx_dir_est

phasor

(FFT)

Figure 5.4: FPGA Module Structure

Samples and a clock signal are received from the AD9361. All of the signal processing is performed in the ADC clock domain. A double-clocked FIFO queue is used to transfer data from the ADC clock domain to the FPGA core clock domain which is running at 100 MHz.

5.2. MATLAB

The first step in the signal processing chain is to calculate the polar representation of the I/Q data. This was implemented in three ways:

• Instantaneous translation using a CORDIC IP block. • Fourier transform of signal before translation. • Görtzel’s algorithm before translation.

The instantaneous translation method was implemented using a CORDIC block.

The Fourier Transform approach was implemented using an FFT IP block with a maxi-mum frame length of N=1024 samples.

The Görtzel algorithm was implemented using complex multiplication IP blocks for the multiplications while the rest of the algorithm was implemented manually. The Görtzel al-gorithm implementation performed well in simulations but did not work on the FPGA. Due to time constraints it was decided that debugging this phasor calculation method was a low-priority.

In order to avoid spurious direction estimates when no signal is transmitted during pulses a thresholding function was implemented at this stage. The simplest form is a direct com-parison to some user defined constant. The final implementation also includes debouncing functionality where the signal amplitude is required to have been stable for some amount of samples before the signal is considered present.

If a signal is considered to be present the direction estimation module receives the phasors and performs the calculations described in the theory chapter. The sin and sin´1functions were implemented using 10 bit wide LUTs implemented using a ’arbitrary function module’ with function values generated from MATLAB.

The direction estimates are finally passed to the data packer module. The data packer module allows user selection of the received data streams.

5.2

MATLAB

In addition to the network client a number of different visualization and analysis tools were written for MATLAB. All testing and measurement was automated and controlled from the MATLAB interface. In order to do this a pulsed sinewave generator was also controlled from MATLAB.

6

Results

6.1

ADC Calibration

The measured values are fairly close to correct. The worst case scenario is off by approxi-mately 20%. 0.7 -60 80 0.8 0.9 -70 60

Correction factor

Input Power [dBm]

Gain [dB]

-80 1.1 40 1.2 -90 -100 20

6.2. Noise Measurements 0 10 20 30 40 50 60 70

Gain [dB]

Correction factor

Figure 6.2: Voltage calibration factor γ at maximum input power.

6.2

Noise Measurements

The system inherent noise spectrum is shown in figure 6.3.

0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 Frequency [GHz] 109 20 21 22 23 24 25 Amplitude [bits]

6.3. Antenna Directivity Measurements

This measurement is what led to the final carrier wave frequency of 5.5GHz being used for this system.

Noise power was measured at maximum gain in a room temperature environment with terminated RF ports and an input bandwidth of 10 MHz. An LO frequency of 5.5 GHz was used to avoid contamination from the inherent noise sources shown in fig. 6.3. With these parameters a noise power of ´39 ˘ 1 dBm was measured for all channels which matches well with the theoretical values derived in chapter 2.

A correlation analysis between the channels and components was performed as well, un-der the same circumstances as the noise power test. These results are shown in table 6.1.

I

Q

I

Q

I

1.000

0.005

-0.009

0.016

Q

-

1.000

0.007

0.011

I

-

-

1.000

-0.004

Q

-

-

-

1.000

Table 6.1: Input channel noise correlation.

These show that the intra-channel noise correlation is not very severe, definitely not enough to induce a substantial bias (as shown in figures 2.16 and 2.18). The inter-channel noise correlation is worse, with I0 and Q1having the strongest noise coupling. The inter-channel correlation is probably worse because of the AD9361 IQ-demodulation error com-pensation modules working to separate the intra-channel components.

6.3

Antenna Directivity Measurements

Two types of matched antenna pairs were tested for this report. One set of narrow-lobed horn antennas and one set of wide-lobed spiral antennas shown in figures 3.2 and 3.3 respectively. The results are shown in figures 6.4 and 6.5.[1] [7]

The primary practical source of bad measurements appears to be signal reflections from surrounding objects as well as the ground. Indeed - the first set of spiral antenna measure-ments collected were completely unusable because of ground reflections. This was mitigated by raising the system higher above the ground while also leaning the transmitting antenna backwards until the ground was outside of the main lobe.

6.3. Antenna Directivity Measurements -100 -80 -60 -40 -20 0 20 40 60 80 100 Angle [Deg] -18 -16 -14 -12 -10 -8 -6 -4 -2

Relative Gain [dB] Left Antenna Right Antenna

-60 -40 -20 20 40 60

Incident signal angle [Deg]

-0.8 -0.6 -0.4 -0.2 0.2 0.4 0.6 0.8

Amplitude Ratio _Measured

Linear Fit Cubic Fit

Figure 6.4: Measured gain and amplitude ratio functions for horn antenna pair.

The wide lobes of the spiral antennas are likely to blame for the bad measurement results in 6.5. The test environment included quite a few nearby reflective surfaces and due to prac-tical reasons it was impossible to measure elsewhere. As such the spiral antenna amplitude ratio looks a lot worse than the horn antenna one.

6.4. Calibration -100 -80 -60 -40 -20 0 20 40 60 80 100 Angle [Deg] -14 -12 -10 -8 -6 -4 -2 Gain [dB] Left Antenna Right Antenna -60 -40 -20 20 40 60

Incident signal angle [Deg]

-0.8 -0.6 -0.4 -0.2 0.2 0.4 0.6 0.8

Amplitude Ratio _Measured

Linear Fit Cubic Fit

Figure 6.5: Measured gain and amplitude ratio functions for spiral antenna pair.

6.4

Calibration

6.4. Calibration

-80 -60 -40 -20 20 40 60

Angle [Deg]

-Phase [Rad]

Figure 6.6: Horn antenna phase difference measurements and separation distance estimate.

-80 -60 -40 -20 20 40 60

Angle [Deg]

-Phase [Rad]

Figure 6.7: Spiral antenna phase difference measurements and separation distance estimate.

The estimated antenna distances for the horn and spiral antennas are 0.29 and 0.11 m respectively. The horn antennas were mounted in holes 0.2 m apart and the spiral antennas were mounted 0.1 m apart. The large discrepancy for the horn antennas is likely because of a failing point source approximation. When squinted the openings of the horn antennas were approximately 0.3 m apart which matches with the estimates. The spiral antenna estimate was much closer to the physically measured distances, likely because they are much smaller and rotating them around their center of mass does not displace their apertures.

6.5. Direction Estimates

6.5

Direction Estimates

The final direction estimates are shown in two forms. The first is a plot of the estimate density and the other is a plot of the mean estimate.

Horn antenna results

-30 -20 -10 10 20 Actual direction -30 -20 -10 10 20 30 Estimated direction

6.5. Direction Estimates -30 -20 -10 10 20 Actual direction -30 -20 -10 10 20 Estimated direction

6.5. Direction Estimates -30 -20 -10 10 20 Actual direction -30 -20 -10 10 20 Estimated direction

6.5. Direction Estimates -30 -20 -10 10 20 Actual direction -30 -20 -10 10 20 30 Estimated direction

6.5. Direction Estimates -30 -20 -10 10 20 Actual direction -30 -20 -10 10 20 Estimated direction

6.5. Direction Estimates

Spiral antenna results

-30 -20 -10 10 20 Actual direction -40 -30 -20 -10 10 20 Estimated direction

6.5. Direction Estimates -30 -20 -10 10 20 Actual direction -40 -30 -20 -10 10 20 30 Estimated direction

6.5. Direction Estimates -30 -20 -10 10 20 Actual direction -30 -20 -10 10 20 Estimated direction

6.5. Direction Estimates -30 -20 -10 10 20 Actual direction -30 -20 -10 10 20 30 Estimated direction

6.5. Direction Estimates -30 -20 -10 10 20 Actual direction -40 -30 -20 -10 10 20 Estimated direction