RL 695 DETERMINATION OF SEM POLES FROM FREQUENCY RESPONSES J. M. Pond and T.B.A. Senior Radiation Laboratory Department of Electrical and Computer Engineering The University of Michigan Ann Arbor, Michigan 48109 Abstract A curve fitting and pole extraction algorithm has been developed and applied to exact frequency domain data for the surface fields on a perfectly conducting sphere. The data are fitted extremely closely and for at least a handful of the lowest order SEM poles, the extracted poles and their residues are in good agreement with their known exact values. Unfortunately, this is not true if the frequency response is degraded in accuracy. In particular, noise effects are explored, and it is found that for noise levels typical of the best experimental data, it is no longer possible to locate more than (at most) the dominant SEM pole to a reasonable degree of accuracy. Acknowledgement This work was supported by Mission Research Corporation under Contract No. F29601-78-C-0082. RL-695 = RL-695

DETERMINATION OF SEM POLES FROM FREQUENCY RESPONSES 1. Introduction The singularity expansion method (SEM) is based on the analytic properties of the electromagnetic response of a body as a function of the complex frequency s. For a passive body the singularities are confined to the left half of the complex s plane, and a knowledge of these singularities can characterize the response to any excitation. If the body is finite and perfectly conducting, the only singularities in the finite part of the plane are simple poles, which occur in complex conjugate pairs [1,2]. It is fundamental to SEM that the poles are a property of the body alone. The pole locations (-but not the residues) are unaffected by a change in illumination, and if a collection of poles is extracted from computed or measured data for the response, the SEM poles can be distinguished from numerical artifacts by their positional invariance. Cataloging the true poles is therefore a simple method of summarizing information about a body, and their extraction from measured data could serve as a means of target identification. Several numerical algorithms exist for determining the SEM poles from frequency domain data. One of these is an iterative method developed by Sharpe and Roussi [3] and based on a technique of Levy [4]. It is essentially a least squares method that fits the data with a rational function. The iteration linearizes the calculation and also reduces the

-2 - excessive weighting of the higher frequencies that a straight least squares computation normally produces. The program was initially applied to measured data for the axial current at several locations on a thick cylinder (for which the SEM poles are unknown) over a frequency range spanning the first five longitudinal modes. In every instance the rational function obtained gave an excellent fit to the measured data curve, but as the illumination changed, all poles except the dominant one showed substantial movement in the complex plane. The extent to which the lack of success was due to the program itself, the selection of such parameters as the sampling interval and the order of the rational function, or to the noise and other inaccuracies in the measured data, was not apparent. In the time domain it is found [5] that pole extraction is quite sensitive to noise. To see if this same sensitivity exists in the frequency domain and gain experience in the application of the program, it is helpful to consider data whose accuracy can be controlled. The only finite body for which the frequency response, poles and residues can be easily obtained to any accuracy is the sphere. In the following sections the determination of the poles and residues from frequency domain data for the surface fields on a perfectly conducting sphere is discussed. After a brief description of the numerical algorithm and the computation of the exact surface fields, poles and residues (Section 2), the extraction of the poles and residues from the frequency response data is discussed (Section 3). In Section 4 we then consider the effect of noise and other data inaccuracies on the pole extraction process.

-3 - 2. Formulation Over any finite frequency range the electromagnetic response of a body can be approximated by a ration function whose poles can be found. It is assumed that a subset of these approximate the SEM poles which are dominant in this frequency range and can be distinguished b; their positional invariance to a change in excitation of the body. Given a (complex) frequency response F(wZ) where wo, Q = 1,2,...,L, are sampled (real) frequencies, the numerical algorithm employed fits this with a rational function m a + a w +...+ amW (4 -+. n a *(m < n). (1) l b + blw. + b n The initial curve fit is obtained when the error L E D(o)F(W) - N(Q) ) = -1 is minimized, subject to the constraint bo = l.O,by solving the simultaneous set of equations aaj 0 j =,,..., aj for the coefficients a. and b.. The square of the resulting 3 J denominator is then used as a weighting factor to improve the rational function fit, giving rise to an iterative procedure. At the kth stage of iteration, the coefficients are obtained by minimizing

-4 - E = j') - < k() Dk() (2) and so on until the error is less than a pre.specified value. A program has been written to implement this curve fitting routine. The program has three parameters: the orders of the numerator and denominator polynomials, M and N respectively, and the maximum allowed error which terminates the iteration, which must be chosen at the outset. At the conclusion of the program the poles and residues of the rational function approximation are computed. The process is then repeated using other (distinct) data for the response of the same body, and those poles which are common to most of the results are identified as SEM poles of the body. To optimize the process, it is helpful to consider data for a frequency response whose poles and residues are known precisely. This is true in particular for the surface field on a perfectly conducting sphere. A sphere of radius a is illuminated by the plane wave = x eiz/c Hi = -yY e1/c where Y is the intrinsic admittance of free space, c is the velocity of light in vacuo, and the time factor e it has been assumed and suppressed. In spherical coordinates the total tangential magnetic fields at r = a are [6]

- 5 - lie YTI c, sin 4) = T2 cost wi.th T (to a 01) nflz' in~1 2n + 1 n(n + I) ___ ___ ___ P~~(l Cos o) + ~ F~J2 )( a) DO n (IDa ) T w 30) C co a ni. n + l 2 n ~ _ _I_ _ _ __II _ 1 n(n + 1) W + (I-) (Cos a ) ( 2 ) (o asinl 0 P '( i)Cos a) n ( 3b) where P(l)(cos a) n is the Legendre function of degree n and order unity -as' defined by Stratton [771 and ( 2)X ~ n (Ix) Al h(2 )() n ( 2) ' (x dx ~h(2 ) (x) where h(2I)(X) is the spherical lankel function of the second kind of 01 order n.

-6 - By appropriate truncation of the infinite series representations, it is a simple matter to compute T and T to any desired accuracy. 1 2 T and T were computed for 6= 0(45)180~ and 0.2 < wa/c < 7.0 to 1 2 - six decimal accuracy. The functions (.2) Cx) and -2)(x) are proportional to polynomials in x of orders n and n+l respectively whose zeros are the SEM poles. In terms of the complex frequency s = iwa/c, all zeros lie in the left half plane, and those which do not lie on the negative real s axis occur in complex conjugate pairs. As shown, for example, by Martinez et al. 18], the zeros can be arranged in layers lying successively further from the imaginary s axis. When ordered from the right, the odd (even) numbered layers are the electric (magnetic) mode resonances produced by the zeros of ( 2) -is) and E 2)(-is) respectively. In general the dominant SEM poles are those in the first (z = 1) layer, and the nth pole numbered up from the negative real s axis is a zero of (2)' (-is). T and T can be expressed as 1 2 Rm(O) R(e) T (-i ) E= -, T (-is,e) = 2 (4) S - Si2 S - S 11-1 m= where the Sm are zeros of eithe (-is) or (-is), and R (o) and Rm(e) are tie residues of T and T respectively at s = s. 1 2 1 2 m

-7 - 3. Exact Data Analysis The curve fitting algorithm was applied to the computed data for T as a function of frequency in an attempt to extract 4 or 5 dominant 2 pole pairs with sufficient accuracy to allow us to determine the e dependence of their residues. Data were available for 0.2 < wa/c < 7.0 in increments of 0.1 and 0.02. Since the first layer poles within this frequency range are dominant in determining the response it was expected that they were the ones most likely to be located accurately. To apply the algorithm there are a number of parameters which must be chosen, some of which relate to the data and others to the curve fitting process. Regarding the data, there are the minimum, maximum and increments of wa/c and, in our case, the choice of phase reference for the frequency response. Since the lower order poles are of prime concern, it was natural to choose min wa/c to be the smallest value for which data was generated, namely 0.2; and to avoid handling more data than was clearly necessary, max wa/c = 4.0 was selected with increments of 0.1. The computed data of Section 2 are phasereferenced to a plane perpendicular to the z axis through the center of the sphere. For all e except ir/2, arg T varies almost linearly 2 as a function of frequency, and this translates into a roughly sinusoidal variation of the real and imaginary parts which are the inputs to the curve fitting process. This variation can be minimized, by choosing the point where the field is computed as the phase reference, which results in a numerically easier curve to fit. Once the curve fitting was accomplished and the poles and residues determined, the phase reference was returned to the original location.

-8 - Three parameters involved in the program itself are the orders M and N of the numerator and denominator polynomials and the maximum allowed error Enax Since the set of SEM poles is infinite in number and the response remains finite as ca/c --, it would seem that the accuracy of curve fit should increase with M and N, and that a logical choice w^ould be M = N. Numerically, however, problems are experienced when M and/or N are large whereas if N is small there are too few poles available to simulate the data. It was therefore anticipated that there would be an optimum range of It and, perhaps, M depending on the frequency span of the data and the particular characteristics of the computer. The error Emax relates to the convergence of the iterative process and is not directly a measure of the curve fit or the accuracy of pole extraction. WJhen running the program, Emax was -8 set at 10 and the iteration was terminated when this value was achieved, or after 20 iterations, whichever came first. In many instances the maximum allow1ed error was not obtained, but the curve fit was still excellent. As a measure of the curve fit, the mean square error L 2 - 1 (5) fit L- K - D.( ) 9J=1

-9 - was computed, (c.f. (2)), where the polynomials are those obtained from the final iteration. Due to the limited precision with which the data were stored, any value of Efit less than 0.25 x 10 7 was essentially zero. Since the curve fit was excellent in most cases, it was not unusual for this to occur. All of the initial runs were carried out for 0 = 0 (for which T = T ). It was found almost immediately that numerical difficulties 2 1 arise if N exceeds (about) 25, and for M = N, if N exceeds (about) 18. Ii, either instance the exponential range of the computer (Amdahl 470/V8) was exceeded. We therefore chose M < N, and because of the restriction on N, limited the frequency span of the data to wa/c < 4.0 to allow for a reasonable number of curve fitting poles in addition to the SEM poles that were sought. Since the frequency range considered contains 4 pole pairs in the first layer, an initial curve was obtained with M = 7 and N = 8. Although the curve fit was quite good and the iteration converged, the extracted poles were closer to the imaginary axis, by about a factor of two. Increasing M and N brought increasing accuracy of both the curve fit and pole location. For M = 11 and N = 12 the modulus (and phase) of the simulated response were graphically indistinguishable from the data within the frequency range spanned by the data, but as seen from Fig. 1, there are discrepancies outside the range. Table 1 compares the locations of the first four SEM poles with those of the extracted poles for 0 =, AOoa/c = 0.1 and four M and N combinations. In each case Efit 0 and the agreement in pole locations improves with increasing order of the polynomials. The best results are for M = 15 and N- 16 in the sense that a further increase in M and/or N gives no improvement. Similar comparisons for a variety

-10 - of polynomial orders from 8 to 18 have shown that the accuracy of the extracted poles is best for M = N-1 and diminishes for N b 20. For given M and N a decrease in the sampling interval from 0.1 to 0.02 has no appreciable effect. Shifting the phase reference of the data to a plane through the center of the sphere decreased the accuracy of the pole extraction as expected. An accurate curve fit does not imply a comparable accuracy in the extracted pole locations, and when the locations of the true poles are unknown, it is necessary to vary the illumination conditions, e.g., change e, and use the positional invariance of the true poles as the criterion of accuracy. We also comment that Efit is unrelated to Emax fit max and for the above cases the specified error 10-8 was never achieved prior to the completion of the allowed 20 iterations. The effect of changing e is shown in Table 2, which gives the extracted pole locations for e = 0(45)180~ with M = 15, N = 16 and Awa/c = 0.1. The accuracy does not vary significantly with e, though it can be seen that at 0 = 90~ where the first and third poles are not excited, the accuracy of the second and fourth poles is better than before. As an example, the residue R3(8) of the third pole for T is 2 2 plotted in Fig.. 2. The fourth pole is fairly close to the upper limit of the frequencies spanned by the data. To improve its accuracy and to locate the next pole or two, it is natural to increase max oa/c to 5 or 6. The best agreement is obtained with M = 17 and N = 18. Although the curve fit is again excellent, as it was for M = 15 and N = 16 with the smaller data set, the first three poles are not quite as accurately

-11 - located, but the fourth through sixth are in reasonable agreement. Unfortunately, to increase the data span still more and extract further poles requires the use of polynomials of higher order. An alternative approach is to retain the same span of data and 'window', i.e., shift the span to encompass those poles which are sought. This is illustrated in Table 3 for three different M and N combinations applied to the data for 2.0 < wa/c < 6.0. In terms of the accuracy of the extracted poles, the case M = 15 and N = 16 is best. The fourth through sixth poles are located more accurately than with the larger frequency span, but the first pole is not picked up at all, and the second is considerably in error. This is hardly surprising since the first two poles are no longer spanned by the data. As a result of the above investigation, the following conclusions can be drawn. The data should fully span the imaginary parts of the poles to be located. If n SEM poles are spanned, N should be in the range 3n to 4n with M = N-l, but N should not exceed (about) 25 to avoid numerical difficulties. This upper limit decreases with increasing max wa/c and is almost certainly machine dependent as well. For a greater span of data and/or to extract more than a handful of SEM poles, it may be necessary to process the data using frequency windows. 4. Effect of Noise In most practical applications of the pole extraction method, the data for the frequency response have been obtained by measurement or by the numerical solution of a less than perfect model of the scatterer. Inevitably such data are subject to noise and other uncertainties, and it is important to see how the accuracy of both the curve fit and the

-12 - SEM pole extraction are affected. For this purpose, two types of noise were considered: numerical inaccuracies in the form of data limited to k decimal places, and added Gaussian white noise of various amplitudes. For the first study, the real and imaginary parts of T2(0) which were originally accurate to six decimal places were rounded to k decimals with k progressively reduced. The data used spanned 0.2 < wa/c < 4.0 in increments of 0.1 and 0.02, and since a rational function with M = 15 and N = 16 had proved to be effective in the absence of noise (i.e., when k = 6), this function was chosen. Because the accuracy of the extracted poles was slightly better for wa/c = 0.02, this sample interval was used in all of the noise studies. As k was reduced down to 1, the extracted poles became increasingly inaccurate as shown in Table 4, and for k = 2 even the dominant pole was substantially in error; however, the curve fit remained good. As k decreases, each pole moves closer to the imaginary s axis, and the general behavior is similar to that found when fitting the exact data with rational functions of progressively lower order. This suggests that by increasing M and N we might be able to overcome some of the noise effects and thereby improve the accuracy of the extracted poles. Using 1M = 23 and N = 24 with data having k = 3 and 2 produces only a slight improvement, primarily for the data with k 2. In the second study Gaussian distributed white noise was added to the exact data for the real 'and imaginary parts of T (0), wa/c = 2 0.2(0.02)4.0. The noise, was produced by a random number generator for which the mean and standard deviation could be specified. In all cases the mean was chosen to be zero and the standard deviation varied

-13 - to change the noise level. For noise with standard deviations 10-5 and 10-4, the curve fitting and pole extraction results were comparable, for similar levels of data uncertainty, to those obtained when the data was rounded. As a final test the curve fitting and pole extraction algorithm was applied to measured data for the field component T at the front 2 (e = 0) of a metallic sphere 6 inches in diameter. The data were obtained in an anechoic chamber over the frequency range 0.118 to 4.4 GHz, corresponding to 0.2 < wa/c < 7.0, but only the data for 0.2 < wa/c < 4.0 were used. The results of pole extraction with a rational function having M = 15 and N = 16 are given in Table 5, and are comparable to those for Gaussian noise with a standard deviation of 10-4. Although only the magnitudes are shown, Fig. 3 demonstrates the fit obtained with measured data. 5. Conclusions Using rational functions a curve fitting and pole extraction algorithm has been developed and applied to exact frequency domain data for the surface fields on a sphere. The data are fitted extremely closely and for at least a handful of the lowest order SEM poles, the extracted poles and their residues are in good agreement witi their known values. It is possible that the method could be further refined

-14 - to yield a few more poles, but the performance is already close to the limits of the computer. Overall, the success is comparable to that achieved by Brittingham et al [9] using a frequency domain Prony's method. Unfortunately, the situation is very different if the frequency response is noisy or degraded in accuracy in a manner typical of measured data. Although the curve fit is still good, even a small amount of noise is sufficient to produce considerable discrepancies between the extracted and true (SEB) poles, and for noise levels characteristic of the best experimental data, it is impossible to locate more than (at most) tie dominant SEM pole to any degree of accuracy. Since these conclusions have been reached using only a single algorithm applied to the frequency response of a sphere alone, it cannot be inferred that all frequency domain methods will be equally affected by noise. Nevertheless, the algorithm is c reasonably sophisticated one and is highly successful in curve fitting the data. Likewise the sphere, though a less resonant structure than (say) a thin wire, has a frequency response which is not unlike that of more complicated targets such as an aircraft. For these reasons, it is probable that the conclusions have general validity.

-15 - References [1] Baum, C. E., "The singularity expansion method," in Transient Electromagnetic Fields (ed. L. B. Felsen), Springer-Verlag, Berlin, 1976. [2] Sancer, M. I., and A. D. Varvatsis, "Toward an increased understanding of the singularity expansion method," AFWL Interaction Note 398, May 1980. [3] Sharpe, C. B., and C. J. Roussi, "Equivalent circuit representation of radiation systems," AFWL Interaction Note 361, April 1979. [4] Levy, E. C., "Complex curve fitting," IRE Trans. Automatic Control, vol. AC-4, pp. 37-44, May 1959. [5] Cho, K. S., and J. T. Cordaro, "Calculation of the SEM parameters from the transient response of a thin wire," IEEE Trans. Antennas Propat., vol. AP-28, pp. 921-924, November 1980. [6] Bowman, J.J., T.B.A. Senior and P.L.E. Uslenghi, "Electromagnetic and acoustic scattering by simple shapes," North-Holland Publ. Co., Amsterdam, 1969. [7] Stratton, J. A., "Electromagnetic theory," McGraw-Hill Book Co., Inc., New York, 1941. [8] Martinez, J.P., Z. L. Pine and F. M. Tesche, "Numerical results of the singularity expansion method as applied to a plane wave incident on a perfectly conducting sphere," AFWL Interaction Note 112, May 1972. [9] Brittingham, J.N., E. K. Miller and J. L. Willows, "Pole extraction from real-frequency information," Proc. IEEE, vol. 68, pp. 263-273, February 1980.

Table 1: Comparison of exact and extracted pole locations for 0 = 0, Aca/c = 0.1 and various M,N Extracted Exact M=12, N=14 M=13, N=14 M=14, N=16 M=15, N=16,. I I -0.500000 -0.5010 -0.5007 -0.5007 -0.5000 ~iO.866025 +i0.8641 +iO.8659 +i0.8657 +iO.8658 -0.701964 -0.6964 -0.7052 -0.7027 -0.7025 +i 1.80740 ~il.808 il.804 +il.803 il.806 -0.842862 -0.8436 -0.8334 -0.8408 -0.8399 ~i2.75786 ~i2.739 ~i2.762 ~i2.766 +i.2.760 -0.954230 -1.056 -0.9309 -0.9199 -0.9439 ~i3.71478 ~i3.575 ~i3.642 i 3.685 ~i3.678

Table 2,: Comparison of exact and extracted pole locations for Aw~a/c = 0.1, M = 15, N = 16 and e = 0(45)1800 I Extracted Exact I0 & = 0 1= 45e = 900 e = 1350 f e = 1800 -0.500000 ~iO.866025 -0. 701 964 ~i I180 740 -0.842862 ~i 2 75 786 -0. 954 230 ~i 3 71 478 - 0. 5000 ~iO.8658 - 0. 7025 ~il1. 806 - 0.8399 ~i 2 760 -0. 9439 ~i 3 678 - 0. 5000 ~iO.8661 not excited -0.8418 ~i2.760 -0.9490 not excited -0. 7020 ~il1.807 not excited -0.9565 +i.3716 - 0.499 7 ~iO.8661 not excited -0.8419 ~i2.759 -0.9341 Zi3. 732 -0 5001 ti 0. 8660 -0. 7008 tili.808 -0.841 0 ~i 2 751 -0.931 0 ~i 3 622

Table 3: Comparison of exact and extracted pole locations for o = 0, and wa/c = 2.0(0.1)6.0 Extracted Exact M = 15, N = 16 M = 17, N = 18 M = 19, N = 20 -1. -0.500000 +i0.866025 -0.701964 +il.80740 -0.842862 ~i2.75786 -0.954230 ~i 3.71478 -1.04764 -i4.67641 -1.12891 -i5.64163 not located -0.6365 ~il.802 -0.8449 +i2.746 -0.9582 fi3.714 -1.041 +i4.656 -1.079 ~i5.467 not located -0.7005 +i 1.730 -0.8558 +i2.765 -0.9422 +i3.728 -1.011 -+i4.658 -1.054 +i5.497 not located -0.6025 il.816 -0.8536 ti2.730 -0.9787 +i3.738 -1.011 +i4.728 -0.9399 ti5.658 - -. I.

Table '; Comparison M = 15 and of exact and N = 16, with extracted poles for e = 0, wa/c = 0.2(0.02)4.0, data sets rounded to k decimal places Extracted Exact..,.._._._,_,_ ----.. k = 6 k = 4 k = 3 k = 2 k = 1 -0.500000 -0.5000 -0.5011 -0.5140 -0.4081 -0.3764 +i0.866025 +iO.8658 +iO.8639 +iO.8577 +i0.8369 +iO.9174 -0.701964 -0.7025 -0.6957 -0.7005 -0.4924 -0.2961 ~il.80740 1il.806 +il.811 +il.861 +il.643 +il.939 -0.842862 -0.8399 -0.8323 -0.7548 -0.4727 -0.1581 ti2.75786 +i2. 760 +i2.732 +i2.842 +i2.516 ~i2.964 -0.954230 -0.9439 -0.9614 -0.6875 -0.3981 -0.04094 +i3.71478 +i3.678 +i3.517 ti3.784 +i3.924 +i3.257 Efit <0.25 x 10-7 0.27 x 10-1 0.27 x 10-1 0.27 x 10-1 0.31 x 10-1 fitIII

Table 5: Comparison of exact and extracted pole locations for measured data at e = 0 with 0.2 < /a/c < 4.0, M = 15 and N = 16 Extracted Exact Experimental -0.500000 -0.3920 ~i0.866025 ~i0.8575 -0.701964 -0.3530 ~il.80740 ~il.927 -0.842862 -0.2764 ~i2.75786 i 3.066 -0.954230 not located i 3.71478 E 0.77 x 10-1 fit

3.0 2.0 E I 1.0 0.0 i 0.0 1.0 2.0 3.0 4.0 5.0 6.0 7.0 wC o c Fig. i: Comparison of T (0 ) ( —) and the curve fit ( ---) obtained with a/c 2(0)4 M 11 and N 12. with ma/c = 0.2(0.1)4.0, M = 11 and N = 12.

1.0 \ \ \! \ \\! R|(0) \! \ \/ \-1.0 I. I \ 0 90 O80 Q (degrees) Fig. 2: Real ( — ) and imaginary ( ---) parts of the residue, R3(e), of the third pole of the first layer and the extracted real (O O ) and imaginary (X X) parts.

3.0 - -- 2.0 - - - LL 1.0 0.0 I ____ 0.0 1.0 2.0 3.0 4.0 wa cQ C Fig. 3 Comparison of experimental data at o = 0 ( —) and the curve fit ( ----) obtained with 0.2 < oa/c < 4.0, M = 15 and N = 16.