Robust Iterative Technique for High-Resolution Spatial Processing and Spectral Estimation

1. Field of the Invention

The present invention relates to data processing techniques, and more particularly to a nonlinear technique for high-resolution data processing to diminish the effect of system corruption.

2. Description of the Prior Art

Prior data processing techniques, such as the maximum entropy and maximum likelihood methods, are described by J. P. Burg in his Ph. d. dissertation entitled "Maximum Entropy Spectral Analysis" at the Department of Geophysics, Stanford University, Stanford, Calif. in May 1975, and by J. Capon in the Proceedings of the IEEE, Vol. 57, pp. 1408-1418, August 1969, entitled "High-resolution Frequency Wave Number Spectrum Analysis", respectively. These approaches use the input time series data and perform mathematical operations on the power levels directly. The results are complex methods having long averaging times and high computational loads.

A general iterative technique is described by R. A. Wagstaff in the Journal of the Acoustic Society of America, Vol. 63, pp. 863-869, 1978, entitled "Iterative Technique for Ambient Noise Horizontal Directionality Estimation from Towed Line-Array Data". The described technique produces a deconvolved spatial spectra for one time period and one heading.

What is desired is a simple, fast technique that can produce high-resolution spectra from the output of a conventional processor, such as a beamformer or spectrum analyzer.

Accordingly, the present invention provides a nonlinear technique for high-resolution data processing which produces high-resolution power spectra from the output of a conventional processor. First, a "guess" of the spectral estimate is formed. Then the spectral estimate is convolved with the system response function to produce estimated output levels. The differences between the estimated output levels and the measured output levels are used to modify the spectral estimate to produce a new "guess". This process is repeated until certain criteria are satisfied.

Therefore, it is an object of the present invention to provide a simple and fast technique for obtaining high-resolution power spectra.

Another object of the present invention is to provide high-resolution power spectra from the output of a conventional processor.

Yet another object of the present invention is to provide a nonlinear technique for producing high-resolution power spectra which uses simple mathematical operations.

Still another object of the present invention is to provide a nonlinear technique for producing high-resolution power spectra which is robust, i.e., relatively insensitive to system and measurement errors.

Other objects, advantages and novel features of the present invention will be apparent from the following detailed description when read in conjunction with the appended claims and attached drawing.

FIG. 1 is a pictorial representation of the problem addressed by the present invention.

FIG. 2 is a graphic representation of the effect of a system upon the received acoustic field.

FIG. 3 is a pictorial representation of the nonlinear technique for high-resolution data processing according to the present invention.

FIG. 4 is a flow chart of the data processing technique according to the present invention.

FIG. 5 is a pictorial representational of the choice of an initial field estimate.

FIG. 6 consisting of (a)-(e) is a graphic illustration of the nonlinear technique for high-resolution data processing according to the present invention.

FIG. 7 is graphical depiction of field estimates at successive iteration steps for four lines of white noise.

FIG. 8 is a detailed graphical depiction of the final field estimate of FIG. 7.

FIG. 9 is a graphical illustration of the final field estimate for four "pairs" of plane waves.

FIG. 10, comprised of FIGS. 10(a) and 10(b), is a graphical illustration of the final field estimate for an idealized spectrum of three sinewaves and colored noise.

The basic philosophy of the nonlinear technique described herein as shown in FIGS. 1 and 2 is that a model of a source signal can be constructed from (1) the output data of a system and (2) knowledge of the system response function, which model will give a high-resolution spectral estimate of the source signal. The technique is described in terms of a sonar system which has a line array with outputs in the form of a set of measured beam levels from a beamformer. The system response function required is a set of beam patterns for the array. These two sets are all that are required by this technique to estimate the spatial spectrum of the acoustic field being observed.

Referring now to FIGS. 3 and 4 the procedure by which the present technique generates a high-resolution spectral estimate is as follows. First, a guess or estimate of the spectrum is made. This "guess" can be nearly any smooth function including a flat spectrum, but is generally a normalized copy of the initial measured response output as shown in FIG. 5. This spectral estimate is convolved with the beam patterns of the array to get estimated beam levels. For most applications it is sufficient to use only the mainlobe beam pattern response. Second, these estimated beam levels are compared with the measured beam levels from the beamformer. Third, the differences in the estimated and measured beam levels are used to modify the spectral estimate within the coverage area (generally only the area covered by the mainlobes) of the beams to generate a new "guess". This process is repeated until certain criteria are satisfied. Then the iterations stop and the result is a high-resolution spectral estimate of the acoustic field.

The iterative process can be viewed as a computerized technique for the gradual molding of a lump of clay (the spectral estimates). By successively using the beam differences, which change with each iteration, the peaks and troughs in the clay can be gradually broadened or narrowed to reduce the differences to zero or to an acceptable level. When this is achieved it can be said that the system cannot tell the difference between the spectral estimate and the measured spectrum.

The spectral estimate is generated from the previous one as explained below where

u=independent variable

i=beam number

I=total number of system outputs (or beams)

j=iteration number

N.sub.j (u)=power level of jth spectral estimate

B.sub.i (u)=beamformer power level response function of ith beam (usually .ltoreq.O with O on beam boresight)

R.sub.i =beamformed output level of ith beam

R.sub.i,j =estimate of ith beamformed output level obtained from jth spectral estimate ##EQU1## .DELTA.R.sub.,j =R.sub.i,j -R.sub.i C.sub.j (u)=accumulation of changes .DELTA.'.sub.i,j (u) in level plus a pedestal (P) to be made to jth spectral estimate N.sub.j (u)

.DELTA..sub.i,j (u)=part of the function N.sub.j +B.sub.i above the threshold t.sub.a or t.sub.s

.DELTA.1.sub.i,j (u)=.DELTA..sub.i,j (u)(-.DELTA.R.sub.i,j)/T where T=t.sub.a or t.sub.s as appropriate

P=a pedestal equal to 0,5

A.sub.j (u)=step function equal to P plus number of times a particular N.sub.j (u) is to be changed (A.sub.j (u).gtoreq.P)

A.sub.j (u)=A.sub.j (u) smoothed over the variable u (A.sub.j (u).gtoreq.P)

t.sub.a =threshold for adding decibels to N.sub.j (u) (6 dB is representative)

t.sub.s =threshold for subtracting decibels from N.sub.j (u) (10 dB is representative)

The algorithm can be summarized by the recursive relation ##EQU2## with stopping criteria that j.ltoreq.K, an arbitrary number such as 25, or where the root means square error .sigma..sub.j is minimum or some arbitrary low value such as 0.001 dB where ##EQU3## In the case where the outputs R.sub.i are from a fast Fourier transform (FFT) beamformer, the calculation of .sigma..sub.j is over only the beams in real acoustic

FIG. 6 contains several plots to illustrate the present technique. Plot (a) contains a solid curve for the jth spectral estimate N.sub.j (j could be zero for the first spectral estimate, or "guess") and a dashed curve for the j+1 spectral estimate. The circles above N.sub.j are the beamformer output levels R.sub.i. The +'s distributed along the u axis about 0 dB are the errors or differences .DELTA.R.sub.i,j which indicate the quality of the estimated spectrum. In this case N.sub.j needs improvement. For this example the modification procedure will be illustrated with beams i and i+1. The modification of the other beams will be omitted in the discussion, but it is understood that they go through the same process during one complete iteration. The order in which beams are considered is not important. They could be considered simultaneously using vector instructions to increase computation speed. It is only necessary that the functions C.sub.j and A.sub.j contain the accumulated results from all beams within the spectral range for which resolution is being improved.

The solid curves of the plots 6(b) and 6(c) are plots of N.sub.J +B.sub.i and N.sub.j +B.sub.i+1, the parts of the spectral estimate which contribute power to the spectral areas covered by beams i and i+1. Since they are of N.sub.j +B.sub.i and not just N.sub.j in the beam coverage areas, they give the relative contribution of N.sub.j to the beam power in the individual beams. The relative contribution at a given u determines the shape of these curves. The highest level contributes the most. This point could occur at almost any location in the beam coverage area, not just near beam boresight.

The first difference considered, .DELTA.R.sub.i,j, is positive. This indicates that the spectral estimate N.sub.j is too high within the beam coverage are a of beam i and needs to be reduced. This is accomplished by normalizing .DELTA..sub.i,j, the shaded area above the threshold t.sub.s, by multiplying it by the ratio-R.sub.i,j /t.sub.s to form .DELTA.'.sub.i,j and accumulating in C.sub.j (plot 6(e)). By this means the modification of N.sub.j is done at the locations in u where N.sub.j +B.sub.i is above the threshold, not where the values of N.sub.j are the greatest within the beam coverage area. This is done because fewer decibels need to be added or removed to affect .DELTA.R.sub.i,j and a better fit is obtained. These are nonlinear (in power) operations that are responsible for some of the attractive features of this technique.

Plot 6(d) contains the functions A.sub.j and A.sub.j. The step function A.sub.j is an accumulator that contains P, which is 0.5, for this illustration, plus the number of times N.sub.j is modified at a given u. The function A.sub.j is simply the function A.sub.j passed through a smoother.

Plot 6(c) contains N.sub.j +B.sub.i+1. In this case the difference R.sub.i=1,j is negative which indicates decibels must be added to N.sub.j within the beam coverage area of beam i+1. Therefore, .DELTA.'i+1j will be stored in the accumulator C.sub.j and the accumulator A.sub.j will be incremented one unit for all u covered by .DELTA..sub.i+1,j.

The values of t.sub.a and t.sub.s have been empirically chosen to be 6 and 10 dB, respectively. The reason for the difference between the two is that spectral lines are narrow and are produced in the present technique by addition. No matter how narrow the lines are, they are still detectable by the sonar system provided they have sufficient level. Hence, the range in u for addition should be narrower than for substraction. However, troughs or holes in the spectral estimate must be broad or they cannot be justified by the data. a very narrow, low level hole in the spectrum is not detectable in the system output and should not be accepted in the spectral estimate. Thus, the subtraction of energy should be carried out over a wider spectral band than the addition to avoid creating "invisible" holes in the spectrum.

The complete spectrum is obtained by applying the above-described procedure to all other beams.

Finally, the step function for A.sub.j is smoothed to produce A.sub.j, shown as the dashed curve in FIG. 6(d). This smoothing is necessary to keep the final correction function C.sub.j /A.sub.j without discontinuities and N.sub.j+1 reasonably smooth. N.sub.j+1 is obtained by adding C.sub.j /A.sub.j to N.sub.j as illustrated by the dashed curve in plot 6(a) which also contains N.sub.j. The differences between the two estimates (N.sub.j and N.sub.j+1) result from one iteration of the present technique.

FIG. 6 is for illustration purposes and the values are not to scale. Further, the functions A.sub.j, A.sub.j, C.sub.j and N.sub.j cannot be obtained for a given range in u without using all B.sub.i which contain any u within that range. For example, if results for only a small angular section are desired, as few as five out of 64 beams could be used in the spectral estimate for the range of u of interest (assuming B.sub.i is limited to the mainlobe response). This permits reducing the computational load and "zooming in" on interesting subsets.

The following algorithm of the present technique has been coded in FORTRAN for a minicomputer to use 64 or 128 beam (I) levels to estimate the acoustic field (N) on either 1024 or 512 points (u). ##SPC1##

The following examples used 64 conventional beam levels obtained by FFT from 40 hydrophones with Hann weighting and zero padding. FIG. 7 illustrates the present technique by showing successive field estimates, N.sub.j, as the iteration process progresses to a final estimate N.sub.14. For this example the true field consists of broadband (or spatially) white noise plus four uniformly spaced plane waves having signal-to-noise ratios of -10, 0, 10 and 30 dB on a single hydrophone. The initial estimate N.sub.o at the top progresses to N.sub.14 at the bottom which is shown in greater detail in FIG. 8.

The dashed vertical lines in FIG. 8 are the original discrete sources and the dashed horizontal line is the original white noise. The dots are the received beam levels R from the spatial FFT. Note the similarity between the shape of the curve which would result from connecting these dots and the initial estimate N.sub.o in FIG. 7. The pulses distributed about 0 dB are the differences .DELTA.R.sub.i,j. The residual error is not zero since the present technique is an approximate method, but the mean-squares error is about 0.25 dB and the largest individual error is slightly above 2 dB. Another characteristic of the present technique is the accuracy in locating peaks in the field, which accuracy increases with the signal-to-noise ratios of those peaks.

FIG. 9 presents an example where the plane waves of FIGS. 7 and 8 have been replaced with "pairs" of plane waves. The spectral resolution limits of the present technique are shown by the separation in angle of each pair which varies with the signal-to-noise ratio. The separations are 1.25 BW, 1 BW, 0.5 BW and 0.375 BW for the respective pairs where BW is the 3 dB beamwidth of the FFT beamformer. The dB levels of the field estimate N and the beam levels R are stored as integer numbers in 0.1 dB steps. As in the previous example, the residual error is due to the approximate character of the present technique with an rms dB error of 0.4 dB.

FIG. 10 compares the present algorithm with an idealized spectrum (FIG. 10(a) composed of three sinewaves and colored noise which corresponds to an infinitely long signal. To obtain the final spectral estimate (FIG. 10(b) the 64 points were first analyzed with a conventional FFT. Sixteen successive overlapping FFTs of 48 points, weighted and padded with 16 zeros, were dB-averaged to obtain the FFT output. This output was input to the present algorithm. The result is a final spectral estimate giving correct positions and levels for the sinewaves, and a reasonably good estimate of the colored noise part of the spectrum, without creating "false peaks" in the noise region.

In many cases the system response function is not perfectly known due to sensitivity change, sensor malfunction without warning, etc. which modify the system response function. The mainlobe response is least affected while sidelobes are drastically affected. But the present techniques use only the mainlobe response and therefore is relatively insensitive to these errors. The present technique can be used with little or no averaging, i.e., on each acquistion, and thus can be linked to the output of a standard beamformer to provide high-resolution spatial spectra.

The present technique requires more elementary operations and more storage than other techniques for small arrays, but is very attractive for large arrays. This is because the operations of the present algorithm are very simple and grow in proportion to the square of the array size, not in proportion to the cube as do most other methods. The memory size grows only in proportion to the number of beams. Even for "small sizes" the large number of the operations is compensated for by the simplicity of the operations; all computations are done on dB levels (except calculating beam levels) and require low precision (10 bits integer additions, subtractions, and comparisions). Even for multidimensional spectral analysis the present technique is still practical, while most other methods need too much storage and calculations to be feasible. The technique can be added to any conventional beamformer system as an add-on "black box" to improve the conventional picture in its totality or locally (zoom), and easily can be expanded to solve problems of acoustic "tomography" and ambiguity resolution.

Thus, the present invention provides a non-linear technique for high-resolution data processing by using a system output and knowledge of the system response function to determine the spectrum of a signal. The mathematics are simplified by operating on decibel quantities in the output domain, consisting of simple additions and subtractions of decibels with the exception of adding the power in N.sub.j +B.sub.i to get the level of the ith beam R.sub.i. The technique is relatively insensitive to errors in the output data and knowledge of the system response functions, and can work on a single acquisition to provide virtually realtime analysis.