RL-948 Time-Frequency Analysis in Radar Backscatter Problems Christopher Joseph McCormack (Doctoral Thesis) December, 1996 RL-948 = RL-948

Time-Frequency Analysis in Radar Backscatter Problems by Christopher Joseph McCormack A dissertation submitted in partial fulfillment of the requirements for the degree of Doctor of Philosophy (Electrical Engineering) in The University of Michigan 1996 Doctoral Committee: Professor William J. Williams, Co-Chair Research Scientist Valdis V. Liepa, Co-Chair Professor Alfred 0. Hero III Associate Professor Kamal Sarabandi Professor Yue Ying Lau Professor Leon Cohen, Hunter College, CUNY

ABSTRACT Time-Frequency Analysis in Radar Backscatter Problems by Christopher Joseph IMcCormack Co-Chairs: William J. Williams, Valdis V. Liepa Time-frequency techniques provide new and unique insights for analyzing electromagnetic scattering problems. These techniques transform functions of time or frequency into two dimensional functions of both time and frequency to reveal non-stationary characteristics of the signal. The theory developed herein justifies applying the frequency-time transform to wide bandwidth signals illuminating stationary targets. The frequency-time representation of the return provides more information about the target and the scattering than regular Fourier analysis. Along with the position of the scattering centers, frequencytime analysis gives insights on the target's composition and configuration. In addition, the performance of these transforms when applied to noise are examined and quantified. The statistics of a Cohen's class time-frequency transformation are derived and verified numerically. Applying the time-frequency techniques to sampled continuous-wave radar data from a dynamic target provides insight into target motion and generate estimates of the target parameters. After considering a figure of merit for evaluating the time-frequency distribution generated, a customized kernel, determined using a genetic algorithm, is used to improve the performance of the standard spectrogram, the Wigner distribution, and the binomial distribution.

~Christopher Joseph McCormack 1996 All Rights Reserved

To Lisa ii

ACKNOWLEDGEMENTS I would like to thank those who made my program at Michigan possible, starting with my advisors, Valdis Liepa and William Williams. They were willing to take the time and effort to deal with a graduate student who didn't fit the regular mold and a research project outside of the scope of their normal work. My thanks go to my fellow graduate students, both in the Radiation Laboratory and in Professor William's time frequency group. They've answered many questions over the last three years, listened to and provided ideas, and made the whole experience more human and interesting. Special thanks go to Joseph Burns of the Environmental Research Institute of Michigan and Randy Haupt of the United States Air Force Academy for providing computer support, chamber support, and a ready technical sounding board, especially in the early stages of my research. Finally, there was no way I could have succeeded at Michigan without the full support of my family. They have made countless sacrifices to allow me to follow this dream, many of which I may never realize. iii

TABLE OF CONTENTS DEDICATION............................................ ii ACKNOWLEDGEMENTS.................................... LIST OF TABLES................................... LIST OF FIGURES.................................. LIST OF APPENDICES............................... iii viii ix xi CHAPTERS 1 INTRODUCTION............. 1.1 Analysis Domains........ 1.1.1 Time Domain...... 1.1.2 Frequency Domain. 1.1.3 Time-Frequency Domain 1.2 Scattering Problems....... 1.2.1 Dispersive Phenomena. 1.2.2 Dynamic Targets... 1.2.3 Noise Corrupted Signals 1.2.4 Customized Kernels... 1.3 Summary............ 2 BACKGROUND.............. 2.1 Traditional Spectral Analysis 2.1.1 Standard Usage..... 2.1.2 Power Spectral Density. 2.1.3 Limitations....... 2.2 Time-Frequency Analysis........................................ 1 1 1 2 2 2 3 3 3 3 4 5 5 5 6 6 7 7 9 9 10 10 10 11........................................................................ 2.2.1 Ambiguity Domain Signal Transformations.. 2.2.2 Autocorrelation Domain Signal Transformations 2.2.3 Sampled Signal Analysis............. 2.3 Desirable Time-Frequency Distribution Properties... 2.3.1 Non-Negative................... 2.3.2 Real Valued.................... 2.3.3 Time Shift Invariant............... 2.3.4 Frequency Shift Invariant.................. 11 iv

2.3.5 Time Marginal............ 2.3.6 Frequency Marginal........... 2.3.7 Instantaneous Frequency........ 2.3.8 Group Delay............... 2.3.9 Time Support.............. 2.3.10 Frequency Support........... 2.3.11 Reduced Interference.......... 2.4 Kernel Design Procedures........... 2.5 Alias-Free Formulation............. 2.6 Dual of transformations............ 2.7 Emphasized Kernels.............. 2.7.1 Spectrogram............... 2.7.2 Wigner Distribution.......... 2.7.3 Binomial Distribution......... 2.7.4 Genetic Kernel............. 2.8 Review of Time-Frequency Developments... 2.9 Time-Frequency Applications......... 2.9.1 Marine Mammal............ 2.9.2 Biomedical................ 2.9.3 Mechanical Engineering........ 2.9.4 Electromagnetics............ 2.10 Summary.................... 3 STATIONARY TARGETS.............. 3.1 Introduction................... 3.2 Background................... 3.3 Duality and Time-Frequency Transformations 3.4 Scattering and Dispersion........... 3.4.1 Reflective Returns........... 3.4.2 Creeping Waves............. 3.4.3 Combined Modes............ 3.5 Transforms Considered............. 3.5.1 Spectrogram............... 3.5.2 Wigner............... 3.5.3 Binomial................. 3.6 Target Geometries............... 3.6.1 Single Sphere.............. 3.6.2 Two Spheres: In-line.......... 3.6.3 Broadside................ 3.7 Simulation Parameters............. 3.8 Results...................... 3.8.1 Single Sphere.............. 3.8.2 Two Spheres: In-Line Orientation... 11 12 12 13 13 13 14 14 14 15 15 15 16 16 17 17 18 18 18 18 18 19:20.20. 20 21.22 23 23 23.23 923.24 94...... 25 2 4 95 26 2)7 28 28.... 32 3.8.3 Two Spheres: Broadside Orientation Perpendicular Polar ization.............................. 3.8.4 Two Spheres: Broadside Orientation - Aligned Polarization 3.9 Summary................................ 36 40 44 v

4 STATISTICS FOR COHEN'S CLASS TFD OF WHITE GAUSSIAN NOISE 15 4.1 Introduction............. 4.2 Continuous Transform....... 4.2.1 Continuous Mean...... 4.2.2 Continuous Variance.... 4.3 Discretized Transform....... 4.3.1 Discretized Mean...... 4.3.2 Discretized Variance.. 4.4 Sampled Transform......... 4.4.1 Data Limitations...... 4.4.2 Sampled Mean....... 4.4.3 Sampled Variance...... 4.5 Specific Distributions........ 4.5.1 Discrete Wigner kernel.. 4.5.2 Binomial kernel....... 4.6 Numerical Results.......... 4.7 Summary.............. 5 DYNAMIC TARGETS........... 5.1 Introduction............. 5.2 Background............. 5.3 Wire Targets............. 5.3.1 Single rotating wire..... 5.3.2 Multiple rotating wires... 5.4 Spherical Targets.......... 5.4.1 Electromagnetics problem 5.4.2 Processing and Estimation 5.4.3 Test results......... 5.4.4 On-axis shifts........ 5.4.5 Off-axis shift......... 5.4.6 Higher order effects..... 5.4.7 Customized kernel..... 5.5 Summary............................. 4 5............... 4 6.............. 7................ 4 8...... 53............... 5 3............... 5 4.......56.......57............... 5 8............... 5 9............... 6 4............... 64............... 65............... 66............... 68 39 39 70 70 71 75 77 79 80 85 86 87 87 88 88 89 89 90 90 93 94 ()4 95 95 97 98 98 99 6 EFFECTS OF ADDING NOISE AND COMPLEXITY SIGNALS......................... 6.1 Background.................... 6.2 Signal M odels................... 6.2.1 Point Targets............... 6.2.2 Distributed Target............ 6.3 Noise Model.................... 6.3.1 Noise Sources............... 6.3.2 Additive White Gaussian Noise..... 6.4 Renyi Entropy................... 6.4.1 Distribution Evaluations......... 6.5 Analysis and Results............... 6.5.1 Receiver Operating Characteristics.. 6.5.2 Entropy Measures............. TO SCATTERED.. *. *. *................................. *................ *... *. *....................................................... vi

6.6 Summary........................... 7 DEVELOPMENT AND USE C 7.1 Background...... 7.1.1 Genetic Algorith 7.1.2 Problem Design System Optimization RE Summary....... )F A CUSTOMIZED KERNEL.................. ms................ ~sults............................... ~ o 7.2 7.3 8 CO NCLUSIONS................................. 8.1 Stationary Targets........................... 8.2 Dynamic Targets............................ 8.3 Noise Statistics............................. 8.4 Effects of Noise and Complexity.................... 8.5 Customized Kernels........................... 8.6 Value of Time-Frequency Techniques................. 8.6.1 Image Processing........................ 8.6.2 Ultrawideband Radar...................... A PPEN D IC ES...................................... BIBLIOGRAPHY.................................... 100 101 101 101 102 103 104 105 105 105 106 106 106 106 106 107 108 122 vii

LIST OF TABLES Table 2.1 Selected TF Kernels.............................. 8 2.2 Desirable Distribution Properties......................... 11 4.1 W igner Statistics................................... 67 4.2 Binomial Statistics.............................................. 68 5.1 Case A Results................................... 85 5.2 Case B Results................................... 86 5.3 Case C Results..................................... 86 5.4 Case D Results................................... 86 5.5 Case E Results.................................. 87 viii

LIST OF FIGURES Figure 3.1 3.2 3.3 3.4 3.5 3.6 3.7 3.8 3.9 3.10 3.11 3.12 3.13 3.14 3.15 3.16 3.17 3.18 3.19 3.20 3.21 3.22 3.23 4.1 5.1 5.2 5.3 5.4 5.5 5.6 5.7 5.8 5.9 Single Sphere................................... Two Spheres: In-Line.............................. Two Spheres: Broadside............................. Single Sphere: Spectral Response........................ Single Sphere: Inverse Fourier Transform................... Single Sphere: Spectrogram........................... Single Sphere: Wigner Distribution....................... Single Sphere: Binomial Distribution...................... In-Line Spheres: Spectral Response....................... In-Line Spheres: Inverse Fourier Transform.................. In-Line Spheres: Spectrogram.......................... In-Line Spheres: Wigner Distribution......................... In-Line Spheres: Binomial Distribution..................... Broadside - Perpendicular Polarization: Spectral Response.......... Broadside - Perpendicular Polarization: Inverse Fourier Transform...... Broadside - Perpendicular Polarization: Spectrogram............. Broadside - Perpendicular Polarization: Wigner................ Broadside - Perpendicular Polarization: Binomial............... Broadside - Aligned Polarization: Spectral Response............. Broadside - Aligned Polarization: Inverse Fourier Transform......... Broadside - Aligned Polarization: Spectrogram................ Broadside - Aligned Polarization: Wigner................... Broadside - Aligned Polarization: Binomial.................. Noise Variance vs Transform Length...................... Dynamic Measurement System Model..................... Single Wire Geometries............................. Single Wire Backscatter Patterns............................ Single Wire Waveforms............................ Single Wire Power Spectral Density....................... Single Wire Binomial Distributions............................ Multiple Wire Geometries............................ Multiple Wire Backscatter Patterns...................... Multiple Wire Waveforms............................ 25 25 26 29 30 30 31 32 33 34 35 35 36 37 38 39 39 40 '41 42 42?1 43 44 66 71 72 72 73 74 74 76 76 77 ix

5.10 NIultiple Wire Spectral Densities. 78 5.11 Multiple Wire Binomial Distributions.................... 78 5.12 Target Dimensions................................ 79 5.13 Centers of Rotation................................ 80 5.14 Orbiting Spheres Waveform......................82 5.15 FFT of Case A Waveform............................ 83 5.16 FFT of Case E Waveform.............................. 83 5.17 128-Point Spectrogram of Case E Waveform.................... 84 5.18 Case A Time-Frequency Representation..................... 84 5.19 Case E Time-Frequency Representation.................. 5.. 6.1 Point Target Geometry................................ 91 6.2 Transformations for Single Point Scatterer..................... 92 6.3 Two-Point Scattering Geometry............................ 92 6.4 Two-Sphere Scattering Geometry........................ '93 6.5 Wire Scattering Geometry............................ 94 6.6 Receiver Operating Characteristics....................... 99 6.7 Entropy Margins......................................... 99 7.1 Genetically Derived Elementary Function....................... 103 7.2 Genetically Derived Time-Frequency Distribution................ 104 x

LIST OF APPENDICES APPENDIX A Notational Conventions............................. 109 B Computer Subroutines............................... 110 xi

CHAPTER 1 INTRODUCTION 1.1 Analysis Domains Many different representations may be used to analyze signals. The transformation of a signal from one form to another may allow simplification of the math (for example, replacing convolutions with multiplications), or it may provide better insight into the behavior of a signal (or system) by relating the measured signal to expected signal characteristics. Some commonly used one-dimensional transforms include the Fourier transform, thie Laplace transform, and the Z-transform. Another type of transform moves from one (iimension to multiple dimensions. These multidimensional transforms also include Cohen's class of time-frequency distributions [12, 32] and the wavelet transform [6, 32]. 1.1.1 Time Domain The time domain representation for a signal typically represents the most; natural form of a signal. Conditions of realizability force any real-world signal to have finite starting and ending times. This domain also relates directly to the direct measurements that can be obtained. The time domain representation provides the simplest view of aperiodic or non-stationary signals. When accurately measured, the time domain of the signal contains all the information available in the signal. Other domains might make certain relationships more clearly visible, but no additional information is generated. 1

1.1.2 Frequency Domain Fourier analysis is normally introduced by starting with a function of time, g(t), and finding an equivalent function of frequency, G(f). This uses the Fourier integral to substitute a dependence on the frequency variable f instead of the time variable t. The frequency domain representation emphasizes the periodic aspects of a signal. For linear system analysis, this provides a valuable analysis tool, both for insight into the nature of a system or signal and for evaluation of a system's response. It is important to remember that the frequency domain representation is a mathematical concept. Theoretically, a finite time signal will have a frequency domain representation that is not band limited. Infinite frequency limits are not a problem when the Fourier domain is treated mathematically. 1.1.3 Time-Frequency Domain While both the time and frequency domain representations of a signal are complete, certain relationships are more apparent in the time domain (signal transitions, starts, and stops), while other effects are easier to identify in the frequency domain (noise level, signal bandwidth). Time frequency analysis attempts to maintain both time and frequency as joint variables to show relationships based on their interactions. 1.2 Scattering Problems Time-frequency representations have been used extensively for analyzing non-stationary systems where the signal spectra varies as a function of time. In some cases., the objective was to characterize individual components of the signal. In other cases, the overall patterns in the time-frequency representation served as the metric for making decisions about the signal [1. 12]. Most of the work using time-frequency techniques has focused on biological and biomedical signals [1, 43]. Very little had been done to address the use of these transforms for problems involving electromagnetic scattering. The remainder of this dissertation explains the application of different time-frequency procedures to several different types of scattering problems. 2

1.2.1 Dispersive Phenomena Traditional signal analysis for signals resulting from electromagnetic scattering assumes stationary. non-dispersive medium. Standard frequency domain techniques cannot properly identify effects which vary as a function of frequency. This is important for most nonidealized scattering targets. The dispersive nature of non-metallic materials and complex scattering modes tends to smear the transformed response. Time-frequency techniques allow the tracking of these effects. 1.2.2 Dynamic Targets The basic Fourier transform assumes signal stationarity in the sense that the frequencies of the signal components do not change with time. For a dynamic target, this is a poor assumption. Using time-frequency techniques to analyze non-stationary targets maintains the signal's relationship to time while displaying how the frequency content of the signal varies with respect to time. This shows how the target movements drive the scattered spectrum. 1.2.3 Noise Corrupted Signals When considering signal processing techniques, the usefulness of these techniques has to be evaluated with respect to contaminated signals. Most radar and remote sensing systems are noise-limited in their performance [36]. It was important to examine how well timefrequency analysis techniques operate in a noise-filled environment. These distributions tend to be very sensitive to noise when trying to track scattering modes containing very low signal levels. 1.2.4 Customized Kernels In most sections of the dissertation, the signals were processed using standard kernels such as the Wigner, the spectrogram, and the binomial. This section looks at the use of customized kernels to generate the time-frequency representation. A genetic algorithm was developed which maximized the third-order Renyi entropy for the time-frequency distribution. This customized kernel provided performance similar to a binomial and consistently exceeded the binomial and the Wigner distribution when evaluated. 3

1.3 Summary Time-frequency techniques clearly demonstrated their value for analyzing radar based signals. For stationary targets examined with a swept frequency signal, these techniques can isolate different scattering centers and scattering modes. For dynamic targets, they provide a view of the target's dynamics unavailable in regular Fourier analysis. 4

CHAPTER 2 BACKGROUND This chapter outlines the basic mathematics developed and used in this research. Starting with the Fourier transform to move between the time and the frequency domains, the idea of changing domains to gain insight for a signal is extended to joint time-frequency and frequency-time domains. The chapter also examines the history of time-frequency analysis and discusses earlier developments and applications of these techniques. 2.1 Traditional Spectral Analysis The Fourier transform and inverse Fourier transform are a set of relationships which allow representing a function in two different domains. The time, or temporal, domain uses time as the independent variable. For the frequency (or spectral) domain, frequency serves as the independent variable. Converting a time function into its frequency domain representation simplifies many operations and provides insight into different aspects of the signal that may be difficult to discern from the time domain representation. 2.1.1 Standard Usage The Fourier transform is given by roo G(f)= g(T)e-j27fTdT (2.1) -00 The corresponding inverse transform is roo g(t) = G(A)ei2AtdA (2.2) -00 Different disciplines employ various permutations of the Fourier transform relationships. By using f to represent the frequency (in hertz), the transform equations above avoid a 5

factor of 1/27r that appears when working with u (in radians/sec) as the frequency variable [15]. For discrete signals, the relationship takes a series format to relate the frequency domain series G(p) with the time domain series g(n) N-1 G(q) = (n)e-jg qn (2.3) n=O where p represents the frequency index and n represents the time index. The corresponding inverse transform is 1 Ni 2r g(t) -= E G(q)ej (2.4) q=0 This discrete form treats the original g(n) as evenly spaced samples of a periodic signal. If the samples are taken AT seconds apart, and there are N samples covering one period of the signal, then in the frequency domain the samples of the spectrum will be spaced 1/AT hertz apart, giving a frequency resolution of AF. The spectrum will also be periodic, with N distinct values. When dealing with two-sided spectra, the frequency range typically is taken from -(N/2)AF to (N/2 - 1)AF. 2.1.2 Power Spectral Density The signal information in the frequency domain is split between the real and imaginary portions of the spectrum. It is often more convenient to work with the power spectral density, a real-valued representation equivalent to the magnitude squared values of the Fourier transform. It can be obtained calculating the spectrum first SX(f) = I_ [x(t)]12 (2.5) or by taking the Fourier transform of the signal's autocorrelation function Sx(f) = F[Rx(T)] (2.6) where the autocorrelation is defined as R(r)( = J X(t 2 )x* (t-2 )dt (2.7) J-00 / 2.1.3 Limitations The Fourier transform has proved itself to be a powerful and useful analysis tool, but it does have limitations when applied to complicated signals. Fourier transforms are poorly 6

suited for non-stationary signals. Although no information is lost when transforming a signal across domains, certain aspects and characteristics can be obscured and difficult to interpret. For instance, a particularly difficult aspect to distinguish with traditional Fourier transforms are changes in the signal's frequency, either due to the signal starting and stopping, or moving from one frequency to another. 2.2 Time-Frequency Analysis Several methods exist for extending the basic ideas of the Fourier transform to adequately represent non-stationary signals. Instead of completely removing the dependence on the original independent variable, it is carried over and combined with a second independent variable. So instead of the transform A[g(x)] = G(y) (2.8) the new transformation gives B[g(x)] = G(x,y) (2.9) There have been many different transformations developed to provide this joint frequency representation. A large family of these, collectively known as Cohen's class, share many characteristics and cover the most popular time-frequency transformations [11, 12]. These transforms can be calculated via the ambiguity domain or via the autocorrelation domain. 2.2.1 Ambiguity Domain Signal Transformations The time-frequency transformations used in this study are members of Cohen's class. These transformations generate a two-dimensional function of time and frequency which represents the signal energy per unit time per unit frequency [11]. Members of this class have the general form [7] 1 00o roo roo / r ^- / Cg(t, WI) = J J J e(-Trw)X( T)g ( + -)g - )dp dTd (2.10) and move from the time-valued function, g(t), to a function representing the distribution of the original signal in a mixed time-frequency domain. The resulting function of time and frequency, Cg(t.w, WI) depends on O(w, r), the transformation kernel chosen. This can be modified slightly to use non-radian frequencies by substituting W -+ 27rf and ( - 27Xv. 7

Applying this to the previous equation gives Cg(t, f, ) eJ(27rv^-27rfT-27rvt) O(, T)g ( + - g ( -- dludTd (2.1 1) A Cohen's class transformation can be viewed as having three steps: moving from g(t) - A(v, Tr), applying a weighting function, and using a double transformation to go from Ag(v T) Cg(t, f). The first step moves the signal from the time domain into a two-dimensional domain, similar to the ambiguity domain used in radar signal analysis [36, 47]. In the radar context, T represents a time shift, and f is a frequency shift. xc X(, f)= J g(t)g*(t+T)ei27rftdt (2.12) -00 The first sub-transformation with Cohen's class gives something slightly different from the definition in equation 2.12. It is called the symmetric ambiguity function and is defined to be Ag9(,T) =- (t + 9 g ( -- e2-dt (2.13) = X*(,)e-3jr7 (2.14) Once the signal is moved to this ambiguity domain, a kernel function, 4(V, T), is applied. The choice of this kernel determines what kind of time-frequency transformation is used. Table 2.1 lists several of the distributions discussed by Cohen [12], along with their transformation kernels and some properties. Distribution Kernel Resolution Cross Name (yv, ) Time Frequency Terms Wigner-Ville 1 High High High Spectrogram f h*(u - ~')e-jvuh(u + ')du Low High Low Exponential e-22/a High High Low Table 2.1: Selected TF Kernels Each of the three distributions listed has desirable properties and liabilities. The spectrogram, also called the short time Fourier transform, involves a direct tradeoff between time and frequency resolution. Depending on the chosen window, h(r), this time-frequency distribution can provide very high frequency resolution but only at the cost of poor temporal resolution. It often provides an unrealistic viewpoint of the time-frequency structure of complex signals [40]. 8

The second distribution, the Wigner, was introduced in the context of quantum mechanics [42] and later adapted and applied to signal processing [39]. The Wigner-Ville distribution in general provides the highest time and frequency resolutions. The drawback when using this form is the presence of significant cross terms between all possible groups of actual signal components. This causes significant clutter and confusion when complicated, multi-component signals are analyzed [13, 14, 19]. The exponential distribution is a Reduced Interference Distributions, or RID. The goal of a RID is to maintain high time and frequency resolution while minimizing the contribution from cross terms [7]. 2.2.2 Autocorrelation Domain Signal Transformations In many cases, it is convenient to work with the transformation in a different form. Instead of moving the signal into the ambiguity domain, the kernel can be applied using a convolution in the autocorrelation domain. In this case, the time frequency distribution is given by [21] roo Cg(t, f; () = Rg(t, T)e-2fTdr (2.15) -00 using Rgr(t, 7) = g + ( u - 4(t - u, )(t- du (2.16) and roo,(tr) = J (v, )e-J2 tdv (2.17) -00 By defining the local autocorrelation function Rg(t,T) = g t + ) g (t- 72 (2.18) this can also be written as roo Cg(t, f; A) = Rg(t, 7) t )(t, 7)e-j2rXfdr (2.19) -00 Knowing the autocorrelation domain representation of the kernel allows one to determine the joint time-frequency representation with only two computationally intensive operations, a convolution and a Fourier transform, instead of the three Fourier transforms required by the original formulation in equation 2.11. 2.2.3 Sampled Signal Analysis Since most work involves sampled data, the transformations take on slightly different discrete forms, with discrete transforms, as shown in equations 2.3 and 2.4, replacing the 9

continuous transforms in equations 2.11 and 2.19. The ambiguity domain based expression becomes C9(,p) E E E 9 (m + ) g* - )- Pl qn ell, 2 qm (2.20) N /=0I q=O m==O 2 2 while the autocorrelation domain expression is given by Cg(n p) g -~ 9 k- - k, I)e-jpl (2.2[) 1=0 k=-oc even I These formulations exhibit the same basic properties as seen in their continuous counterparts [21]. A major difference is the sampled signal is assumed periodic, with the period length matching the signal vector g(n). The overall effect is to generate a two-dimensional time-frequency distribution which is periodic in both directions. Properly defining the summation limits is essential to avoid problems of aliasing. 2.3 Desirable Time-Frequency Distribution Properties When considering different time-frequency distributions, there are several desirable properties we would like the distribution to exhibit. They range in importance, with different kernels satisfying different properties. Ideally, this joint distribution would behave as a power or energy density function and would share properties found in two-dimensional probability density functions. Table 2.2 lists these properties. 2.3.1 Non-Negative Since the time-frequency distribution can be viewed as an energy density distribution, it does not make sense for any values in the time-frequency plane to be negative. If these negative values can be accepted as a numerical artifact, loosening this requirement makes meeting the other properties possible. The only common time-frequency kernel to support this property is the spectrogram [22]. 2.3.2 Real Valued Since the time-frequency distribution is related to the energy spectral density, all values should be real. To achieve this, the time-frequency kernel possesses Hermitian symmetry with respect to v and T in the ambiguity domain 0(iT, -) = 0(-', -T) (2.2:2) 10

PO Non-Negative P1 Real Valued P2 Time Shift Invariant P3 Frequency Shift Invariant P4 Time Marginal P5 Frequency Marginal P6 Instantaneous Frequency P7 Group Delay P8 Time Support P9 Frequency Support P10 Reduced Interference Table 2.2: Desirable Distribution Properties and possesses Hermitian symmetry with respect to T in the autocorrelation domain )(t, T) = *(t,-T) (2.23) 2.3.3 Time Shift Invariant If the original signal is shifted in time, the time-frequency distribution should exhibit an equivalent shift along the time axis. This condition will be met with any kernel in a Cohen's class transformation. 2.3.4 Frequency Shift Invariant If the original signal is shifted in frequency, the time-frequency distribution should exhibit an equivalent shift along the frequency axis. The Cohen's class transformations automatically satisfy this condition regardless of the kernel selected. 2.3.5 Time Marginal Paralleling the concepts of a two-dimensional probability density function, the timefrequency distribution should give the correct value for the instantaneous signal power at any time by integrating across the frequency axis g(t)l= J Cg(t, f)df (2.24) 00For this condition to be met, the ambiguity kernel must be equal to 1 along the frequency For this condition to be met, the ambiguity kernel must be equal to 1 along the frequency 11

axis, and the autocorrelation kernel must be a delta function along the time axis. 0(,O) = 1 (2.25) 0(t, 0) = 6(t) (2.26) 2.3.6 Frequency Marginal As just mentioned for the time marginal, the power spectral density should result after integrating out the dependence on t. roo G(f)12 = Cg(t, f)dt (2.27) -o00 This condition can be met using ambiguity kernels having the value of 1 along the lag axis (v = 0) or equivalently, an autocorrelation domain kernel which when integrated with respect to t will equal 1 for all values of r. 0(0, r) = 1 (2.28) roo J/ (t, T)dr = 1 (2.29) -00 2.3.7 Instantaneous Frequency The average frequency at any time should equal the instantaneous frequency. Treating the time-frequency distribution like a probability distribution allows expressing the average frequency as favg(t) =- J f(tf)df(230) ff_ Cg(t, f)df and the instantaneous frequency 1 &0(t) fi(t) 2 (2.31) where 0(t) is the signal phase at time t. To achieve this property, the kernel must satisfy several conditions. It must provide a proper time marginal (equations 2.25 or 2.26). Also, the partial derivative of the ambiguity domain kernel taken with respect to the lag variable, r, and evaluated at 7 = 0 must be zero for all M 9(v, ir) = 0 (2.32) O7 r=0 12

2.3.8 Group Delay For a given frequency, the average time across the time-frequency distribution should equal the group delay of the original signal. The average time is given by f, tCg(t, f)dt ag f Lug (2 3:3) fu c,(t, f)dt and the group delay t (f) 1 0(f) (2.34) - 27r Of with 0(f) giving the signal phase as a function of frequency. To achieve this property, the kernel must satisfy several conditions. It must provide a proper frequency marginal (equations 2.28 or 2.29). Also, the partial derivative of the ambiguity domain kernel taken with respect to the frequency shift variable, v, and evaluated at v = 0 must be zero for all r 0 (2.35) 9V V=0 2.3.9 Time Support The time support property states if a signal has non-zero values only over the range Itl < tc, then the time-frequency distribution should only have non-zero values over the same range. This property corresponds to an ambiguity domain kernel satisfying ((, r)e-j27vtdf = 0 (2.36) -00 or an autocorrelation kernel where (t, T) =0 (2.37) for 1T| < 21tl. 2.3.10 Frequency Support Frequency support means if the spectrum G(f) is zero for all If > f,c then the timefrequency distribution Cg(t, f) must also be zero when Ifl > fc. The kernel required for this property has an ambiguity domain kernel where -oo J /(v, a)eJ2rf dr = 0 (2.38) for all < 00 for all IvI < 2| |1. 13

2.3.11 Reduced Interference The final desirable property for a time-frequency distribution is to have small interference terms between the signal components. While not precisely defined as the previous properties, this is very important since interference terms can seriously hinder interpretation of the time-frequency distribution. To get reduced interference terms, the ambiguity domain representation of the kernel must behave as a low-pass filter with respect to both the frequency shift, v, and the time lag, T, axes. 2.4 Kernel Design Procedures Although the properties listed in table 2.2 look formidable, it is possible to meet all the properties except PO (non-negativity) through a straightforward design procedure [22]. The basic steps involve taking a smooth elementary function, h(t), defined over the region -1/2 < t < 1/2. The function should enclose an area of 1, be symmetric about t = 0, and go smoothly to zero as Itl approaches 1/2. Once this elementary function is defined, the ambiguity domain kernel is given by the Fourier transform, with the product vr substituted for the frequency variable ( T) = h(t)eJ2f dt (2.39) f f=VT Moving back into the autocorrelation domain, using the properties of Fourier and inverse Fourier transforms, the kernel becomes $(t, T) = (v, T)e-J27'tdv (2.40) | T| ( ) (2.4:L) = h (2.41) 2.5 Alias-Free Formulation When evaluating a discrete autocorrelation function, a problem arises due to the arguments m + 1/2 and m - 1/2. While they cause no problems for continuous signals, a discrete signal is typically treated as undefined (or zero) for non-integer arguments. This limitation has the effect of zeroing out the autocorrelation for all odd lag values. When viewed along the lag axis (1), instead of having P points separated by AT, there are only P/2 points separated by 2AT. The end result is undersampling the distribution by a factor of 2, allowing a frequency span of only -(N/4)AF to (N1/4 - 1)AF) before aliasing occurs. 14

To avoid this problem, it is possible to re-formulate the definition of the autocorrelation function to provide a half step shift to place the values back on the sample points. Phasing problems are avoided by incorporating a corresponding negative half step into the kernel definition. 2.6 Dual of transformations A quick examination of equations 2.1 and 2.2 shows the forward and inverse Fourier transforms differ only in the sign of the exponential term. This means any properties or relationships exhibited by the forward transform will also apply to the inverse transform simply by substituting g(-t) for G(f). With the discrete transforms given by equations 2.3 and 2.4, the same relationship almost holds. The only difference is a scaling factor 1/N that applies to the inverse transform. This similarity between forward and inverse transforms, or duality [15], allows recasting the time-frequency transformations into frequency-time transformations. Given a function of frequency, one may find the temporal response as a function of frequency. Problems of this type are typical in measured electromagnetic scattering data. 2.7 Emphasized Kernels As presented, Cohen's class covers a broad range of time-frequency distributions. Emnphasis has been placed on four different distributions: the spectrogram, the Wigner distribution, the binomial distribution, and a genetically derived distribution. 2.7.1 Spectrogram The spectrogram, also known as the short time Fourier transform, or STFT, provides the most intuitive time frequency distribution. While it can be cast in the standard Cohern's class form with the kernel shown in table 2.1, the usual approach centers a window on the time of interest and takes the Fourier transform of the windowed data. Finding the magnitude squared of the result completes the transformation. The resolution provided by this transform depends on the length of the window. A long window provides good frequency resolution, but it degrades the transform's ability to track the temporal position of the signal. A narrow window improves the ability to determine the time of a signal component accurately, but this comes at the expense of coarse frequency 15

resolution. The spectrogram is best suited for signals with a slowly varying frequency content. While the spectrogram provides a non-negative time-frequency representation, it, does not satisfy the time or frequency marginals. 2.7.2 Wigner Distribution The Wigner distribution, also referred to as the Wigner-Ville distribution, was originally used in the study of quantum statistical mechanics, working with position and momentum instead of time and frequency [42]. Later work by Ville independently derived this transformation for signal processing, using time and frequency [39]. The ambiguity domain representation for the Wigner is a delta function along the time axis, t, and is independent of the time lag, T. The delta functional form of the Wigner generates two properties for the distribution. Having a delta function in the kernel makes it unnecessary to carry out the convolution with the kernel, making the transform quick to calculate. This also provides the highest possible resolution in the time frequency domain. The Wigner distribution possesses nine of the eleven desirable properties listed in table 2.2. The only two it fails to meet are non-negativity, PO, and it is not a reduced interference distribution, P10. A drawback of the Wigner distribution is the lack of any smoothing in the time-frequency domain. Large cross terms exist even when the true comrnponents are clearly separated. While this is tolerable when few signal components exist, the number of cross terms, N, generated from N actual components is given by N(N - 1)(2.42) NV = (2.42) 2 2.7.3 Binomial Distribution The binomial distribution is typically generated using a discrete autocorrelation domain kernel which is closely related to the continuous exponential distribution in table 2.1. The binomial distribution sacrifices some of the resolution provided by the Wigner to reduce the cross terms, similar to how a tapered antenna array trades off pattern beamwidth for reduced sidelobes. The binomial kernel gives a reduced interference distribution which satisfies all the desired properties except non-negativity. 16

2.7.4 Genetic Kernel Chapter 7 discusses a genetic algorithm used to design a time-frequency distribution. The algorithm selected coefficients for a sinusoidal representation of the elementary function h(t) as outlined in section 2.4. The distributions found via this method were very similar to those provided by the binomial. Like the binomial, they met all the desirable properties except non-negativity. This kernel exceeded the binomial performance when evaluating the entropy of the resulting distribution. Chapter 7 gives more information on this kernel. 2.8 Review of Time-Frequency Developments The basics of time-frequency techniques trace back to Fourier analysis through the short time Fourier transform. The spectrogram began seeing extensive use for speech signals in the 1940's [8, 31] due to the inadequacy of simple Fourier transforms for analyzing nonstationary signals. In 1932, Wigner introduced a means of estimating momentum and position in quantum mechanics [42]. The basic techniques were applied in a signal processing context by Ville in 1948 [39], leading this approach to often be referred to as the Wigner-Ville distribution [32]. This had the advantage of much higher resolution than the spectrogram. Alternative representations for the joint time-frequency domain were developed through the 1950's and 1960's with researchers such as Page, Levin, Rhihaczek, and Born-Jordan. The widely varying characteristics of these different distributions were consolidated in 1966 by Cohen into a general class of transformation for joint distributions [11]. With this representation, an arbitrary kernel q(0, T) differentiated the different distributions. This formulation allowed analysis and comparisons between kernels and led to experimentation on new kernels [12]. In 1980, Claasen and Mecklenbrauker brought the Wigner distribution to the forefront with a series of articles discussing the strengths, weaknesses, and applications for this analysis approach [8, 9, 10]. A lot of recent work has led to a better understanding of the fundamental limits of performance for Cohen's class distribution and finding kernel and computational techniques to provide the clearest time-frequency representations [2, 18, 22, 21, 45]. 17

2.9 Time-Frequency Applications Time-frequency representations have been used extensively for analyzing non-stationary systems where the signal spectra varies as a function of time. In some cases, the objective was to characterize individual components of the signal. In other cases, the overall patterns in the time-frequency representation served as the metric for making decisions about the signal. 2.9.1 Marine Mammal A large amount of work has been conducted analyzing the whistles and other sounds generated by marine mammals. Decomposing these signals in the time-frequency domain shows dolphins and whales using frequency modulated sonar signals employing sophisticated techniques [1]. 2.9.2 Biomedical Another area where time-frequency techniques have been applied heavily is the analysis of other biological signals. Some of these signals have been acoustic, such as differentiating benign jaw clicks from those requiring corrective treatment [1]. Other medical uses have examined the body's own electrical signal, such as trying to characterize the electroencephalogram (EEG) signals of brainwave activities associated with epileptic seizures [1]. The heart's electrical signals, measured with an electrocardiogram (EKG) have been analyzed to try and differentiate the signals of a healthy heart from a diseased heart [43]. 2.9.3 Mechanical Engineering Time-frequency techniques are also finding uses in disciplines such as mechanical engineering. Researchers have monitored vibrations from cutting tools, trying to determine characteristics in the tool's chatter that may indicate imminent failure [1]. 2.9.4 Electromagnetics Cohen's class of time-frequency distributions have not been used heavily in the field of electromagnetics. Some work has employed spectrograms and wavelet analysis [23, 25, 281, but high resolution time-frequency techniques like the Wigner and the binomial have been largely ignored. The typical application has been the analysis of scattering modes when analyzing a target illuminated by a wideband system. 18

2.10 Summary Time-frequency analysis has been largely ignored in the study of electromagnetic scattering. While some related work has been done, the emphasis has been on using the wavelet transform instead of Cohen's class distributions such as the Wigner or binomial distributions. The typically frequency dependent nature of scattering and the time-varying spectra associated with dynamic targets make time-frequency techniques an attractive tool for working with scattering problems. 19

CHAPTER 3 STATIONARY TARGETS 3.1 Introduction Anechoic chambers can be used to measure both the spatial pattern and the frequency response of a target. As long as the measurement system provides a stable, adjustable frequency source, the chamber can also provide the information necessary to evaluate the impulse response of the target and chamber. 3.2 Background Anechoic chamber measurements to determine the pulse response of a target normally do not use an actual pulse. Rather than deal n with the physical difficulties of generating high power, short duration pulses, a measurement system illuminates the target with a stepped continuous wave (CW) signal. Instead of trying to accurately record and analyze a pulse return in real time, the frequency is selected, the system is allowed to reach steady state, and a complex-valued data point (magnitude and phase), G(fn) is collected before moving on to the next frequency [24, 41]. Using the inverse Fourier transform gives the target and chamber response as a function of time g(tn) = F-1 [G(fn)] (3.1) which approximates the impulse response. If N samples are taken, each separated by Af hertz, these samples can be transformed to the time domain where they will provide N samples of the impulse response over a total temporal span of I/Af seconds. Assuming the wave propagates at a velocity equal to the speed of light in free space, 20

= c = 3 x 108m/s. This means the time response can be expressed as a range response g(r) = 2g(t) (3.2) This same approach to finding the time response can also be applied to numerically simulated results. Finding the spectral response provides the sampled impulse response of a virtual anechoic chamber with a maximum unambiguous range Rmax = 2\f (3.3) Measurements within this virtual chamber have a range resolution determined by the maximum and minimum frequencies used in calculating the transform. For N frequencies spaced Af apart, the range resolution is 2(N- )Af (3.4) 3.3 Duality and Time-Frequency Transformations The typical time-frequency transformations move a function of time into a time-frequency space using the Cohen's class transformation C9(t, f,) /_) /_ eJ27r(u-T-^(t)^)(,r T)g(u + -) g* (u - -)dudrd (3.5) This can also be viewed as a Fourier transform to convert the signal's local autocorrelation into an ambiguity function, applying a weighting function, or kernel, then doing a twodimensional Fourier transform to move from the ambiguity domain to the time-frequency domain. Cg(t f) = Fr1 (F- {Fuv) g (+ ) ( - ( )}) (3.6) Typically, instead of transforming the local autocorrelation into the ambiguity domain, the kernel is applied by performing a convolution along the time axis, t, in the autocorrelation domain. After convolving, taking the Fourier transform of the lag variable, Tr, yields the time frequency distribution. Cg(t, f) = FTf [R9(t, r) (t b)(t, T)] (3.7) with Rg(t, T) as the local autocorrelation and Ot representing the convolution with respect to t. When beginning with a frequency domain signal, G(f), the spectral autocorrelation is RG(f P) = G f + 2 G* f- ) (3 8) 21

Taking the inverse Fourier transform with respect to f converts the spectral autocorrelation into the same ambiguity domain representation that would have resulted if starting with g(t). This means the time-frequency representation for G(f) is given by C-(t f) = Srf (-\ [^\ [G (p+ ) G -)] T}) (3.9) This result can also be obtained by applying an inverse Fourier transform on the frequency shift variable, v, from the spectral autocorrelation after convolving with the spectral autocorrelation kernel, P(Lv,f) CG (t, f) = Ft' [RG(V, f) ~f '(V, f)] (3.10) Working through these transformations, the duality between forward and inverse Fourier transforms provides a simplification. Substituting G(f) for g(t) and working through the temporal time-frequency distribution steps will generate the result CG(f,t). This simple switching of axis works provided b(f I,) = 1(jf) (3.11) This rotational symmetry will happen whenever the ambiguity domain kernel is a product kernel such that (V, T) = q(VT) (3.12) All the kernels used in this chapter satisfy that condition. 3.4 Scattering and Dispersion When considering the scattered radar return from a target, a variety of modes may contribute. Based on the target composition and physical arrangement it is possible to separate the returns from scattering centers which are in different range bins on the target. Scattering modes may be divided into dispersive and non-dispersive modes. A nondispersive mode will behave identically at all frequencies. A dispersive mode will depend on a mechanism which is frequency dependent [33]. When processed using a standard Fourier transform, sharp resolution is possible on nondispersive modes. The energy returned through the dispersive modes will appear spread across a band of range bins since the signal travels at a non-constant velocity which changes as a function of frequency. However, using a frequency-time transformation clearly shows the difference between the dispersive and non-dispersive modes, and can also differentiate various dispersive modes. 22

3.4.1 Reflective Returns A specular return follows from geometric optics and is found by tracing simple reflections off the target surface with the angle of reflection being equivalent to the angle of incidence relative to the local normal on the target. These returns can involve single reflections directly back to the radar or multiple bounce reflections. Over the range of frequencies considered these modes are non-dispersive. 3.4.2 Creeping Waves Another scattering mode considered involves creeping waves which attach themselves to the target and travel around the target before shedding energy back toward the radar. As these waves travel around the target, the velocity of propagation as well as the energy shed toward the radar will vary as a function of frequency. At low frequencies, the waves travel closely attached to the surface, while at higher frequencies the waves are offset slightly from the surface. Using a frequency-time transform, the position of these scattering centers will move with frequency. The slope will be determined by how far the energy travels while in a dispersive mode. These modes are difficult to track accurately because they carry much less energy than the specular returns for simple targets. 3.4.3 Combined Modes Finally, there are combined modes to consider. These modes have energy travel a path that is sometimes specular and sometimes dispersive. They can also be categorized by how much of the energy path is along a dispersive mode. 3.5 Transforms Considered For this study, three different kernels were used to perform the frequency-time transformation. The kernels were chosen because of their widespread use, ease of implementation, and the variety in their distribution characteristics. 3.5.1 Spectrogram This transform is also referred to as the short time Fourier transform, or STFT. It involves moving a window in time across the signal and taking the transform for each 23

window. This transform gives a non-negative distribution with good noise immunity arid low cross terms. The distribution satisfies neither marginal characteristic and generally provides poor resolution. 3.5.2 Wigner While the spectrogram is one of the most intuitive time-frequency distributions, the Wigner employs the simplest kernel, a delta function. The result is a very high resolution technique which suffers from high interference terms. 3.5.3 Binomial The binomial distribution was developed as a general purpose kernel to provide high resolution with low interference terms. While the binomial cannot match the resolution of the Wigner, the lack of large interference cross terms typically results in a clearer joint domain representation, especially with complicated signals 3.6 Target Geometries For this study, the basic target was a 15 centimeter radius perfectly conducting sphere. It was chosen due to the canonical nature of the target. It was simple to numerically model, and the scattering from the sphere is relatively simple and well understood. The sphere scatters energy using two basic mechanisms: a simple reflection governed by geometrical optics, and dispersive scattering based on creeping waves traveling around the surface of the sphere. 3.6.1 Single Sphere The simplest case, shown in figure 3.1, consists of a single perfectly conducting sphere in free space illuminated by an incident plane wave. The sphere modeled had a radius, r, of 15 centimeters. The solution for the scattered field was obtained by evaluating the Mie series [46] at the frequencies of interest. The expected returns consist of a reflective return off the front of the sphere which appears a distance r ahead of the phase reference at the sphere's center. Most of the reflected energy comes from this mode. The next contributor is the first creeping wave which travels an additional 7rr beyond the center of the sphere. Since the range values are based on two-way propagation times, this scattering mode appears 7rr/2 behind the sphere center. The other detectable mode involves a second creeping wave, one 24

15 cm 7 radar incident scattered Figure 3.1: Single Sphere which travels one-and-a-half times around the sphere before shedding energy back towards the radar. Since the traveling wave continually sheds energy as it travels around the sphere, the actual impact of this component is negligible. 3.6.2 Two Spheres: In-line The second geometry, shown in figure 3.2 involved two spheres, each with a radius, r, 1 meter O radar -- -- incident scattered 15 cm Figure 3.2: Two Spheres: In-Line of 15 centimeters, placed with a distance, d, of 1 meter between their centers. Adding the interactions between the spheres leads to several more modes making contributions to the scattered signal. The primary scattering centers are the front face of each sphere, located r in front of the phase center (the center of the leading sphere) and d - r behind the phase center. The next modes involve the first creeping wave for each sphere. These appear rr/'2 behind each sphere center, or at ranges of 7rr/2 and d - irr/2 behind the phase center. When examining the returned signal, the scattering modes which involve creeping waves 25

will not have a stationary range in the frequency-time distribution. The energy in these dispersive modes will also drop off as the frequency increases. 3.6.3 Broadside In the third sphere setup, the line between the spheres is oriented perpendicular to the direction of wave propagation. For this geometry, the polarization of the illuminating wave becomes significant when considering the interactions between the spheres due to the behavior of creeping waves. Figure 3.3 shows the two spheres at a broadside orientation to the source..As seen.n s /radar. incident scattered 15 cm Figure 3.3: Two Spheres: Broadside the previous two models, the specular reflection off the front edge of the sphere will be the first scattering mode. In this case, the energy off both spheres will reinforce each other and provide a single large return a distance r in front of the phase reference. The combining effect will also give a large creeping wave return. Both these modes are essentially single sphere modes happening independently with each sphere. They do not depend on the incident polarization. The closest scattering mode to involve interactions between the spheres is the reflection of energy off one sphere towards the second, then off the second back in the direction of the source. This mode has an apparent position (d/2) - r behind the phase reference. The remaining modes are summarized in figure 3.3. Dule to the sensitivity of some modes to the polarization of the illuminating field, two 26

trials were computed. The first trial oriented the E-field parallel to the line connecting the spheres. This allowed sphere interactions based on creeping waves. The second trial oriented the E-field perpendicular to the line between the spheres. This essentially drove the creeping waves' contributions to cross-sphere effects to zero. 3.7 Simulation Parameters The targets outlined above were simulated numerically using two different packages: a Mie series code for the single sphere and a body of rotation code for the two-sphere cases. In each case, the monostatic backscatter measured included the magnitude and phase of the electric field. For the two-sphere targets, the spectrum simulated covered the band from 1.00 gigahertz through 7.55 gigahertz in 50 megahertz steps. This provided 132 samples for processing. The 50 megahertz sampling interval provided a maximum unambiguous range scale of 3 meters using the relationships 1 Tmaz = (3.13) Af 1 (3.14) 50 MHz = 20 nsec (3.15) and c Rmax = Tmax (3.16) 2 = 1.5 x 108- 20 nsec (3.17) sec = 3m (3.18) The range resolution, Ar, depends on the span of frequencies covered by the spectral data. With spectral information covering the range fmin to fmax, the inverse Fourier transform can provide range resolution based on Ar = - (3.19) 2 (fmax - fmin) assuming the signal propagates at the speed of light. 27

For the two-sphere targets, when fmax = 7.55 GHz and fmin = 1.00 GHz, the spatial resolution, Ar, is 2.29 centimeters. In general, the resolution in the frequency-time domain will depend on the length of the transform used after applying the selected kernel in the autocorrelation domain. With a frequency sampling interval of Af, an N point transform equates to Ar = 2(N - )f (3.20) The single sphere was simulated over a band of frequencies starting at 1.005 gigahertz and extending to 8.865 gigahertz with a sampling interval, Af, of 60 megahertz. These frequency values provide a maximum unambiguous range, Rmax, of 2.5 meters with a range resolution, Ar, of 1.91 centimeters. 3.8 Results 3.8.1 Single Sphere The single sphere did not provide a large amount of information. While there are a number of scattering modes to examine, the dynamic range between the specular return and the creeping waves makes even the first creeping wave difficult to see, and the second time around creeping wave is not discernable without special processing. This second time wave was only visible using double precision math for all calculation along with increasing the frequency sampling interval and the band of frequencies swept. Spectrum The initial spectrum for a sphere starts low and oscillates through the resonance region as it settles towards the final value (figure 3.4). These oscillations are due to constructive and destructive interference between the reflection from the front of the sphere and the creeping wave travelling around the sphere. The oscillations decay since the creeping wave intensity decreases with frequency. At first glance, the phase response is essentially linear as the number of wavelengths to the front of the sphere increases with frequency. However, close examination shows small deviations due to contributions from the creeping wave. Inverse Fourier Transform Applying an inverse Fourier transform and scaling the horizontal axis by the two-way propagation time provides a range profile that clearly shows the leading edge of the sphere 28

Single Sphere Spectrum 151 14 - -2 - 2 4 6 8 2 4 6 8 Frequency (GHz) Frequency (GHz) Figure 3.4: Single Sphere: Spectral Response 0.15 meters in front of the phase reference (see figure 3.5). This plot also clearly shows the creeping wave at its expected positions of -rr/2 = -0.236 meters. The second-time around creeping wave should occupy a position at -37rr/2 = -0.707 meters, but it is not visible in this plot due the extremely small contribution it makes over this frequency range. Spectrogram The spectrogram for the single sphere is shown in figure 3.6. Here the specular return shows as a very clear band located 0.15 meters ahead of the phase reference, while the creeping wave is located at an approximate position of 0.24 meters behind the phase reference. The creeping wave component is fairly clear at the lowest frequencies and becomes less pronounced as the frequency increases. Once again, the second-time around creeping wave lacks enough energy to be clearly discernible, but the distribution does exhibit a very slight line at the correct range. Wigner Distribution The Wigner distribution for the single sphere is shown in figure 3.7. Like the spectrogram, the Wigner distribution clearly shows the specular and the creeping wave re 29

Single Sphere - IFT -10H-20 F- -30 -E -50K -60 F -70 - Figure 3.5: Single Sphere: Inverse Fourier Transform Single Sphere Spectrogram 2 () E 4a 2 j 4 D 6 / Frequency (GHz) Figure 3.6: Single Sphere: Spectrogram 30

05 0.5 11 - - i j1 2 3 4 5 6 7 8 Frequency (GHz) Figure 3.7: Single Sphere: Wigner Distribution turn. The major difference is the appearance of very large cross terms between the actual signal components. The actual signal component should be the specular reflection at +0.15 meters, the creeping wave at -0.235 meters, and the second-time creeping wave at -0.707 meters. These components give rise to cross terms at (0.15 - 0.235),/2 = -0.0423, (0.15 - 0.707)/2 = -0.278, and (-0.235 - 0.707)/2 = -0.471 meters. This transform is periodic along the frequency axis, so additional cross terms appear due to the repetitions of the three signal components. The specular and the specular's periodic copy cause the large cross term which appears at -1.175 meters. The creeping wave, although smaller than a nearby cross term, shows three expected features. In addition to appearing at the correct range, the component has a decrease in intensity as the signal frequency rises, and the creeping wave line also has a slight upward slope, indicating the propagation velocity for that mode varies with frequency, slowing down slightly as the frequency increases. Binomial Distribution The binomial distribution, shown in figure 3.8 shares some characteristics with the Wigner, but contains significantly reduced cross terms. While the cross terms between the 31

Single Sphere Binomial 0.5... -.:. -. " i,- - - ~~{:i-::: —::.W,,^:^ ^^ i. 2 3 4 5 6 7 8 Frequency (GHz) Figure 3.8: Single Sphere: Binomial Distribution specular and the creeping wave which appears at -.043 meters still exists, there are no visible cross terms due to the periodic images of the components. This representation also shows the changing intensity and position of the creeping wave as seen with the Wigner. Missing from this distribution is any visible indication of the second-time creeping wave. Single Sphere Summary Since the variation in return with respect to frequency was small, the spectrogram actually provided a relatively sharp view of the target in the joint frequency-time domain. It was unable to discern the changing propagation velocity in the creeping wave mode which was evident in the Wigner and binomial distributions. 3.8.2 Two Spheres: In-Line Orientation When positioned one behind the other, two spheres should provide two specular scattering modes and two readily discernible creeping wave modes. These spheres will also scatter based on the interaction of the spheres. The primary interaction between the spheres is the wave that reflects off the back sphere, travels back to the front sphere where it reflects off the rear surface of the sphere, then reflects off the back sphere a second time before returning 32

to the radar. Placing the phase reference at the midpoint between the two spheres means the specular returns should appear at distances of +0.65 and -0.35 meters, the creeping wave returns at +0.265 and -0.735 meters, and the triple reflection return at -1.05 meters. Spectrum Figure 3.9 gives the spectral return. The spectral response shows rapid variations in the In-Line Spheres Spectrum c m 'a 0) 2 a. MU Frequency (GHz) Figure 3.9: In-Line Spheres: Spectral Response scattered signal as the components combine in and out of phase with changes in wavelength of the illuminating plane wave. With the single sphere, these oscillations settled out as the frequency increased because only one significant mode existed. With the two spheres, the creeping wave contribution also drops off with increasing frequency but the specular returns continue to make significant contributions throughout the frequency sweep. Inverse Fourier Transform The inverse Fourier transform for this signal shows five clearly discernible peaks (figure 3.10). The first four peaks correspond to the speculars off the front of each sphere (at 0.65 meters and -0.35 meters), the creeping wave around each sphere (at.264 meters and -.736 meters). The return at -1.05 meters is due to part of the specular return off the 33

!n-Line Spheres - IFFT -51 '- I I I,I 10 I I —.5 -1 -0.5 0 0.5 1 1.5 Range (meters) Figure 3.10: In-Line Spheres: Inverse Fourier Transform second sphere being reflected by the first sphere before scattering back to the source. Spectrogram Viewed with the spectrogram (figure 3.11), the two speculars and two creeping waves are clearly visible. The doubly reflected mode also is visible. As the frequency increases, there is a slight drop in intensity for the creeping wave modes. The apparent position for the creeping wave also moves slightly over the frequency range since the creeping wave travels around the spheres faster at lower frequencies. Wigner Distribution The Wigner distribution for this signal, shown in figure 3.12, suffers from many large interference terms. All five components, two speculars, two creeping waves, and the triple bounce, along with their periodic copies combine to generate cross terms. All this interference makes it difficult to directly view the lower energy scattering components. Binomial Distribution The binomial distribution shown in figure 3.13 illustrates how much clearer a reduced interference distribution can be. In this case, most of the cross terms that cluttered the 34

In-Line Spectrogramn E -LU a) E a) m c co T Frequency (GHz) Figure 3.11: In-Line Spheres: Spectrogram In-LUne Wigner -0 5 4 Frequency (GHz) 6 7 Figure 3.12: In-Line Spheres: Wigner Distribution 35

In-Line Binomial... I:. ==.= = ==,c-=....:x: x::X,> Y x x:::X::c::::::::::.::::::. ' = =...:_.::::: -r-: -:c x i - o -:.:.::: ---: -: -—:; ':_.'i:'.~ ' '.:-.9.:c*-e- 89.......... ~......... _,,; _~;;;;:::zc: *.<:4.....*..:.:.:.::.:..::~.:.:.' c~ - <*:.: *~ * *S#.*~ S a <., c* * ^V i":...:.c:>.: cr.::Xi:-.d -:: k * 8 is: a -* - a:* & >:.*~ 41 >% >w:~ *: i:: -: -:.::. *.~: 9.;:5e- F:~ s: &- * Q q.s id a -r -ar.~~ ~~ -i:5,~ ~^ ^ - -: *:: *: fi.:i e Z:r~ ~R *:~ s; 9. 9 ~-9 *.# a ~r js a xS.-.a -.0 ~.s asA r~ V~ 'X: I.' '. 56+l B^<. > - -w 8. PC~ X~ a a L-l -I ~~ * - ^.-XZ - <*. ^.~ *B *ao 0 4* ~Z~ i; x *:: * $~~~~~~~:Q ~:Cx. e~:eX:stc*:z0:S:4. Xo;=o 'X c;:'I 'R. t ~~:. m zee p::-a-X f Y:~:-.::~-)*: z: ~ 9: Rn )10::11: 99.......................... ~ ~ DlRI,- 1~ -9 a: k a: ar —.,Figure 3.13,:..In-Line Spheres: Binomial Distribution:' 0.5 1 2 3 4 5 6 7 Frequency (GHz) Figure 3.13: In-Line Spheres: Binomial Distribution Wigner distribution have been suppressed. The specular reflections are clearly shown, and the creeping waves off each sphere can also be seen. The triple bounce component at -1.05 meters is too small to show up in this distribution. In-Line Sphere Summary With this geometry, the main contribution comes from the front of each sphere. The creeping waves of the two spheres provide a small return best seen in the spectrogram. The spectrogram also shows the doubly reflected component at -1.05 meters. 3.8.3 Two Spheres: Broadside Orientation - Perpendicular Polarization Positioning the two spheres broadside to the radar changes the number and types of modes expected in the returned signal. In addition to the returns from a single sphere, several modes exist to transfer energy between the spheres before sending it back towards the radar. The signal polarization also starts to play a major role in this geometry. Unlike the single sphere and the in-line spheres, this problem is no longer axially symmetric with respect to the direction of propagation. The first case considered will orient the electric 36

field perpendicular to the line connecting the sphere centers. This means creeping waves will not scatter energy from one sphere in the direction of the other. All sphere interactions with this polarization are based on reflections as governed by geometric optics. Each sphere will act independently to provide the creeping wave component as would be seen for a single sphere. Spectrum As shown in figure 3.14, the spectral response is similar in some ways to the single sphere Broadside - Perpendicular Pol Spectrum 0.71 / 1 c 0.55 0.5 0.45 0.4.s 0 - Figure 3.14: Broadside - Perpendicular Polarization: Spectral Response response. The reflections off the front and the creeping waves off each sphere reinforce each other to provide a reinforced single sphere type return. The reflections travelling between the spheres cause this spectrum to deviate from the single sphere response. Inverse Fourier Transform Figure 3.15 shows the inverse Fourier transform for the broadside spheres. The main peaks correspond to the front of the spheres at +0.15 meters, the creeping waves at -0.24, and a reflected mode that bounces once off each sphere before returning to the radar at -0.64 meters. 37

Broadside - Perpendicular Pol - IFT -10_I, -,,l, - -20 -30 -2,-40 -C -50 - -60 -70 -80, ----1.5 -1 -0.5 0 0.5 1 1.5 Range (meters) Figure 3.15: Broadside - Perpendicular Polarization: Inverse Fourier Transform Spectrogram The spectrogram in figure 3.16 shows the same three scattering modes discussed in the previous section due to the front of the spheres, the creeping waves; and a double reflection return. In addition to showing the position for these returns, the drop in creeping wave intensity as the signal frequency increases is also evident. Wigner Distribution Once again, the Wigner distribution, show in figure 3.17 shows the strong return off the front of the spheres very clearly, but it has trouble with losing the weaker components in the cross terms. The creeping wave shows up at the low end of the frequency sweep, but it drops off quickly. The cross terms between the front of the sphere and the creeping wave returns are evident across the entire frequency range. Binomial Distribution The reduced cross term levels of the binomial distribution, shown in figure 3.18, make it much easier to view the creeping wave behavior. In particular, the drop in intensity and the change in apparent position are both evident in the binomial distribution. 38

Broadside Pemrpendicular Spectroaram - * 0.5 1 2 3 4 5 6 7 Frequency (GHz) Figure 3.16: Broadside - Perpendicular Polarization: Spectrogram Broadside Perpendicular Wigner -1 -0.51 E g, nmn a 0l 0.51 Figure 3.17: Broadside - Perpendicular Polarization: Wigner 39

Broadside Perpendicular Binomial.......:: j::......................... -o:,, -......,.,,.,,.,<.,:'"'Z'"::.:.,.,.,.,:....................... 1F 8.::.:.:.:.:.:............................ -::.:`::.:.:.:.:..;.~.:.:.~.... -O.05 -~ _ _ x< _ _ _ _X _ _- _ -- _ _-__ _ __.'__"E.'.-.' i, co n1 2 3 4 5 6 7 Frequency (GHz) Figure 3.18: Broadside - Perpendicular Polarization: Binomial Broadside With Perpendicular Polarization Summary When illuminated with this polarization, the broadside spheres acted very similarly to a single sphere. Although modes existed which representated interactions between the two spheres, their contribution was relatively small, especially when viewed in the timefrequency domain. 3.8.4 Two Spheres: Broadside Orientation - Aligned Polarization Keeping the same target position and orientation, the radar system was rotated ninety degrees so the incident electric field was in alignment with the line between the spheres' centers. This polarization will support all the scattering modes seen in previous section but will generate additional modes using creeping waves to transfer energy between the spheres. Spectrum The spectral response in figure 3.19 is very similar for the broadside spheres, with the phase response virtually identical. Once again, the magnitude of the spectrum shares some characteristics with the single sphere but behaves differently as the frequency increases. 3.. 4 T wo Spheres- B r a s d e O i n a t o - A i g e:oarzto:'pe<'.t<:.::. 40

Broadside - Aligned Pol Spectrum g a) I V1 CI.'D a).rCL Figure 3.19: Broadside - Aligned Polarization: Spectral Response Inverse Fourier Transform The inverse Fourier transform show in figure 3.20 contains the same peaks seen with the other polarization. It shows the front of the spheres, the creeping wave, and the reflection from one sphere, off the second, and back to the radar. Three additional modes appear in this plot. The return at -0.53 meters represents energy reflecting from the first sphere toward the second, then traveling around the back of the second sphere before returning to the radar. The second new mode at -0.74 meters follows a creeping wave 90 degrees around to the back of the first sphere, travels to the second sphere, and reattaches as a creeping wave for another 90 degrees before heading back to the radar. The final new mode employs a creeping wave around the first sphere, then reflects off the second sphere back towards the first where it reflects back towards the radar. This mode appears at a range of -0.88 meters. Spectrogram Figure 3.21 shows the spectrogram for the broadside spheres with the incident wave aligned with the line between the spheres. Several components appear in this distribution at ranges of 0.15, -0.24, -0.35, -0.53, and -0.63 meters. 41

Broadside - Aligned Pol - IFT 1 n -). -20 -30-. -40 Range (meters) Figure 3.20: Broadside - Aligned Polarization: Inverse Fourier Transform Broadside Aligned Spectrogram -1.5 -80 -1.5 -1 -0.5 0 0. 5 1.5 Frequencyge (Hzmeters) Figure 3.20: Broadside - Aligned Polarization: Inverse Fourier TSpectrogransform - 1 -- 1 2 3 4 5 6 7 Frequency (GHz) Figure 3.21: Broadside- Aligned Polarization: Spectrogram 42

Wigner Distribution The large number of scattering modes leads to many cross terms in the Wigner distribution. shown in figure 3.22 The only clearly illustrated modes are the specular reflection off Broadside Aligned Wigner -1.5 __ _ _- _ _ _ __ _, _X:X.} -IM -:.,:.:.::-.-:-~4 0X4O~. -:-: -0.5 - -------------------- 12 3 4 5 6 7 Frequency (GHz) Figure 3.22: Broadside - Aligned Polarization: Wigner the front of the spheres and the creeping waves. There is a slight indication of the multiple bounce reflective component at -0.64 meters. Binomial Distribution The binomial distribution in figure 3.23 does not have large cross terms, but only three modes show clearly: the sphere fronts at 0.15 meters, the creeping waves, initially at -0.24 meters, and the combined specular/creeping wave mode at -0.53 meters. The binomial does clearly illustrate the dispersive nature of the creeping wave mode. Broadside Aligned Polarization Summary As expected, more interactions were evident between the two spheres using this polarization. This was due to the direction of travel for the creeping waves around the spheres. 43

Broadside Aligned Binomial:I-,.......... -.....:........:....:.:..:. -0.5 aj 2 cr 0.5,...:,...::::..,,~~,........:.:...:......... l- WAWA-.... -'.:.....::'...':.,...::., r.:..-a......s-.. F I it & Mr= - 8 aka P II M -M. -M or '-&W -'...*A -,.: I I... I 1..F... I. !.. I I .: m-'.... I I -%.:,..I::' 1 I 2 3 4 5 Frequency (GHz) 6 7 Figure 3.23: Broadside - Aligned Polarization: Binomial 3.9 Summary Applying Cohen's class distributions to generate time-frequency distributions using swept frequency signals was successful and informative. This technique provided information on the targets and their scattering modes not readily seen using either the frequency or the time domain representations of the signal. By using time-frequency techniques, dispersive effects were highlighted, even though the dispersion involved with these targets was small. 44

CHAPTER 4 STATISTICS FOR COHEN'S CLASS TFD OF WHITE GAUSSIAN NOISE 4.1 Introduction Time-frequency analysis is based on transforming a signal from either a function of time, g(t), or frequency, G(f), into a function of both time and frequency, Cg(t, f). This representation has the advantage of capturing the features of non-stationary signals and has been heavily used with acoustic signals[l]. Two frequently used time-frequency transformations generate the Wigner distribution and the spectrogram [37]. These distributions belong to a larger class of distributions referred to as Cohen's class [12]. The derivation here will focus on this general case. This chapter examines the mean and variance for the time-frequency distribution of a noise signal. Both analytical and numerical values are presented for a sampled signal processed with discrete transformations. This analysis assumes a complex noise only signal, g(t), comprised of two independent, identically distributed random variables, gr(t) and gi(t), which are normally distributed with a white spectra. This means g(t) is independent of g(r) unless t = r. The mean, variance, autocorrelation, and spectrum for the noise signal are defined as: E[g(t)] = 0 (4.1) var[g(t)] = 2 (4.2) Rg () = a26() (4.3) G(f) 12 = ir[g(t)]12 = 2 (4.4) This white noise assumption creates problems with the continuous signal since it implies infinite signal bandwidth and infinite power. A finite result is possible when using a kernel 45

with finite extent in the ambiguity domain, but many transforms, such as the Wigner, binomial and spectrogram, do not satisfy this constraint. However, the band-limited nature of a discrete signal will allow this for any transformation kernel. Throughout this paper, all integrations will have limits going from -oc to +,0o unless otherwise specified. With the discrete formulas, the indices run from 0 to N -- 1, where N is the number of data points considered or 0 to P - 1, where P is the length of the transform used. 4.2 Continuous Transform A continuous time signal, g(t), can be represented in the joint time-frequency domain using a transformation belonging to Cohen's class [12], expressed as Cg(t, f) = g (u+ 2 g ( - 2 (, T)eji2('u-t-f)dudvdr (4.5) The three integrations in the preceding equation represent two forivard and one inverse Fourier transforms. Using the relationships Ftf[x(t)] = Jx(T)e-i27rfrd (4.6) = X(f) (4.7) ftl[X(f)] = fX(A)eJ2AtdA (4.8) = x(t) (4.9) allows re-writing the transformation as Cg(t, f) = Ff ( t {u(-.v) [9 (u + ) 9* (u - ) (, Tv)}) (4.10) Cohen's class transformations cover a wide range of time-frequency distributions, including the spectrogram, the Wigner-Ville, and the reduced interference distribution (RID) [29]. From the view of implementation, this form involves three Fourier transformations and a simple multiply. Since multiplies in one domain correspond to convolutions in the other domain, equation 4.10 can be re-written as Cg(t, f) = (u + u - g u (t - u, T)e-J2 fdudT (4.11) or, in a more compact form Cg(t, f) = Frf [Rg(t, r) Ot )(t, 7)] (4.12) 46

with Rg(t, T) as the Local autocorrelation, Ot representing the convolution with respect to t, and p(t, T) as the autocorrelation domain kernel. 4.2.1 Continuous Mean Performing a Cohen's class transformation on this signal gives a time-frequency distribution Cg(t, f). Start with the expected value of this distribution E(t,f)[Cq(t. )] = E(tf) [J g (u + ) - ) (V, T)e2(u-t-f)dudvdr] (4.13) where the operator E(tf)[Cg(t, f)] represents the expectation of Cg(t, f) and can be written explicitly as E(t,f)[Cg(t,f)] = f C (t, f)fg(t, f)dtdf (4.14) using fg(t, f) to represent the probability distribution function for the noise signal. The expectation operation can move inside the integrals over r and v since they represent linear transformations between time and frequency domains. This gives E(t,f) [Cg (t, f)] = J (, ) g u + 9g -u- e) du(T) e-j2(t+ dvdr (4.15) The expectation operation will return zero unless r = 0. For the case of T = 0, the expectation is a2. This allows rewriting the equation as E[C' (t, f)] = f (26(T)(v, T)ej2(Uv-tfT)du dvdT (4.16) Moving the constant a2 outside the integral and nesting the three integrals yields E[Cg(t, f)] = U2 J [f (J eJ2udu) (v, Tr)e-j217rtdv] 6(T)e-J2rfrdT (4.17) The inner integration gives the Fourier transform of a constant which leads to a delta function in the frequency domain, 65(v). E[C (t, f)] = fc2 (v( v, rT)e-i2vtdv(T)e fdT (4.18) Both integrals contain a delta function of the integration variable. This means the integrand is sampled at the point where the argument of the delta function equals zero. Cg(t,f) = E[Cg(t, )] = 20(o,0) (4.19) 47

The expected value is the noise power (the variance of the noise signal, or a2) multiplied by the transformation kernel, ((v, r), evaluated at the origin. For the distributions of concern, 0 has a value of 1 at the origin,[7] giving the final result E [Cg(t, f)] = (4.20) This can also be evaluated using the autocorrelation expression of the transform, which begins E[C9(tj)] = E g u + g* ( - ) (t - u,) e-TJ2fdudud (4.21) Once again, the expectation returns a non-zero value only if T = 0 in which case it becomes a2. Substituting this relationship and reversing the order of integration means E[Cg(t, f)] = J J a2(r)(t - u, T)e-J2fdudr (4.22) = a2 f/ (t- u,0)du (4.23) = a2 (4.24) The last step in the two preceding derivations assume the kernel generates a distribution which satisfies the time marginal property. To satisfy the time-marginal property, (f, 0) must be zero [22], which means the kernel in the autocorrelation domain must have the property ~(t, 0) = 6(t). 4.2.2 Continuous Variance Finding the variance for the time-frequency distribution is more complex. The derivation begins with var [Cg(t, f)] = E [Cg(t, f) - E[Cg(t, f)]12] (4.25) Since the expected value of Cg(t, f) is a constant, the argument of the expectation can be expanded and simplified. var [Cg(t, f)] = E [Cg(tf)C(t, f) - Cg(t, f)Cg(t, f) - C (tf)Cg(t, f) + C*(t, f) C*(t, f)] (4.26) =E [C(t, f)C (t, f)] - E [Cg(t, f)] C(t f) -E [C(t, f)] C(t, f+ Cg(t, f) C(t, f) (4.27) = E[Cg(t, f)C;(t, )] — C(t, f) C(t,f) (4.28) = E[Cg(t, f)C(t, f)] -- o410(0 0)12 (4.29) 4:8

The first term can be expanded using the definition for the transformation C9(t( ) = J (u+ ) * (u 2)'( T) T T)ei27(vu-tfT)_ d/UddT (4.30) which means E [IC(t, f)2] = E [ J ( g 9 )g (u + li) 9 (uvt -f r - )d uddr (, ) 27r(\, v+t\ i g9 (v + -) g (v- - *(Ae27-dvddst (4.31) Reordering the integrals and moving the expectation to include only the nondeterministic elements E [\Cg (tf)\2] ~ J J [(u (, 2 U) ( v-j —9(v + ) 9v )] *0 (v, T) * (A, S)eji27(lvu-t-frT-A,+At+f s) dudvdTdvdAds (4.32) Evaluate this integral using the expectation properties of white gaussian noise. The values of g(tl) and g(t2) are independent unless tl = t2. Independence allows separating the expectation of their product into the product of their individual expectations. For the case t l 7 t2 E [g(tl)g(t2)] = E [g(ti)] E[g(t2)] (4.33) When applied to equation 4.32, if tl $ t2: t3 $ t4 then the expectation is E[g(tl)g*(t2)g*(t3)g(t4)] = E[g(t)] E[g(t2)]E [g(t3)]E [g(t4)] (4.34) = 0 0 0 0 (4.35) = 0 (4.36) since the individual distributions are zero-mean gaussian random variables. A non-zero value results unless all four times are matched in pairs, that is, only two time instances t1 and t2 (which may be the same), need to be considered E [g( tl)g* (tl)g* (tl)g(t)] E [g(tl)g*(tl)g* (t2)g(t2)] E [g(tl)g*(t2)9*(tl)g(t2)] E [g(tl)g*(t2)9* (t2)g(tl)] = E[[g(tl)l4] = E [lg(tl)l2] E [g(t2)12] == E [lg(ti) 2] E (2) = E [{g(tl)}2] E [{g*(t2)}2] (4.37) (4.38) (4.39) (4.40) 49

Equation 4.37 requires u = v with the additional condition that r and s are both zero. Expanding this for the complex value grl + jgil gives E[g(ti)l4] = E [gr, + j 4ii1 ] (4.41) = [1 + 47 + 2g419 ] (4.42) 4 4 2 2 = 3<r4 + 3'4 + 2cr-i (4.43) = 3 +34 +2 (4.44) = 2 4 (4.45) since gr and gi have the same distribution the total noise power u2 is split evenly between the real (in-phase) and imaginary (quadrature) components, making ar2 = a/2 and a2 = a2/2. Equation 4.38 corresponds to u: v and r = s = 0. Expand using grl + jgil and gr2 + g9i2 E [g(tl)l] E [g(t2)2] = E [g9l + jill]] E [gr2 + gi22] (4.46) = E [g21 + g] [r22+ 2] (4.47) = 2 2+ a2 + 22 2 (4.48) Cr 0 +UaJ7 + 2r(z (4.48) = a4 (4.49) The same result occurs for equation 4.39, used when u = v and T = s - 0. Equation 4.40 is valid when u = v and T = -s: 0. Here the values at each time instance do not get conjugated and multiplied to give the magnitude squared. E [{g(tl)}2] E [*(t2)}2] = E [(grl +jgil)2] [(r2- jgi2)2] (4.50) = E [9r - + 2igrigi] E [Ei -g - 2igr29i2] (4.51) = ( + 0) ( - + o) (4.52) = U + 0 - 2r202 (4.53) = 0 (4.54) with r = oi = a// allowing the move from equation 4.53 to equation 4.54. The three non-zero cases can be accounted for by the substitution E 9 ( + 9 (- ) g ( + ) (- )= [6(r)6(s) + 6(u - v)6(T - s)] a4 (4.55) into equation 4.32 E [Cg(t, f)2] = 50

////// [6(T)6(S) + 6(t( - ')6((T - S)] 4 O(, 7)0* (A, s)e27r(vu-lt-fr-AXv+t+ y s) dududTdvdAds (4.56) Rewrite this integral as a sum of two terms, one accounting for the case T = s = 0 and the other T = s: 0. E [ICg(t, f)2] = a4(A _+ B) (4.57) where A = JJJJJ 6(T)6(s)(G,( T)&* (A, S)ej27r(vu-v't - T-Av+t+fs)dudvdTdvdA ds (4.58) B = JJ//Jj(U V)6(T) - s)4(v, T)(A, s)ei2(vu"-t-fr-Av+At+fs)dudvddvdAds (4.59) In A, the expression 6(s)6(T) allows removing two levels of integration resulting in A = JJJ f (), 0)*(A, 0)ej27(vu-lt-Av+At) dudvdvdA (4.60) Regrouping to isolate the u and v dependencies A = ( eJ2u du) (/ e-i2rAdv) (v, 0)*(A, 0)eji2(At-t)d dA (4.61) The inner integrals act as Fourier transforms of constant arguments and can be replaced with the delta functions, 6(u) and 6(A), leaving A = / J 6(v(A)()0(,, 0) (A, 0)ji2(A-t)dvdA (4.62) which simplifies to A = ~(0, 0)~* (0, 0) = 10(0,0)12 (4.63) In the expression for B, 6(u - v) and 6(r - s) allow the substitution of v = u and s = T and removal of the v and s integrals, leaving B = / J /fc)(v,T)c*(A,T)eJ2r(-tu t+-AU t)dudvdTdA (4.64) Evaluating the integration over u results in Je-j27r(-)udu = u(,-_){(1 (4.65) = 6(u - A) (4.66) 51

When substituted into the expression for B, this allows removing the A integration. B becomes B = J / O(, T)*(V, T)eJ2(O)dvdr (4.67) = J IQ(VT)2dvdr (4.68) Taking the expression for variance in equation 4.29 and substituting equations 4.63 and 4.68 (A and B) var[Cg(t, f)] = E[lCg(t, f)12] -- E[Cg(t, f)]2 (4.69) = 4[A4 + B] - 4 10(0, 0) 2 (4.70) 41 (0, o)12 + ~,4 J / C(' T) 2dvdT -_ 4(o, 0)|2 (4.71) = 04 J/ l(v, T)12dvdT (4.72) Unfortunately, this integral cannot be evaluated for all possible kernels. A popular class of kernel functions are product kernels -that is, O(v, r) = -(vrT) -with the value of 1 at the origin (necessary for desirable marginal properties in the distribution [21]. For these kernels, the integrals are unbounded. Viewed from the autocorrelation domain, the mean square value is represented by E[lCg(t, f)2] = E [ g(u + )g*(u - )(t - u, T)e-J2rfdud * J g(v + 2)g*(v - ) -(t-v, s)e- 2sdvds] (4.73) Moving the expectation to cover only the random elements and reordering the integrals E[ICg(t, f)[2] E[C f2E g(u + )g* (u - ')g(v + )g*(v 2 2 2 2 ] ~ (t - u, r)h)(t - v, s)e-j2fr e-ji2~7fdudTdvds (4.74) As described earlier, the expectation can be represented inside this integral as two delta functions. E[lCg(t, f)\2] = f j[6(r)6(s) + 6(u - v)6(T-s)] a4 * (t - u, T)'(t - v, s)e-J2If(T+S)dudTdvds (4.75) Writing this as the sum of two integrals and applying the properties of the delta functions gives the following result E[[Cg(t, f) 2] = o4 Jl u(t - U,0)'(t - v, O)dudv + 4 f (t - u, T)(t - u,-r)dudT (4.76) 52

To simplify this, the standard kernel property 0(t, 0) = 6(t) [22] allows removal of the first integral. Using a change of variables, integration limits of infinity, and symmetry properties of standard kernels allows writing the t - u arguments as u. E[ICg(t, f)l2] -= 0 +4 JJ / '(ur)' )(a, -r)dudr (4.77) The variance, obtained from subtracting the square of equation 4.24 from equation 4.77, has the same form as found earlier when working in the ambiguity domain. var[Cg(t, )] = '4 J / (U, T))(U,- -r)dudT (4.78) As seen before, these integrals with infinite limits are unbounded. 4.3 Discretized Transform Results form the continuous transformation shows an ill-behaved expression for variance because of the infinite power required for a truly white noise signal. For discrete signals, a white signal is physically realizable because discrete signals are inherently band-limited. Returning to the continuous definition of the time-frequency distribution Cg (t, f) 9 u+ g ) 9 (u- ) ~(v, T)eJ2 (u f)dudvdT (4.79) To discretize this, nAT will represent time and pAf the frequency where AT is the sampling interval and,Af the corresponding frequency resolution. Cg(nAT, pAf) = EEE g mAT+ 2 )* (mT - (qf, lAT) I q m 2 2 lj 27r(qA fmAT-qAfnAT-pAflAT) X f AT (4.80) The sampling interval, AT, can be set to 1 second without loss of generality. Since the frequency resolution, Af, is given by the reciprocal of the window length, NzAT, the summations can be re-written as Cg(n,p) = g (m + * (m- ( ) N m-n (4.81) 4.3.1 Discretized Mean To find the expected value for the discrete distribution, the non-deterministic g(n) terms need to be kept inside the E[.] operator E[Cg(n,p)] = 53

IE g [ 2 ( 2 +1 ( - I e, -(qm-qn-pl) E (m + (m - ) (4.82) I q m The expectation for the gg* term is a251, so the equation becomes 1 2 q 2,1,,q l (qm-qnpp) E [C, (n, p)] = - 5 5 a2 S ( ') e) I::(qm- q -p ) (4.83) I q mn Resolving the I summation E[Cg(n,p)] = -E o) - (4.84) q m 0\2 2'*-q 7rmq = ( e-j Eej N (4.85) q m The second summation, which involves the complex exponential with respect to m will evaluate to zero unless q = 0, in which case the summation becomes N. The mean value may be re-written as E [Cg(n, p)] = o2 bqq, ( o 0) e N (4.86) q N Performing the last summation and applying the property 0(0, 0) = 1 to obtain E [Cg(n,p)] = a2 (4.87) which matches the expected value found for the continuous case. 4.3.2 Discretized Variance Turning to the variance calculation, we begin with the discrete Cohen's class transformation from equation 4.81. As previously explained, the new term needed is E [\C,(np)\2]= E[ SS {g (m + L) g* (m -A) t (~,) <(qm-qn-pI) } E { (c g* (c - ( d) ej2 (bcbn ) (4.88) d b C J Rearranging the summation symbols and collecting random variables under the expectation operator E [lCg(n p)2 = N2 EEEl~d E [9 2 ) 2} 2)\ ( q ( d b c (q I 2 m ~1 0 -,I,d el'jN (qm-bc)+(b-q)n+p(d-1l (4-89) (N.89 54

As with the continuous variance, the expectation operation has three cases with non-zero results. These are accounted for by substituting E [g (m + g (m - * (c + d (c - = [66d + (m)(d)] (4.90) The variance becomes E[Cg (n,p)12] = 2 4 [d + (m-c)(1-d))]! q m d b c. ( ) (,d) eI N (qm-bc+(b-q)n+p(d-l)) (4.91) As with the continuous case, breaking this expression into two pieces and incorporating each 6 function relationship results in a4 E [ICg(n,p)12] -[A B] (4.92) E^^ A +B} (4.92) where A = ESE 6d(q l0*(bd) I q m d b c j 2 (qm-bc+(b-q)n+p(d-l)) (4. N T(4.93) B = SESEESEE(m-c)(-d) (Ql) (b ( d) I q m d b c j 2N (qm-bc+(b-q)n+p(d-l)) (4.94) N.gJ^(^ +(&(?),+(.-0 (4.94) The expression for A simplifies when summed over I and d because of the properties of the 6 functions ( N 0 ~, 0 e Nj (qm-bc+(b-q)n) - [= N2'b ] (, 0 ( 0 ejb- (qm+(b-q)n) (4.96) q m b C The bracketed summation has two possible values. If b $: 0, the summation goes to zero. When b = 0, it equals N. In that case - qm, 0*0 A=[v N [ (q ( N 0) * (, O) e-J Nqn (4.97) q L m The summation over m will also become N when q = 0 and is zero for all other q. This simplifies the expression for A to A = VN2(0, 0)4* (0,0) = NV21(0,0)12 (4.98) 55

Applying the 6 functions to simplify B B E — q) eJ I (qm-bm+bn-qn) (4.99) I q m b e EE E [ej(qb)m ( I) ( e) (bn (4.100) I q b m The bracketed summation is non-zero only when q = b, when it evaluates to N. This leaves the result B=N E N ',) (4.101) I q Substituting for A and B in equation 4.92 and using the result in the expression for variance var {C(n,p)} E [Cg(n,p)2] - Cg (n,p)2 (4.102) 4 = [A + B]- (a2~(0 0))2 (4.103) = E (,i) ] (4.104) N= aE (,I) (4.105) These results parallel the continuous expressions for A and B obtained in the preceding section. In contrast to the infinite integral limits used for the continuous case, the discrete variance calculation involves finite summations. This allows numerically evaluating the variance. 4.4 Sampled Transform A parallel approach considers evenly spaced gaussian noise samples, g(n), where g(n) = g(nAT). Transform the g(n) samples using the discrete counterparts to the Fourier transforms N-1 np[x(n)] = E x(n)e-j pn (4.106) n=O = X(p) (4.107) [X(p)] = X(p)e pn (4.108) (n (4.109) = x~n) (4.109) 56

These allows expressing the transformation of the time signal to the sampled time-frequency domain as Cg(np) = (Z {Fm(-q) [g (m + ) (m - g)] (q, 1)}) (4.110) The Fourier transforms in equations 4.106 and 4.108 involve a periodicity, either real or implied, over a period of N.At in the time domain which provides a periodic spectra with frequency resolution of Af = 1/(NAt), with repeats every 1/At hertz. To express the transform using the autocorrelation kernel, convolve the kernel with the autocorrelation and take the Fourier transform with respect to the lag variable, 1. C9(np) = Tip [Rg(n, 1) O V)l (n, 1)] (4.111) where Rg(n,l 1) in the local autocorrelation and i(n, 1) is the autocorrelation kernel. The Fourier transform of the autocorrelation kernel gives the ambiguity domain kernel, O(q, 1) Ynq[l(n, 1)] = l(q, 1) (4.112) 4.4.1 Data Limitations There is an important limitation to keep in mind when using this sampled data approach. The continuous definition for Cohen's class implies time shifts of ~1/2 when determining the discrete autocorrelation function Rg(n, 1). Odd values of I this requires values at points between the samples. Two techniques for dealing with non-integer sample points yield quite different results [21]. One approach treats these intermediate values as zero. Unfortunately this traditional approach discards information provided by odd lag values and halves the maximum unaliased frequency. An alternative approach uses separate methods for odd and even lag values. By slightly modifying the definition for the kernel and the autocorrelation function, it is possible to avoid undersampling the signal. This alias-free formulation allows an unambiguous frequency range of half the sampling frequency in comparison with the traditional formulation, which has a range of one-fourth the sampling frequency. The derivations below follow the alias-free formulation in both ambiguity and autocorrelation domains. The traditional methodology differs by eliminating terms with odd lag values. 57

4.4.2 Sampled Mean Calculating the expected value for the time-frequency distribution computed from sampled data requires evaluating E[Cg9(np)] =E p (F n {Fm(-q) [g (m + ) (m - (, })] (4.113) When expanded, this becomes E[Cg(np)] =E - E g (m + g* (m - ) )( I)- P e 1=0 N q=0 m=O / (4.114) XMoving the deterministic portions of this expression outside the expectation operator P1 NP IV2 Nr 2 7r E[Cg(n,p)] = - E E E gE [ + -) m -e- e- e 1=0 q=0 m=0 (4.115) The expectation will be non-zero only when I = 0 and the arguments for g and g* are equivalent. When I $ 0, the independence of the samples drives the expectation to zero. Accounting for this by replacing the E[gg*] term with a261, makes the expected value P N N E[Cg(np)] = Z 2 (q )e-Jl e-qm (4.116) l=0 q=0 m=O Using the delta, function to resolve the summation over I (72 N N E[Cg(n,p)] = E E )(q, )qm (4.117) q=O m=0 In the summation over m make the substitution N ej =m N6q (4.118) m=0 The result is or2 N E[C (n,p)] = ~ N6q(q, O)e-j Pn (4.119) q=0 = C2)(0, 0) (4.120) To satisfy the marginal properties, the kernel equals zero at the origin [18], with the final result E[Cg(n,p)] =2 (4.121) 58

4.4.3 Sampled Variance Finding the variance of a distribution requires both the mean value and the mean square vue. The mean square value calculated using the ambiguity domain form of the timefrequency distribution given in equation 4.110 is E [lCg(n p) 2] = E {lp (F-ln {m-) [g (m+ 2 - (m -)] (q, 1}) = E{p ( — qn {m(-q) [g (rm+ 2- - (m-2)] 0(q,)}). -dp (1bn {Fcb [9* (C + - g (c- 2)] *(b,d)})} (4.123) Expanding the Fourier and inverse Fourier transforms E [ Cg(n,p)\2] = E Fp g m + )g ( (+m-2 g- -(q, )e e i P:=0 N =0 m=0 d\ ( \ 2d d0 b.cE g (c + ) 9 (c- - N bnej b (4.124) d=2 2b= c=O Expanoving constants outside the summations and deterministic portions outside the expectation E[ Cg(np) 2] 2/ r 2/ d i= N, EEE9 (m 29 (c+^2)r c (4.125) =0 q=0 m=O d=0 =0 c=0 g(q c ) c (b, d)e ep p(d- )ei N (e-N)nejbc -j"m (4.125) To avoid an expectation of zero, the four samples under consideration must be matched in conjugate pairs. A non-zero value occurs when m = c and 1 = d, or when l = d = 0. m =c ld = 0 24 m=c l =dm O a4g m c = = d= 0 a4 [ (m + ) * (m - ) * (c + (c -) = cr46mc^ d + 41d (4.126) 2[!, 2/ 2/ 2 - 59

The mean square value becomes E[ Cg (n,p) 12] = 4 P N NA P N N 2 (EE E E5 E 5 E m-cS-d~(q, 1)* (b, d)ej p(d1=0 q=O m==O d=O b=O c=0 P N N P N N + E Ei ~ 616bd(q, 1)* (b, d)eJ p (d-)e 2 (q1=0 q=0 m=0 d-=O b=O c=0 -1) e j' (q-b)nej - bc 2j qm el N el N C - N -b)n ej N bce -jNqm 2re NI (4.127) Resolve the summation over the delta functions E [|C(n,p)|2] = 4 P N N N N2 5E E E (q, l)4*(b, )eN (q-b)nej (b-q)m 1=0 q=0 m=0 b=O N N N N + E E (q, 0)O* (b, )eJ N (q-b)n e- J (4.128) q=O m=O b=O c=O In the first set of summations, the only m dependence occurs in the complex exponential term. This term becomes N when q = b and is zero for all other cases. E[ C,(n,p) 2] 4 P N N2 NEE ](q ) 12 1=0 q=O N N N N + E E E 0(q, O)* (b, O)e N (- be-j (4.129) q=0 m=0 b=0 c=0 For the second set of summations, c and m dependencies occur only in the exponential terms and have a non-zero result only when q = 0 and b = 0. E [cn p) 2] =c (N This reduces to E[ICg(n,p)2 ]= P N \ E E [o(q, 1)I2 + N2 1(0, 0)t2) 1=0 q=0 4t P N ) N =0 l(q, )2 + =4 1=0 q=O (4.130) (4.131) The variance is the mean squared value minus the square of the mean value found in equation 4.121 var [Ig(t, f) 2].4 P AN = 0 q ( )l2 l=O q=O (4.132) 60

Determining variance in the autocorrelation domain requires evaluating var [ (n,p) 2] = E [|(, (n,p)2] - (E [C(n,p)])2 (4.133) Evaluate this expression by separating the odd and even lag values, 1, to allow an aliasfree formulation for Cg(n,p)[21]. 9C(n.p) n g (1 2P = n1=0 even I P-l 0 / 1 27 ------ PI gy(k+ - -k++ (k+ - '(n-k, )e ko 2 2/ 2 2/ 1=0 odd I (4.134) where /'(n, 1) in the modified ambiguity free-kernel for odd values of 1. Mlodifying the lag index covers positive and negative values, ranging from -P/2 + 1 to P/2 - 1. The local autocorrelation for a lag of -P/2 must be zero to satisfy symmetry properties of the transform. even I -pP-1 + Z E ^g (k + )g - +2 +'(n+-_ - 2 +1 1=:-P+-1 k=-oo odd I (4.135) This expectation becomes P 2 00\ j 2jr PI E[\Cg(nr,p)2] = E ~ k= g k g* (k- - ((n-k,l)e-J P 1=-P+1 k=-oo / even I 2 -J2-1 + E ( 2 + 2 9 (k + -- ) 1'(n-,)e-P 1=-P+1 kc=-oo 2 2 / odd I (4.136) Breaking up this incredibly cumbersome mess into only slightly more workable nightmares E[lCg(n,p)22]= E[A+ B + C +D] (4.137) 61

where = E[A] + E[B] + E[C] + E[D] k=-c m= —g 2 2 (4.138) P-1 2 A += even I — 1 2 d= P -+ 1 even d g9 (r+ - + g* (m- - k, l)i(n - m, d)e-J P p(l+d) 2 2 (4.139) 2:L.E B -- 1 even I '(m+ 2 00 00 I 1 E E 9g k+ )9*(k-) d=-P+1 k —m-oo m= — oo 2 2 odd d - + ) * m + - ) V;(n - k, 1)'(n - m, d)e-i p(l+d) (4.140) 2 / 2.-2 2 00 00 / 1I\ / 1 c = 1 E g(k++l 9g*' k++ =-:p+1 d=-P+1 k=-owm= —oc 2\ 2 odd I even d g9 (m+ g-) * (m - '(n - k, l)(n - m, d)e-pP(l+d) 2 2 -, 1 (4.141) pI -1. P-1 D = 5 E ( 1 + )g (k+ - + ) l=2-=+1 d=2P+1 k=-oom=-o 2 odd l odd d 9 (m + + 9g* (m+ - )'(n- k, l)'(n- m, d)e-i'p(l+d) (4.142) 2 2 2+ 2, + - The expected value of A simplifies using the relationship E (k+ g* (k - 9* g (m + d) (m - = a4 (616d + 6k-m6l-d) L 2 2/ 2 2/ \ / \ (4.143) so E[A] P-1 2 == _a4 E 1= P + even I - -1 d= P+ even d 00 00 E: E~ [t61d + 6k-mbl-d] k=-oo m=-0o * (n-k, I)~((n - m, d)e-j p (l+d) 00P-1 00 + S 1= P +21 k=-oo even I (4.144) (4.145) 62

Following the same basic steps used to find E[A] will lead to the expected value of D. Since I and d must be odd in the expression for D, they cannot be zero. This means there will not be the 616d term as in equation 4.144. The expected value of D simplifies to P 1 2 OC E[D] = 4 " (n - k,l(n- k, -1) (4.146) =- P +1 k=-oo odd I The expected values of B and C are zero. In equation 4.140 the g(m +.5 + d/2) and g* (m +.5 - d/2) terms can combine as a complex conjugate pair only if d = 0, which is not possible since the summation limits d to odd values. The g(m+-.5+d/2) and g*(m+.5-d/2) cannot match up with the g(k + l/2) and g*(k -l/2) terms. The spans, \dl and Ill, cannot be equal since I must be even and d must be odd. The same limits on the summation variables force the value of C to zero in equation 4.141. The variance of the distribution now becomes var[C,(n,p)] = E[lCg(n,p)2] - IE[Cg(np)]12 (4.147) = E[A] + E[D] - (a2)2 (4.148) 00 00= - = L o4 E m, O(n - k, O))(n-m,0) k=-oo TM=-oO P 1 21 0 + E E: (n — k, l)>(n-ki-l) 1=-P+1 k=-oo even l P2 00 + E: a 4'(n -- k, l)'(n - k, -1) (4.149) 1=-P+1 Ak=-oo odd I Since all the k based summations have infinite limits, the specific value of n does not matter, allowing n = 0. var[Cg (n,p)] P 'X) 00 2 00 = (4 -1 + E E (-,0)~ (-k,1)) +-k,-l) k=-oo m=-0 1=-P +1 k=-oo00 even I P 2 oo + E V '(-k, )'(-k,-1 (4.150) 1=-P+t kc=-oo odd I 63

If we assume the kernel provides a time-frequency distribution with the correct time marginal, the kernel at zero lag can be treated as a delta function, 4'(n, 0) = 6n, forcing the first summation to a value of 1 var[C (n, p)] = 4 P -1 E E (-kI)V(-k,-1) + I —=P+1 k= —o even. I P1 2-1 00 E E 0'(-k, 1)'(-k, -1) 1=-2P+1 k=-oo odd I (4.151) 4.5 Specific Distributions Finding actual values for the mean and variance requires a specific transformation kernel. This section numerically examines the mean and standard deviation obtained when using two members of Cohen's class, the Wigner distribution and the binomial distribution. The analytic results derived include both the alias-free and the traditional formulations of these distributions. 4.5.1 Discrete Wigner kernel The discrete Wigner kernel, or perhaps more properly the quasi-Wigner kernel [30], has the alias-free form (n,) = 6n 1'(n, 1) = +,5n-1 which makes the alias-free variance which makes the alias-free variance (4.152) (4.153) var[Cg(n, p)] — 1 = 4 E 1== -P+1 z even I P = 4 E 1=:- P+1 even I = 4 (3_-1) 4 00 Z 6-k8-k + k=-oo P 1 2 ~ /I 1 2 E E k + -k1 (4.154) 1=-P+1 k=-o 2 odd odd I P-1 2 odd odd I 1 2 (4.155) (4.156) 64

As mentioned in section 4.4.1, the traditional formulation sets the function 0'(n, I) to zero. This forces all the rows in the ambiguity domain associated with odd lags to zero. When applied to the variance, this gives PI — 1 2 var[Cg(n,p)] = a4 E 6 -k (4.157) l=- P+1 k=-oo even I P1 2 7-i = (4 1 (4.158) 1= -P+1 even I 4 (P -1) (4.159) 4.5.2 Binomial kernel Kernels which generate a reduced interference distribution (RID), such as the binomial and the spectrogram, provide lower values for variance in the time-frequency domain than the Wigner. The RID kernels replace the uniform ambiguity domain weighting of the Wigner kernal with a tapered window before transforming the values into the time-frequency domain. This provides a measure of smoothing to the final distribution at the cost of some frequency and time resolution. Determining the variance for the binomial kernel is more involved since unlike the Wigner, the kernel does not have a simple delta function form in the autocorrelation domain. Exploiting the symmetry in the kernel allows expressing the variance as P 1P-1 2 1 00 2 00 var[Cg(n,p)] = 4 2 [i(k,)]2 + E E [I (k)] (4.160) 1=-P+1 k=-oo 1=-+1 l k=-oo even I odd I Each row of the binomial kernel is given by the binomial coefficients, scaled so the sum across the row is unity. Both 4O and 4' have the same form, allowing the two summations to be combined. P-1 var[Cg(n,p)] =4 a 5 [4(k,1)]2 (4.161) I- -P+1 k= —oo — 2 Row I of the kernel has I + 1 non-zero elements. These elements have the unscaled values based on the binomial coefficients (1 ) =m( )! (4.162) 7 -m!(/ m,)' 65

with m = 0... 1. Each element in row I must be scaled by 2-1 to force the sum of the values across the row to equal 1. Replacing the k summation with a summation over m gives -P 1 l 1 2 2 i 1 var[C9(nrp) ]= 4 Z T( 1 )' (4.163) /= ---iTh m=0r te a -e f This is symmetric about 1 = 0 reducing the variance for the alias-free form to var[Cg(n,p)] - 4 1 + 2 E 21 ]} /=1 m:-O ' ' (4.164) In the traditional form, the even lag rows have the same values, but the odd lag rows contain only zeros.. Since summing across the zero rows is not necessary, the variance can be re-written as var[C,(np)] = 41 +2 Z 221 b (l ) (4.165) /2:= 2 even I Figure 4.1 shows how the variance of the Wigner and binomial distributions vary as the Binomial Variance 200 150 0.Q sO > c4 0 100 d) 0 0 Cd.P4; IL4 col ItA/ Traditional 0 50 100 150 200 250 300 Transform Length CL- -- I - I -I I 1 f I. I 0 50 100 150 200 250 300 Transform Length Figure 4.1: Noise Variance vs Transform Length length of the transform increases. 4.6 Numerical Results The numerical values were derived from a vector containing 4096 independent, normally distributed point with mean 0 and variance 1. The time-frequency distribution was found for 66

this vector, using traditional and alias-free formulations with varying frequency resolution. The values for the time-frequency distribution filled a matrix Cg(n,p). The trial mean was given by mean = E [Cg] = N Cg(n,p) (4.166) n p where N is the number of time samples, in this case 4096, and P is the transform length and the number of frequency bins which ranged from 4 to 256. The variance was estimated by squaring each element in the time-frequency domain, finding the average of these squared values, and subtracting the square of the mean value. This can be expressed as 1 var - Cg(n, p)2 - E [Cg]2 (4.167) n p Tables 4.1 and 4.2 show how well the analytically obtained values for the distribution Traditional Alias-Free Window Mean Variance Variance Length ideal trial ideal trial ideal trial 4 1 0.9978 1 1.0461 2 1.9973 8 1 0.9978 3 3.0689 5 7.0689 16 1 0.9978 7 7.0730 11 11.0606 32 1 0.9978 15 14.8377 23 22.7618 64 1 0.9978 31 30.7465 47 46.5828 128 1 0.9978 63 62.1522 95 93.6993 256 1 0.9978 127 124.6105 191 187.3551 Table 4.1: Wigner Statistics statistics match the values obtained through a computer trial. These tables show excellent agreement between the theoretical and experimental statistics, with only two trial variance values differing from the expected analytical values by more than four percent. These results show a growth in variance similar to that displayed when computing a periodogram [31]. A more complete understanding would be provided by considering the case of a desired signal, s(t), added to the noise considered in this chapter. This would allow evaluation the performance of Cohen's class transformation as a signal estimatorman in terms of bias an asymptotic performance. 67

Window Length Mean ideal trial Traditional Variance ideal trial Alias-Free Variance ideal trial I, 4 8 16 32 64 128 256 1 1 1 1 1 1 1 0.9978 0.9978 0.9978 0.9978 0.9978 0.9978 0.9978 1.0000 1.7500 2.7480 4.1145 6.0152 8.6812 12.4538 1.0461 1.8054 2.8194 4.1613 6.0542 8.7077 12.3190 2.0000 3.3750 5.2842 7.9568 11.7164 17.0188 24.5074 1.9973 3.4049 5.3348 7.9594 11.6903 16.9346 24.1765 Table 4.2: Binomial Statistics 4.7 Summary Noise driving a Cohen's class transformation provides a time-frequency distribution with characteristics based on the kernel used and the length of the transform. Comparing the quasi-Wigner and the binomial kernels demonstrated a significant edge in performance for the binomial as the frequency resolution (and the transform length) increased. The quasi-Wigner displayed slightly worse than linear growth in the distribution variance as the transform length increased, 0(N), while the binomial's variance was slightly greater than O(VN). 68

CHAPTER 5 DYNAMIC TARGETS Studies of long wires and pairs of orbiting spheres demonstrate the usefulness of timefrequency distributions for the analysis of radar signals. They served as simplified models for moving blades of an engine propeller or a helicopter rotor. In general, the backscatter signal was difficult to interpret in either the time or the frequency domains. Applying the binomial distribution, a discrete time-frequency distribution, allows clearly associating each sphere with its corresponding doppler return. The binomial distribution provided a detailed view of the target dynamics, opening the way for target classification and identification. The structures and details available in the time-frequency domain were not readily exploited in the original signal representation. 5.1 Introduction In addition to looking at dispersive targets, time-frequency techniques are also directly applicable to dynamic targets. With the target in motion, the scattered signal will have non-stationary characteristics. Focusing on these changing characteristics allows estimation of parameters relating to the target's geometry and insight into the target dynamics. Radar systems combine elements of electromagnetic theory with signal processing to estimate target parameters. A rudimentary system may only try to determine whether or not a target is present in a region. More elaborate systems provide estimates of the target's location and velocity and direction. Some systems extend the parameter estimation to provide target identification [36]. Simple radar analysis can determine position and radial velocity. The target location is resolved in polar coordinates by using time delay, AT, in the returned signal to measure range and the pointing direction of a narrow-beam antenna to provide the necessary angular 69

coordinates (azimuth and elevation). Transforming the returned signal into the frequency domain gives the doppler shift imposed on the signal by the target. Knowing this shift and the original carrier frequency gives the radial velocity of the target (how quickly it is moving directly towards or away from the radar). While these techniques work reasonably well for a simple return, multiple targets, or multiple scatterers on a single target can obscure the analysis. Interaction terms make it difficult to determine the contribution from individual scatterers or scattering modes. Advanced methods like time-frequency analysis can isolate components from such multiple returns. 5.2 Background The basic motion selected for this work was rotation about a point. The targets consisted of wires and spheres. The objective for the time-frequency analysis was to determine instantaneous target velocities based on the doppler shift of the return signal. For a carrier frequency of J, the doppler shift for a point target moving with a radial velocity of v is fd 2vf (5.1) using c as the velocity of propagation. Figure 5.1 shows the basic measurement system modeled to measure these doppler shifts. The target was illuminated with a 1 GHz plane wave, then the scattered signal from the target was converted down to the baseband in-phase, I, and quadrature, Q, channels. For an object following a circular path, the radial velocity will follow a sinusoid with the maximum directly proportional to the rotational rate, frot, and the distance from the center of rotation, r. The maximum radial velocity is 27rfrotr, so the maximum doppler shift is 47r frot r fc fdmax = f f (5.2) The rotational rate and distance were selected in conjunction with the sampling frequency to avoid any ambiguities in the instantaneous doppler frequency. 5.3 Wire Targets Thin, straight pieces of wire provided a relatively simple yet potentially important target for analysis. The insights provided from a single wire carry over into multiple wire targets. 70

I Q Figure 5.1: Dynamic Measurement System Model Rotating wires serve as a simplified representation of rotating structures such as airplane propellers and helicopter rotors. 5.3.1 Single rotating wire The first set of dynamic targets examined were long, perfectly conducting wires, modeled using the Numerical Electromagnetics Code (NEC) [5]. These simple targets illustrate how time-frequency distributions can be applied -to a radar scattering problem. Since the wires used were electrically long, the backscatter pattern has several expected peaks. In addition to a sharp peak when the wire is broadside to the incident field, two other broad peaks result due to traveling waves set up on the wire. These peaks appear at [24] Op = 49.35 A/i (5.3) For the wires used, 2.5 and 5 wavelengths in length, these scattering wave peaks occur at 22.1 and 31.2 degrees. These long wires demonstrated a second scattering property. The primary scattering center was located at the tip of the wire furthest from the radar. As the wire rotates, the primary doppler contributor will follow the movement of the outer tip of the wire. 71

Three different single wire cases were tested. In two of the trials, the wire was rotated around one end. In the third trial, the rotational center was placed 2.5 wavelengths beyond the end of the wire. This allowed the end of the 2.5 wavelength wire to reach the same radial velocity, and therefore generate the same adoppler shift, as the 5 wavelength wire. Figure 5.2 illustrates these cases. 4 i Y 4Y A' A 5 2.5 a 5 2.5 L y N X X x Figure 5.2: Single Wire Geometries The corresponding backscatter wires are shown in figure 5.3 2.5 Wavelengths patterns for the 2.5 wavelength and the 5 wavelength 5 Wavelengths 180 270 270 Figure 5.3: Single Wire Backscatter Patterns 72

Figure 5.4 shows the magnitude of the time response for these three cases. A significant 2.5 Wavelength 150 200 2.5 Wavelengft with Offset 0 ____ 35 150 200 5 Wa~vlPength 300 _____ 350 ' 0.2 h 150 200 250 300 350 Tim2.5 Wavelengtdex with Oset 'E 0.2 -150 200 250a~et 300 350 0.6 - 20.4 -150 200 250 300 350 Time Index Figure 5.4: Single Wire Waveforms feature is how both 2.5 wavelength targets generate the same received signal power as a function of time. The difference in rotational centers only appears in the phase of the received signal. Taking the Fourier transform of the signals provides the power spectral densities shown in figure 5.5. Because the scattering centers on the wire have constantly changing velocities, the doppler of the scattered signal has a smeared nature, covering a range of values. The differences between the 2.5 and 5 wavelength wires can be seen in the maximum doppler achieved, with the 5 wavelength wire having twice the frequency spread. Another interesting feature in figure 5.5 is the difference between the 2.5 wavelength wires. As expected, the wire with the offset rotational center reaches the same maximum frequency as the longer wire. Both ends of the shorter wire are in motion, so neither tip stays stationary at the rotational center to make a significant contributions at the lower doppler frequencies. The binomial distribution for the single wires is shown in figure 5.6. In the joint timefrequency domain several characteristics of the signal become obvious. The rotational aspect of the wire dynamics show up in the sinusoidal envelope formed in the binomial distribution. The maximum doppler frequency occurs when the wire is perpendicular to the direction of 73

2.5 Wavelength PSD 150 -.F 100 - 50 - I I I I - T -- - I -1 I d A — -250 -200 -150 -100 -50 0 50 100 150 200 250 2.5 Wavelength with Offset PSID 1502 -(." I I,AAAAA.- - I 61 -250 -200 -150 -100 -50 0 50 100 150 200 250 5 Wavelength PSD 200 a) 150 1000 50 7- H I -1, I. I An- I -- -AA- I AAA I. I2,, L 0 L-L - ----- ---- -250 -200 -150 -100 -50 0 Frequency Index 50 100 150 200 250 Figure 5.5: Single Wire Power Spectral Density 2.5 Wavelength Wire la UC-) CL a' U 150 200 2 'v~50~i Wr 300 350 Figure 05.6: Single Wire Binomial Distributions 74

propagation, so it corresponds to one of the peak amplitudes in the return. The broadside orientation also allows each point along the wire to contribute energy to the return, creating a broadband frequency return as all dopplers between 0 and fmax contribute. The largest returns in the backscatter signal are due to the traveling wave returns. These show up four times in each cycle. The extreme size of this component obscures the details of the time-frequency distribution around these points in addition to creating cross terms. Another characteristic of the scattered wave shown by the time-frequency distribution is how the scattering favors the end of the wire furthest from the radar. On these plots, the sinusoidal pattern generated is asymmetric through a cycle. The trace is darker from the negative maximum doppler point through the positive maximum doppler one half cycle later. This corresponds to the time the outer tip is further from the radar than the inner point. The other half cycle has a lighter component and corresponds to the time when the fast moving outer point is closer to the radar. For the wires rotated about their end, the zero frequency line is stronger during this half cycle. For the offset wire, there are no significant contributions from the spectrum around zero hertz. However, there is still a jump in emphasis following the rear point. In this case, neither end of the wire is stationary, so the zero doppler contribution is very small. The offset wire has an interesting overall pattern. It has the same basic features as the other 2.5 wavelength wire, but a null band appears at the low doppler frequencies. At the maximum doppler points, the offset wire reaches the same peak frequency as the longer wire, as expected, but it cannot support modes with zero doppler shifts since all the points on the offset wire are traveling fast enough to generate a doppler shift of half the peak value. 5.3.2 Multiple rotating wires The scattering off a 2.5 wavelength wire was extended to consider multiple wires. In this set of trials, the 2.5 wavelength section served as the initial target. Then a second, identical, wire was attached at the rotation point, positioned 120 degrees away from the first wire. A third wire attached to the rotation point and positioned 240 degrees away from the first wire completed the symmetric three-wire set. The geometries involved are shown in figure 5.7. The backscatter for the single wire was already shown in figure 5.3. The corresponding backscatter patterns for the double and triple wire targets are given in figure 5.8. The targets were rotated at while illuminated with a planar incident field with the electric field polarized parallel to the plane of the wires. Figure 5.9 shows the magnitude of 75

y y y 2.5 >4 >4 x x x Figure 5.7: Multiple Wire Geometries Double Wire Target Triple Wire Target 0 0 270 270 Figure 5.8: Multiple Wire Backscatter Patterns 76

Single Wire 0.42 0 " 200 210 220 230 240 Douge0Wire260 270 280 290 300 0.6 -'- 1Il '., 0.4 0 -- 200 210 220 230 240 250 260 270 280 290 300 Figure 5.9: Multiple Wire Waveforms the received backscatter signal. The corresponding power spectral density, obtained using the Fourier transform, is shown in figure 5.10. As with the single wire trials, the dominant returns are due to traveling waves. With careful processing of the scattered signal, it is possible to discern portions of two sinusoidal patterns in the two-wire target. This is clearest in the vicinity of the maximum doppler shifts. Figure 5.11 shows how increasing the complexity of the target return rapidly fills the time-frequency response with components and cross terms that make interpretation difficult. The rapid changes in the power of the backscatter signal contribute to this difficulty. There is one additional problem encountered when trying to interpret the three-wire scatterer compared to the single and double wire cases. The three-wire target is the only one in the set that exhibits complete rotational symmetry - each third of a rotation is indistinguishable from the other two thirds. This makes it extremely difficult to make analysis without a priori knowledge of the target makeup. 5.4 Spherical Targets This study used two spheres traveling in simple circular orbits. Although the spheres were illuminated by a planar, continuous-wave signal, the return from the two moving 77

Single Wire PSD 200J.. Q) C 100 o.I.... L.ULL JILX AA A A,,, -250 -200 -150 -100 -5Double IirePS 0 100 150 200 250 400 - -- E' 200 -01I -250 -200 -150 -100 -59riple VYire Ps130 100 150 200 250 400 3 1200 --250 -200 -150 -100 -50 0 50 100 150 200 250 Frequency Index Figure 5.10: Multiple Wire Spectral Densities Single Wire TFD -200 c O' ---,I- 200 200 210 220 230 2, ouble2 6r T o 270 280 290 300, 200 __ 200 210 220 230 240 250 260 270 280 290 300 5.1iplW le irnoe Triuin;:-200........... 200 210 220 230 240 250 260 270 280 290 300 78

z X radar 4 incident scattered 15 cm Figure 5.12: Target Dimensions scatterers undergoes rapid changes in both amplitude and phase. Using the scattered signal allows determining the number of scatterers, their rotational frequency, and position of the spheres with respect to the center of rotation. By using time-frequency techniques, the typical estimates of these values had errors less than five percent. 5.4.1 Electromagnetics problem The basic problem examined involved two metallic spheres as sketched in Figure 5.12. The target consisted of two perfectly conducting spheres with their centers separated by one meter. Each sphere had a radius of 15 centimeters. When illuminated by a 1 GHz incident plane wave, the multiple spheres created several different modes of reflected electromagnetic energy. These modes included specular reflections, creeping waves, and multiple bounce reflections. The original data for this experiment was generated by running a body of revolution simulation using two spheres as targets. The returned signal values were based on monostatic, far-field calculations. The simulation covered a frequency range of 1.0 to 7.55 GHz and considered observation angles over a span of 90 degrees. Viewed in polar coordinates, the spheres were placed with one center at (0.0, 0.0, 0.5) and the other center at (0.0, 0.0, -0.5). The angle 0 ranged from 0 to 90 degrees, sampled every half degree. The symmetry of the two spheres allowed this data set to be extended to cover view angles from 0 to 360 degrees 79

1 meter D C B A 15 cm E Figure 5.13: Centers of Rotation with respect to the positive z-axis, measured in the x-z plane. The carrier frequency of 1 GHz insured the doppler shifts imposed on the signal by the target rotation would be small enough to avoid undersampling in the frequency domain. 5.4.2 Processing and Estimation To create a dynamic problem, numerical simulation placed the spheres in circular orbits around five different rotational centers, referred to as cases A through E. Figure 5.13 shows the positions chosen. Each of the five cases provided different return signal characteristics. Each trial used two seconds of data sampled at a rate of data sampled at a rate of 512 hertz. After applying a 512 -point transformation, this provided resolution of 1 hertz in the frequency domain. The rotational rate was chosen as 3.5 hertz so the maximum expected doppler shift would not give an undersampled signal. Calculations for each sampling instant determined the view angle, then interpolated between the values in the simulated data set. In all cases, each sphere maintained a constant distance from the rotational center. Both spheres traveled at the same angular rate. For each of the five cases, estimates were generated of four parameters: the period (To), the distance from the rotational center to the first sphere (r1i), the distance from the rotational center to the second sphere (r2), and the distance between the spheres (ro). The period was determined using the autocorrelation, R(T), of the time domain samples. The number of autocorrelation peaks, N, and the time shifts for the first and last peaks, Tr 80

and TN, provided the estimate of the period TO 'TN - 71 /- A) TN - T1 To- 1- (5.4) The radii of the sphere orbits were based on the measured peak doppler shift, fdmaxz the estimated rotational period, To, and the known illumination frequency, fi. Using equation 3.2, the only unknown value after measuring the peak doppler shift was the radius, r. The estimate for the radius r1 or r2 becomes rl,2= 4fC (5.5) The final value, ro, was found using r1, r2, and the phase shift, 0, between their peak doppler points. The estimate of 0 depends on T1 and T2, the times each sphere generates their maximum doppler shift T2 - T1 0 = 360 2 (5.6) To Using this phase shift in the law of cosines gives the sphere separation, ro as ro= r + T 2 - 2rlr2 cos 9 (5.7) When generating these estimates, the period was calculated autonomously by a computer. The maximum doppler shift frequencies and times values were manually extracted from plots of the time-frequency distribution to allow separating the doppler components associated wit-: each sphere. Since the spheres traveled in a circular path, their radial velocity with respect to a stationary observer varied sinusoidally. The velocity passed through zero at the spheres' closest and furthest position with respect to the observer, that is, when the sphere, the rotational center, and the observer were in line. The maximum velocity points were offset by ~90 degrees from the zero velocity points. The doppler shift imposed on the scattered signal is directly proportional to the sphere velocity, and is given by the relationship 2v fd = - (5.8) Since the two targets cannot occupy the same position of the same orbit simultaneously, each will have a different doppler frequency over the course of a cycle. In the cases we considered, the doppler returns will all have the same period, but will vary in their peak frequency deviations and their relative phases. Figure 5.14 shows a segment of the backscatter waveform. Since the test signals are 81

300 350 400 450 500 550 600 Time (sample) Figure 5.14: Orbiting Spheres Waveform distinguishable only by their phasing, this figure represents the magnitude of the complex value signal for all five cases. The test signals show dramatic differences when viewed in the frequency domain. Two examples of the fast Fourier transform on the signal gives the spectral magnitude shown in Figures 5.15 and 5.16. Unfortunately, the continuously changing doppler frequencies cannot be discerned from these graphs. These spectral graphs can also used to estimate the period and the maximum doppler frequency for the signal. Figure 5.17 shows one of the possible spectrograms for the Case E signal. This was calculated using a 128 point window, giving a frequency resolution of 4 hertz (compared with the 1 hertz resolution provided by a 512 point binomial transformation). As the figure shows, this form makes it difficult to associate specific frequencies with particular times. The binomial distribution gives a different perspective by taking the signal into the time frequency domain, shown in Figures 5.18 and 5.19. These graphs show two sinusoidal components, one created by each sphere, distinguished by their magnitudes and phases. These differences translate into distance from each sphere to the center, and from the spheres to each other. 82

ULL Ua) 'a.F Figure 5.15: FFT of Case A Waveform I I I - T - I I I l 70 60.. tlI ULL LL a) 'o c (1 50 40 30 20 10 -250 -200 -150 -100 -50 0 50 100 150 200 250 Frequency (Hz) Figure 5.16: FFT of Case E Waveform 83

urlt t I.I.. I 200 150o 100 ' 50 - > c 0 i- -5C) < < -K — I 1 7! IF I P -100 cz - 150 - -200 - -4ou r- I I I I I 300 350 400 450 Time (sample) 500 550 600 Figure 5.17: 128-Point Spectrogram of Case E Waveform 250'. 200 150 -100 -I 50 i a50 300 350 400 450 500 550 600 Time (sample) Figure 5.18: Case A Time-Frequency Representation 84

-150 To rl r2 r0 Actual 146.29 0.50 0.50 1.00 % Error 0.03 0.98 0.98 1.97 Table 5.1: Case A Results 5.4.3 Test results CasTe A symmetric The first rotating sphere target had the rotational center at the midpoint between the spheres, indicated by position A in Figure 5.19: Case E3. This required no additional manipulation of the signal since the rotational center matches the original simulation. Since the spheres were equidistant from the center and traveling at the same angular rate, the velocity of one sphere was always the negative of the other. This means the two spheres impose equal and opposite doppler shifts on the returned signal. Table 5: Cse 5.1 summarizes es to the high level of symmetry in this case, a priori knowledge that the value of To corresponded to half the actual period was used when generating the estimates. used when generating the estimates. 85

T r r2 ro Actual 146.29 0.25 0.75 1.00 Estimate 146.25 0.27 0.76 1.05 % Error -0.02 6.44 0.98 4.75 Table 5.2: Case B Results To rl r2 ro Actual 146.29 0.00 1.00 1.00 Estimate 146.33 0.00 1.00 1.00 % Error 0.03 0.00 0.30 0.30 Table 5.3: Case C Results 5.4.4 On-axis shifts The next three cases share a common difference with the original data set. Case B moved the rotational center closer to one sphere. In case C the rotational center was colocated with the center of one sphere. Case D positioned the rotational center on the extension of the line connecting the sphere centers. In these three cases, the return signal value from each angle requires an additional phase factor, to account for the apparent movement of the original rotational center and the wavelength, A, of the illuminating signal. To move the origin p meters along the axis between the spheres used a phase offset a given by = Ap cos(O) (5.9) The larger errors in table 5.2 are due to difficulties in estimating small lengths due to the small doppler shifts they produce. This same effect also contributes to the large errors encountered in Case E below. The results for Case C, summarized in table 5.3, and Case D, summarized in table 5.4, had lower errors, corresponding to the larger radii being estimated. To r r r2 ro Actual 146.29 0.50 1.50 1.00 Estimate 146.33 0.50 1.51 1.02 % Error 0.03 -0.38 0.53 2.14 Table 5.4: Case D Results 86

To rl r2 ro Actual 146.29 0.38 0.92 1.00 Estimate 146.25 0.41 0.92 1.13 % Error -0.02 6.92 — 0.35 12.58 Table 5.5: Case E Results 5.4.5 Off-axis shift The rotational center is placed at an arbitrary position, specified as an angle and distance from the symmetric center. This case requires two parameters to determine the phase corrections. Like the previous case, the distance p factors in, but it also involves the angle Ooff between the z-axis and the line connecting the original and new centers of rotation. This yields a phase correction of 4rrp = l cos(0 + 0off) (5.10) As shown in Table 5.5, Case E exhibited the largest errors, in part because it was the only case where requiring estimates of the angular offset between the two time-frequency domain traces. In the previous four cases, we place the doppler components either exactly in phase, or 180 degrees out of phase. For the arbitrary rotational center placement, the error estimating the phase difference combines with any errors estimating r1 and r2 to give the total error in the ro value. 5.4.6 Higher order effects There are several higher order effects that can arise in a dual sphere problem. The spherical targets will support a creeping wave, where part of the impinging energy follows the surface of the sphere, shedding energy as it travels around, some of it directed back towards the source. Energy following this path can be distinguished from the primary returns due to an additional time delay from their longer pathat. They do not show up with this time-frequency analysis because the primary (or specular) return and the creeping wave return undergo the same doppler shift. Note that at extremely high sampling rates, two doppler lines could arise, offset in time by the additional travel time of the creeping wave. For a 15 centimeter radius sphere, the wave travels an additional 0.15(7r+ 2) = 0.77 meters, or 2.57 microseconds in free space. This would require a sampling rate over 388 MHz. Another observable effect involved waves reflecting off both spheres in succession before

return to the source. Due to the motion of the reflection points on the surface of the sphere, the path length traveled by the wave between the spheres varies with the observation angle, giving a roughly linear value over each half cycle of a rotation. The energy involved in this mechanism is very low compared to the specular peaks, and does not show up in most of the sample runs. 5.4.7 Customized kernel Since this problem was constrained to simple rotational motion, it is possible to design a kernel optimized to extract the expected doppler shapes. I did not explore that option since the rotating spheres were only serving as an example to justify research using more complicated targets. Customizing the kernel to the problem may provide dramatic improvements in the resolution and fidelity of the results, particularly when the signal is buried in noise [45]. 5.5 Summary Time-frequency techniques worked very well on this problem. Examining the binomial distribution allowed us to extract information not readily available in other representations of the signal. From this information, we could estimate the physical parameters of the target with good accuracy, with only two estimates out of the group of twenty exceeding a five percent error. The techniques used here are applicable to other types of problems. Passive systems which measure emitted energy could also use this type of time-frequency analysis. In astronomy, an analogous problem might be the signals obtained from a binary star system. 88

CHAPTER 6 EFFECTS OF ADDING NOISE AND COMPLEXITY TO SCATTERED SIGNALS The previous chapters presented the application of time-frequency techniques to several electromagnetic problems. While the results have been promising, their applicability is limited due to the idealized conditions employed. Using simulations for the experiments provided signals free of the noise components with which all physical systems must contend. Evaluating this aspect is essential since other promising techniques for scattered signal analysis, such as the singularity expansion method (SEM), have shown a loss of effectiveness once a realistic level of noise was added to the system [4, 26]. For this chapter, noise signals were inserted at various power levels to determine the robustness of the time-frequency techniques. The signal degradation was monitored using the Renyi third order entropy [44] and the effects of increasing the target's complexity examined. 6.1 Background An electronic system, especially one based on receiving scattered electromagnetic waves, must deal with a variety of noise sources. These sources may include thermal noise, shot noise, atmospheric noise, solar noise, cosmic noise, and urban noise. In some instances, these contributions can overpower the desired signal, making it impossible to extract any useful information. While some systems can increase their signal power to make the noise contribution insignificant, many radars cannot do this due to size, cost, or power source limitations. Instead, the signal processing subsystems in the radar must compensate to extract useful information from the imperfect received signal. To investigate the impact of noise, a simplified scattering model using a point target 89

was first investigated. This scattered signal was then corrupted by a noise source, in this case modeled using band-limited additive white Gaussian noise (AWGN). The third-order Renyi entropy provided a quantitative measure of the noise's impact. Various amounts of noise were then added to corrupt the signal to see how the received information degrades as the signal to noise ratios (SNR) decreases. Finally, these techniques were repeated to corrupt the signals from the two orbiting spheres discussed in chapter 5. 6.2 Signal Models This study examined four target sets, two based on ideal scatterers and two based on realistic physical models. The data sets represent the scattered field from a target positioned along the y-axis and rotated around the z-axis. The targets were rotated at a rate of 3 revolutions per second and illuminated with a 2 GHz carrier. The scattered field was sampled for one second at 512 hertz. Placing the target pieces within 0.75 meters of the origin yields a maximum expected doppler shift of 2v fd (6.1) 47r3(0.75) (6.2) 0.15 = 188.5 hertz (6.3) 6.2.1 Point Targets To focus on the noise effects, the simplest possible data signal was used. Therefore, instead of the accurately modeled sphere scattering discussed earlier, a point target was placed in the far field of the radar antenna, traveling in a circular orbit of radius r as illustrated in figure 6.1. Since the target is in the far field, the return amplitude remains constant as the range changes slightly with target motion. The target motion is significant with respect to illuminating signal wavelength so the phase will depend on the target-radar separation. The target's angular position at time t depends on the frequency of rotation (frot) and is given by 27frott. For the carrier frequency (fc), the signal will have a wavelength A 300, 000km/s (6.4) J c This means the signal will undergo a round-trip phase shift (0) of 0(t) = 47rcos(27frott) (6.5) Ac

L.75 meter r Z " radar incident scattered Figure 6.1: Point Target Geometry radians. The returned signal is g(t) = ei(t) (6.6) With the target traveling in a circular path, the range to the target is given by a sinusoidal function. This means the time-frequency distribution, based on the velocity, or the derivative of the range, is also given by a sinusoid. The doppler shift imposed on the signal is based on the instantaneous radial velocity Vrad(t) = rsin(21rfrott) (6.7) and the carrier frequency wavelength 2Vrad(t) fdop(t) = - ---- (6.8) Figure 6.2 shows the difference between a traditional Fourier transform and a time frequency transformation for an orbiting point target. The binomial distribution clearly shows how the target's doppler changes as a function of time. The first two targets consisted of idealized point scatterers. The one point target placed the scatterer at the point y = 0.75 (0.75 meters from the center of rotation). The scattered signal from the single point target was unique since it exhibited no power fluctuations, only variations in phase. The two-point target, shown in figure 6.3, used this scatterer and added a second scatterer at y = -0.25 (0.25 meters on the opposite side of the center of rotation). Due to constructive and destructive interference between the scatter fields from the two points, variations occur in both the magnitude and phase of the signals. The simplified nature of these signals made them suitable for evaluation using the receiver operating characteristic. Calculating the receiver operating characteristic requires 91

Fourier Transform 80 -70 -60 i 50 40 20 Binomial Distribution -250 -200 -100 1 -50 i I I,,,, S, 0 ( I 250 I 2/ \JI 250 --- —— I 0 0.2 0.4 0.6 0.8 Time (sec) for Single Point Scatterer -ZUO -10u U 1WU 0UU Frequency (Hz) Figure 6.2: Transformations [.75 meter - 1 t I L.25 meter incident radar scattered Figure 6.3: Two-Point Scattering Geometry 92

.75 meter z radar incident scattered.15 meter Figure 6.4: Two-Sphere Scattering Geometry an apriori time-frequency distribution to determine where the targets components should appear. This definitive template for distinguishing false alarms from legitimate detection is only possible for the simplest targets. 6.2.2 Distributed Target The other two targets were more difficult to analyze since the ideal point scatterers were replaced by numerical simulations modeling physical targets which generated more complicated scattered fields. The first realistic target, shown in figure 6.4, consisted of two metallic spheres, each 0.15 meters in radius. The centers of the spheres were positioned at y = 0.75 and y = -0.25 meters. The final target, shown in figure 6.5, was a 1 meter long wire placed along the y-axis between the points y = 0.75 and y = -0.25 meters. The scattered signals from these targets were not suitable for evaluation using the receiver operating characteristic since one could not dictate a priori what the time frequency representations should be for these returns. Each target signal was normalized to have a mean square value of 1.0. To create a noise-corrupted version, the target signal was added to an appropriately scaled noise signal. The noise signal was modeled as complex valued, with both the in-phase and quadrature channels having uncorrelated, normally distributed values with mean values of zero and equal variances. 93

.75 meter z radar CIA incident scattered Figure 6.5: Wire Scattering Geometry 6.3 Noise Model Adding noise to the signal required a computational method of modeling the noise sources. My goal was to approximate the type of signal degradation that would be found in a physical radar system. The first step in developing the noise model was to identify the expected noise sources and their impact on the received signal. The second step was to determine a mathematical model for simulating these sources of noise to add to the desired signals. 6.3.1 Noise Sources All electrical systems are subject to thermal noise due to the random movements of electrons within the system components. These random movements are directly related to the temperature of the components since temperature on the macroscopic level is due to kinetic energy on the atomic level. Any lossy component (including the antenna) will generate small noise signals at temperatures above absolute zero. A special case of thermal noise called shot noise affects systems which deal with small signals and large amounts of amplification (like a radar). This comes from the random creation and loss of carriers in the amplifying components of the radar. In addition to its own thermal noise, an antenna can also receive external noise signals, such as atmospheric noise, solar noise, cosmic noise and urban (man-made) noise. In most instances, thermal noise has the greatest impact on (94

system performance [27]. 6.3.2 Additive White Gaussian Noise While removing noise from the random movements of electrons and charge carriers is impossible, modeling their impact on a system employs a straight-forward technique. Concern here was with the combined effect of an uncounted number of individual contributors, all following the same statistical model. Combining all these independent contributors and applying the central limit theorem leads to an approximately gaussian distribution for the noise voltage. Since the signals of interest are complex valued, the noise model must provide complex valued samples. The received signal can be viewed as two orthogonal channels: the in-phase channel, I. corresponding to the real portion of the complex signal; and the quadrature channel, Q, corresponding to the imaginary portion of the complex signal. For a narrowband additive white Gaussian noise source, the noise signal is given by n(t) = ncr(t) cos 27rfot + ns (t) sin 27rfot (6.9) where fo is the operating frequency of the radar, and the functions nc(t) and rins(t) represent a Gaussian distributed, zero-mean random amplitude for the in-phase and quadrature channels. With orthogonal I and Q channels, the values nc(t) and rins(t) are independent and identically distributed. The noise power is given by Pn = E[n2(t)] = E[n2(t)] + E[n2(t)] = 2a2 (6.10) While the noise in each channel has a Gaussian distribution, the magnitude of the noise, In(t)l = x/lncl2 + ns2I, follows a Rayleigh distribution, and the phase is uniformly distributed between 0 and 27r radians. This representation of additive white noise has assumed a band-limited receiver. The receiver band-limits the noise using a low-pass filter after mixing the scattered return with the carrier frequency. 6.4 Renyi Entropy With all the processing performed on the scattered signals, one needed a method to evaluate the resulting time-frequency distributions. At first viewed the processed time-frequency distributions were studied and qualitatively judged how well expected features appeared in 95

the processed return. While this qualitative approach could distinguish broad differences when processing simple signals, it was inadequate for differentiating the noise corrupted cases. These cases required a more objective, quantitative measure of performance. This figure of merit would also allow the development of automatic, iterative techniques. One measure of the complexity of a time-frequency distribution is the entropy. For a two-dimensional signal Cs(t, f), the Shannon entropy [17] would be H(Cs) =-J C(t f)log2Cs (t f)dtdf (6.11) There are some major problems with this entropy measure. Due to the log2 operation, this definition will only apply if the distribution Cs is greater than zero for all values of t and f. Zero values often occur in time-frequency distributions due to the limited support of time-frequency components, and negative values are a result in many distributions which allow negative values in order to satisfy the marginal distributions. To avoid zero and negative values in the entropy calculation, a more generalized entropy definition was used. For a bivariate density, P(x, y), the Renyi entropy is given by 1-a f f P(x, y)dxdy Ha(P) = 1og2 (6.12) When a = 1, the Renyi entropy matches the Shannon definition. By allowing the summation of values to occur before taking the logarithm, the Renyi entropy can avoid the zero and negative values that cause the logarithm to fail. Converting the double integrals with respect to time and frequency (f f dxdy) into a double summation ( E 6t6f) yields the Renyi entropy for a discrete distribution 1 / s(n, k H, (C [n, k]) = 1 I EZ C ([n'l k'])) +log2 6tf (6.13) This form is easily coded in a computer. To compare the relative entropy between two distributions, the log2 6t6f term is constant and can be dropped. The Renyi entropy exhibits five key properties when analyzing multi-component signals [3]: * Ha attempts to measure the number of components * Ha does not count cross components for odd values of a > 1 * Ha has an upper and lower bound * Ha values are invariant to signal time and frequency shifts 96

* Ha exhibits extreme sensitivity to the phases of closely spaced components. The first four properties are beneficial for this application, but the last property can cause some problems. This phase sensitivity manifests itself as misleadingly high or low entropy values for a distribution due to the contribution of cross-term components. Using reduced interference distributions reduce the cross terms and help avoid this sensitivity. In general, a can take on any positive value greater than or equal to one. This restriction guarantees the existence of the entropy measure [3]. However, for non-integer values of a, the entropy will yield complex numbers of limited use. When a takes on an even integer value, the entropy expression loses the ability to distinguish positive and negative values in the time-frequency distribution. This has led to the suggestion of using a = 3 for the entropy calculations [44]. In this study, the Renyi entropy values provided two different types of information. First, the entropy represents how much uncertainty exists in the distribution. Noise will increase the uncertainty, so a good processing technique will reduce the noise components in the final distribution, and therefore reduce the entropy. The second use for the Renyi value is to identify the number of significant signal components. The larger the entropy value, the more scatterers can be placed in the target field. With the proper selection of 6t and 6f, the number of scatterers can be estimated using [3] Nscat = 2Ha(C) (6.14) Unfortunately this value is valid only for idealized time-frequency distributions. Typical results have the components of the target signal spread in both the time and frequency directions, resulting in a reduction of the signal uncertainty, and an underestimated value for Nscat. 6.4.1 Distribution Evaluations After finding the time frequency distribution, a numerical metric objectively evaluated the result. The first measure used was the receiver operating characteristic [38]. This involved generating an ideal distribution for comparison with the signal's time-frequency distribution. Applying a threshold to the realized signal should give a matching distribution. Errors can be expressed in radar terms, with a probability of detection (Pd) and probability of false alarm (Pf a) for each threshold. Plotting Pfa versus Pd provides a way to compare different distributions and the impact of adding different amounts of noise. Ideally, the curve would allow a -Pd approaching 1.0 while keeping Pfa approximately 0. Graphically, 97

the closer the curve gets to the upper left corner of the plot, the better the performance. Figure 6.6 shows somne representative receiver operating characteristic curves. The other measure used was the Renyi entropy of the distribution [3]. This metric attempts to measure a distribution's information content. For a discrete time frequency distribution, Cg[n, m], the Renyi entropy is given by H3(C9[ in m]) =- log2 E (: 92 / 9, m ] ) log (6.15) 2 n m \Zn1 LEM'/ C91n m where 6t and 6f are the time and frequency resolution. Since entropy measures the amount of disorder in a system, the maximum entropy occurs for the distribution of a noise only signal. As the signal to noise ratio improves, the entropy decreases. Using different time-frequency transform kernels generated distribution with widely varying entropy values. This made it difficult to compare values from different distributions. For example, the high resolution nature of the Wigner distribution consistently provided lower entropy values than the binomial distribution for the same signal. Instead of the absolute entropy value from different transforms, the entropy margin, or the entropy improvement (in bits) over a pure noise signal was emphasised. 6.5 Analysis and Results The analysis for each target considers the scatter signals with signal to noise ratios (SNRs) between 10 down to -10 decibels, with the special cases of zero noise (SNR = oo) and zero signal (SNR = -oo). They were transformed to the time-frequency domain via the short-time Fourier transform (spectrogram), the binomial transform, and the Wigner transform. 6.5.1 Receiver Operating Characteristics Figure 6.6 shows some representative receiver operating characteristic (ROC) curves. Each plot shows the ROC for the spectrogram (lowest curve in each set), the Wigner distribution (center curve with circled points) and the binomial distribution (top curve). The binomial and Wigner distributions significantly outperformed the spectrogram, with the binomial typically holding a slight edge. The differences were most apparent with the two-point scatterer since the additional signal component led to more cross terms (and more errors) in the Wigner distribution. 98

0,. P.-nt T.,grt 0. 9 0. a I 0 0.7 I 02:4 0.4 P-ob f F.I.. Al.rrr (Rf.) T-. P..,,t T.1Qrt a.4.3.2 0 0.2 0 '4 0 a O a I Fl,.b.1 F.I.. Al*,- (Pf.) Figure 6.6: Receiver Operating Characteristics 1.8 1.6 1.4 1.2 1 0.8 0.6 0.4 0.2 0 - -0.2 0 One Point 2 4 6 Case (inf: 10:3:0:-3:-10 dB) Wire 1.8 1.6 1.4 -0.8 0.6 0.4 0.2 0 - -0.2 -0 2 4 6 C:ase (inf:10:3:0:-3:-10 dB) Figure 6.7: Entropy Margins 6.5.2 Entropy Measures The performance advantage of the binomial distribution over the Wigner for complicated signals was also apparent from the entropy calculations. Figure 6.7 shows the entropy improvement over pure noise for the simple single point scatterer and for the wire. With the single point, the Wigner (shown with circled points) was the top performer, beating the binomial at high SNRs. But as the SNR decreased, the Wigner's performance dropped to roughly match the binomial. As the SNR dropped below -3 dB, the smoothing properties of the spectrogram gave it a slight advantage, although all three distributions performed poorly. For the complicated scattered signal off the wire, the results were a little different. Once again the spectrogram yielded the smallest entropy improvement over pure noise. 99

The existence of significant interfering cross terms dropped the Wigner's entropy margin (marked with circles) almost down to the spectrogram level. The binomial, a reduced interference distribution, performed significantly better than either the spectrogram or the Wigner. This advantage was evident until the SNR dropped to -10 dB, at which point the target return was indistinguishable from the noise. 6.6 Summary The characteristics of the signal, as well as the signal's quality, as measured by the signal to noise ratio, have a significant impact on how well different time-frequency distributions will perform. The receiver operating characteristic can give a clear idea of the transformation performance but requires an idealized distribution for comparison. The inability to determine an ideal distribution makes the receiver operating characteristic difficult to apply except for simple canonical targets. The Renyi entropy values did not require any additional information about the signal and could be applied to any signal, giving results consistent with qualitative evaluations of the distributions. The work presented in this chapter also demonstrated how a reduced interference distribution like the binomial can outperform the Wigner. Although the Wigner can provide higher resolution with respect to both time and frequency, this resolution comes at the expense of large interfering cross terms that can degrade the time frequency representation. 100

CHAPTER 7 DEVELOPMENT AND USE OF A CUSTOMIZED KERNEL The Wigner and binomial distributions represent general purpose time-frequency kernels. Although they may provide excellent time-frequency representations for many signals, they make no effort to exploit characteristics such as periodicity which may exist in a signal. This section examines the use of customized time-frequency kernels which improve on the performance of the general kernels. To examine how customized kernels might be applied, a genetic algorithm was developed to maximize the Renyi entropy of the time-frequency distribution. The kernel design was constrained to provide a kernel with reasonable properties, as outlined in chapter 2. 7.1 Background Due to the generality of Cohen's class of transformation, essentially an infinite number of different time frequency representations can exist for a given signal [11]. Based on the results in chapter 6, the goal was to find a kernel which provided time-frequency distributions with higher Renyi entropy values than the Wigner or the binomial. 7.1.1 Genetic Algorithms Genetic algorithms are an optimization technique that models itself on the process used in evolution and natural selection [20]. Parameters related to the value to be optimized are coded into genes, and sets of genes are put together to make chromosomes. Allowing the computer to use survival of the fittest to update the available gene pool serves as an optimization technique. This has been applied to a wide variety of problems, ranging from 101

economic prediction [20] to pattern classification [35] to array antenna design [16]. The genetic algorithm was chosen for this problem due to the complicated mapping between the parameters driving the kernel and the resulting cost function. It also provided an opportunity to apply an optimization technique that has not been used in this area. Basic Procedure After deciding how to encode the necessary parameters genetically, the procedure followed is simple. An initial random population is generated. These genes are rank ordered based on their fitness to solve the selected problem. Those scoring the highest are carried into the next generation or iteration. The poorest performers are discarded. The next step is to fill the vacancies in the population. This is done by performing a genetic crossover operation between two parents selected from the previous generation's fittest chromosomes. Once this new generation is filled, the fitness of the new members is evaluated, the chromosomes are rank ordered, and the process of discarding members and creating a new generation are repeated. Terminating the Algorithm The genetic algorithm is a global search technique that typically performs well, but it does not have guaranteed convergence. The code implemented used two criteria for stopping. The first was based on the top performer not changing over a given number of generations. The second condition limited the total number of iterations. 7.1.2 Problem Design The approach followed to select the elementary function, h(t), needed for the kernel was to use a genetic algorithm to select candidates. According to the design rules [22], h(t) should be a smooth function defined over the range -0.5 < t < 0.5 which goes smoothly to zero at the ends when Itl = 0.5 and has an area of 1. The elementary function was defined using a cosine based series N h(t) = Z an cos ((2n - 1)27rt) (7.1) n=l where there are weight values, at, which are assigned values in the range 0 < a < 1. These weights are applied to the odd harmonics of the cosine, so each term will have a maximum at t = 0 and go to zero at Itl = 0.5. 102

The target used for these trials was a point target traveling in a circular path, similar to the case shown in figure 6.1. The point was rotated at 5 hertz a distance of 1 meter from the center of rotation while illuminated with a 1 gigahertz plane wave. The scattered signal was sampled at a rate of 512 hertz. These values were selected to provide several cycles over a 1 second record without doppler frequency ambiguities. 7.2 System Optimization Results Once the basic encoding was selected, trials were conducted to see how well the genetic algorithm could select a kernel. The parameters chosen were 5 weighting coefficients, each coded using 5 bits. The population size was set to 30, with a maximum of 15 generations before terminating the genetic algorithm. While using more weights, more bits, a larger population, and a longer running algorithm should allow greater performance, the scope was selected to show the feasibility and to keep the problem scaled for execution on a 486 computer system. To generate a time-frequency distribution for the test signal, the genetic algorithm selected the elementary function shown in figure 7.1. This function was selected over a Elementary Function O ' s \ - 0.8 - 0.7 - 0.7 -0.4 1 0.3 - /\ 0.2 - 0.1 K -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 Figure 7.1: Genetically Derived Elementary Function period of 15 generations, during which time the third-order Renyi entropy of the resulting time-frequency distribution increased from an initial value of 11.2000 to the final value of 103

11.3174. For comparison, the third order Renyi values for the Wigner and the binomial distributions of this signal are 10.1904 and 11.0409. Figure 7.2 shows the resulting time-frequency distribution for the point scatterer. This Genetic Derived TFD -30; -20 -__ ___ 2: ' ~: t 30i 50 100 150 200 250 300 Time Index I.350 400 450 500 Figure 7.2: Genetically Derived Time-Frequency Distribution is very similar to the distribution obtained using the binomial kernel. While the genetically derived distribution shares the desirable distribution properties of the binomial discussed in chapter 3, the genetic algorithm did provide slightly lower cross terms near the maximum doppler points. 7.3 Summary The results of using the genetic algorithm for this problem were very promising. The results were very close to the binomial distribution, showing the technique could meet the standard results. This approach needs to be taken further to try and exploit characteristics such as periodicities in the signal. It should be possible to have the kernel tuned to certain signal parameters such as the general shape and period of the signal. Used with a genetic algorithm, this would provide an adaptive processing technique giving clearer time-frequency distributions than the standard Wigner and binomial kernels. 104

CHAPTER 8 CONCLUSIONS This work was interdisciplinary, mixing elements and concepts from electromagnetics with those of signal processing. The purpose was to find some new applications for some popular and powerful signal processing techniques. 8.1 Stationary Targets Examining stationary targets, using a sweep of frequencies, allowed estimating the range to scattering centers on the target. By restating the transformations used in time-frequency analysis for use in frequency-time analysis, it was easier to determine what scattering mechanisms were contributing to the scattered signal. The transformations not only returned information on the position of the scattering centers, but also provided insight on the dispersion, if any, exhibited by the scattering mode. The understanding provided through this type of analysis can be applied to help control the scattering from a target in either a monostatic or bistatic situation. 8.2 Dynamic Targets The situation with analyzing a dynamic target was very different. Here the radar generated a continuous wave signal. The information was contained in the doppler shift imposed on the signal when scattered by the target. Performing time-frequency analysis on the received signal generated a distribution showing the target's doppler components as a function of time. This doppler history can be used to estimate different target parameters. 105

8.3 Noise Statistics Looking at how noise is processed by a transformation belonging to Cohen's class provided an understanding on the impact of noise and how the noise in the time-frequency distribution relates to the transform selected. The results show a significant advantage for the reduced interference distributions, such as the binomial, over the standard Wigner. 8.4 Effects of Noise and Complexity Time-frequency distributions are easiest to interpret when there are few signal components and a very high signal to noise ratio. As more components appear in the timefrequency distribution, they allow the formation of cross terms, typically degrading the clarity of the distribution. This degradation occurs whether the additional component was actually generated by the target or was simply caused by noise. The third order Renyi entropy proved to be a valuable metric for numerically evaluating the quality of a time-frequency distribution. These entropy values corresponded to the subjective quality ratings an experienced viewer would assign to the different distributions. The results showed the superiority of the reduced interference distribution except when dealing with very simple targets with almost no noise corrupting the signal. 8.5 Customized Kernels The work done using a genetic algorithm provided performance exceeding the Wigner and binomial distribution for a simple target. Before exploring this area more thoroughly, careful consideration must be given to the fitness function and the genetic encoding used. 8.6 Value of Time-Frequency Techniques This research has shown time-frequency analysis as a valuable tool for analysis of radar and other electromagnetic based signals. The ability to identify and evaluate different scattering modes provides insight not readily apparent in other signal representation. 8.6.1 Image Processing Once a signal has been processed using a time-frequency transformation, the resulting function of two variables is similar in many aspects to an image. It appears that applying 106

image processing techniques, such as noise reduction or edge enhancement, may provide the next logical step for processing scattered signals. The image processing approach also leads into the area of target classification. 8.6.2 Ultrawideband Radar The processing of ultrawideband radar signals is a natural extension to the work in this dissertation. Since an ultrawideband radar cannot apply narrowband assumptions, the ability to process a signal in a joint time-frequency domain should prove beneficial. 107

APPENDICES 108

APPENDIX A Notational Conventions Continuous Discrete q(v, T) q(q, I) ambiguity domain kernel,4(t, 7) (n, I) autocorrelation domain kernel Cg(t, f) Cg(n, p) time frequency distribution g(t) g(n) time domain signal gr (t) gr(n) real portion of signal gi (t) gi(n) imaginary portion of signal G(f) G(q) frequency domain signal f q frequency of interest v, A q, b dummy frequency variable t n time of interest u, v m, c dummy time variables T, s l, d lag variables Ot convolution along the t axis E[x] expected value of x (also written x).tf Fourier transform from g(t) to G(f) Ftl1 inverse Fourier transform from G(f) to g(t) 109

APPENDIX B Computer Subroutines This appendix lists the major Matlab routines used in this research. The routines assume the signal will be provided as a one-dimensional vector. The values may be real or complex. The values may also be pre-processed to make the input signal analytic which reduces the number of cross terms in the final distribution. 110

B.1 Function myspct As mentioned in the text, a spectrogram is also known as a short term Fourier transform and is generated by centering a window function at the time of interest and taking a regular Fourier transform of the windowed data. function spct=myspct(a,hl); % spct=myspct(a,hl) Calculate the two-sided spectrogram % for the given vector (real or complex). %. If the transform's half length (to be % consistant with mybtfd and mywtfd) it % defaults to half the vector length. % Values returned are arranged so the rows %. correspond to frequencies from -hl to hl-1, and % columns correspond to time bins (0 to N-l). % cjm - 31jan96 [n,m]=size(a); % n=#samples if n*m > max([n m]) % only true for error spct='myspct error - only works with a vector' else if n==l; n=m; m=l; a=a.';end % make input a column vector tl=n; if nargin~=1; tl=hl*2;end % and set transform length tmp=[zeros(tl/2,1); a; zeros(tl/2,1)]; spct=zeros(tl,n); for i=l:n, spct(:,i)=fftshift(fft(tmp(i:i+tl-1))); end spct=abs(spct).^2; end 111

B.2 Function myfwig This routine generates the Wigner distribution for the input vector. The frequency resolution depends on the maximum time lag applied when calculating the local autocorrelation. This routine works with the alias-free formulation with calculations performed in the autocorrelation domain [21]. function wtfd=myfwig(a,l); % wtfd=myfwig(a,l) Calculate the two-sided wigner time-frequency 7, distribution for the given signal vector. % If the maximum lag (1) is not specified, it %, defaults to the vector length. % Values returned are arranged so the columns % correspond to frequencies from -1 to 1-1, and % rows correspond to time bins (0 to N-l). I. Z/, this routine combines code from mylac.m and mywig.m - also parallels mybtfd.m % cjm - 27aug95 - modified 2apr96 to save memory [m,n]=size(a); % if n*m > max([n m]). tfd='mybtfd error - only works else if n==l; n=m; m=l; a=a.';end % ml=n; if nargin~=l; ml=l;end % n=#samples, m will be used for lag only true for error with a vector' make input a row vector and set maximum lag value. first find the local autocorrelation maxd=min([n,ml]); % limit loop to range of delays ac=conj(a); % conjugated a a=[a zeros(l,maxd)]; % padded a lacm=zeros(maxd,n); for row=l:maxd; lacm(row,:)=a(row:row-l+n).*ac; 112

end. next apply the wigner kernel ', we know [ml,n]=size(lacm); for r=l:maxd/2; row=r*2; lacm(row-l,:)=[zeros(1,r-1) lacm(row-l,1:n-r+l)]; lacm(row,:)=0.5*([zeros(1,r-1) lacm(row,1:n-r+1)]... +[zeros(1,r) lacm(row,1:n-r)]); end; %, finally, convert from autocorrelation domain to time-frequency domain wtfd=real(fft([lacm; zeros(l,n); conj(flipud(lacm(2:ml,:)))])); 7. (put negative frequencies first) wtfd=[wtfd(ml+1:2*ml,:); wtfd(l:ml,:)]; end 113

B.3 Function myfbin This routine generates the binomial distribution for the input vector. The frequency resolution depends on the maximum time lag applied when calculating the local autocorrelation. This routine works with the alias-free formulation with calculations performed in the autocorrelation domain [21]. function btfd=myfbin(a,l); % btfd=myfbin(a,l) Calculate the two-sided binomial time-frequency 7, distribution for the given signal vector. % If the maximum lag (1) is not specified, it % defaults to the vector length. % Values returned are arranged so the columns % correspond to frequencies from -1 to 1-1, and % rows correspond to time bins (0 to N-l). % this routine combines code from mylac.m and mybin.m % - also parallels mywtfd.m % cjm - 23aug95 updated to reduce memory use 2apr96 [m,n]=size(a); % n=#samples, m will be used for lag if n*m > max([n m]) % only true for error tfd='mybtfd error - only works with a vector' else if n==l; n=m; m=l; a=a.';end % make input a row vector ml=n; if nargin-=l; ml=l;end % and set maximum lag value % first find the local autocorrelation maxd=min([n,ml]); % limit loop to range of delays ac=conj(a); % conjugated a a=[a zeros(l,maxd)]; % padded a lacm=zeros(maxd,n); for row=l:maxd; lacm(row,:)=a(row:row-l+n).*ac; 114

end % At this point, lacm holds half the autocorrelation matrix, with all % values slid to the left. This means the second row should be shifted % one HALF column to the right, the second row should shift one full % column to the right, etc. This shifting will be accounted for when % the kernel is applied below. /.'%%%% % next apply the binomial kernel % we know [ml,n]=size(lacm); krn=[.5.5]; for row=2:maxd;; lacm(row,:)=conv(lacm(row,l:n-row+l),krn); krn=conv(krn,[.5.5]); end; % The matrix lacm now contains values ready for transforming via the % FFT into the time/frequency domain. The values should correspond % to the upper half of the matrix generated via use of dosauto and % bin.m in the FTP distributed files. %%%%% % finally, convert from autocorrelation domain to time-frequency domain btfd=real(fft([lacm; zeros(l,n); conj(flipud(lacm(2:ml,:)))])); % (put negative frequencies first) btfd=[btfd(ml+1:2*ml,:); btfd(l:ml,:)]; end 115

B.4 Function renyi This function finds the Renyi entropy for a two-dimensional array [44]. These function served as the fitness function for the genetic algorithm kernel. function h=renyi(cs,a,dt,df) % h=renyi(cs,a,dt,df) Calculate the Renyi entropy for the time% frequency distribution cs (matrix) using 7, default order (a) of 3, and default delta % values (dt, df) of 1 (for 0 bits) % cjm - 10oct95 if nargin<4; df=l; end % handle default values if nargin<3; dt=l; end if nargin<2; a=3; end h = log( sum(sum( (cs/sum(sum(cs))).^a )) * dt * df) / (log(2) * (1-a)); 116

B.5 Function genalg This genetic algorithm tries to find a Cohen's class distribution that satisfies most of the desirable properties outlined in chapter 2 and maximizes the Renyi entropy of the timefrequency distribution. The routine uses integer coding of the genes and keeps the top half of each generation to serve as the parents for the next generation [20, 34]. function [mg,mx]=genalg(sig,ml,nw,bw,pop,mr,mi) %. [mg,mx]=genalg(sig,ml,nw,bw,pop,mr) %/ genetic algorithm to find the 'best' kernel for a Cohen's %7 class TFD. The signal (sig) goes through a transform with %7 max lag ml and a kernel using a cosine series with nw odd %7 harmonics, with each parameter weight coded in bw bits.,% The gene pool is pop members long. The algorithm will %. iterate until the maximum performer is constant for mr %7 cycles or mi cycles are completed. 7. Returned values are the 'best' gene per iteration and the scores. 7. cjm - 30jun96 iter=0; %7 iteration count rpts=1; %7 number of repeats gl = nw*bw; %7 gene length mx = []; mg = []; []; ]; 7. zero out the 'best' arrays pool = floor(rand(pop,gl)+.5); %7 randomize the gene pool score = zeros(pop,1); %7 clear out scores; for i=pop/2+1:pop %7 after first time, we'll know these score(i)=real(renyi(myfgen(sig,ml,pool(i,:),nw,bw))); end while( (iter<mi) & (rpts<mr) ), for i=pop/2:-1:1 score(i)=0; for pt=i+l:pop if prod(pool(i,:)==pool(pt,:))==1, score(i)=score(pt);end end 117

if score(i)==O, score(i)=real(renyi(myfgen(sig,ml,pool(i,:),nw,bw))); end end [score,i]=sort(score); % arrange by lowest to highest entropy pool=pool(i,:); mx=[mx; score(pop)]; % keep info on best mg=[mg; pool(pop,:)]; rpts=rpts+1; if prod(lb==pool(pop,:))==0, rpts=l; end lb=pool(pop,:); iter=iter+1; for i=1:2:pop/2 % propagate next generation cp=floor(rand(1)*(gl-2)+2); % pick crossover point pool(i,:) = [pool(pop+l-i,l:cp-1) pool(pop/2+i,cp:gl)]; pool(i+l,:) = [pool(pop/2+i,l:cp-1) pool(pop+l-i,cp:gl)]; end cp=floor(rand(l)*gl+l); i=floor(rand(1)*pop/2+1); pool(i,cp) = abs(pool(i,cp)-l); [iter rpts mx'] end % pick random mutation point % limit to new offspring % flip the bit % optional feedback to user 118

B.6 Function myfgen This function generates a time-frequency distribution using the kernel as represented by the genetic algorithm. The returned matrix exactly matches the layout generated by myspct.m, myfwig.m, and myfbin.m. function tfd=myfgen(a,l,code,nw,bw); % tfd=myfgen(a,l,code,nw,bw) % Calculate the time-frequency distribution using % a weighted cosine kernel. The n weights are coded % values between zero and 1 using bw bits. The % weights apply to an expansion of odd cosine % harmonics. % Values returned are arranged so the columns % correspond to frequencies from -1 to 1-1, and % rows correspond to time bins (0 to N-1). % this routine combines code from mylac.m and mybin.m - also parallels % mywtfd.m - derived from myfbin.m % cjm - 23aug95 updated to reduce memory use 2apr96 % cjm - 29jun96 recoded for genetic decoding and weighting [m,n]=size(a); % n=#samples, m will be used later for lag if n*m > max([n m]) % only true for error tfd='mybtfd error - only works with a vector' else if n==l; n=m; m=l; a=a.';end % make input a row vector ml=l; % first find the local autocorrelation maxd=min([n,ml]); % limit loop to interested range of delays ac=conj(a); % conjugated a a=[a zeros(1,maxd)]; % padded a lacm=zeros(maxd,n); 119

for row=l:maxd; lacm(row,:)=a(row:row-l+n).*ac; end % At this point, lacm holds half the autocorrelation matrix, with all % values slid to the left. This means the second row should be shifted % one HALF column to the right, the second row should shift one full % column to the right, etc. This shifting will be accounted for when % the kernel is applied below. %%%%% % new code for genetic coding y. bwts = [-1 2.^(-1:-1:-bw-1)]; % numeric weights for each code bit mf = 1:2:nw*2; % multiplier for odd harmonics wts = zeros(nw,l); % find weight values whos;pause for i=l:nw wts(i)=code((i-l)*bw+*b+l:ibw)*bwts; % code must be row vector end for row=2:maxd % now apply kernel to autocorrelation krn=cos( ( (l:row)'*pi/(row+l)-pi/2 ) * mf ) wts; krn=krn/sum(krn); % make sure row sums to one lacm(row,:) = conv(lacm(row,1:n-row+1),krn); end y. % end of new code T m xf'/.'.'.m'. y. % The matrix lacm now contains values ready for transforming via the 120

% FFT into the time/frequency domain. The values should correspond % to the upper half of the matrix generated via use of dos.auto and % bin.m in the FTP distributed files. % finally, convert from autocorrelation domain to time-frequency domain tfd=real(fft([lacm; zeros(l,n); conj(flipud(lacm(2:ml,:)))]));. (put negative frequencies first) tfd=[tfd(ml+1:2*ml,:); tfd(l:ml,::)]; end 121

BIBLIOGRAPHY 122

BIBLIOGRAPHY [1] NM. Amin, L. Cohen, and W. J. Williams. Methods and Applications for Time Frequency Analysis, August 1993. [2] NM. G. Amin, "Minimum variance time-frequency distribution kernels for signals in additive noise," IEEE Transactions on Signal Processing, vol. 44, no. 9,, September 1996. [3] R. Baraniuk, P. Flandrin, and O. Michel, "Measuring time-frequency information content using the Renyi entropies," Technical Report 9512, Rice University, July 1995. [4] C. E. Baum, E. J. Rothwell, K.-M. Chen, and D. P. Nyquist, "The singularity expansion method and its application to target identification," Proceedings of the IEEE, vol. 79, no. 10, pp. 1481 —1492, October 1991. [5] G. J. Burke, "Numerical electromagnetics code -- NEC-4 method of moments, Part I: User's manual," Technical Report UCRL-MA-109338, Lawrence Livermore National Laboratory, January 1992. [6] L. Carin and L. B. Felsen, "Wave-oriented data processing for frequency and timedomain scattering by nonuniform truncated arrays," IEEE Antennas and Propagation Magazine, vol. 36, no. 3, pp. 29-43, June 1994. [7] H. I. Choi and W. J. Williams, "Improved time-frequency representation of multicomponent signals using exponential kernels," IEEE Transactions on Acoustics, Speech and Signal Processing, vol. ASSP-37, no. 6, pp. 862-871, June 1989. [8] T. A. C. M. Claasen and W. F. G. Mecklenbrauker, "The wigner distribution - a tool for time-frequency analysis, part i: Continuous-time signals," Philips Journal of Research, vol. 35, no. 3, pp. 217-250, 1980. [9] T. A. C. M. Claasen and W. F. G. Mecklenbrauker, "The wigner distribution - a tool for time-frequency analysis, part ii: Discrete-time signals," Philips Journal of Research, vol. 35, no. 4/5, pp. 276-300, 1980. [10] T. A. C. M. Claasen and W. F. G. Mecklenbrauker, "The wigner distribution - a tool for time-frequency analysis, part iii: Relations with other time-frequency signal transformations," Philips Journal of Research, vol. 35, no. 6, pp. 372-389, 1980. [11] L. Cohen, "Generalized phase-space distribution functions," Journal of Mathematical Physics, vol. 7, pp. 781-786, 1966. 123

[12] L. Cohen, "Time-frequency distributions - a review," Proceedings of the IEEE, vol. 77, no. 7, pp. 941-981, July 1989. [13] P. Flandrin, "Some features of time-frequency representations of multi-component signals," in IEEE International Conference on Acoustics, Speech and Signal Processing, volume 41B, pp. 4.1-4.4, 1984. [14] P. Flandrin and F. Hlawatsch, "Signal representations geometry and catastrophes in signal processing," in Mathematics in Signal Processing, T. Durrani, J. Abbiss, J. Hudson, R. Madan, J. Mcwhirter, and T. Moore, editors, pp. 3-14, Clarendon Press, Oxford, 1987. [15] J. D. Gaskill, Linear Systems, Fourier Transform, and Optics, John Wiley and Sons, 1978. [16] R. L. Haupt, "An introduction to genetic algorithms for electromagnetics," IEEE Antennas and Propagation Magazine, vol. 37, no. 2, pp. 7-15, April 1995. [17] S. Haykin, Communication Systems, John Wiley and Sons, second edition, 1983. [18] S. B. Hearon and M. G. Amin, "Minimum-variance time-frequency distribution kernels," IEEE Transactions on Signal Processing, vol. 43, no. 5, pp. 1258-1262, May 1995. [19] F. Hlawatsch, "Interference terms in the Wigner distribution," in Digital Signal Processing 84, V. Cappellini and A. Constantinides, editors, pp. 363-367, North Holland, 1984. [20] J. H. Holland, "Genetic algorithms," Scientific American, pp. 66-72, July 1992. [21] J. Jeong and W. J. Williams, "Alias-free generalized discrete-time time-frequency distributions," IEEE Transactions on Signal Processing, vol. 40, no. 7, pp. 2757-2765, November 1992. [22] J. Jeong and W. J. Williams, "Kernel design for reduced interference distributions," IEEE Transactions on Signal Processing, vol. 40, no. 2, pp. 402-412, February 1992. [23] H. Kim and H. Ling, "Wavelet analysis of radar echo from finite-size targets," IEEE Transactions on Antennas and Propagation, vol. 41, no. 2, pp. 200-207, February 1993. [24] E. F. Knott, J. F. Schaeffer, and M. T. Tuley, Radar Cross Section, Artech House, 1985. [25] H. Ling, J. Moore, D. Bouche, and V. Saavedra, "Time-frequency analysis of backscattered data from a coated strip with a gap," IEEE Transactions on Antennas and Propagation, vol. 41, no. 8, pp. 1147-1150, August; 1993. [26] E. K. Miller, "A study of target identification using electromagnetic poles," Technical Report UCRL-52685, Lawrence Livermore Laboratory, March 1979. [27] J. Minkoff, Signals, Noise, and Active Sensors: Radar, Sonar, Laser Radar, Wiley Interscience, 1992. 124

[28] A. Moghaddar and E. K. Walton, "Time-frequency distribution analysis of scattering from waveguide cavities," IEEE Transactions on Antennas and Propagation, vol. AP41, no. 5, pp. 677-680, May 1993. [29] C. J. MCCormack, V. V. Liepa, and W. J. Williams, "Time-frequency analysis of radar target backscatter," in Advanced Signal Processing Algorithms, F. T. Luk, editor, volume 2563, July 1995. [30] J. C. O'Neill and W. J. Williams, "New properties for discrete, bilinear time-frequency distributions," in Proceedings of the IEEE-SP International Symposium on TimeFrequency and Time-Scale Analysis, pp. 505-508, 1996. [31] A. V. Oppenheim and R. W. Schafer, Discrete-Time Signal Processing, Prentice-Hall, 1989. [32] B. Porat, Digital Processing of Random Signals, Theory and Methods, Prentice-Hall, 1994. [33] D. M. Pozar, Microwave Engineering, Addison-Wesley, 1990. [34] R. L. Riolo, "Survival of the fittest bits," Scientific American, pp. 114-116, July 1992. [35] R. L. Riolo. Introduction to Evolutionary Algorithms, March 1996. [36] M. I. Skolnik, Introduction to Radar Systems, McGraw-Hill, second edition, 1980. [37] J. P. Stephens, "Advances in signal processing technology for electronic warfare," Journal of Electronic Defense, pp. 41-50, September 1995. [38] H. L. Van Trees, Detection, Estimation, and Modulation Theory - Part I, John Wiley and Sons, 1968. [39] J. Ville, "Theorie et applications de la notion de signal analytique," Cables et Transmissions, vol. 20A, pp. 61-74, 1948. [40] W. A. Watkins, "The harmonic interval: Fact or artifact in spectral analysis of pulse trains," Marine Bio-acoustics, vol. 2, pp. 15-43, 1966. [41] M. W. Whitt, F. T. Ulaby, and K. Sarabandi, "Polarimetric scatterometer systems and measurements," in Radar Polarimetry for Geoscience Applications, F. T. Ulaby and C. Elachi, editors, chapter 5, pp. 191-217, Artech House, 1990. [42] E. Wigner, "On the quantum correction for thermodynamic equilibrium," Physics Review, vol. 40, pp. 749-759, 1932. [43] W. J. Williams, "Biological applications and interpetations of time-frequency signal analysis," Proceedings of the IEEE, September 1996. [44] W. J. Williams, M. L. Brown, and A. C. H. III, "Uncertainty, information and timefrequency distributions," in Advanced Signal Processing Algorithms, Architectures and Implementations II, volume 1566, pp. 144-156, 1991. [45] W. J. Williams and T. Sang, "Adaptive RID kernels which minimize uncertainty," in Proceedings of the IEEE-SP International Symposium on Time-Frequency and TimeScale Analysis, pp. 96-99, 1994. 125

[46] W. J. Wiscornbe. MIEV Program and Documentation, August 1992. [47] P. M. Woodward, Probability and Information Theory, with Applications to Radar, -McGraw-Hill., 1953. 126