NUM"1II CA\L METHI1OD)S 01 No SI, FORI(1I r NKIk(lINV 1)t)>1\IN K David A. Ksienski and Thomas M. Willis* Radiation Labl)oratory Department of Electrical Engineering and Computer Science The University of Michigan Ann Arbor, Michigan 48109 ABSTRACT This paper addresses the numerical problem of performing a pole residue expIansion on noisy data. A method of combining dissimilar data sets to iaclicve an effective increase in signal-to-noise ratio is tested on computer Sc encrated data simulating scattering from a sphere and from a thin wire. The dat.l is corruI)ted with wlite Caussian noise 1and Lthe allgorithlim is tested for signal-to-noise ratios rangingl from -10 dB to 60 dB. An iterative version of ttli algorithm is also tested. Finally, tthe problemi of filtering noise from relatively quiet data is discussed and a novel filtering algorithm is presented..I NTR(OD)UCT ION This paper addresses the numerical problem of performing a pole residue expansion on noisy data. Previous investigators [1,2] have shown that even relatively low noise levels severely diminish the accuracy of the pole residue expansion. In a recent paper by Ksienski [2], a method of combining dissimilar daita sets to achieve an effective increase in signal-to-noise ratio was presented and shown to be effective for sphere data at relatively low noise levels. In the' present paper the method is tested on sphere data at much hlilth r noise levels as well as on data associated with scattering from a thin wire. An iterative technique suggested in Ksienski [2] is also tested. Finlillv, the problem of. filtering noise from relatively quiet data is discussed and a novel filtering algorithm is presented. 2. [MULTIPILE DATA SETS 2.1 Formulation. Before presenting the numerical results associated with comltining dissimilar data sets, a brief summary of the theoretical foundation is presented. A more detailed presentation is given by Ksienski [2]. The datai will be represented using the standard pole residue expansion M = km ak F(j)- ++ N(l) (1) k ( -s + - s m=1 AV& Bell Laboratories Scholar. AT.',XI' Bell Labor.ntories Sclolarr. RL-781 = RL-781

where k is an index to the various lmel'asurem-ents and Ni (j( ) represents noise and clutter. The representation in terms of conjugate pole pairs is used to explicitly show that Fk(jo) is an even function of the circular frequency u alnd is associated with a real function of time. We wish to combine K such data sets in a manner so as to have a maximum signal-to-noise ratio in the composite data set, where the noise is assumed to have zero mean and be independently and identically distributed. Assuming a linear combination with arbitrary complex weighting coefficients, the composite data set is K = w kFk k I k F (ji) c Ci)Ilp K >Zw = Wk k=l M /I f ak. m J 1 - s m= l a km j1 - s m ) + Nk.(i)) (2) ~and by interchanging the order of summationi and noting the poles s with measurement number k, m invariance of the F (j,.) comp M c m m = > + ---- - Sm - m= 1 ) + N(jk) (3) s m wlet re K m k= K k akm wkak k kin (4) In general, it will only be possible to maximize one b or c at a time. This m m is because for different choices of m the akm will not vary in unison with k. lThus, thle increased signal-to-noise ratio, and hence increased accuracy, will only be obtained for the pole s or s. However, since the specification of m Il the pole is arbitrary, several poles may be obtained through successively emIphasizing different b or c. For the present, we restrict ourselves to m m obtaining an improved estimate of s, which necessitates maximizing lb 2. Froml (4), we may write 1 1 I? 1 t r Iw a11 ' (6 ) T T whlere =w [w,,... ] and a = [a,a,...,a 1 2 I 11 12 K1 igenerality, the weighting vector may be constrained to noti ing the (lot product formulation of (6), thie lmaxilnuil Without loss in -T satisfy w w = 1. b will occur for 1 Then., W = 777PT (7)

;and for tilis Loptimal Clcoit o1 tl' wei!ll i, ve'tor, 1) v = '.; I'o I evaluate the energy associated with the noise N(jo), we take the expectation of the noise squared K 2 E| N(j) 2 = E wNk(j) k=l land assuming the variance of the noise is equal to o2 in each of the ite. ti s ucle n t s F (j F ) K EIN(ji) 2 = <- E w 2 = o2 k=l Tlthus the energy associated with the noise will not increase with w normalized as in (7). If K data sets are combined each containing an approximately equal excitlation of the mode associated with the desired pole, the signal-to-noise ratio will increase by about a factor of K. For any set of residues, the resulting weighting coefficients will produce the maximum possible increase in signal-to-noise ratio. However, this is dependent upon an accurate knowledge of the residues which is generally not available. Two solutions to this problem have been devised, and both focus upon increasing the accuracy of thle weighting vector w. First, the error due to variation in the estimate of the real part of the pole (which often constitutes a majority of the error) nlma be reduced by modifying (7). Assuming that the primary effect of a pole Iand its associated residue is in thie immediate vicinity of the pole, it is possible to get a first estimate of the residue wihich would have resulted had til pole been constrained to hlave a different real part. Since for a fixed residue the component of the signal at a point on the j axis closest to the pole is inversely proportional to the real part of the pole, an improved weilgting vector may be obtained by normalizing each element of the weighting vector (i.e., the residue) by the real part of the pole associated with each residue. Thus (7) is used with a redefined as a) ' Re(s ) Re(SK 11 21 1 where si refers to the estimate of the pole s obtained from the ith data set. Th'l second method of improving w requires additional computation. Since the app)lication of the algorithm should result in improved residues (using the teclhnique described below), these improved residues mav be used to construct a 1new and presumably more accurate weighting vector. Thus, improvements may t.be made throuIh thle iterative applicatiol of the algorithm. 'lTe. results of licle pplication of the original algorittm as well as the iterative version are pre.sented in the next section.

After the composite data set is formed, it mutist be subjected to a pole and residue extraction procedure. The algorilthm due to Levy 3 and Sanathanan and Koerner [4] assumes that the poles and residues exhibit conjugate symmetry. This is generally appropriate as it is equivalent to requiring the frequency data to correspond to a real function of time. Untlfortunately, F (j) does not have conjugate symmetry. This may be seen comp from (3), (4), and (5), where for an arbitrary complex weighting coefficient wk, the resulting b and c will not be equal. This problem may be k' m m circumvented by using the algorithm derived in Ksienski [2], which does not force the residues to occur in conjugate pairs. The implementation of this algorithm requires negative frequency information which is obtained from the dat.a sets Fk(jw) by employing conjugate symmetry. After the accuracy of the pole location is improved, the next step is to compute the residues relative to the new pole estimate. Constraining a pole to a particular location is not an easy task [1], but the algorithm presented in Ksienski [2], by constraining the poles to occur in conjugate pairs without similarly restricting the residues, provides a simple alternative. Since the data set F (jo,:) corresponds to a real function of time, its negative frequency values are given by k(j ) = F k(-j). A second set of data, call it I(j,), is constructed so tha1lt TL( j) = -l(-j) and thus corresponds to an imaginary function of time. I(j ) may then be added to F k(jui) and the sum expanded as a series of poles which are then constrained to occur in conjugate pairs. Noting that Fk(ja) hlis an even real part and an odd imaginary part while l(j;() has an odd real part and an even imaginary part, the pole residue expansion of F (jpj) is immiiediately separable into components due to F (Jyi) and I(j(). Specifically, the pole residue expansion of F k(j,) may be written as M (b + c )/2 (b + c )/2 < m m m m j I 2? 5tm- s j (0- s ^J J' m - s m m while the pole residue expansion of L(j',) is M (b - c )/2 (b - c )/2,m m m m mLO -- S -m j ( - m- sI Thus, the second data set may be used to constrain the poles of F (j:') witiyout corrupting the residues. k 2.2 Results. The procedure was tested on two objects, a sphere and a thin wire with length-to-diameter ratio of 50 to 1. For the case of the wire, poles were extracted from computer generated data associated with current measurement on the surface of a wire. The wire was divided along its length into 50 sections to facilitate analysis by the NEC program [5]. iThe wire was iltluminated with a plane wave electromagnetic field, with the axis of the wire par>llel to the direction of the electric field. The frequency of the il llinination was stepp)ed from ((aL/c'1) = 0.2 to 4.0 in steps of 0,1. The computttr simulated current mieasurements of the 5th, 10th, 15th, _Oth, and 25th Sc t'Lions were then individually subjected to a pole residue expansion. A pole

coiimolliiio Lo o a L l ive d ta slts: Was l ocILteLd atL (:-;IcI) - -0.08/ + iO0.8855. Tlle t's ti l; m t. of tlIe I)ole;s 1 (ltali'ned fr)in tll' I-ve tt' tIl 'ts g'1-re'e t o tLhins va;lltI to within 0.2 percent. Tllis value is also in relative agreement with that obtained by Tesche [6]. For the sphere, the pole residue expansion was performed on computer generated data corresponding to current probe measuremenits at 0 = 0, 10, 20, and 30 degrees, witlhi a/c = 0.2 to 4 by 0.1. The location of the primary pole for the sphere may be determined analytically [1], and lies at sa/c = -0.5 + i.866. Since the benefit of the algorithm may be attributable to.increasing the siginal-to-noise ratio, the amount of improvement which might be reasonably exptected can be determined by first examining tlle increase il accura:cy which occurs when the signal-to-noise ratio is directly varied. The data was corrupted with computer generated white, Gaussian noise at levels necessary to produce signal-to-noise ratios of -19 dB to 60 dB in steps of 2.0 dB. The signal-to-noise ratio is defined as the total energy in the original uncorrupted data sets, Fk(jy(), divided by the total expected noise energy contained in these data sets. The range of errors associated with the pole as extracted from the Fk(jwo) is shown in Fig. 1 for each noise level as falling between the two horizontal marks. The fact that the pole for the wire is much more, resonant than the pole for the sphere permitted the SEM to be performed at much higher noise levels for the wire than for the sphere. The errors plotted in Fig. 1 for signal-to-noise ratios less than 15 dB are those associated with the wire data, while those for signal-to-noise ratios greater than 10 dB are associated with the sphere data. For each noise level a composite data set, Fk(-j) was created using the above procedure, and the specific case of wire data with signal-to-noise ratio of 0 dB is shown in Fig. 2 as the dashed line. The pole was extracted from F (jo)) and the error comp associated witl this pole is marked with a triangle. A second composite data set was created using the estimates of thie poles and residues as obtained via thel first F (j1i) and the error associated withl the pole extracted from this comp second F ( jw) is marked with a circle. Finally, the error associated with coomp the pole as extracted from the original data set Fk(j.)), with the strongest excitation of the dominant mode (the 25th section for the wire and the 0 = 0 case for the sphere) is marked with an X. There are several important features which are exhibited in Fig. 1. The most noticeable is the apparent effect of tile pole resonance on the accuracy of extraction of the desired pole; the expansion performed on the wire data wittl signal-to-noise ratio of 10 dB produced more accurate results than sphere data with signal-to-noise ratio of 50 dB. Additionally, the rate of implrovement in accuracy of pole extraction with increasing signal-to-noise ratio is certainly not a smoothi function. For tile wire, there is a large jump at -2 dB and for the sphere there is a 15 dB region centered at 40 CBI where increasing the signal-to-noise ratio has a negligible effect upon the accuracy of tlie pole extraction. Given these peculiarities, the improvement resulting fromn extracting the pole from F (j0 ) produces fairly reasonable results. comp In regions where the slope is steep the improvement resulting frqm extracting thle pole from the composite data set produces very nice results (remn-mbering tlhlt tihe percent error is graphed on a logarithmic scale), while in regions

200 100 10 *n - i * 0 iT I 1 i I I ) SPHERE I I 0@ 0 oS 0 UI EL 0 0 I C I 0 WIRE i 0.1 I I I I I I I -_ -t0 0 10 20 30 S/N (dB) 40 50 60 Fi. 1: Vertical bars delimit range of errors for poles extracted from Fk(js ); X indicates error for pole extracted from Fk(j(i) with greatest S/N; triangle indicates error for pole extracted from first compositL data set; circle indicates error for pole extracted from second composite data set.

,t II % 9 L I I I I, I I I 3 1 I* I I I I! I 3 - I _C1T er mode nd te often is bttr than any of the oriial etiate. is with t exception of the region where ll of te inistin; l edstimates t omposle r in error by greset ter t 5 percet, s til genrally mkes twhre asscit slope is efi effectively z ero, the improvlm ent is nevliibr, wlen t llst sme wole, the polt e estiates o t racted from the composite data set is almost lwas hottr theanmple, tat ined drm te dta st wir dtl n stimate reaf te stole wexcittion ofis the desired mode, a"nd quite often is better thlan any of the orina estimates. iaccu'lralte' to 7 percent anid ano t1her estilmate which hias a greater than3 80 percent error. Not only does the algorithm work in this case, but the iterative

I))t 1 it It i l on 1tl ll it,( i, t1it. tlll i; 'ali tO ts e; l' t i III ip I) V (I ( es i lll e tal tilt residues to obtain an even better estimate for the desired plole. 3. FILTERING OF DATA 3.1 Background. While experimenting with the method of Levy, Sanlthanan, and Koerner [4], [5] (LSK algorithm) on sphere data, sparsely sampled data sets (39 points) and densely sampled data sets (189 points) with identical signal-to-noise ratios were found to yield poles and residues of comparable a;ccuracies. This observation is somnewhat disconcerting since thel dcnser data set hlas almost five times the information of the sparse one and yet the LSK algorithm can not extract any better results. This fact has encouraged research into the possibilities of filtering the data before using it to find poles and residues. Unfortunately the use of filtering to improve the alccluralc of the poles found by tile ISK algorithm is somewhat more difficult thL:.;l first expected. Indeed several clif-ferentlt tI.ctoLrs tend to make it so. First, even if a filter is successful in improving the signal-to-noise ratio by, for example, what would otherwise be a respectable 10 dB, there is no assurance that the locations of the poles will improve. The authors found that sphere data with signal-to-noise ratios of 20 dB, 30 dB, and 40 dB all Ihad comparable errors in tle locltion ao te ples. On)ly when the s ignal-t — noi-se ratio increased to 50 dB and 60 dB did the locations of the poles tend to improve. Another problem is that for sphere data very little noise is acccptable, as stated before the noise must be some 50 dB or 60 dB down from the signal for good results. This means that what ever filting scheme is used it must work on the noise without corrupting the signal. Ideally the filter should not alter the data at all in the limit as the noise tends to This places several limitations on the choice of the nature of the filter. For example, filters which exhibit a "phase shift" cannot be used since they will move the location of the resonant peaks in the data thus shifting the ima.inary part of the pole. Also standard zero phase low pass filters cannot be used' since the't round out peaks in the data thius changing the real part of the pole to be less resonant. Additionally various types of transform filters which suffer from aliasing problems are unacceptable since they also corrupt the signal. Despite these difficulties the authors have developed a filtering scheme which has had some success in improving the location of the poles for sphere datal. This scheme is a type of finite window moving aperture filter. However, whereas the standard linear moving aperture filter [7] has a fixed weighting function, this scheme is a nonlinear filter wliere the weighting function is allowed to change from data point to data point. The filtering process begins with discrete frequency domain data, H(k), consisting of a signal plus noise, S(k) and N(k). The real part and imaginary part of the data will be filtered separately so that the analysis will concern only real quantities. A weighting function, W (k) where integer k ranges from -M to M, is used to generate filtered data, H(k) as followsi: HI(m) = W (mk)t(m (+ *k) k = -,,

whiler e i tLe order o0 tLhe Iilte r. Thie ftil terle d liat is tlln used w i tll the LSK algorithm to find poles and residues. The following section is a derivation of the weigIhting function W (k). m 3.2 Filter Derivation. The assumptions about the noise that will be made are as follows: E[N(k)] = 0 E[N(k)N(j])] = 0 if k # j F[N -(k) ] = -2. ( -) (2) (;3) Further, thle following restrictions will apply to thle weighting function M k=-M 1 (4) W (k) = W (-k) m m (5) Thet weighting function will be an even function so that the peaks in the data will not shift in position. It is also assumed that the signal is varying slowly enough compared to the width of the window so that it can be accurately represented by a quadratic expansion, that is: S(m + k) = a + a k + a k; for -M < k < M. 0 1 2 - - (6) 1Of course the values of a, a, and a will depend on the choice of m. 0 1 2 With this filter the estimate of S(m) will be M a = W (k)[S(m + k) + N(m + k)] 0 k=- m k=-M ( 7) and so an error measure can be defined as: = (a - a )Y 0o M = a - W (k)[S([m + k) + N(m + k)] o -M k ---il I (8) Now under tlhe previously stated assumptions tlle expected value of menasure can be found the error M M FI 1 -a a - Wl (k<)S(m + k) + o o r = k- -k=-1 W (k )S ( + k) M + o > w i (k) Z-j m i k=-:

>1 E() = 2 L k.k' 4 (k)J 2 1 m -- k=-H J M + a > W' (k) mk k=-H = 4a2 > k2W (,k) 2 in k=l - M + 2c. k=l r M W2 (k) + 2 I - 2 E W (k) m m k=l (9) By the proper choice of the W (k) values it is desirable to minimize the m exteccted value of the error measure. Therefore the values of the weighting function must satisfy,: ( J ),w (- J) II1 8 I2.I k-W (k) + 4a2W (J) 4c 1 - 2 W (k) q m m m k=l - k= - M = 4oC[ (J) - 1] + w (k)[ 8o + SJa -k' m m 2 k=l for J = 1,2,...,M 0; (10) This leads to M linear equations which can be solved for the M unknowns W (1) m tlthoulh W (M). The values of W (0) can be found from relation (4). For exacllple, when M = 1. 2a2 + 02 2 ' 2 w (0) = m 2a2 3 2a -+- 3o M = 2 34 a2 + (2 W (0) = - m 5(14 a2 + 2) 2 W (1) m = _____ — 2 + 32 2a'- + 3o w (1) m 24 a2 + o -2 5(14 a2 + - ) -6 a2 + c ' -_ 2 5(14 a + o ) 2 w (2) m Notice that if a = 0, that is the signal has no curvature, the weightting function is uniform so that the most noise possible can be eliminated. Notice also if j" = 0 then the estimate of S(m) is correct to the extent that the daltai can be represented by a quadratic expansion. Therefore, to find tlhe optimum weightinl, funct ionl or a given size window the variance of the noise and the a coefficient for eachdata point.must be known. Usuallny the noise variance is at least approximately known,

however the a coefficients require some knowledge of the signal without noise. 2 Hence the values of these a coefficients can only be estimated from the signal 2 plus noise data. For this reason the filter will only work on data with a relatively high signal-to-noise ratio (> = 20 dB). The a coefficients are found by a best curve fit with the data over the 2 range of the window plus one point. The error between the curve' fit and data is defined as: M+l1 n = [H(m + k) - (a + a k + a k2 + a k3)]2 (11) 0 2 3 k=-M-l The coefficients a through a are then chosen to minimize n. This leads to 03 the following equation for the a coefficients 2 k2 * H(m + k) - 1 Z k2H(m + k) a = (12) 2 2 (2k)-21*2 k4 Using this relation and knowing the noise variance o2 weighting functions for each data point can be generated and the data filtered. 3.3 Results of Filtering. Figure 3 shows noisy sphere data (S/N = 40 dB), and the same data filtered with M = 1 and M = 2 filters. The filter was tested on various sphere data with noise. On sparsely sampled data (39 points) no improvement was made in the location of the first dominant pole for any signal-to-noisd ratio. On more densely sampled data (189 points) improvement was made when the signal-to-noise ratio was between 50 dB and 60 dB. With signal-to-noise ratios greater than 60 dB the unfiltered data yielded the correct result to one percent accuracy anyway and filtering did not improve it. However, with a 55 dB signal-to-noise ratio the average error of the first dominant pole was 7 percent. After filtering with M = 1 this average error fell to 3.5 percent and with M = 2 fell to 2.5 percent. It should be noted that similar improvements can be obtained with more traditional filters, such as a uniformly weighted moving aperture, for this data. However a uniform weighted filter actually increases errors with other types of data sets, whereas our filter will not. For example, the unfiltered error of the 39 point data set with 55 dB signal-to-noise ratio averaged 7 percent, but after filtering with a uniform aperture with M = 1 the error went to 8 percent and with M = 2 climbed to 31 percent. Our filter however, had an average error of still 7 percent for both M = 1 and M = 2 filters. In conclusion the presented filter can improve the location of poles found with the LSK algorithm if both the density of the sampled data and the initial signal-tonoise ratio are sufficient. In addition this filter will not statistically increase errors of the poles and residues where more traditional'filters may.

25 2.0 -150 1 -162 o U, rq Lu 0 3 Ia X 1.5 ID1 -175 '187 0.5 0 0 2 3 4 FREQUENCY Fig. 3(a): 2.5 2.0 Sphere data filtering. (Unfiltered data) - -150 i -IIO - -162 - -175 I - -187: t! 1.5 w a I-.J a. ]! 1.0 0.5 FREQUENCY Fig. 3(b): Sphere data filtering. I (Filtered data M =;1)

2.s r -ISO - -150 2.0.5 _ -162 ~ / \ i J 1.0 ' a - '175 0.5 1-e7 0 I 0 I 2 3 4 FREQUENCY Fig. 3(c): Sphere data filtering. (Filtered data M = 2) ACKNOWLEDGEMENT This work was supported by the Air Force Weapons Laboratory under Contract No. F29601-82-K-0024.. REFERENCES 1. Pond, J. M. and T.B.A. Senior (1982), "Determination of SEM poles from frequency responses", Electromagnetics, Vol. 2, 55-67. 2. Ksienski, D. A. (1984), "Pole and residue extraction from measured data in the frequency domain using multiple data sets", Radio Science, Vol. 20, 13-19. 3. Levy, E. C. (1959), "Complex curve fitting", IRE Trans. on Automatic Control, Vol. AC-4, 37-44. 4. Sanathanan, C. K. and J. Koerner (1963), "Transfer function synthesis as a ratio of two complex polynomials", IEEE Trans. on Automatic Control, Vol. AC-8, 56-58. 5. Burke, G. J. and A. J. Poggio, (1977), Interaction Note 363, "Numerical Electromagnetic Code (NEC)", Air Force Weapons Laboratory, Albuquerque, N.M. 6. Tesche, F. M. (1972), "On the singularity expansion method ah applied to Electromagnetic scattering from thin-wires", Interaction Note 102, Air Force Weapons Laboratory, Albuquerque, N.M. 7. Oppenheim, A. B. and R. W. Schafer (1975), Digital Signal Processing, pp. 237-269, Prentice-Hall, Englewood Cliffs, N.Y.