Application of l1 Inversion to Underwater Acoustic Signal Processing Charles F. Gaumond1 and Geoffrey F. Edelmann2 1

Acoustics Division, Naval Research Laboratory, Washington DC, USA. [email protected] 2 Acoustics Division, Naval Research Laboratory, Washington DC, USA. [email protected] Abstract—A novel method to solve ill-posed problems is presented with a brief mathematical introduction. The method generates l1 solutions that have few strong components in contrast to l2 solutions with minimum energy. The method is compared with l2 method using numerical and at-sea data with various numbers of samples and with background noise levels. [Work was funded by ONR.] Index Terms—Sonar Signal Processing, Inverse Problems, Fourier Analysis

1. Introduction One type of radiated acoustic signal that is potentially detectable is a tone, or single frequency component. In the ocean, the detection of this tone is complicated by fading due to propagation, changes in the mechanical operating state of the source, and Doppler effects due to motion of the source, receiver and acoustic channel [1]. The detection of the tone is limited by the above modulations that cause frequency spreading as well as oceanic noise from shipping, weather and biological sources that create the background above which the tone must be detected. Our mathematical understanding of frequency is based on two approaches to analysis. The first is based on the Fourier integral of functions derived using the L2 norm in a Hilbert space [2]. The second is based on the finite Fourier series, or discrete Fourier transform, that represents a function as a finite number of samples [3]. This letter presents the

application of compressive sampling to finding the presence of a tone in background noise [4]. Tones can arise from rotating machines in ships and the background noise from more distance ships, wind rain and biologics. This application demonstrates that compressive sensing techniques work on real-world problems despite low signal-to-noise ratio (SNR), modest number of measurements, and temporal variability found in underdetermined, inverse problems of underwater acoustics. First, the basic equations are presented and two examples, one involving ideal data and the other at-sea data. The ideal data is numerically generated and the at-sea data is taken from a single element on a towed array. The paper concludes with a summary and discussion.

2. Method Compressive sampling method is applied to the inversion of the following signal model

PROCEEDINGS OF SYMPOL-2011

and

N

y(m) = ∑ W (n , m)Y(n ) , n =1

(1) N

where W (n, m ) = exp(i 2πf n t m ) , using at-sea or data numerically generating with the formula

Y 1 = ∑ Y (n ) .

The minimization of ⎛ 2πm f ⎞ ⎟ + υ(m ) y(m) = A sin ⎜⎜ ⎟ ⎝ N fs ⎠

(2)

where the amplitude A and frequency f are arbitrary and unknown. The second term υ is assumed to be Gaussian, zero-mean additive noise. If we assume that there are M < N time samples. Thus there are fewer samples of y than the number of frequency samples in the search space. Note that in the discrete Fourier transform the number of frequency and time samples are the same and that the samples are equally spaced. In the discrete Fourier transform, t m = m f s and f n = nf s N where f s is the sampling frequency. However, in the case considered here, the matrix W is not square and the inversion of Eqn. (1) is illposed. In this case, N > M and the inversion generates finer resolution of the spectrum than appears to be supported by the limited data. However, the ability to generate finer resolution depends on the overall duration of the time sample with M samples taken from that interval.

(4)

n =1

Y 1 generates the

solution that has the fewest significant components. On the other hand, the minimization of Y 2 generates the solution that has the least energy. Thus compressive sensing method generates another class of solution from the minimum energy method. The solutions generated by the compressive sensing method are called sparse solutions. The use of minimum energy or compressive sensing therefore depends on prior assumptions about the solution. If the appropriate solution has minimum energy then that method should be used. However, if the sparse solution is more appropriate, then compressive sensing should be used. The data in Eqn. (2) have a single frequency component and random, additive noise. Clearly the assumed solution is sparse rather than minimum energy. Note that more than one tone could be present and the exact number need not be known a priori.

The numerical solution of this ill-posed system of equations requires an additional condition on the solution. The Moore Penrose inverse is often used to solve an underdetermined set of equations [5]. This inverse invokes the condition that the solution is a minimum l2 norm Y 2 .

The set of M sample points should be chosen at random, namely with equal probability form the N samples determined from the Nyquist sampling theorem. This method of choosing sample points is called the uniform uncertainty principle (UUP) [4]. The number of samples required to reconstruct of spectrum with fewer than K samples is proportional to K log 2 N for the general case where UUP holds [4].

In contrast to that approach, compressive sampling assumes the minimization of the l1 norm Y 1 so that the solution is

3. Numerical Result

yˆ(m) = arg min Y 1 , such that y = WY , (3)

The testing of the algorithm is done by performing the inversion of Eqn. 1 using data generated with Eqn. 2. The sampling frequency

Gaumond et al.: Application of Compressive Sampling to Underwater Acoustic Signal Processing

is 6125 Hz, the frequency samples are equally spaced {f n = nf s N , 0 ≤ n ≤ N − 1} and the time samples are chosen from {t m = m f s , 0 ≤ m ≤ M − 1} . In the top panel of Fig. 1, all M time samples are used. Note that the inversion is perfect. The circle shows the amplitude A and frequency f used in Eqn. 2 to generate the data vector y .

Figure 1. Numerical examples of inversion of Eqn. 1. The top panels show the inverse of full W . The middle panels show the MoorePenrose inverse and the bottom panels show the l1 inversion. The in-phase and quadrature components are displayed to show the phase. For the middle and lower panels, 32 time samples were chosen at random from the 2048 equally spaced time samples. The inversion was performed using the YALL1 algorithm [6]. The YALL1 algorithm requires Eqn. 1 to be regularized so that it is self-adjoint. In this paper we used the QR method, with W' = QR so that

R'\y = R'\R' Q' Y ,

(5)

and therefore the regularized data is R'\y and the regularized W matrix is R'\R' Q' . After this regularization, the middle panels show the resulting inversion using the Moore

Penrose inverse. Note that the solution has a definite peak at the correct frequency, but there is a significant amount signal energy spread throughout the spectrum. The bottom panels show the result of using the YALL1 algorithm for minimum l1 inversion. The spectrum is perfectly reconstructed because 32 is a sufficient number of sample points to reconstruct the sparse spectrum that contains one non-zero component. Another numerical demonstration shows further degradation resulting from a limited number of samples. In the previous example, a data-deficient, non-square matrix was regularized to enable the minimum energy inversion of Eqn 1. The data deficiency arose from randomly choosing 32 time samples from a set of 128 equally spaced samples. This example shows the degradation due to choosing the first 128 and using the inverse Fourier transform. This is compared with the use of l1 minimization with 128 randomly chosen samples. Figure 2 shows the results in a similar format as in Fig. 1. In this case the inversion of the square W matrix using 128 equally spaced time samples and 128 sample frequencies produces a poor reconstruction because the frequency f of the data lies exactly between two of the coarsely sampled frequencies. In the lower two panels a 128 by 2048 W matrix, with finer frequency resolution and randomly spaced time samples, is used. The time samples are, as before, randomly chosen from 2048 samples that are equally spaced in time. The difference between the lower and upper is remarkable, with the l1 method achieving the correct result while the l2 method completely fails. This situation was set up for the failure of the l2 method is expected, of course, because the frequency was chosen using sampling theory to be pathological. The example does demonstrate the ability to achieve finer PROCEEDINGS OF SYMPOL-2011

frequency resolution with coarser time sampling. Note that the time samples are drawn from the time interval 2048 / f s . The ability to achieve finer frequency resolution appears to depends on the overall duration of the time series.

the other frequencies and over different noise and sets of time samples.

Figure 3. The output SNR as a function of input SNR with the number of frequency samples held constant at 2048 and the number of time samples at 128.

Figure 2. Numerical examples of the inversion of Eqn. 1. The top panels show the inversion of the 128 by 128 W matrix using equally spaced samples. The lower panels show the l1 inversion of the 128 by 2048 W matrix with randomly chosen time samples. The last case will show the output SNR as a function of input SNR with the numbers of frequency samples and time samples held constant. The input SNR is defined as SNR input =

A std(ν )

Note that there is linear relationship between the input and output SNR. The inversion fails for input SNR below 8, on a linear scale. Note also that at higher values of input SNR, there is a moderate increase in SNR generated by the inversion. The input and output SNRs are equal with an input SNR of approximately 15 dB.

(2)

And the output SNR is defined as SNR output =

mean( Y(f f s ) )

std (Y(m ≠ m 0 ))

(2)

with m 0 = f f s and the mean of Y at the frequency of the data is taken over different noise and sets of time samples. Similarly, the standard deviation of Y is averaged over all

Figure 3. Comparison of output to input signal to noise ratios as a function of the number of random data samples.

Gaumond et al.: Application of Compressive Sampling to Underwater Acoustic Signal Processing

4. Experimental Results Data were acquired from the OREX05 experiment that was performed off the coast of Oregon in 2005. The received signal from one element of a towed array is shown below in Fig. 4. Note that there are very loud signal components below 500 Hz; these are active transmissions and will be ignored here. There is one frequency line near 1300 Hz that was radiated from a nearby ship. This data is bandpassed filtered, using a 10% Tukey window from 1 kHz to 2 kHz. The first 2048 samples of the resulting data containing this frequency line are analyzed using l2 and l1 methods.

new W matrix for each set of t m . The coherent average appears to lower the background noise, namely the signal other than the tone line. The l1 spectrum of background noise is in general lower and also has a different appearance. The l1 spectrum tends to be composed of spikes rather compared with the l2 spectrum.

Figure 5. Comparison of the l2 inversion using 2048 data points with the coherent average of four l1 inversions using 512 data points. The vertical is linear.

Figure 4. Spectrogram of a signal from one element of a towed array. The intensity is proportional to the absolute value of the spectrogram. The comparison in Fig. 5 used 512 randomly chosen from 2048 time samples with l1 minimization of Eqn. 1. There are four groups of 512 distinct time samples, and the resulting inversions were coherently averaged to generate the results shown in the lower panels of Fig. 5. Ordinarily a power spectrum is generated by incoherently averaging overlapping spectra. In the case of l1 inversion presented here, the four sets of randomly chosen samples are coherent with each other. Each of the four sets is taken coherently from the same 2048 set. This coherence requires a

Figure 6. Comparison of the results shown in Fig. 5. The signals are compared with a dB scale. The peaks of each spectrum is normalized to one. The l2 spectrum is a dotted line and the l1 spectrum is a solid line.

PROCEEDINGS OF SYMPOL-2011

The above figure shows the comparison of the two spectra in terms of decibels. Note that l1 method suppresses all of the background noise. Note that the same set of data was used to generate each spectrum. The suppression is not uniform; there are peaks present in the l1 results.

5. Discussion and Conclusions A recent numerical method of inverting a matrix is presented. This compressive sensing method generates a class of solution that differs from that generated with the minimum energy, or l2, method. The solutions generated by the compressive sensing method are called sparse solutions. This method uses a minimization of the l1 norm when the problem is underdetermined. When the data is caused by a sparse process, the underdetermined inversion can be determined using the l1 technique. Sparseness can occur in spectral analysis because several mechanical processes, like gear noise or other rotating machinery, produce repetitious signals. Sparseness can also occur in beamforming when a few strong signal sources are placed at discrete bearing angles from an array. Sparseness does not always occur. Random background noise is not sparse in frequency. If the signal is dominated by a sparse signal, the signal can still be detected and reconstructed . The success of the inversion requires the SNR in the frequency domain to be sufficiently high. This is shown in Fig. 2, the sufficient level of SNR is approximately 8 dB for the specific, numerical case treated. The success of the inversion also requires a sufficient number of samples. This situation is shown in Fig. 3, again for the specific case considered. The method is also demonstrated with data taken at sea where the signal frequency is not known and the background noise is not Gaussian and white. The ability to coherently average results from independent sparse

samples is also demonstrated in Fig. 5. The suppression of real world noise is shown in Fig. 6. In our opinion, there are many potential uses of this method.

ACKNOWLEDGMENTS The authors acknowledge support with data from Laurie Fialkowski. This work was supported by the Office of Naval Research.

References [1] A. Gerlach, ‘Impact of the Ocean Acoustic Transfer Function on the Coherence of Undersea Propagations,’ IEEE Trans. Acoust. Speech and Sig. Proc, Vol. ASSP28, 1980, pp. 145-158. [2] D.C. Champeney, ‘A Handbook of Fourier Theorems,’ Cambridge University Press, Cambridge, 1987, pp. 1-43. [3] ibid., pp. 155-166 [4] E. Candes, J. Romberg and T. Tao, ‘Robust Uncertainty Principles: Exact Signal Reconstruction from Highly Incomplete Frequency Information,’ IEEE Trans. Inf. Theory, Vol. 52, 2006, pp. 489-509. [5] G. Golub and C.F. Van Loan, ‘Matrix Compuations,’ Johns Hopkins University Press, 1989, pp. 246-249. [6] J.F. Yang and Y Zhang, “Alternating direction algorithms for l1 problems in compressive sensing,” CAAM TR09-37, 2009, arXiv:0912.1185v1; software available as download at http://yall1.blogs.rice.edu/ (accessed June 1, 2011).

Gaumond et al.: Application of Compressive Sampling to Underwater Acoustic Signal Processing