Method for separation of data into narrowband and broadband time series components

[0003] (1) Field of the Invention

[0004] The present invention relates to a method for detecting signals in noise. More specifically, the present invention relates to the separation of known or unknown time series data into its narrowband and broadband components.

[0005] (2) Description of the Prior Art

[0006] Many applications require the identification of a desired signal within undesired random signals (noise) which is received with, and often interferes with, the desired signal. For example, in sonar systems randomly generated sounds from both natural and man-made sources give rise to a noise that interferes with desired acoustic signals. The detection and identification of a specific target, such as an underwater vehicle, requires a system which can detect a signal corresponding to the target within received data containing both the signal and noise.

[0007] Underwater acoustic signals are often complicated, consisting of a superposition: $\begin{array}{cc}x\left(t\right)=\sum _{k=1}^p{n}_k\left(t\right)+b\left(t\right)t_o\le t\le t_o+T& \left(1\right)\end{array}$

[0008] of non-stationary narrowband (hereinafter designated as NB) n_k(t) and broadband (hereinafter designated as BB) b(t) components. Such signals arise from a variety of sources, such as ship machinery, marine mammals, drilling platforms, and active sonars. FIG. 1 illustrates what a typical underwater acoustic spectrum might look like, consisting of a superposition of narrowband line-like components plus a broadband component. The BB component itself can be colored, with many local spectral peaks and valleys.

[0009] There is a need for the separation of the signal x(t) into the constituent NB:n_k(t)=S^p_k=1n_k(t) and BB:b(t) time series components when little or nothing is known about the NB and BB components and only a short data record is available. This is of great importance in sonar, especially passive sonar, where the desirable signal (for detection, classification, and localization) is often either the NB or BB component and the other is regarded as interference.

[0010] Existing approaches primarily deal with power spectrum estimation of the composite data, rather than recovering the constituent narrowband and broadband time series. Having the separated narrowband and broadband time series available is very useful since it allows many additional forms of processing (processing which is impossible or difficult to do using only the power spectrum), such as extraction of time series statistics, improved wavelet and Wigner analysis, pattern recognition, and parametric modeling.

[0011] The separation of the time series data into the NB and BB time series components is difficult. Wiener filtering (see S. Haykin, “Adaptive Filter Theory, Third Edition”, Prentice Hall, 1996) is not practical since the covariance or spectral densities of the NB and BB components are not known. Parametric methods, e.g., MA (Moving Average), AR (Auto Regressive), ARMA (Auto Regressive Moving Average) modeling (see P. Stoica et al., “Introduction to Spectral Analysis”, Prentice Hall, 1997), require choosing a model type for the underlying broadband component and for each of the NB components present. This is difficult since nothing is known about the NB and BB components. Adaptive methods applied directly to the time series, such as adaptive notch filters and line enhancers (see S. Haykin, supra) and Principal Component Inverse (PCI) method (see D. Tufts et al., “Data Adaptive Estimation by Singular-Value Decomposition of a Data Matrix” Proc. IEEE, Vol. 7, pp. 684-685, 1982; I. P. Kirsteins et al., “Adaptive Detection Using Low Rank Approximation to a Data Matrix” IEEE, Trans. Aerospace and Elect. Sys., Vol. 30, No. 1, pp. 55-57, 1994), tend to perform poorly when the broadband spectrum has a large dynamic range. That is, if a weak NB component is present in a “valley” of the BB spectrum, the notch filter tracker or PCI method might lock onto a nearby peak of the BB spectrum and filter it as the NB component, rather than the true NB component.

[0012] Accordingly, an object of the subject invention is to provide a new method, based on Thomson's multiple taper spectral estimation (hereinafter designated as MTSE) technique (see D. J. Thompson “Spectral Estimation and Harmonic Analysis”, Proc. of IEEE, Vol. 70, No. 9, pp. 1055-1096, 1992; “An Overview of Multi-Window and Quadratic-Inverse Spectrum Estimation Method” (ICASSP-94, Vol. VI, pp. 185-194, 1994; P. Stoica, supra), and PCI method (see D. Tufts, supra; I. P. Kirsteins, supra), which alleviates the problems of false tracking inherent in adaptive notch filters and adaptation to the background BB and NB components.

[0013] Thomson proposed Multiple Taper Spectral Estimation (MTSE) (see D. J. Thompson, infra; P. Stoica, supra) for spectral analysis of complicated non-stationary data consisting of lines plus a background component with continuous spectrum in which the physical processes generating the data were poorly understood. A multiple taper spectral estimate, as shown in FIG. 2, is given by: $\begin{array}{cc}|{\stackrel{\_}{X}\left(f_k\right)}^2=\frac{1}{K}{\Vert U^{H}D\left(f_k\right)x\Vert }^{2}F& \left(2\right)\end{array}$

[0014] where U is a N×K matrix whose K columns are the principal Discrete Prolate Spherioidal Sequence's (DPSS's) u^k_n(N,W), which are the eigenvectors of the N×N matrix, $\begin{array}{cc}\left[R\right]m,n=\frac{\sin 2\pi W\left(n-m\right)}{\pi \left(n-m\right)},& \left(3\right)\end{array}$

[0015] arranged to correspond to the eigenvalues in descending order, D(f)=diag(1, e^−i2xf, . . . , e^{−i(N−1)2xf}), W is the analysis bandwidth parameter, and K^˜2NW.

[0016] The matrix U in equation 2 acts as a lowpass filter with bandwidth W. Thus, the projection of the data vector x onto V_k=D(f_k)U^Hbecomes a bandpass filtering of the data to [f_k−W, f_k+W]. X (f) is approximately the average energy in the band [f_k−W, f_k+W]. For a given set of frequency points, f₁, . . , f_M, the spectrum estimate of equation 2 is analogous to filtering x into the subbands {[f₀−W, f₀+W], [f₁−W], f₁+W, . . . , [f_M−W, f_M+W]}, as shown in FIG. 2, and then calculating the average energy in each band.

[0017] The multiple taper method provides a simple and effective way of locating and removing tonals from the underlying continuous BB spectrum component. To avoid mistaking local BB spectral peaks as NB components and for avoiding interference from adjacent tonals, Thomson proposed that the estimation and removal of tonals be done separately within each subband, as depicted in FIG. 2. The main idea is that the projection of x onto V_keffectively isolates the frequency band [f_k−W, f_k+W] from out of band tonals, and if W is properly chosen, the background noise spectrum is approximately locally flat or white. Thus, effects from out of band tonals and problems of locking onto local spectral peaks of the broadband component are minimized.

[0018] An object of the invention is to use the methodology of MTSE to recover the NB:n(t)=Σ^p_k=1n_k(t) and BB:b(t) time series components. In concatenation of the matrix filter banks V_kinto matrix [V₁|V₂| . . . |V_h], the total filter bank output can be written as z=V^H_x. If rank [V]=N (where N is the number of time series samples in x), then×can be reconstructed from z using x=(V^H)#z, where operator # means pseudo-inverse. This suggests that an estimate of the broadband time series might be x_BB=(V^H)#_Zcleanwhere

[0019] Z_clean^T=[clean(V₁^Hx)^T|clean(V₂^Hx)^T. . . |clean(V_m^Hx)^T]  (4)

[0020] and the operator “clean” depicts the tonal removal procedure depicted in FIG. 2, applied to the subbands where tonals are present.

[0021] Two problems with the above procedure are: (1) the NB components may not necessarily be pure tones (they could be from processes with bandwidth), and (2) leakage of NB components into adjacent subbands is not dealt with, since it is not possible to construct ideal passband filters.

[0022] A further object of the invention is to alleviate these two concerns by estimating and removing the NB components using the PCI method and then developing a simple method to estimate and remove leakage components.

[0023] A still further object of the invention is to provide a new procedure for removing the BB time series from the NB components.

[0024] With the above and other objects in view, as will hereinafter appear, a feature of the present invention is the provision of a method for separation of acoustical data into narrowband and broadband time series components, the method comprising the steps of:

[0025] (1) constructing a data matrix in which time series samples of acoustical data are arranged in a forward-backward linear predictor matrix;

[0026] (2) decomposing a subband of the matrix into an in-band component and an out-of-band component;

[0027] (3) estimating narrowband components;

[0028] (4) estimating out-of-band leakage;

[0029] (5) removing the in-band narrowband component and the out-of-band narrowband component;

[0030] (6) repeating steps (2)-(5) until narrowband components from all of the subbands of the matrix are removed from the acoustical data; and

[0031] (7) reconstructing the broadband time series;

[0032] whereby to provide separately a narrowband or broadband time series estimate indicative of the presence of a signal in the acoustical data.

[0033] In accordance with a further feature of the invention, there is provided a method for detecting a target in an ocean environment, the method comprising the steps of:

[0034] (1) providing acoustical data by operation of a sonar system, the acoustical data including narrowband components indicative of presence of the target and broadband components resulting from randomly generated sounds;

[0035] (2) constructing a data matrix in which time series samples of acoustical data are arranged in a forward-backward linear predictor matrix;

[0036] (3) decomposing a subband of the matrix into an in-band component and an out-of-band component;

[0037] (4) estimating narrowband components;

[0038] (5) estimating out-of-band leakage;

[0039] (6) removing the in-band narrowband component and the out-of-band narrowband component;

[0040] (7) repeating steps (3)-(6) until narrowband components from all of the subbands of the matrix are removed from the acoustical data; and

[0041] (8) reconstructing the broadband time series;

[0042] whereby to provide a narrowband time series estimate indicative of the presence of the target in the ocean environment.

[0043] The above and other features of the invention, including various novel details and combinations of steps, will now be more particularly described with reference to the accompanying drawings and pointed out in the claims. It will be understood that the particular method embodying the invention is shown by way of illustration only and not as a limitation of the invention. The principles and features of this invention may be employed in various and numerous embodiments without departing from the scope of the invention.

[0044] A more complete understanding of the invention and many of the attendant advantages thereto will be readily appreciated as the same becomes better understood by reference to the following detailed description when considered in conjunction with the accompanying drawings wherein like reference numerals and symbols designate identical or corresponding parts throughout the several views and wherein:

[0045] FIG. 1 is a graph which depicts a typical underwater signal spectrum;

[0046] FIG. 2 is a graph which depicts a multiple taper method and tonal removal;

[0047] FIG. 3 is a graph depicting division into local subbands and leakage;

[0048] FIG. 4A is a graph depicting a signal envelope;

[0049] FIG. 4B is a graph depicting interference and signal spectrum;

[0050] FIG. 4C is a graph depicting a Principal Component Inverse residual spectrum;

[0051] FIG. 4D is a graph depicting a residual spectrum as obtained in accordance with the invention herein;

[0052] FIG. 5 is a graph depicting a power spectrum of simulated interference and signal division into subbands; and FIG. 6 is a graph in which normalized mean-square estimation error of the Principal Component Inverse method is compared with the method described herein.

[0053] In order to exploit the low rank structure of the NB components in the signal and noise data, the time series samples are arranged into a forward-backward linear predictor matrix, the data matrix so constructed appearing as: $\begin{array}{cc}X=\left[\begin{array}{cccccc}{x}_1& \dots & x_{N-L+1}& {\stackrel{\_}{x}}_L& \dots & {\stackrel{\_}{x}}_N\\ x_2& \dots & x_{N-L+2}& {\stackrel{\_}{x}}_{L-1}& \dots & {\stackrel{\_}{x}}_{N-1}\\ ⋮& \dots & ⋮& ⋮& ⋮& ⋮\\ x_L& \dots & x_N& {\stackrel{\_}{x}}_1& \dots & {\stackrel{\_}{x}}_{N-L+1}\end{array}\right]& \left(5\right)\end{array}$

[0054] Assuming that there are NB components present in a kth subband, columns of X are decomposed into an in-band component:

[0055] Z_k^j=V_k^HX  (6)

[0056] and an out-of-band component:

[0057] Z_k^tj=V_k^+HX=[V₁|V₂| . . . |V_k−1|V_k+1| . . . |V_h]^HX  (7)

[0058] as illustrated in FIG. 3. The kth matrix filter V_kis a passband filter designed for the band [f_k−W, f_k+W]. If desired, a more sophisticated filter design can be used, as is known in the art. The filters V_kdo not have to be orthogonal and can be overlapping.

[0059] The Principal Component Inverse (PCI) method of reduced-rank interference nulling is then used to estimate the NB component in the kth subband by solving: $\begin{array}{cc}{s}_k\mathrm{subj}.\stackrel{\min }{\mathrm{to}}\mathrm{rank}\left[s_k\right]=r^{{\Vert Z_k-S_k\Vert }^2}F& \left(8\right)\end{array}$

[0060] A solution can be obtained from the Singular Value Decomposition (SVD) of Z_kand is given by: $\begin{array}{cc}{\hat{S}}_k=\sum _{k=1}^T{\hat{\sigma }}_k{\hat{u}}_k{\hat{v}}_k^H& \left(9\right)\end{array}$

[0061] where the û_kand {circumflex over (v)}_kare the r principal left and right singular vectors, respectively, of Z_k, and the {circumflex over (σ)}_kare the r respective principal singular values. If r tonals are present in the kth subband, then the rank of V_k^HX will also be r, and the matrix corresponding to the tonals can be estimated by finding the best rank r approximation filtered data matrix. Non-tonal NB components can also be effectively modeled as low rank (see L. L. Scharf et al., “Rank Reduction For Modeling Stationary Signals”, IEEE Trans. on Acoustics, Speech, and Signal Proc., Vol. ASSP. 35, No. 3, pp. 350-354, 1987).

[0062] Inasmuch as it is not possible to design perfect passband filters V_k, NB components in the kth subband leak into the adjacent subbands (FIG. 3). This is not desirable. However, one can estimate or extrapolate the NB out-of-band leakage from the in-band NB estimate S_kand then remove it.

[0063] Denoting the BB component in X as N, and the NB component as S, the in-band and out-of-band components can be written as: $\begin{array}{cc}{Z}_k^H=X^H{V}_k=S_k^H\underset{\mathrm{NB}}{V_k}+\underset{\mathrm{BB}}{N^H{V}_k}& \left(10\right)\end{array}$

[0064] and: $\begin{array}{cc}{Z}_k^{†H}·=X^H{V}_k^{†}=S_k^H\underset{\mathrm{NB}}{V_k^{†}}+\underset{\mathrm{BB}}{N^H{V}_k^{†}}& \left(11\right)\end{array}$

[0065] respectively. The formula immediately above (11) is rewritten as:

[0066] Z_k^tH=B^HQ+N^HV_k^t  (12)

[0067] where B is an r×2(N−L+1) matrix whose rows span the row space of S, and Q is an r×K (K is the number of columns of concatenated matrix filter bank V) matrix such that S^HV_k^t=B^HQ. It is then necessary to estimate the leakage V_k^tHS, which might appear impossible since it requires knowledge of S, which is not known. However, estimating the NB components, above, has provided an estimate of the row space of S_kfrom the SVD of Z_kas shown in equation 9, B=[û₁]. . . [û_r]. Therefore, to estimate S^HV_k^t, one substitutes the estimate of B in place of B in (12), and solves for Q by least squares fitting to Z_k^tHusing equation 12. The final estimate of the out-of-band leakage is given by:

[0068] Z_k^t=Z_k^tjB^H(BB^H)⁻¹B  (13)

[0069] where j is the iteration number.

[0070] The in-band NB component

[0071] Z_k^j+1=Z_k^j−S_k, is removed,  (14)

[0072] and the out-of-band NB component

[0073] Z_k^j(j+1)=Z_k^tj−Z_k^t, is removed  (15)

[0074] and the counter is incremented j=j+1.

[0075] The above-described steps of subband decomposition, estimation of NB components, and estimation of out-of-band leakage and removal of NB components is repeated until NB components from all subbands are removed.

[0076] The BB time series is obtained from the cleaned Z, denoted as Z_clean, by first reconstructing the BB time series:

[0077] X_BB=(V^H)#Z_clean  (16)

[0078] present in the columns of X, using the pseudo-inverse of V^Hand then arithmetically averaging all elements in X_BBwhich correspond to the same time series sample using the method developed in D. Tufts et al., supra, for reduced-rank signal enhancement, which is denoted by the operation:

[0079] X_BB=Average [X_BB]  (17)

[0080] The NB time series estimate is then:

[0081] X_NB=X−X_BB  (18)

[0082] Weighted averaging proposed in A. A. Shah et al., “Estimation of the Signal Component of a Data Vector”, Proc. of ICASSP 92, 1992, can be used to obtain even better performance.

[0083] A computer generated, FIG. 4A, a BB transient-like signal 130 samples long that had a triangular shaped spectrum, and then added three interfering complex tones. The signal plus interference spectrum is plotted in FIG. 4B. The new procedure described above to remove the three interfering tones from the BB transient was applied. Briefly, the Nyquist band was divided up into ten non-overlapping subbands using matrix filters based on DPSS's sequences, as described earlier. The spectrum of the estimated BB component for the new method is plotted in FIG. 4D and the PCI method is plotted in FIG. 4C. The new procedure removed the tonals with little distortion to the BB transient, while for comparison, the PCI method applied directly to the time series caused considerable distortion to the transient by locking onto the transient spectral peak and missed one of the tonals.

[0084] Next, a comparison was made of the new method against the Wiener filter method (assuming full knowledge of the NB and BB covariance) and the PCI method for estimating three complex tones scaled by complex Gaussian random variables embedded in colored Gaussian BB noise. The power spectrum of the three tonal signals and broadband noise for Signal to Noise Ratio (hereinafter designated as SNR) 0 db, is measured relative to the power of the three tones and BB component and division into subbands are plotted in FIG. 5. One hundred independent simulation trials were performed using computer generated 130 sample segments of signal-plus-noise to measure the mean-square estimation error (MSE) of the signal estimate for the new method and PCI method, assuming rank 3. The Wiener MSE was calculated theoretically. The MSE normalized by the total signal power is plotted for all three methods in FIG. 6 as a function of SNR (the BB noise component power was varied while the signal power was kept fixed).

[0085] At high SNR, the new method and the PCI method both performed well, being close to the Wiener MSE (see FIG. 6). However, at moderate to low values of SNR starting around 15 db and lower, the new method significantly outperformed the PCI method, while the MSE is still close to the Wiener filter.

[0086] This is because the PCI method began to mistake the local spectral peaks of the BB component for the NB components, while the new method using local subband processing avoided this problem.

[0087] From practice of the above-described method, there is desired a new, robust procedure for automatically separating time series data into its constituent narrowband and broadband time series components. The new procedure avoids the pitfalls of adaptive notch filters, PCI method, or other similar algorithms, of mistaking and filtering local spectral peaks of the broadband component as narrowband components by decomposing the data vector into local subbands by a bank of matrix filters. In piecewise fashion, the narrowband components are estimated and then filtered from each subband using the PCI method. Finally, the filtered components are coherently recombined to obtain the narrowband and broadband time series estimates. Computer simulation results have shown that the new procedure works well and can have performance close to the clairvoyant Wiener filter.

[0088] There is thus provided a method for detecting a contact in an ocean environment and, more generally, a method for separation of acoustic data into a narrowband and broadband time series component, indicative of the contact, and a broadband time series component resulting from randomly generated sounds, i.e., background noise.

[0089] It will be understood that various changes in the details, and method steps, and arrangement of steps, which have been herein described and illustrated in order to explain the nature of the invention, may be made by those skilled in the art within the principle and scope of the invention as expressed in the appended claims.