# Estimation and mitigation of phase noise in a coherent distributed passive radar system

TABLE OF CONTENTS Page List of Figures iv List of Tables ix Chapter 1: Introduction 1 Chapter 2: Manastash Ridge Radar 4 2.1 Hardware 5 2.2 FM Broadcasts 7 2.3 Software 9 2.4 MRR Targets 13 Chapter 3: Time Series Analysis 14 3.1 ARMA Models 14 3.2 Wiener Process 19 3.3 ARIMA Models 21 3.4 ARFIMA Models 21 Chapter 4: Oscillator Phase Noise 23 4.1 PSD Specifications 24 4.2 Previous Work 27 Chapter 5: Phase Noise Analysis 30 5.1 Propagation Through MRR Processing 31 5.2 Effect on Range-Doppler 35 5.3 Presence of Additive Noise 35 5.4 Estimation of the Time Series 39 5.5 Detection of an Interferer 41 l

5.6 Structure of the Phase Noise Time Series 45 5.7 Performance of Time Series Analysis 45 5.8 Consistency of Models 55 5.9 Wiener Model Analysis 59 5.10 Summary 63 Chapter 6: Correction Techniques 64 6.1 The Sharpness Metric 64 6.2 Autofocus Algorithm 65 6.3 Voltage Level Correction 66 6.4 Autofocus Performance 67 6.5 Mitigation by Deconvolution 68 6.6 Simulation Study 70 6.7 Summary 74 Chapter 7: Performance Analysis 75 7.1 Summary Statistics 75 7.2 Mapping the Performance 80 7.3 Linear Regression 80 7.4 Local Linear Regression 82 7.5 Gaussian Process Regression 82 7.6 K-Means Classification 83 7.7 Cross-Validation 83 7.8 Regression Results 83 7.9 Summary 85 Chapter 8: Conclusion 90 8.1 Summary 91 8.2 Future Work 92 Bibliography 94 Appendix A: Identity Proofs 101 A.l Re[E[g(X)}]=E[Re[g(X)}} 101 n

A.2 E[X[n}*Y{n]] = E[X[n})*E[Y[n]} 102 Appendix B: Choice of Sample Rate 103 Appendix C: Stationarity Tests 105 C.l Priestly and Subba Rao 105 C.2 Kwiatkowski, Phillips, Schmidt, and Shin 105 Appendix D: Particle Swarm Optimization 107 Appendix E: Cramer-Rao Bound 109 i n

LIST OF FIGURES Figure Number Page 2.1 The Manastash Ridge Radar System 5 2.2 A block diagram of the MRR receivers. Both the PLL and the digitizer card receive control signals from the processor 6 2.3 A block diagram of the MRR processing algorithm 10 2.4 An example range-Doppler output from MRR. The color shows a nor malized power scale in dB. The large target that spans from 800 - 1100 km is ionospheric turbulence from auroral echoes. The strong target around 50 km, 0 m/s is the reflection from Mt. Rainier. The weak target around 70 km, 100 m/s is most likely an airplane. Meteors are also occasionally spotted on MRR at widely varying ranges above planes and below aurora. The vertical streaking that occurs around targets are ambiguity sidelobes associated with the FM waveform be ing transmitted 12 3.1 Example ACF and PACF functions for MA and AR processes. The dashed lines show significance levels where 95% of values that are truly 0 will appear between them 18 3.2 An example path of a Wiener process with a2 = 0.8 and a timestep of 0.01 20 4.1 A block diagram of a PLL 24 4.2 A PSD plot of an ideal sinusoid affected by a sample phase noise process. A 10 kHz sample rate signal, x(t), was created as x(t) = exp(j27r/ot+0(t)), where /o =2 kHz and

4.4 Sample power distribution curves of phase noise processes. When the noise power level of the DUT is being measured using a reference with significantly less noise power over all frequencies (such as Refl), we can assume that the reference is not contributing to the noise power spectrum. If a reference such as Ref2 is used, where the reference noise power is on par with or exceeds that of the DUT, we must recognize that the measured phase noise at that frequency is including significant contribution from the reference. If the noise characteristics of such a reference are well known, we can still infer the noise characteristics of the DUT at these frequencies 28 4.5 A block diagram of a phase discriminator 28 5.1 Three plots of a sample Doppler spectrum for a single range containing a large clutter return and two relatively smaller target returns. The progression from left to right shows increasing amounts of phase noise existent in the system. Notice how the ability to separate the two targets becomes lost. With more phase noise added to the system, even a single target would become impossible to detect in the presence of the clutter: either widening clutter peaks may envelop targets of interest or the rising noise floor may grow high enough to make targets undetectable. By the third plot, a potential target has also appeared around 95 Hz even though it does not appear in the first plot 36 5.2 A plot of power at the zero Doppler frequency vs. range lags. Maximum power occurs at lag 41, corresponding to 410 //s and 61.5 km 40 5.3 A sample realization of

5.9 A sample realization of 4>e[n] computed over the same 10 seconds worth of data as in Figure 5.3 but using a decimation factor of 1000 47 5.10 A sample realization of (j)e[n] computed over the same 10 seconds worth of data as in Figure 5.3 but using a decimation factor of 4000 47 5.11 A sample realization of e[n] computed over the same 10 seconds worth of data as in Figure 5.3 but using a decimation factor of 8000 48 5.12 A sample realization of 4>e[n] computed over the same 10 seconds worth of data as in Figure 5.3 but using a decimation factor of 16000. ... 48 5.13 A sample realization of

5.25 A comparison of PSDs estimates for eJ<^tn' using a Wiener model for (j)t[n\. The 3 points plotted show the expected values obtained from Table 4.1 assuming the 2 PLLs are the sole contributors to system phase noise. The solid line shows the PSD computed from Equation 5.66 using the average a2 from Table 5.1 (a2 = 0.0091 at 6.25 Hz). The dashed line shows the PSD computed using the estimate of a2 that minimizes the mean-squared error (MSE) between the PSD and the 3 points (a2 = 0.344 at a 6.25 Hz rate) 62 6.1 A comparison of a true Doppler with the result of the deconvolution algorithm. The data used is an actual reference signal processed with a simulated scatter signal without any decimation 71 6.2 A comparison of the deconvolution algorithm using simulated data with the true phase noise free spectrum and the pre-processed spectrum. We have zoomed in on the target of interest because a view of the entire 2 kHz Doppler spectrum would not show enough detail to distinguish the curves 72 6.3 Plots of the absolute error between the curves in Figure 6.2 73 6.4 The same curves from Figure 6.3 but focused in to the range of Doppler shifts that contain the target of interest 73 7.1 An example target Doppler profile illustrating the measurement of SNR. The noise floor was estimated using the mean value of a target- free range. The figure is shown in dB for ease of reading although the actual computations of signal power, noise power, and their ratio are performed using linear values. Once the SNR ratio is computed, it is then converted into dB 76 7.2 An example target Doppler profile illustrating the measurement of SNR. The noise floor was estimated using the mean value of a target- free range. The figure is shown in dB for ease of reading although the actual computations of peak power, sidelobe power, and their ratio are performed using linear values. Once the PSR ratio is computed, it is then converted into dB 77 vn

7.3 A comparison of target Doppler profiles computed from a data set from Table 7.1 at 97.3 MHz. The target, which contains Doppler shifts ranging from approximately 100-300 Hz, sees an increase of about 2 to 3 dB/Hz per Doppler bin. The power levels shown are arbitrary - only the relative powers are relevant. We note that there has been some amplification of high frequency content, evidenced by the mitigated curve being more jaggy 79 7.4 Plot of the residuals from fitted SNR values using the linear regression parameters from Table 7.6. The 2 groups correspond to the different illuminating signal frequencies 86 7.5 Plot of the residuals from fitted PSR values using the linear regression parameters from Table 7.6. The residuals from PSR do not show the same grouping effect we see in Figure 7.4 from the different illuminating signals 87 7.6 A terrain map of the local Seattle area, showing the location of the 96.5 MHz FM and 97.3 MHz FM transmitters. Map made using Google Maps (http://maps.google.com) 89 B.l Aliasing the FM band by 56 MHz and 72 MHz sample rates 104 via

LIST OF TABLES Number Page A phase noise comparison (in dBc/Hz) between the Nova Engineering NovaSource NS2-0045-025 and the EM Research SPC-PLXO-56-06 ac cording to their specification sheets, "n/a" is used if no information is provided for that frequency offset. The Nova PLL has a tunable output frequency whereas the EM PLL is fixed 25 Estimates of a2 for the consecutive sets of 10 second data in Figure 5.13. 56 Estimates of a2 for the first 17 consecutive sets of 10 second data in Figure 5.16 57 A summary of the effect of the mitigation procedure over 120 targets. Mean AMetric refers to the mean change in the given metric using 120 targets. Std. Dev. AMetric refers to the standard deviation of the change in the given metric over those same 120 targets 78 ANOVA table for the data using Equation 7.3. As all three factors are found to be significant at the 5% level, we choose to include all three effects in predicting SNR 81 ANOVA table for the data using Equation 7.4. As all three factors are found to be significant at the 5% level, we choose to include all three effects in predicting PSR 81 Cross-validation results showing the optimal values of K for the regres sion/classification techniques. In the case of Gaussian process regres sion, a vector is shown representing K as discussed in Section 7.5. . 84 A comparison of the performance of the function mapping techniques. Each value reported is the MSE for the performance of the given regres sion/classification scheme, after cross-validation, attempting to predict the given metric 84 Estimated linear regression parameters 85 IX

ACKNOWLEDGMENTS I would like to thank my advisor, John Sahr, for the opportunity and support to work on this research, and I would like to thank my committee members for their help in revising and improving this work. I would also like to thank my labmates Melissa Meyer, Zac Berkowitz, and Laura Vertatschitsch for their help not only in developing this work but also in blurring the. line between work and play. x

DEDICATION / would like to dedicate this work to all my friends and family. XI

1 Chapter 1 INTRODUCTION "Well, I believe in the soul... the small of a woman's back, the hanging curve ball, high fiber, good scotch, that the novels of Susan Sontag are self-indulgent, overrated crap. I believe Lee Harvey Oswald acted alone. I believe there ought to be a constitutional amendment outlawing Astroturf and the designated hitter. I believe in the sweet spot, soft-core pornog raphy, opening your presents Christmas morning rather than Christmas Eve, and I believe in long, slow, deep, soft, wet kisses that last three days." — Crash Davis, Bull Durham Range and velocity estimation is a critical function for many radar remote sensing systems. Accurate Doppler spectral estimates place stringent requirements on the quality of the collected data. The local oscillators providing the reference for the radar receivers are largely responsible for the resolution of the target time and frequency in each range Doppler cell. The phase noise of the oscillators effectively distorts the range Doppler signature of targets. Additional phase noise in the system from sources such as atmospheric effects also contributes to deviation in Doppler bin estimates from a phase noise free estimate. Using real data from a coherent distributed passive radar system, we will extract estimates of the phase noise and determine what can be learned from its structure. By performing time series analysis, we will examine potential stochastic models for these phase noise processes. Comparing various stochastic models, we will see which models appropriately describe the nature of the phase noise process. These models include parameters that can be computed regularly to monitor system performance.

2 We will also offer potential mitigation techniques to apply to the radar system's output in an attempt to improve the range-Doppler output dat a products. The clas sic autofocus algorithm, an optimization algorithm adapted from other phase noise correction scenarios, is discussed in our context. We also put forward an alternate approach not based on solving an optimization, which we found to produce superior results. In Chapter 2, we describe the radar system that generates the data we will be using for analysis. Both hardware and software level descriptions are detailed. The hardware perspective relates one source of phase noise in the system: that which is contributed by the phase-locked loop. The software processing is integral to our analysis as it is the basis for our estimation procedure. Chapter 3 is a primer on the stochastic processes we will refer to work in this work. ARMA models and Wiener processes are commonly used and should be familiar to most readers who deal with stochastic processes on some level. ARIMA and ARFIMA models, which are also presented, are less commonly studied but are very useful for our discussion of phase noise. In Chapter 4, we discuss phase noise as it relates to oscillators. Oscillator phase noise is very relevant to our discussion since we know t hat our radar system's phase noise contains at the very least a contribution from the phase-locked loops. Most pre vious work in phase noise analysis and simulation has been performed in the context of oscillators. Chapter 5 tracks phase noise as it propagates through the radar system. We derive the mathematics that shows the effect phase noise has on the radar's output and illustrates our technique for estimating the phase noise time series. In the process, we discover an interesting experimental result, where we are able to identify a likely interfering transmitter by examining the phase noise. We are also able to examine the spectral structure of the phase noise at very low frequencies, corresponding to frequencies very close in to the carrier in the context of oscillator phase noise. This

3 is a view of phase noise power close to carrier that is not usually seen. Using estimates of the phase noise, we then attempt to fit the stochastic models discussed earlier to the data, and we conclude that the Wiener model is most appro priate. The resulting models provide us with parameter estimates for each model, and we consider whether the model parameters have any interpretation as far as system performance. In Chapter 6, we first discuss the autofocus algorithm, which is a phase noise mitigation technique derived in other applications. Then we propose an alternate mitigation technique based on deconvolution, which produces better results and is also faster to execute. We use simulation results to demonstrate the technique. Chapter 7 examines the performance of our mitigation technique. We find that our technique achieves an average improvement of 2.5 dB in signal-to-noise ratio (SNR) and an average improvement of 1 dB' in peak sidelobe ratio (PSR) over our collected observations. This boost in SNR translates to approximately doubling the sensitivity of our radar, purely through processing techniques. In addition to collecting summary statistics describing performance, we attempt to fit a mapping to the performance data. The mapping is used in order to predict the improvement attainable from our mitigation algorithm given the quality of the data going into it, and we find that linear regression provides a useful linear model in predicting performance. In Chapter 8, we summarize this work, noting the main results of our analysis. We also suggest topics of research that, if pursued, could continue and/or complement this work.

4 Chapter 2 MANASTASH RIDGE RADAR "There are some places that the road doesn't go in a circle. There are some places where the road keeps going." — David, Pleasantville The Manastash Ridge Radar (MRR) is a coherent distributed passive radar system based at the University of Washington. Leveraging the power of commercial FM broadcasts, it detects scatter from a variety of targets in the atmosphere. The radar was designed to detect meter scale plasma density irregularities in the ionosphere [17,57]. In addition it is sensitive to meteor trails (out to about 850 km) and aircraft (about 200 km). We will analyze the presence and effect of phase noise using the data and processor of MRR. Interest in passive radar systems has grown considerably in recent years as a variety of technological advances have reduced the expense of computation, making passive radars more accessible. As they rely upon existing high power transmitters, passive systems are inexpensive, safe, and avoid introduction of new interference problems. Existing TV and FM broadcasts in the VHF and UHF bands provide substantial radio illumination at wavelengths that are useful for geoscience as well as aerospace applications. Lockheed-Martin1 and Roke Manor2 both have developed projects using passive radar technology [51]. 1 http: //www. lockheedmartin. com 2http://www. roke. co.uk

Targat Y FM Transmit! Scattered Energy :tw X % hS tntenwt Connection University o( Washington Receiver Manastash RMge Observatory Receiver Figure 2.1: The Manastash Ridge Radar System 2.1 Hardware Several reports have described MRR and its observations [34,37,58,59]. Its primary data outputs, range-resolved Doppler power spectra and interferometry [36], are on par with other VHF coherent scatter radars used for ionospheric research. Figure 2.1 shows the fundamental composition of the system. A receiver at the Univer sity of Washington records FM broadcasts as a reference signal. A second receiver located at the Manastash Ridge Observatory (MRO), near Ellensburg WA, collects the signals scattered from targets in the atmosphere such as airplanes, meteors, and ionospheric turbulence associated with the aurora borealis. MRR has the added geo graphic advantage that the Cascade Mountain Range, which contains Mt. Rainier as the largest peak, suppresses most of the direct path transmitter power that reaches MRO. Mt. Rainier, however, does reflect a large amount of power to the MRO re ceiver and routinely appears as a large, fixed target. Its presence in the range-Doppler outputs is used to indicate whether MRR is functioning properly. The reference and

6 Y GPSAn»mi GP5 Receiver \ •10 MHz- Wwt t ocked loop ^i V Bf Antenna Bandpass Flltori -100 MHz- KMHst Dlgtesr Cart ^ 1! •100 MHz. Figure 2.2: A block diagram of the MRR receivers. Both the PLL and the digitizer card receive control signals from the processor. scatter signals are time identified with GPS receivers for use in processing. The scat ter signal is sent across an Internet connection to the University of Washington to be processed with the reference signal. Figure 2.2 shows a block diagram of the MRR receivers. The RF antenna col lects and passes signals through bandpass filters, preserving the FM frequency range (centered around 100 MHz), then sends it to the digital receiver in a PC. The GPS receiver provides a high quality 10 MHz sinusoid to the phase-locked loop (PLL). The PLL uses the 10 MHz signal to generate the 56 MHz3 signal required by the digitizer card. The PLL's 56 MHz signal is not as "clean" as the 10 MHz, carrying a larger amount of phase noise in it (see Chapter 4). The digital receiver selects the specific FM channel of interest (e.g. 97.3 MHz) by bringing that frequency to baseband and filtering out the other channels. The digital receiver also downsamples the signal, 3See Appendix B for a discussion on the choice of sample rate for the digitizer card.

7 providing 100 kHz I/Q data samples for use in processing. 2.2 FM Broadcasts Passive radar relies upon the fact that the correlation properties of commercial FM broadcasts resemble those of bandlimited white noise, which is desirable as true white noise has a two dimensional delta function for its selfambiguity4. Empirical studies [23] show that broadcast music (of all styles) works better than speech. The pauses are periods of low modulation and thus long correlation time, which destroys the range resolution. The form for a monaural FM broadcast is taken to be x(t) = Re AeJucteja Sl a{t')dt' (2.1) where LOC is the carrier frequency and a(t) is the content signal5. The lack of a lower limit in the integral in Equation 2.1 indicates that the initial start time is arbitrary. It only serves to set the initial phase of the phase accumulator. As a simple technique for demodulating an FM signal involves evaluating the difference of the phase signal, it does not matter what that initial phase value is. We will model a(t) as a random process with mean zero, unit variance, and band width characteristic of audio signals. The parameter, a, relates the modulation in dex6. The factor Ae>Wct is nonrandom high frequency content associated with the transmitter carrier frequency, so we will only consider the random component in ex- 4The ambiguity function is discussed in Section 2.3 Estimation of the parameters for signals of this form is an interesting problem and strategies are highly contingent on the scenario involved [39,43]. 6For more information on FM modulation see [74].

8 amining the statistics of x(t): E Re ejafla(t')dt' = Re E > pa{t')dtr = Re[E[ei ay]] , Y -iV(0,t) = Re[$y(a)] = Re a. = e '• r a 2 t l e 2 2 t (2.2) (2.3) (2.4) (2.5) (2.6) In Equation 2.2 we used the relationship proved in Appendix A.l. In Equation 2.3 are making the assumption that samples of a(t) are independent and identically dis tributed and appeal to the Central Limit Theorem [63] to replace the integral with Y, and we then use the definition of the characteristic function of a Gaussian random variable in Equation 2.4. As the modulation process continues, the phase of the con tent signal continues to add to prior values, and in practice that accumulation has been processing for a long time. This means that we are concerned with t —> oo in which case E [x(t)] —> 0, and we can calculate the autocovariance function (which we define in Chapter 3) of x(t) as ~A*e~j"ote~jafta{t')dt'Aej"o{t+T)ejaft+Ta(t")dt" Ja{- f* a(t')dt'+ft+T cc(t")dt") eja ftt+T a(t')dt'~ = |A| V-W0TE [ejaY] , Y = N(0, r) = \A\2e^0T^Y{a) E[x*(t)x{t + r)] - E = |A|2e^0TE = \A\2ejU0TE = \Afe 2OJU0T, (2.7) (2-8) (2.9) (2.10) (2.11) (2.12) This shows that the FM signal model decorrelates as r increases. The typical band width of FM broadcasts is of the order 100 kHz, suggesting that the correlation time is about 10 microseconds. This has been shown by Meyer [36], who notes that the 100 kHz sampling rate that MRR employs provides nearly uncorrelated samples of

9 x[n], although the sample bandwidth is about 67 kHz. 2.3 Software MRR has several operational modes [38], but we will focus on the most basic. The signal processing algorithm computes the cross-ambiguity of the reference and scatter signals. Equation 2.13 shows the continuous, idealized form of the cross-ambiguity function, where x is the reference signal, y is the scatter signal, r is the range dimen sion, and v is the Doppler dimension. \A(r, u)\2 represents the scattering power of the targets at range r and Doppler shift v. Although Equation 2.13 is not what the MRR processor directly implements, it motivates the algorithms that we actually apply. This equation can be interpreted as the magnitude squared of the Fourier transform of a matched-filtering operation between x and a shifted y: A(r,u) = (2.13) x*(t)y(t + r)e~j27n/tdt . Figure 2.3 shows a block diagram of the actual computation that MRR performs in estimating the cross-ambiguity function. For each range of interest, the scatter signal is shifted by a corresponding number of samples. In routine operation, the sample rates of x and y are 100 kHz and thus each sample corresponds to 10 [is of range delay, or about 1.5 km of monostatic physical range. The processor then performs a point-by-point multiplication between the conjugated x signal and the time-shifted y signal to create what is called the detected signal, a: a[n] = x*[n]y[n + r]. (2.14) Next the processor performs a coherent integration7 on this new series by summing together a number of these points. MRR usually uses a coherent integration factor of d = 50, which reduces the sample rate from 100 kHz to 2 kHz: 49 ad[n] = J^a[50n + fc]. (2.15) fe=0 In this work we may also refer to this coherent integration process as decimation.