Deconvolution for transient classification using fourth order statistics

Considered are numerical experiments involving an exponentially damped harmonic source signature, and underwater acoustics Green's functions (GFs) representing bottom-limited propagation in a wave guide. It has been shown that the fourth order blind deconvolution functional proposed by Wiggins (1978) can be used to generate a class of signals from the propagation distorted signature such that one or more members of the class are "close" to the original signature. There are obvious applications to improving sonar target classification in the presence of multipath distortion. In this paper it is demonstrated that the algorithm based on the Cabrelli (1984) functional (D-norm) can be implemented in a similar fashion, and that it performs significantly better in the presence of additive noise than does the Wiggins V-norm.


I. INTRODUCTION
An area of interest in several fields of underwater acoustics signal processing is the accounting for environmental distortion introduced into a propagating pulse. One useful mathematical model representing this process is that of convolution, which leads to the consideration of the deconvolution problem [2]. When only the received signal is available, and at best, perhaps statistical information about the system transfer function or input is known, then the problem under consideration is called a "blind" deconvolution problem.
In previous work we have discussed the blind deconvolution technique introduced by Wiggins [ 11, and its application to underwater acoustics problems [3,4]. Wiggins' intuitive notion was to use the Varimax or V-norm from Factor Analysis as a measure of sparsenessa statistical property his signal model possessed (reflection seismogram). He termed this method "Minimum Entropy." The idea was to develop a filter that would drive the received signal to a sparser representation and, hence, recover the underlying Green's function (impulse response function (IRF) in engineering terminology). Such a filter gives information about the source signature as well.
By the time of Walden's review paper [SI, it was understood that Wiggins method could be viewed as a higher order statistics method that exploited non-Gaussianity in the IRE; since it essentially seeks to maximize kurtosis. This is appropriate, since positive kurtosis (lepto-kurtosis) is a measure of the degree of departure from a Gaussian distribution. A lepto-kurtic distribution has a longer tail, and long-tailed distributions are appropriate for representing the sparseness (or simplicity) in the structure of the IRFs that was being assumed by Wiggins. In 1985, Cabrelli [6] introduced a method similar to Wiggins' that is apparently little known to the blind deconvolution community. He developed a measure of simplicity (or, sparseness) that he referred to as the D-norm, and arrive at a significantly different algorithm, but similar solutions to the V-norm for high SNR test cases. However, when SNR decreased, Cabrelli noted that his algorithm was more robust.
Wiggins and Cabrelli's blind deconvolution methods were implemented and applied in an underwater acoustics setting, as opposed to the reflection seismology problem, where they were developed. We show that this is appropriate when the Green's functions (or IRF's) have sufficient positive kurtosis, which is the case in typical shallow water propagation scenarios for ranges that are not too great. It is shown that the algorithms can be used to generate a solution class that could be exploited by an automatic classification system. Also, the V-norm and Dnorm functional's relative performance is studied as SNR decreases.
In the present work, II. SIGNAL MODEL When a sound pulse travels through the ocean, especially in shallow water, it will scatter back and forth between the sea surface and the ocean bottom. Further interaction with strata in the ocean bottom can also occur, conspiring altogether to produce a received pulse that can look quite different from the original source signature. An automatic pattern classifier will often be confused by this environmental distortion, and fail to properly classify the source signature.
To account for the effects of environmental distortion in signal processing schemes, a signal model is needed. In this paper, we do not explicitly consider the noise in the model, but rather, seek algorithms that are stable in the presence of noise. We concentrate, then, on accounting for wave propagation effects in the signal model (primarily, multipath), by using the wave equation Let us define the Green's function G(.), such that and assume a separable source Then, it can be shown that D&r,tIr',t') = 6(r-r')6(t-t') F ( r , t ) = F(r)S(t).
Using causality and assuming time invariance, the convolution relation is obtained: The convolution equation provides the basis for our signal model, where S(t) is the source signature initially propagated, G( r, t , t') encodes the multipath environmental distortion effects, and v(r, t) is the received signal at range r , and depth z . The source is at the origin at depth z', and r = ( r , z ) (using cylindrical coordinates, and assuming azimuthal symmetry). In a briefer form, where it is understood that G(t) is associated with some fixed, specified source-receiver geometry.

III. BLIND DECONVOLUTION PROBLEM
We will, in general, assume knowledge of two things: (1) the received signal, and (2) some appropriate, exploitable statistical characterization of G ( t ) (but, not G(t) itself). Under these "blind" conditions, we may, at best, only be able to produce a finite class of candidate solutions. We would then define success as there being one or more members of this solution class that are sufficiently close to S(t) to be useful.
This reduced goal may not be adequate in general problems, but for classification, the solution class can be presented to the classifier to see if there is any recognition. Then, in a sense, the classifier and its training signature set become part of the information exploited in our statistical approach. Of course, additional information or assumptions about the signal model invoked will help reduce ambiguity in the solution class. Also, when multiple received signals, over significantly different transmission paths, are available, significant ambiguity reduction can probably be accomplished. We will address this issue elsewhere.

IV. STATISTICAL SIGNAL MODEL
Wiggins [l] gives the key idea that we employ here in a slightly generalized form. We consider G ( t ) to be a realization of a non-Gaussian random process, input to a linear system whose impulse response is S ( t ) . From various (unpublished) calculations, we know that non-Gaussianity is appropriate because of the high, positive kurtosis of G ( t ) . This is due to its sparseness (or simple structure, or lower entropy) relative to a Gaussian random signal.

How, then, do we accomplish the estimation of S(t)?
The following approach is a generalization, and in some sense, an extension of Wiggins' approach. We wish to design a filter f (t) that increases the sparseness measure of the received signal. The hope is to "drive" the received signal towards G(t) . The rationale is that the convolution process is a smearing process that moves the signal to less sparseness (more Gaussian). Increasing the non-Gaussianity should then drive yr(t) towards G(t).
After deconvolution is used to estimate the filter, an estimate of S(t) is easily provided by the inverse filter f -' ( t ) . Let sp{.} mean some sparse measure. We seek f ( t ) such that To do this, we can use a gradient operator, for example. If Comparing this with our model Wiggins used the idea of "simplicity" or "sparseness" borrowed from factor analysisnamely, the V-norm v(Y)=CY: CY, ' . /I: 1 Here we give Cabrelli's statement of Wiggins' algorithm (also see Walden [ 5 ] ) . Consider N observed signals x I . . . x N .
for each i (i=1, ..., N ) , and let xi berepresented by xi = w*qi + q i , where w is the source signal. Suppose that each signal xi is convolved with the same filter f in order to obtain an output y, = f * X i =(f*w)*q;+f*q;. To determine the filter, the Varimax criterion is then applied to the outputs yi in order to maximize V(y) over all filters The filters are applied to the input to generate the YCij), and the algorithm terminates by computing max( hl p I/IIY 1 1 ) .  Fig. I (a) The short pulse signal and its energy spectrum. (b) The exponentially-damped sinusoid signal and its energy spectrum.

VI. DECONVOLUTION TEST RESULTS
Two test transient signals are used to compare the performance of the Wiggins and Cabrelli deconvolution methods. The first transient is a short duration pulse shown in Fig. l(a). The second is a longer duration exponentiallydamped sinusoid, shown in Fig. l(b). Multipath propagation distorted versions of the signals are generated by convolution of the two source transients with three model Green's functions. The Green's functions were created using a parabolic equation model and geoacoustic parameters from an experimental site in the Atlantic Ocean at the southern end of the Blake Plateau [2]. The receiver depth is 250 m in a 915 m water column, and the Green's functions are calculated at 600 m, 4300 m, and 7900 m upslope from the receiver.
Thus, the Green's functions, shown at the three source-toreceiver ranges in Fig. 2, are realistic representations of those that occur in the ocean. Note that the spacing of consecutive path arrival times is significantly smaller than the time duration of either test source signature, indicating that the multipath arrivals in the propagated signals will not be visually separable.
The two blind deconvolution methods evaluated in this paper require sparse Green's functions, or high kurtosis, to work successfully. Each of the three Green's functions has a significantly large value of kurtosis, at 124.63, 37.98, and 27.33 for the 600 m, 4300 m, and 7900 m ranges, respectively.
To evaluate the performance of the deconvolution algorithms, the source estimate generated by the algorithm is compared to the known source signal using the absolute value of the correlation coefficient, which is one when the two signals are equal. This is a simple, but well-accepted measure of data similarity, and although classifiers that would be used in practice will generally be much more sophisticated and involve many signal features, the correlation coefficient is a convenient numerical measure of performance that lends itself to straightforward comparison between the deconvolution methods, and is thus appropriate for use in this analysis. The Wiggins method requires some level of prewhitening during the autocorrelation inversion for some of the examples. A value of 0.01% is found to be sufficient for stabilization in the cases where pre-whitening is required, and is thus used for all cases with the Wiggins method to preserve consistency and reflect a realistic application scenario where instabilities cannot be predicted. Prewhitening is not required for the Cabrelli algorithm, and in some cases, leads to much poorer results. Filter lengths in each method range from 1 to 50 points. The best result from the source estimates produced with the various filter lengths is chosen.
Blind deconvolution methods would generally be used to improve a received source prior to classification. Thus, it is reasonable to expect that an estimate of the signal passband would be available for use in the deconvolution routines. As such, the source estimates produced by the two deconvolution methods are subjected to a bandpass filter before classification, i.e., before the correlation coefficient is calculated. For the pulse, the passband filter is defined over 25-150 Hz, and for the damped sinusoid, over 0-100 Hz.
The results from the two methods are given in Table 1

VIII. SIMULATIONS WITH NOISE
To assess the performance of the Wiggins and Cabrelli methods in noise, simulated independent Gaussian noise is generated and added to each received signal before application of the deconvolution algorithms. The S N R is defined as the ratio of signal and noise standard deviations and converted to decibels (dB) for display. In these simulations, the bandpass filter is applied twice: first, it is applied to the received signal before input to the deconvolution algorithm, and second, it is applied to the source estimate generated produced by the algorithm.
The results for the pulse signal and S N R s from -10 dB to 30 dB are shown in Fig. 3. At all three ranges, both the Wiggins and Cabrelli methods produce source estimates that are more accurate than before deconvolution. Note that the solid curves depicting correlation coefficients before processing include the initial application of the bandpass filter. At the 600 m range, the Wiggins method produces slightly better results than the Cabrelli method at SNRs above 22 dB, but the Cabrelli method produces better results at all the lower SNRs.
In the 4300 m range case, the Wiggins method produces better results at S N R s above 14 dB. In the 7900 m case, the Wiggins method performs best at SNRs between 14 dB and 20 dB, but the Cabrelli method performs best at all other SNRs. Below about 8 dB at all three ranges, the best Wiggins methods source estimate occurs at a filter length of 1 point (correlation coefficient = 0.667). While this result is significantly better than no processing, it is an artificial result due to the pulse-like nature of the source signal, and the particular biases of the correlation coefficient. Thus, this moderately high correlation coefficient of 0.667 is strictly signal dependent and cannot generally be expected to occur in general. In fact, it does not for the second test signal used in this paper. Fig. 4 depicts the deconvolution results for the damped sinusoid signal. For this signal, the Wiggins method has a significant problem producing acceptable results below an SNRs of 10-14 dB. While the Wiggins method performs slightly better than the Cabrelli method at S N R s above 28 dB for the 600 m and 4300 m ranges, the Cabrelli method shows significantly superior results at lower values of S N R . In fact, the Cabrelli method produces source estimates that are significantly better than before deconvolution at SNRs as low as -10 dB.  Fig. 4 Correlation coefficients between the known damped sinusoid signal and source estimates produced by the Wiggins (x) and Cabrelli ( 0 ) methods in Gaussian noise. The solid lines indicate correlation coefficients before processing.

IX. CONCLUSIONS
Using numerical experiments involving an exponentially damped harmonic source signature and underwater acoustics Green's functions representing bottom-limited propagation in a wave guide, it was shown that deconvolution methods based on fourth order statistics can be used to generate a class of source signature estimates from the propagation distorted signature. A unique solution is not obtained. However, if the IRF kurtosis is sufficiently large, at least for the simple signatures tested so far, it tends to be true that several members of the class are "close" to the original signature. This should be exploitable in improving sonar classification performance in the presence of certain types of environmental distortion.
It is also demonstrated that the algorithm based on the Cabrelli functional (D-norm) performs significantly better in the presence of additive noise than does the Wiggins Vnorm. For the damped sinusoid signal, the Wiggins method fails at even moderate levels of SNR, producing source estimates that are more distorted than the input signals. The Cabrelli method, however, performs well for this signal at SNRs as low as -10 dB.
Techniques to exploit redundancy in multichannel systems to reduce the ambiguity in the solution class (nonuniqueness) are currently being investigated. Preliminary results appear promising.