Abstract:
A method and apparatus is provided for suppressing noise which includes receiving a plurality of signals from an array of sensors and transforming each of these signals to the frequency domain. The transformed signals are beamformed so that noise sources can be identified by bearing and frequency range. A planewave-fit noise-suppression routine is then used to remove identified noise sources from the transformed signals and to provide a noise-suppressed transformed signal having signals from said identified noise sources suppressed.

Description:
STATEMENT OF GOVERNMENT INTEREST 
   The invention described herein may be manufactured and used by or for the Government of the United States of America for Governmental purposes without the payment of any royalties thereon or therefore. 

   BACKGROUND OF THE INVENTION 
   (1) Field of the Invention 
   The present invention generally relates to sonar systems and more specifically to sonar systems particularly adapted for identifying the location of an underwater object(s). 
   (2) Description of the Prior Art 
   Conventional passive sonar systems detect acoustic signals emanating from an underwater object; that is, any device that moves through the water while emitting acoustic signals that sonar can detect. Torpedoes and submarines are examples of such underwater objects. 
   As modern, very quiet submarine platforms become operational in large numbers, new methods of detecting very low level signals from these quiet submarine platforms are desired, especially in the presence of high noise levels from surface shipping, wind biologics, and other sources of ambient noise. 
   Ambient noise from all these sources can be very loud, especially in coastal waters where surface ship densities are largest, and these loud noise signals can prevent the detection of weaker signals of interest. It is therefore, highly desirable to suppress the unwanted noises in order to better detect and track the signal(s) of interest. Noise must be suppressed coherently by estimating both the phases and amplitudes of the noise sources desired to be suppressed. This can be accomplished prior to the beamforming process, during the beamforming process, or in the post-beamforming process. Prior to beamforming, digital time series data sets from the hydrophone array are decimated, filtered, and heterodyned to prepare the multiple time series data sets to be transformed to the frequency domain utilizing a Fast Fourier Transform routine. Fourier transformed frequency domain data sets, commonly referred to as “hydrophone FFTs,” are then beamformed by using any one of numerous beamforming algorithms. A particular beamforming algorithm, called the Fourier Integral Method (FIM), was previously patented by the inventors in U.S. Pat. No. 5,481,505 which is incorporated by reference herein. 
   There have been several prior art methods developed to solve the sonar problem of reducing noise from a loud, near-surface noise source while maintaining the signal level of signals produced by the target of interest (TOI). As used herein, the phrases “near-surface noise source” or “near-surface source” refer to an object (e.g., ship) that is primarily located on or near the ocean surface. An intensive effort has been directed to the area of adaptive beamforming, as evidenced by the development of the well known minimum variance distortion response (MVDR) algorithms. For ideal ocean conditions, when the spatial coherence of the acoustic field is known exactly, MVDR algorithms are optimum in minimizing the total noise field while maintaining the target of interest&#39;s signal level constant. However, there is only a finite time to estimate the acoustic field spatial coherence. Furthermore, errors between the actual and estimated acoustic field spatial coherence degrade the performance of MVDR algorithms rapidly because MVDR algorithms are highly non-linear. MVDR algorithms require the calculation of the inverse matrix of the acoustic field spatial coherence spectral matrix (CSM). Small errors in the estimate of CSM can propagate to become very large errors in the estimate of the inverse matrix of CSM. The CSM is defined as the matrix of all cross-product pairs of individual hydrophone time series Fast Fourier Transforms (FFTs). The CSM is described in detail in commonly owned U.S. Patent No. 5,481,505. Therefore, MVDR algorithms are not robust in realistic open ocean environments, and are severely degraded when short averaging times must be used in tactical sonar systems. 
   A second class of prior art algorithms developed to address the aforementioned problem is referred to as the WHISPR family of processing algorithms. The WHISPR-related algorithms are relatively large. These algorithms rely on one physical principle: the acoustic time series of a near-surface noise source has a significantly greater time variance than the acoustic time series from a submerged target of interest due to the Lloyd&#39;s Mirror effect and several other causes. Although WHISPR has shown some promise on selected acoustic data sets, it has never been developed into a real time system because it is not robust in real ocean environments. 
   U.S. Pat. No. 6,424,596 is directed to a method and system for significantly reducing the acoustic noise from near-surface sources using an array processing technique that utilizes multiple signal classification (MUSIC) beamforming and the Lloyd&#39;s Mirror interference pattern at very low frequencies. Noise from nearby near-surface sources, such as merchant ships, super tankers, fishing trawlers, seismic profiling platforms, or other sources near the ocean surface can significantly interfere with the detection and tracking of a quiet target-of-interest (TOI) located well below the ocean surface. The &#39;596 invention reduces the noise of the near-surface sources, without degrading the signal level and quality of the TOI, by utilizing a unique application of the MUSIC beamforming process to separate the noise and signal subspace. Next, eigenvalue beamforming is used to reduce narrowband energy in selected frequency bins wherein the near-surface noise is radiating. Next, predetermined frequency and magnitude variance parameters are used to eliminate broadband noise emanating from the near-surface sources. 
   However, these prior inventions do not allow a sonar system to selectively reduce sound emanating from an object at a known azimuth and bearing in order to enhance sound emanating from an object of interest. 
   SUMMARY OF THE INVENTION 
   Therefore, it is an object of the present invention to provide a system that suppresses the noise from underwater noise sources of all types. 
   Another object of the present invention is to improve detection of signals present in noisy backgrounds. 
   In addition, it is an object of the present invention to provide a system having improved detection of signals of interest emanating from submerged objects. 
   Accordingly, the current invention provides a method and apparatus for suppressing identified noise sources. This includes receiving a plurality of signals from an array of sensors and transforming each of these signals to the frequency domain. The transformed signals are beamformed so that noise sources can be identified by bearing and frequency range. A planewave-fit noise-suppression routine is then used to remove identified noise sources from the transformed signals and to provide a noise-suppressed transformed signal having signals from said identified noise sources suppressed. 

   
     BRIEF DESCRIPTION OF THE DRAWINGS 
     The appended claims particularly point out and distinctly claim the subject matter of this invention. The various objects advantages and novel features of this invention will be more fully apparent from a reading of the following detailed description in conjunction with the accompanying drawings in which like reference numerals refer to like parts, and in which: 
       FIG. 1  is a block diagram of the noise suppression system of the current system; and 
       FIG. 2  is a graph showing maximization of the response function for a given bearing vector. 
   

   DESCRIPTION OF THE PREFERRED EMBODIMENT 
   The present invention is an extension of the invention described in U.S. Pat. No. 5,481,505, which is herein incorporated by reference. The present invention utilizes tracker inputs from the “Eight Nearest Neighbor Peak Picker (ENNPP)” which is completely described in U.S. Patent 5,481,505. For the sake of brevity, the details of the operation of the system of U.S. Pat. No. 5,481,505, although applicable to the present invention, are not repeated herein, but rather are referenced as needed. While the current invention was conceived as an improvement to U.S. Pat. No. 5,481,505, it should be evident that it is applicable to any system providing bearing and azimuth inputs to the algorithm. 
   In general, the technique of the present invention determines the phases and amplitudes of the noise sources to be suppressed by modeling the noise sources as planewave time series whose amplitudes and phases are determined from the measured time series input and the bearing-time and frequency time history of the noise sources to be suppressed. Detection is enhanced for passive sonar systems after noise suppression, especially if the Fourier Integral Method (FIM) algorithm is used as the beamforming method. The present invention may be furthered described with reference to  FIG. 1 . 
   The apparatus  10  shown in  FIG. 1  includes a horizontal hydrophone line array  12  having M hydrophones that receives acoustic signals in the water for all potential sources including any underwater noise sources of other objects. OBJ 1  and OBJ 2  represent two objects that produce acoustic signals that radiate as multiple planewaves PW 1  and PW 2  respectively. Fast Fourier Transform (FFT) processors  14 , shown as individual processors FFT( 1 ), FFT( 2 ), FFT( 3 ) . . . FFT (M) inclusively, process signals from the corresponding waveforms of M spaced hydrophones in the array  12 . Planewave-fit noise reduction routine  16  receives frequency domain output signals from the FFT processors  14 . Routine  16  provides modified frequency domain output signals FFT′( 1 ), FFT′( 2 ), FFT′( 3 ) . . . FFT′(M). A conventional measured covariance matrix processor  18  receives the modified output signals from routine  16 . The covariance matrix processor  18  output can be weighted by a weighting processor  20 . The weighted covariance processor output is provided to an inverse-beamforming planewave beamformer processor  22  for producing an estimated bearing to a selected object, OBJ 1  or OBJ 2 . 
   The remaining portions of the apparatus  10  utilize the estimated bearing signal from the inverse beamforming planewave beamformer  22  covariance matrix data supplied by the measured covariance matrix processor  18  to produce beam output values. 
   Inverse beamforming processor  22  uses the output of the measured covariance matrix processor  18  in its original or weighted form, thereby forming beam levels in all frequency bins and at a selected number of bearings every FFT averaging time. A peak selection circuit  26  selects those incremental locations that exhibit a maximum with respect to adjacent incremental locations. The foregoing processors operate iteratively over time. 
   An “M of N” tracker circuit  28  comprises a processor that utilizes the succession of signals from the peak selection circuit  26  during each iteration to eliminate false targets and enable a target classifier  30  to select classification  32  of a source object as a target or a noise source. This can also be performed by a user. When an object is identified as a noise source, bearing and frequency range data  34  is provided to planewave-fit noise reduction routine  16 . A target display  36  provides the track of the bearing, range and depth of each target over time. 
   The peak selection circuit  26  can be realized as an eight nearest-neighbor-peak picker as fully described in the U.S. Pat. No. 5,581,505. Application of this circuit results in the detection of all the peaks with relative maxima in beamformed levels on the beamformed FRAZ surface for a given time epoch, also more fully described in U.S. Pat. No. 5,481,505. A peak or relative maxima, beam level can be described by the following parameters: level; frequency; azimuth angle; azimuthal angle width; elevation angle; elevation angle width; and time location. 
   Beam Level on the FRAZ surface as a function of time is input to the ENNPP and tracked by the Inverse Beamforming M of N tracker circuit  20  in a manner more fully described in U.S. Pat. No. 5,481,505. Since sources of interest in detection are assumed to be distant point sources, azimuthal angle width and elevation angle width are not used in the peak picking process of the present invention. 
   The noise-suppression algorithm performs the noise suppression every FFT time epoch using the following mathematical relationship:
 
   FFT   ( M )= FFT ( M )− X ( m ,( k   i   a   k   i   b ),( l   a   ,l   b ))
 
Where:
 
m is the hydrophone index, 1≦m≦M,
 
(k i   a k i   b ) are the I frequency indices over which noise suppression is desired,
 
(l a ,l b ) are the bearing indices over which noise suppression is desired,
 
FFT(m) is the original hydrophone FFT,
   FFT (m) is the residual hydrophone FFT after noise suppression,
 
X( ) is equal to X I ( ) for a one planewave-fit and equals X II  (°) for a two planewave-fit.
 
   The X I ( ) is determined adaptively each FFT time epoch by fitting a planewave model with arbitrary complex amplitude, a k , noise source arrival direction u=cos(θ) (where θ is the arrival angle relative to forward end fire (θ=0) on a line array) to measured FFT data. The best fit for the model is defined as the {a k } and u that minimize the error between the model output and the data. The search space consists of all possible complex amplitudes in frequency intervals (k i   a ,k i   b ) and over all possible bearings in the angular sector of interest. 
   The desired noise source bearing indices (l a ,l b ) are set equal to the planewave model 
                   X   I     =       a   k     ⁢     exp   ⁡     [     -       i2π   ⁢           ⁢     kx   ⁡     (   m   )       ⁢   u       Nc   ⁢           ⁢   Δ         ]                 (   2   )               
where a k  and u have been determined by the search described below. This procedure is repeated every FFT time epoch. In equation 2, k is the frequency bin number with (k i   a ,k i   b ). (Note that  FFT (m)≡FFT(m) when k is not within the fitted frequency bin intervals (k i   a ,k i   b ).) The incoming sound is a continuous pressure field  p (t,x) at time t and location x. N is the number of time series points in the FFT time epoch (FFT size). Δ is the sampling interval in seconds between two adjacent points in the input time series. For each time epoch, t, N samples are collected, with each sample n being taken at a time Δ from the previous sample. Thus, time epoch, t=nΔ for n=0:N−1.
 
   Samples are also taken for each hydrophone along a linear array of hydrophones. x(m) is the position of the mth hydrophone in the linear array. c is the speed of sound in water. Thus, spatial samples are taken at locations x=x(m) for m=0:M−1 where x is the location which is a function of the hydrophone index m, for a total of M hydrophones. In these equations, the colon symbol J:K is used to denote {J,J+1, . . . , K} with J≦K. The total discrete data available is p(n,m)≡  p (nΔ,x(m)) for
 
 n= 0:N−1 ,m =0 :M −1.
 
   In further detail, amplitudes {a k } are determined in accordance with the following principles. If a single planewave arrives from angle u 1 =sinθ 1 , comprised of frequencies {f 1 (k)} k     b     k     a   , the observed complex pressure field  p   1 (t,x) is modeled as 
                         p   _     1     ⁡     (     t   ,   x     )       =       ∑     k   =     k   a         k   b       ⁢         a   1     ⁡     (   k   )       ⁢     exp   ⁡     [     i2   ⁢           ⁢   π   ⁢           ⁢       f   1     ⁡     (   k   )       ⁢     (     t   -       x   c     ⁢     u   1         )       ]             ,           (   3   )               
where the amplitudes {a 1 (k)} k     b     k     a    are complex. For given frequencies {f 1 ,(k)} k     b     k     c   , the amplitudes are chosen so that the weighted fitting error e is minimized. Here, error e is defined as
 
                             e   =       ⁢       ∑     n   =   0       N   -   1       ⁢       ∑     m   =   0       M   -   1       ⁢         w   t     ⁡     (   n   )       ⁢       w   x     ⁡     (   m   )       ⁢              p   ⁡     (     n   ,   m     )       -         p   _     1     ⁡     (       n   ⁢           ⁢   Δ     ,     x   ⁡     (   m   )         )              2                         =       ⁢       ∑     n   ,   m       ⁢         w   t     ⁡     (   n   )       ⁢       w   x     ⁡     (   m   )       ⁢              p   ⁡     (     n   ,   m     )       -       ∑   k     ⁢         a   1     ⁡     (   k   )       ⁢   exp   ⁢     {     i   ⁢           ⁢         α   1     ⁡     (   k   )       ⁡     [     n   -       β   ⁡     (   m   )       ⁢     u   1         ]         }                2           ,                   (   4   )               
where p(n,m) represents samples of the actual data,  p   1 (nΔ,x(m)) represents the fitted data samples, and the temporal and spatial weightings {w 1 (n)} and {w x (m)} are real and positive, and the known dimensionless parameters
 
                       α   1     ⁡     (   k   )       ≡     2   ⁢   π   ⁢           ⁢       f   1     ⁡     (   k   )       ⁢   Δ       ,           ⁢       β   ⁡     (   m   )       ≡         x   ⁡     (   m   )         c   ⁢           ⁢   Δ       .               (   5   )               
The temporal and spatial windows are given as:
 
                                 W   t     ⁡     (   α   )       =       ∑     n   =   0       N   -   1       ⁢         w   t     ⁡     (   n   )       ⁢     exp   ⁡     (       -   i     ⁢           ⁢   α   ⁢           ⁢   n     )                     for   ⁢           ⁢   all   ⁢           ⁢   α     ,         W   t     ⁡     (   0   )       =   1     ,                 (   6   )                         W   x     ⁡     (   γ   )       =       ∑     m   =   0       M   -   1       ⁢         w   x     ⁡     (   m   )       ⁢     exp   ⁡     (     i   ⁢           ⁢     β   ⁡     (   m   )       ⁢   γ     )                     for   ⁢           ⁢   all   ⁢           ⁢   γ     ,         W   x     ⁡     (   0   )       =   1     ,                 (   7   )                 
and the two-dimensional data spectrum is defined as
 
   
     
       
         
           
             
               
                 
                   
                     P 
                     ⁡ 
                     
                       ( 
                       
                         α 
                         , 
                         γ 
                       
                       ) 
                     
                   
                   = 
                   
                     
                       ∑ 
                       
                         n 
                         , 
                         m 
                       
                     
                     ⁢ 
                     
                       
                         
                           w 
                           t 
                         
                         ⁡ 
                         
                           ( 
                           n 
                           ) 
                         
                       
                       ⁢ 
                       
                         
                           w 
                           x 
                         
                         ⁡ 
                         
                           ( 
                           m 
                           ) 
                         
                       
                       ⁢ 
                       
                         p 
                         ⁡ 
                         
                           ( 
                           
                             n 
                             , 
                             m 
                           
                           ) 
                         
                       
                       ⁢ 
                       
                         exp 
                         ⁡ 
                         
                           [ 
                           
                             
                               
                                 - 
                                 i 
                               
                               ⁢ 
                               
                                   
                               
                               ⁢ 
                               α 
                               ⁢ 
                               
                                   
                               
                               ⁢ 
                               n 
                             
                             + 
                             
                               i 
                               ⁢ 
                               
                                   
                               
                               ⁢ 
                               
                                 β 
                                 ⁡ 
                                 
                                   ( 
                                   m 
                                   ) 
                                 
                               
                               ⁢ 
                               γ 
                             
                           
                           ] 
                         
                       
                       ⁢ 
                       
                           
                       
                       ⁢ 
                       for 
                       ⁢ 
                       
                           
                       
                       ⁢ 
                       all 
                       ⁢ 
                       
                           
                       
                       ⁢ 
                       α 
                     
                   
                 
                 , 
                 
                   γ 
                   . 
                 
               
             
             
               
                 ( 
                 8 
                 ) 
               
             
           
         
       
     
   
   The optimal complex amplitudes { a   1 (k)} k     b     k     a    for equation 4 must satisfy the simultaneous linear equations 
                       ∑     k   =     k   a         k   b       ⁢           a   _     1     ⁡     (   k   )       ⁢       W   1     ⁡     (         α   1     ⁡     (     k   _     )       -       α   1     ⁡     (     k   _     )         )       ⁢       W   x     ⁡     (       [         α   1     ⁡     (     k   _     )       -       α   1     ⁡     (   k   )         ]     ⁢     u   1       )           =     P   ⁡     (         α   1     ⁡     (     k   _     )       ,         α   1     ⁡     (     k   _     )       ⁢     u   1         )         ,           ⁢       for   ⁢           ⁢     k   _       =       k   a     :       k   b     .                 (   9   )               
For flat temporal weighting:
 
                       w   t     ⁡     (   n   )       =         1   N     ⁢           ⁢   for   ⁢           ⁢   n     =     0   :     N   -   1           ,           (   10   )               
and fitting frequencies
 
                       f   1     ⁡     (   k   )       =         k     N   ⁢           ⁢   Δ       ⁢           ⁢   for   ⁢           ⁢   k     =       k   a     :     k   b           ,         α   1     ⁡     (   k   )       =     2   ⁢   π   ⁢           ⁢     k   /   N         ,           (   11   )               
we can write an explicit solution to equation 9 for the optimal amplitudes, namely,
     a     1 ( k )= P (α 1 ( k ),α 1 ( k ) u   1 ) for  k=k   a   :k   b .  (12) 
The corresponding optimal (minimum) error is then expressible as
 
                   e   _     =         ∑     n   ,   m       ⁢         w   t     ⁡     (   n   )       ⁢       w   x     ⁡     (   m   )       ⁢            p   ⁡     (     n   ,   m     )            2         -       ∑   k     ⁢          P   (         α   1     ⁡     (   k   )       ,         α   1     ⁡     (   k   )       ⁢     u   1                            (   13   )               
The last term of this quantity  e  is redefined as follows:
 
   
     
       
         
           
             
               
                 
                   ∑ 
                   k 
                 
                 ⁢ 
                 
                    
                   
                     
                       P 
                       ( 
                       
                         
                           
                             α 
                             1 
                           
                           ⁡ 
                           
                             ( 
                             k 
                             ) 
                           
                         
                         , 
                         
                           
                             
                               α 
                               1 
                             
                             ⁡ 
                             
                               ( 
                               k 
                               ) 
                             
                           
                           ⁢ 
                           
                             u 
                             1 
                           
                         
                       
                        
                     
                     ≡ 
                     
                       
                         r 
                         ⁡ 
                         
                           ( 
                           
                             u 
                             1 
                           
                           ) 
                         
                       
                       / 
                       
                         N 
                         2 
                       
                     
                   
                 
               
             
             
               
                 ( 
                 
                   13 
                   ⁢ 
                   a 
                 
                 ) 
               
             
           
         
       
     
   
   We now further minimize the error e by choosing the planewave arrival angle u 1  that maximizes r(u 1 ), resulting in 
                     r   ⁡     (     u   1     )       =       ∑     k   =     k   a         k   b       ⁢              ∑     m   =   0       M   -   1       ⁢           w   x     ⁡     (   m   )       q     ⁢     (     k   ,   m     )     ⁢     exp   ⁡     [     i   ⁢       2   ⁢   π   ⁢           ⁢   k     N     ⁢     β   ⁡     (   m   )       ⁢     u   1       ]                2         ,           (   14   )               
where
 
                     q   ⁡     (     k   ,   m     )       =         1   N     ⁢       ∑     n   =   0       N   -   1       ⁢       p   ⁡     (     n   ,   m     )       ⁢     exp   ⁡     (       -   i2     ⁢           ⁢   π   ⁢           ⁢     kn   /   N       )       ⁢           ⁢   for   ⁢           ⁢   k         =     0   :     N   -   1           ,     m   =     0   :     M   -   1                 (   15   )               
is the temporal discrete Fourier transform of the m-th element data. The maximization of r(u 1 ) by choice of u 1  is depicted in  FIG. 2 . Equation 14 provides guidance as to the preferred processing. The time-space data {p(n,m)} is first transformed into the frequency-space domain {q(k,m)}, and then for the hypothesized arrival angle u 1 , the m-th element component in frequency bin k is scaled and phase-shifted by (2πk/N)β(m)u 1 . This phase shift exactly compensates for that of a single frequency planewave arriving at angle u 1 :
 
                   phaseshift   ⁡     (     k   ,   m     )       =         -   2     ⁢   π   ⁢     k     N   ⁢           ⁢   Δ       ⁢       x   ⁡     (   m   )       c     ⁢     u   1       =         2   ⁢           ⁢   π   ⁢           ⁢   k     N     ⁢     β   ⁡     (   m   )       ⁢       u   1     .                 (   16   )               
Thus, the inner complex sum over m in equation 14 is a coherent one for any planewave arriving at angle u 1 . Finally, the outer sum over k in equation 14 is an incoherent sum over frequencies in the band of interest. This incoherent sum is necessary because no interrelations have been assumed for individual frequency components in the planewave arrivals. After the sum over k is complete, then r(u 1 ) is plotted for all u 1  in the angular sector of interest and its maximum is located at û 1 . The complex sum on element number m in equation 14 cannot be accomplished by a Fast Fourier Transform (FFT) because, in general, β(m)=x(m)/(cΔ) is not linear in m for an unequally spaced line array, but it is linear for a sparse equispaced line array. This complex sum must be carried out by brute force. However, one shortcut available uses the recursion
 
                   exp   ⁡     [     i   ⁢       2   ⁢           ⁢   π   ⁢           ⁢   k     N     ⁢     β   ⁡     (   m   )       ⁢     u   1       ]       =       exp   ⁡     [     i   ⁢       2   ⁢           ⁢     π   ⁡     (     k   -   1     )         n     ⁢     β   ⁡     (   m   )       ⁢     u   1       ]       ⁢     exp   ⁡     [     i   ⁢       2   ⁢           ⁢   π     N     ⁢     β   ⁡     (   m   )       ⁢     u     1   ⁢                 ]                 (   17   )               
for each m and u 1  to generate the k-values needed for the exponentials.
 
   Finding the best spatial weights {w x (m)} is not trivial; they should not simply be taken as flat but should reflect known element locations {x(m)}. A Hann shading can be used for this. 
   After the conditionally best coefficients { a   1 (k)} are found for a specified u 1 , the minimal conditional time-space residual is, from equations 3 and 4, 
                           p   _     ⁡     (     n   ,   m     )       ≡               p   ~     ⁡     (       n   ⁢           ⁢   Δ     ,     x   ⁡     (   m   )         )       -         p   ~     1     ⁡     (       n   ⁢           ⁢   Δ     ,     x   ⁡     (   m   )         )                 =           p   ⁡     (     n   ,   m     )       -       ∑   k     ⁢           ⁢           a   _     1     ⁡     (   k   )       ⁢       exp   ⁡     [       -   ⅈ     ⁢       2   ⁢   π   ⁢           ⁢   k     N     ⁢     (     n   -       β   ⁡     (   m   )       ⁢     u   1         )       ]       .                       (   18   )               
The corresponding FFT of this residual is, from equation 15,
 
                         ⁢               q   _     ⁡     (     k   ,   m     )       =       1   N     ⁢       ∑     n   =   o       N   -   1       ⁢           ⁢         p   _     ⁡     (     n   ,   m     )       ⁢     exp   ⁡     (       -   ⅈ2π     ⁢           ⁢     k   _     ⁢     n   /   N       )                         =       q   ⁡     (       k   _     ,   m     )       -       1   N     ⁢       ∑     k   =     k   a         k   b       ⁢           ⁢           a   _     1     ⁡     (   k   )       ⁢     exp   ⁡     [       -     
     ⁢   ⅈ     ⁢       2   ⁢   π   ⁢           ⁢   k     N     ⁢     β   ⁡     (   m   )       ⁢     u   1       ]       ⁢                     ∑               ⁢     n   =   0                 ⁢     N   -   1         ⁢           ⁢     exp   ⁡     [       -     ⅈ2π   ⁡     (       k   _     -   k     )         ⁢     n   /   N       ]         ︸       Nδ   ⁡     (       k   _     -   k     )                               =       q   ⁡     (     k   ,   m     )       -           a   _     1     ⁡     (     k   _     )       ⁢     exp   ⁡     [       -   ⅈ     ⁢       2   ⁢   π   ⁢           ⁢     k   _       N     ⁢     β   ⁡     (   m   )       ⁢     u   1       ]             ,       for   ⁢           ⁢     k   _       =       k   a     :       k   b     .                         (   19   )               
The conditionally optimal amplitudes are, from equations 8, 10, 11, 12 and 15,
 
                             a   _     1     ⁡     (   k   )       =     P   ⁡     (         α   1     ⁡     (   k   )       ,         α   1     ⁡     (   k   )       ⁢     u   1         )                   =         ∑     M   =   0       M   -   1       ⁢           ⁢         w   x     ⁡     (   m   )       ⁢     q   ⁡     (     k   ,   m     )       ⁢     exp   ⁡     [     ⅈ   ⁢       2   ⁢   π   ⁢           ⁢   k     N     ⁢     β   ⁡     (   m   )       ⁢     u   1       ]       ⁢           ⁢   for   ⁢           ⁢   k       =       k   a     :       k   b     .                       (   20   )               
These complex amplitudes should be evaluated only after the best arrival angle u 1  (namely û 1 ) has been determined. Then, equation 23 should be evaluated, using û 1  in place of u 1 . That is, the unconditionally optimum amplitudes are
 
                             a   ^     1     ⁡     (   k   )       =     P   ⁡     (         α   1     ⁡     (   k   )       ,         α   1     ⁡     (   k   )       ⁢       u   ^     1         )                     =         ∑     M   =   0       M   -   1       ⁢           ⁢         w   x     ⁡     (   m   )       ⁢     q   ⁡     (     k   ,   m     )       ⁢     exp   ⁡     [     ⅈ   ⁢       2   ⁢   π   ⁢           ⁢   k     N     ⁢     β   ⁡     (   m   )       ⁢       u   ^     1       ]       ⁢           ⁢   for   ⁢           ⁢   k       =       k   a     :     k   b           ,                 (   21   )               
and the minimal residual in the frequency-space domain is
 
                   q   ⁡     (     k   ,   m     )       =     {             q   ⁡     (     k   ,   m     )       -           a   ^     1     ⁡     (   k   )       ⁢     exp   ⁡     [       -   ⅈ     ⁢       2   ⁢   π   ⁢           ⁢   k     N     ⁢     β   ⁡     (   m   )       ⁢       u   ^     1       ]                     for   ⁢           ⁢   k     =       k   a     :     k   b         ,               q   ⁡     (     k   ,   m     )               for   ⁢           ⁢   k     ∉       k   a     :       k   b     .                         (   22   )               
Letting the spatially weighted FFT at frequency bin k be defined as
   q   w ( k,m )= q ( k,m ) w   x ( m ),  (23) 
and its spatial autocorrelation for frequency bin k be defined as
 
                     ϕ   ⁡     (     k   ,   j     )       =       ∑   m     ⁢           ⁢         q   w     ⁡     (     k   ,   m     )       ⁢       q   w   *     ⁡     (     k   ,     m   -   j       )             ,           (   24   )               
then, from equation 12, if x(m)=dm (equally spaced line array),
 
                         r   ⁡     (     u   1     )       =       ∑     k   =     k   a         k   b       ⁢           ⁢              ∑     m   =   0       M   -   1       ⁢           ⁢         q   w     ⁡     (     k   ,   m     )       ⁢     exp   ⁡     [     ⅈ   ⁢       2   ⁢   π   ⁢           ⁢   k     N     ⁢     dm     c   ⁢           ⁢   Δ       ⁢     u   1       ]                2                   =       ∑     k   =     k   a         k   b       ⁢       ∑     j   =     1   -   M         M   -   1       ⁢           ⁢       ϕ   ⁡     (     k   ,   j     )       ⁢       exp   ⁡     (     ⅈ   ⁢       2   ⁢   π   ⁢           ⁢   k     N     ⁢     d     c   ⁢           ⁢   Δ       ⁢     ju   1       )       .                         (   25   )               
A FIM-like generalization would be to weight the spatial sum over j to obtain the modification
 
                       r   s     ⁡     (     u   1     )       =       ∑     k   =     k   a         k   b       ⁢           ⁢       ∑     j   =     1   -   M         M   -   1       ⁢           ⁢         w   s     ⁡     (   j   )       ⁢     ϕ   ⁡     (     k   ,   j     )       ⁢     exp   ⁡     (     ⅈ   ⁢       2   ⁢   π   ⁢           ⁢   k     N     ⁢     d     c   ⁢           ⁢   Δ       ⁢   j   ⁢           ⁢     u   1       )               ,           (   26   )               
where the {w(j)} M−1   1−M  are separation weights (real and even about j=0).
 
   The two planewave-fit function, X II , is determined in a similar way to the one planewave-fit function, except that there are assumed to be two coherent planewaves interfering with within the frequency interval (k i   a ,k i   b ) and bearing interval (l a ,l b ) Therefore, X II  is given by: 
                   X   II     =         a   k     (   1   )       ⁢     exp   ⁡     [     -       ⅈ2π   ⁢           ⁢     kx   ⁡     (   m   )       ⁢     u     (   1   )           Nc   ⁢           ⁢   Δ         ]         +       a   k     (   2   )       ⁢     exp   ⁡     [     -       ⅈ2π   ⁢           ⁢     kx   ⁡     (   m   )       ⁢     u     (   2   )           Nc   ⁢           ⁢   Δ         ]                   (   27   )               
and the best fit to the data requires a coherent four-parameter search in two complex amplitude number sets {a (1)    k } and {a (2)   k } and in two real angles, u (1)  and u (2) , for the two planewaves arrivals. In this way, the two planewaves are allowed to interfere with one another coherently. The process to determine X II  is repeated every FFT time epoch.
 
   In deriving this search, the data available is again {p(n,m)} for n=0:N−1 and m=0:M−1, and the fitting frequencies are 
   
     
       
         
           
             
               
                 
                   f 
                   k 
                 
                 = 
                 
                   
                     k 
                     T 
                   
                   = 
                   
                     
                       
                         k 
                         
                           N 
                           ⁢ 
                           
                               
                           
                           ⁢ 
                           Δ 
                         
                       
                       ⁢ 
                       
                           
                       
                       ⁢ 
                       for 
                       ⁢ 
                       
                           
                       
                       ⁢ 
                       k 
                     
                     = 
                     
                       
                         k 
                         a 
                       
                       : 
                       
                         
                           k 
                           b 
                         
                         . 
                       
                     
                   
                 
               
             
             
               
                 ( 
                 28 
                 ) 
               
             
           
         
       
     
   
   The range (k a ,k b ) represents the common frequencies for the two planewave-fits. It is possible to generalize to the case where only some of the frequencies {f k } are common to the two planewaves; the corresponding complex amplitudes {a 1 (k)} and{a 2 (k)} would then be set to zero in the appropriate frequency bins. 
   If two planewaves arrive from angles u 1  and u 2 , the observed complex pressure field over observation interval T=NΔ is modeled as 
                       p   ^     ⁡     (     t   ,   x     )       ≡         ∑     k   =     k   a         k   b       ⁢           ⁢         a   1     ⁡     (   k   )       ⁢     exp   ⁡     [     ⅈ2π   ⁢     k   T     ⁢     (     t   -       x   c     ⁢     u   1         )       ]           +       ∑     k   =     k   a         k   b       ⁢           ⁢         a   2     ⁡     (   k   )       ⁢     exp   ⁡     [     ⅈ2   ⁢           ⁢   π   ⁢     k   T     ⁢     (     t   -       x   c     ⁢     u   2         )       ]               ,           (   29   )               
where the amplitudes {a 1 (k)} and {a 2 (k)} are complex. The two sources are assumed to be uncorrelated with each other. The sampled pressure field is
 
                     p   ^     ⁡     (       n   ⁢           ⁢   Δ     ,     x   ⁡     (   m   )         )       =         ∑     k   =     k   a         k   b       ⁢           ⁢         a   1     ⁡     (   k   )       ⁢     exp   ⁡     [     ⅈ   ⁢           ⁢       α   k     ⁡     (     n   -       β   ⁡     (   m   )       ⁢     u   1         )         ]           +       ∑     k   =     k   a         k   b       ⁢         a   2     ⁡     (   k   )       ⁢       exp   ⁡     [       ⅈα   k     ⁡     (     n   -       β   ⁡     (   m   )       ⁢     u   2         )       ]       .                   (   30   )               
An average squared error is defined between the data and the model as:
 
                   ⅇ   =       1   N     ⁢       ∑     n   ,   m       ⁢           ⁢         w   x     ⁡     (   m   )       ⁢                p   ^     ⁡     (       n   ⁢           ⁢   Δ     ,     x   ⁡     (   m   )         )       -       p   ^     ⁡     (       n   ⁢           ⁢   Δ     ,     x   ⁡     (   m   )         )              2             ,           (   31   )               
where we utilize flat temporal weighting w t (n)=1/N, and allow a spatial error-weighting function w x (m), which is real and positive. Then the error e can be expressed in the form:
 
                       ⅇ   =       1   N     ⁢           ⁢       ∑     n   ,   m       ⁢           ⁢         w   x     (         ⁢           ⁢   m   ⁢           ⁢         )     ⁢           ⁢            p   ⁡     (     n   ,   m     )       -       ∑     k   =     k   a         k   b       ⁢           ⁢         a   1     ⁡     (   k   )       ⁢     exp   ⁡     [     ⅈ   ⁢           ⁢       α   k     ⁡     (     n   -       β   ⁡     (   m   )       ⁢     u   1         )         ]           -           ⁢       (   32   )     ⁢           ⁢     
     ⁢               ∑               ⁢     k   =     k   a                   ⁢       k   b     2         ⁢       a   2     ⁡     (   k   )       ⁢     exp   ⁡     [     ⅈ   ⁢           ⁢       α   k     ⁡     (     n   -       β   ⁡     (   m   )       ⁢     u   2         )         ]                                  =         1   N     ⁢       ∑     n   ,   m       ⁢         w   x     ⁡     (   m   )       ⁢            d   ⁡     (     n   ,   m     )            2           =       ∑   m     ⁢         w   x     ⁡     (   m   )       ⁢     1   N     ⁢       ∑   n     ⁢              d   ⁡     (     n   ,   m     )            2     .                         □             
Letting the FFT be a discrete FFT:
 
                   D   ⁡     (     l   ,   m     )       =         ∑     n   =   0       N   -   1       ⁢           ⁢       exp   ⁡     (       -   ⅈ2π     ⁢           ⁢     nl   /   N       )       ⁢     d   ⁡     (     n   ,   m     )       ⁢           ⁢   for   ⁢           ⁢   l       =     0   :     N   -   1                 (   33   )               
it follows that
 
                           1     N   2       ⁢       ∑     l   =   0       N   -   1       ⁢                  ⁢     D   ⁡     (     l   ,   m     )            2         =       1   N     ⁢       ∑     n   =   0       N   -   1       ⁢                ⁢     d   ⁡     (     n   ,   m     )       )                2     ,           (   34   )               
which is a discrete form of Parseval&#39;s Theorem. The error e follows from equation 32 as
 
                 ⅇ   =       ∑   m     ⁢           ⁢         w   x     ⁡     (   m   )       ⁢     1     N   2       ⁢       ∑   l     ⁢           ⁢              D   ⁡     (     l   ,   m     )            2     .                   (   35   )               
But from equations 32 and 33,
 
                     1   N     ⁢     D   ⁡     (     l   ,   m     )         =       q   ⁡     (     l   ,   m     )       -         a   1     ⁡     (   l   )       ⁢     exp   ⁡     [       -   i     ⁢           ⁢     α   1     ⁢     β   ⁡     (   m   )       ⁢     u   1       ]         -         a   2     ⁡     (   l   )       ⁢     exp   ⁡     [       -   i     ⁢           ⁢     α   1     ⁢     β   ⁡     (   m   )       ⁢     u   2       ]                   (   36   )               
for l=0:N−1, m=0:M−1. Therefore the error e in equation 35 becomes
 
                 ⅇ   =       ∑   m     ⁢           ⁢         w   x     ⁡     (   m   )       ⁢       ∑   l     ⁢                q   ⁡     (     l   ,   m     )       -         a   1     ⁡     (   l   )       ⁢     exp   ⁡     [       -   ⅈ     ⁢           ⁢     α   1     ⁢     β   ⁡     (   m   )       ⁢     u   1       ]         -         a   2     ⁡     (   l   )       ⁢     exp   ⁡     [       -   ⅈ     ⁢           ⁢     α   1     ⁢     β   ⁡     (   m   )       ⁢     u   2       ]                2     ⁢           .                   (   36   )               
Define functions
 
   
     
       
         
           
             
               
                 
                   
                     
                       X 
                       ⁡ 
                       
                         ( 
                         
                           k 
                           , 
                           γ 
                         
                         ) 
                       
                     
                     ≡ 
                     
                       
                         ∑ 
                         m 
                       
                       ⁢ 
                       
                           
                       
                       ⁢ 
                       
                         
                           
                             w 
                             x 
                           
                           ⁡ 
                           
                             ( 
                             m 
                             ) 
                           
                         
                         ⁢ 
                         
                           q 
                           ⁡ 
                           
                             ( 
                             
                               k 
                               , 
                               m 
                             
                             ) 
                           
                         
                         ⁢ 
                         
                           exp 
                           ⁡ 
                           
                             [ 
                             
                               ⅈ 
                               ⁢ 
                               
                                   
                               
                               ⁢ 
                               
                                 γβ 
                                 ⁡ 
                                 
                                   ( 
                                   m 
                                   ) 
                                 
                               
                             
                             ] 
                           
                         
                         ⁢ 
                         
                             
                         
                         ⁢ 
                         for 
                         ⁢ 
                         
                             
                         
                         ⁢ 
                         k 
                       
                     
                   
                   = 
                   
                     
                       k 
                       a 
                     
                     : 
                     
                       k 
                       b 
                     
                   
                 
                 , 
                 
                   all 
                   ⁢ 
                   
                       
                   
                   ⁢ 
                   γ 
                 
                 , 
               
             
             
               
                 ( 
                 37 
                 ) 
               
             
           
           
             
               
                 
                   
                     W 
                     x 
                   
                   ⁡ 
                   
                     ( 
                     γ 
                     ) 
                   
                 
                 ≡ 
                 
                   
                     ∑ 
                     m 
                   
                   ⁢ 
                   
                     
                       
                         w 
                         x 
                       
                       ⁡ 
                       
                         ( 
                         m 
                         ) 
                       
                     
                     ⁢ 
                     
                       exp 
                       ⁡ 
                       
                         [ 
                         
                           ⅈ 
                           ⁢ 
                           
                               
                           
                           ⁢ 
                           
                             γβ 
                             ⁡ 
                             
                               ( 
                               m 
                               ) 
                             
                           
                         
                         ] 
                       
                     
                     ⁢ 
                     
                         
                     
                     ⁢ 
                     for 
                     ⁢ 
                     
                         
                     
                     ⁢ 
                     all 
                     ⁢ 
                     
                         
                     
                     ⁢ 
                     
                       γ 
                       . 
                     
                   
                 
               
             
             
               
                 ( 
                 38 
                 ) 
               
             
           
         
       
     
   
   Minimization of error e results in simultaneous equations for the best fitting coefficients { a   1 (k)} and { a   2 (k)}, namely, 
                                       a   _     1     ⁡     (   k   )       +       W   k     ⁢         a   _     2     ⁡     (   k   )           =     X   (     k   ,       α   k     ⁢     u   1                             W   k   *     ⁢         a   _     1     ⁡     (   k   )         +         a   _     2     ⁡     (   k   )         =     X   (     k   ,       α   k     ⁢     u   2                   }     ⁢   for   ⁢           ⁢   k     =       k   a     :     k   b         ,           (   39   )               
where W k ≡W x (u 1 −u 2 )). The solutions for the optimal amplitudes are, for k=k a :k b ,
 
                         a   _     1     ⁡     (   k   )       =         X   ⁡     (     k   ,       α   k     ⁢     u   1         )       -       W   k     ⁢     X   ⁡     (     k   ,       α   k     ⁢     u   2         )             1   -            W   k          2           ,           (   40   )                     a   _     2     ⁡     (   k   )       =           X   ⁡     (     k   ,       α   k     ⁢     u   2         )       -       W   k   *     ⁢     X   ⁡     (     k   ,       α   k     ⁢     u   1         )             1   -            W   k          2         .             (   41   )               
Decoupling in frequency index k occurs due to the flat temporal weighting and the fact that frequency increment Δf=1T. The optimum value of error e is now, from equation 36,
 
                         ⅇ   _     =       ∑   m     ⁢           ⁢         w   x     ⁡     (   m   )       ⁢       ∑   k     ⁢       {       q   ⁡     (     k   ,   m     )       -           a   _     1     ⁡     (   k   )       ⁢     exp   ⁡     [       -   ⅈ     ⁢           ⁢     α   k     ⁢     β   ⁡     (   m   )       ⁢     u   1       ]         -           ⁢       (   42   )     ⁢     
     ⁢         a   _     2     ⁡     (   k   )       ⁢     exp   ⁡     [       -   ⅈ     ⁢           ⁢     α   k     ⁢     β   ⁡     (   m   )       ⁢     u   2       ]           }     ⁢   q   *     (     k   ,   m     )                         =         ∑     m   ,   k       ⁢           ⁢         w   x     ⁡     (   m   )       ⁢            q   ⁡     (     k   ,   m     )            2         -       ∑   k     ⁢           a   _     1     ⁡     (   k   )       ⁢   X   *     (     k   ,       α   k     ⁢     u   1         )         -       ∑   k     ⁢           a   _     2     ⁡     (   k   )       ⁢   X   *     
     ⁢     (     k   ,       α   k     ⁢     u   2         )                         =         ∑     m   ,   k       ⁢         w   x     ⁡     (   m   )       ⁢            q   ⁡     (     k   ,   m     )            2         -     τ   ⁡     (       u   1     ,     u   2       )           ,               □             
where we have used the solutions for  a   1 (k)and  a   2 (k) from equations 40 and 41, the fact that the initial sum over m and k is independent of u 1  and u 2 , and defined function
 
                   r   ⁡     (       u   1     ,     u   2       )       =           ⁢       ∑     k   =     k   a         k   b       ⁢           ⁢       {                    X   ⁡     (     k   ,         2   ⁢   π   ⁢           ⁢   k     N     ⁢     u   1         )            2     +            X   ⁡     (     k   ,         2   ⁢   π   ⁢           ⁢   k     N     ⁢     u   2         )            2     -               2   ⁢   ℜ   ⁢     {         W   x     ⁡     (         2   ⁢           ⁢   π   ⁢           ⁢   k     N     ⁢     (       u   1     -     u   2       )       )       ⁢   X   *     (     k   ,         2   ⁢   π   ⁢           ⁢   k     N     ⁢     u   1         )     ⁢     X   ⁡     (     k   ,       2   ⁢   π   ⁢           ⁢   k     N       )       ⁢     u   2       }             }       1   -              W   x     ⁡     (         2   ⁢   π   ⁢           ⁢   k     N     ⁢     (       u   1     -     u   2       )       )            2                   (   43   )               
Here             denotes the real part, and from equations 37 and 38 we have function values

                     X   ⁡     (     k   ,         2   ⁢   π   ⁢           ⁢   k     N     ⁢   u       )       =       ∑     m   =   0       M   -   1       ⁢           ⁢         w   x     ⁡     (   m   )       ⁢     q   ⁡     (     k   ,   m     )       ⁢     exp   ⁡     [     ⅈ   ⁢       2   ⁢   π   ⁢           ⁢   k     N     ⁢       x   ⁡     (   m   )         c   ⁢           ⁢   Δ       ⁢   u     ]             ,           (   44   )                     W   x     ⁡     (         2   ⁢   π   ⁢           ⁢   k     N     ⁢   u     )       =       ∑     m   =   0       M   -   1       ⁢           ⁢         w   x     ⁡     (   m   )       ⁢     exp   ⁡     [     ⅈ   ⁢       2   ⁢   π   ⁢           ⁢   k     N     ⁢       x   ⁡     (   m   )         c   ⁢           ⁢   Δ       ⁢   u     ]             ,           (   45   )                 q   ⁡     (     k   ,   m     )       =         1   N     ⁢       ∑     n   =   0       N   -   1       ⁢           ⁢       p   ⁡     (     n   ,   m     )       ⁢     exp   ⁡     (       -   ⅈ2     ⁢           ⁢   π   ⁢           ⁢   k   ⁢           ⁢     n   /   N       )       ⁢           ⁢   for   ⁢           ⁢   m         =     0   :     M   -   1                 (   46   )               
which are defined for k=0:N−1. However, only the values for k=k a :k b  are needed. Equation 46 consists of M N-point temporal discrete Fourier transforms.
 
   To further minimize e in equation 42, we must maximize the quantity r(u 1 ,u 2 ) given by equation 43 by choice of both arrival angles u 1  and u 2 . We may take u 1 &lt;u 2  without loss of generality. Let û 1 ,û 2  be the location of the maximum of r(u 1 ,u 2 ). 
   The factor 
               2   ⁢   π   ⁢           ⁢   k     N     ⁢       x   ⁡     (   m   )         c   ⁢           ⁢   Δ       ⁢   u         
in equation 44 is exactly the phase
 
compensation required at frequency
 
           k     N   ⁢           ⁢   Δ           
and array element m, in order to coherently “line up” all components arriving at angle u (u=0 is broadside). Thus, equation 44 is a “coherent” sum done prior to the “incoherent” sum over k (frequency) in equation 43. The sum in equation 45 cannot be done a priori in closed form because of the irregular element locations {x(m)}; also W x (γ) is complex and will remain so for general array element locations {x(m)}.
 
   Notice, from equation 43, that for given arrival angles u 1 ,u 2  and index k, only three complex quantities need be computed. However, they must be combined according to equation 43, yielding a purely real quantity, which is then summed over the frequency band of interest, k=k a :k b . Equations 44 and 45 cannot be done with the fast Fourier transform because the element locations {x(m)} are presumed unevenly spaced. 
   A measure of the total power in arrival 1 is available, by reference to equation 18 as 
                     P   1     =       ∑     k   =     k   a         k   b       ⁢           ⁢                  a   ^     _     1     ⁡     (   k   )            2         ,           (   47   )               
where { â   1 (k)} are the optimal coefficients obtained from equations 40 and 41, using û 1  and û 2 , after the best angles, namely, û 1  and û 2 , have been determined from the maximization of r(u 1 ,u 2 ) given by equation 43. If P 1  is large, then source 1 could be coherently subtracted from the total input. If P 1  is small, then the coefficients {{circumflex over ( a )} 1 (k)} could be discarded and considered as noise.
 
   With the optimal coefficients {{circumflex over ( a )} 1 (k)} and {{circumflex over ( a )} 2 (k)} for specified u 1  and u 2  the minimal residual is, using equation 30, 
                           p   _     ⁡     (     n   ,   m     )       ≡       p   ⁡     (     n   ,   m     )       -       p   ^     ⁡     (       n   ⁢           ⁢   Δ     ,     x   ⁡     (   m   )         )                             ⁢     =       p   ⁡     (     n   ,   m     )       -       ∑   k     ⁢           ⁢           a   _     1     ⁡     (   k   )       ⁢     exp   ⁡     [     ⅈ   ⁢       2   ⁢   π   ⁢           ⁢   k     N     ⁢     (     n   -       β   ⁡     (   m   )       ⁢     u   1         )       ]           -                               ∑               ⁢   k       ⁢         a   _     2     ⁡     (   k   )       ⁢       exp   ⁡     [     ⅈ   ⁢       2   ⁢   π   ⁢           ⁢   k     N     ⁢     (     n   -       β   ⁡     (   m   )       ⁢     u   2         )       ]       .                   (   48   )               
The corresponding residual in the frequency-space domain, using equation 46 and the best angles û 1  and û 2 , is
 
                       q   ^     _     ⁡     (     k   ,   m     )       =         q   ⁡     (     k   ,   m     )       -             a   ^     _     1     ⁡     (   k   )       ⁢     exp   ⁡     [       -   ⅈ     ⁢       2   ⁢   π   ⁢           ⁢   k     N     ⁢     β   ⁡     (   m   )       ⁢       u   ^     1       ]         -             a   _     ^     2     ⁡     (   k   )       ⁢     exp   ⁡     [       -   ⅈ     ⁢       2   ⁢   π   ⁢           ⁢   k     N     ⁢     β   ⁡     (   m   )       ⁢       u   ^     2       ]       ⁢           ⁢   for   ⁢           ⁢   k       =       k   a     :       k   b     .                 (   49   )               
This residual is devoid of the two strongest planewave arrivals and can now be processed further in order to detect additional weak arrivals. Without this coherent subtraction, the two strong arrivals could override a weak arrival in close angular proximity.
 
   The processing after noise suppression proceeds with the new  FFT (m) determined from equations 1, 2, and 3 being input to the measured covariance matrix processor  16  in  FIG. 1 . 
   Processing can be performed over any frequency interval of interest, but utilizes a minimum frequency bin size for digital processors since most data processors are digital and the frequency spectrum of a beam output is generated by a fast Fourier Transform (FFT) (For analog processors, bandwidth can be any arbitrary value.) 
   Peaks found by the above algorithm are processed with the M-of-N tracker circuits  22  more fully described in U.S. Pat. No. 5,481,505. 
   The advantage of noise suppression is to improve detection of signals in noisy background by reducing noise levels coherently without decreasing the desired signal levels. Signals of interest in this case are submarines operating submerged and producing signals of generally low levels. Noise comes from various sources including surface shipping, wind, waves, marine life, seismic activity, and seismic profilers. The noise sources originating from the sea surface, shipping, wind, waves, and seismic profiling will be suppressed significantly in level, by the application of the algorithms detailed herein. This use of noise suppression in the practice of the present invention, therefore, improves detection of submerged signals of interest. 
   It will be understood that various changes in the detail, steps and arrangement of parts, which have been herein described and illustrated, in order to explain the nature of the invention, may be made to those skilled in the art with the principle and scope of the invention as expressed in the independent claims.